多元统计笔记

多元统计笔记

部分笔记，其他有时间整理，聊作记录

2025-5-12 实验5讲解及第十章

1、显著性的适用是需要检验过的，需要p值（或者是可以直接很明显的看出来显著）
2、多个总体的均值是否存在显著差异？（可检验，有差异，则可以作判别分析，如此，若无差异，则不能适用判别检验，或者判别检验效果不行）
3、异常值，箱型图，分位数...
4、判断在拥有多个典型判别函数时，是不是都用？检验。
5、判别结果一样，但是后验概率会发生改变

第十章 多元回归与多元方差分析

多元回归

1、名称基于响应变量，即因变量为多元的回归分析；多元回归与多重回归（自变量是多个但是响应变量是一个），二者是不同的，注意术语的区别。
2、多元回归更有效，估计的标准误更小。多元回归的到的收益依赖于correlation，即相关性。但其缺点，参数更多，统计推断更复杂，对自由度的要求更高。如果相关性够强，多元，否则，一个一个作回归也是没有问题的。
3、协变量，
4、预测变量是同一组，那怕某组关联性很弱，原因在于线性代数上表述的清晰，适用矩阵代数知识。
5、参数估计与推断，最小二乘与最大似然；估计过程与一元相似。
6、多元\(Y = ZB + E\)，一元回归\(y = x \beta + \sigma\)
7、多元回归退化为Manova的条件？方差分析可以看做特殊的线性回归，一元anova中，其也是两样本t检验的推广，比较多个组时，使用anova。
目标，不同的组之间是否存在显著的差异，或者说，均值是否存在显著差异。
8、MANOVA，多元方差分析，一元方差分析的推广。【整合了响应变量之间的相关性】
9、anova的推断，依赖于方差分解，将响应变量中的总的变异分解为组间效应与组内效应(护着叫做响应误差)。
检验统计量构建组间与组内的比值（服从F分布）；要求响应变量服从正态，并要求方差齐性。（方差齐性假设）
10、MANOVA，数学表达式更复杂，检验统计量，以近似分布代替。

代码示例

##########---------Example 1--------------------
options(digits=2)
quiz.dat <- read.table("class.dat4", head=T)
head(quiz.dat)
quiz.dat1 <- na.omit(quiz.dat)
dim(quiz.dat1)
cor(quiz.dat1[, 4:8])
attach(quiz.dat1)
# MANOVA is done via R function manova
fit <- manova(as.matrix(quiz.dat1[, 4:8]) ~ factor(gender)*factor(major))
summary(fit)
summary.aov(fit) # 输出每个变量的一元方差分析

##########----------------------------------------
# multivariate linear model with 5 quizzes as responses
fit.m <- lm(cbind(quiz1, quiz2, quiz3, quiz4, quiz5)~gender*major, data = quiz.dat1)
summary(fit.m)
# linear regression with (quiz4, quiz5) as a multivariate response
fit2.m <- lm(cbind(quiz4, quiz5) ~ quiz1+quiz2+quiz3+gender+major, data = quiz.dat1)
summary(fit2.m)

##########----------------------------------------
# hypothesis test
library(car)
# test if there is a difference between male and female students in the quiz scores
linearHypothesis(fit.m, "genderM")
# test if the averages of Math major students and Stat major students
# are significantly different from Comp major students
linearHypothesis(fit.m, "0.5 * majorMath + 0.5 * majorStat")

例二

##########---------Example 2-------------------
consum.1 <- read.table("consum2007.txt", head = T)
[consum.1[1:3, ]]
consum.east <- subset(consum.1, East==1) [, c(1, 2, 3, 6,7)]
consum.west <- subset(consum.1, west==1) [, c(1, 2, 3, 6,7)]
consum.central <- subset(consum.1, (East!=1)&(west!=1)) [,c(1, 2 ,3 ,6 ,7)]
# compute the mean vectors for three areas
(consum.east.mean <- colMeans(consum.east))
(consum.west.mean <- colMeans(consum.west))
(consum.central.mean <- colMeans(consum.central))
consum.east$Area <- 1
consum.central$Area <- 2
consum.west$Area <- 3
consum.s <- rbind(consum.east, consum.central, consum.west) # 按行合并
Y <- as.matrix(consum.s[ , 1:5]) # 5 response variables
manova.con <- manova(Y ~ factor(Area), consum.s) # 这的factor(Area)是必须的，不然会将其认为连续的
summary(manova.con)
summary.aov(manova.con)
# testing equality of covariance matrix
(S1 <- cov(consum.east[, -6]))
(S2 <- cov(consum.west[, -6]))
(S3 <- cov(consum.central[, -6]))
(n1 <- dim(consum.east)[1])
(n2 <- dim(consum.west)[1])
(n3 <- dim(consum.central)[1])
(p <- length(consum.east.mean))
(Sp <- ((n1-1)*S1 + (n2-1)*S2 + (n3-1)*S3)/(n1+n2+n3-3))
(lRM <- (n1-1)/2*log(det(S1))+(n2-1)/2*log(det(S2))+(n3-1)/2*log(det(S3))
      -(n1+n2+n3-3)/2*log(det(Sp)))

k <- 3
c1 <- ((n1-1)+(n2-1)+(n3-1)-(n1+n2+n3-3)*(2*p^2+3*p-1)/(6*(p+1)*(p-1)))
c1
(df=1/2*(k-1)*p*(p+1))
pvalue <- 1-pchisq(-2*(1-c1)*lnM,df)
pvalue
## 利用Box'M 检验协方差齐性
library(heplots)
boxM(consum.s[,1:5],consum.s[,6])

附录

数据

class.txt

2025-5-19 实验6讲解及第九章

广义线性模型

非线性模型允许响应变量与协变量之间是一个非线性的关系，但是要求响应变量服从正态分布
在实际中Y不一定是连续的，不一定是服从正态分布的。

\[LM:Y = x \beta+\epsilon\\ NLM:Y = g(x,\beta)+\epsilon\\ GLM:Y:N(\mu,\sigma^2)B(1,p)Poisson(\lambda)\\ \]

指数族分布 The Exponential Family

不是指指数分布，只是写出来是指数形式
指数族分布包含正态分布，泊松分布，二项分布，包含很多普通分布，指数分布是其广义形式。

定义：

\[f(y|\theta, \phi) = \exp\left\{\frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi)\right\}, \quad (1) \]

其中a,b,c与其参数唯一依赖。
\(\theta是位置参数，\phi是扩散参数。\)

作业

利用密度函数的正则性，积分为1；

\[\int f(y|\theta, \phi)dy = 1\\ EY_i = \int y f(y|\theta, \phi)dy = b(\theta)\\ \]

推导：

\[E(y_i) = \mu = \frac{\partial b(\theta)}{\partial \theta}, \quad \text{Var}(y_i) = a(\phi) \frac{\partial^2 b(\theta)}{\partial \theta^2}. \]

正态分布，泊松分布，二项分布验证上述结论

二项分布时的\(\theta\)是关于\(\mu\)的一个logistic变换。

\(\theta\)决定广义线性模型中的连接函数怎么取
GLM中方差和期望是有关系的，方差是期望的函数。

非常重要：回归模型的类型是由谁决定的？响应变量，而不是预测变量

GLM的一般形式

线性预测器：

\[g(E(Y_i)) = {x_i}^{T} \beta\\ \]

连接函数：依赖于\(\mu_i\)

\[g(E(Y_i)) = \eta_i\\ \]

对响应变量取对数不是广义线性模型，如对收入和GDP取对数，不是广义线性模型。值域对应到负无穷到正无穷。
广义线性模型（GLM）可以写成：

\[g(\mu_i) = x_i^T \beta, \quad i = 1, 2, \cdots, n, \quad (2) \]

其中 (g(\cdot)) 是单调可微函数，称为链接函数。因此，GLM有两个组成部分：

• 响应变量 (y_i) 服从指数族分布；
• 链接函数 (g(\cdot)) 描述了均值响应 (\mu_i) 如何与线性预测子 (\eta_i) 相关。
  评论：
• 链接函数转换 (E(y_i))。
• 链接函数不转换 (y_i) 本身。

估计

似然法

\[l(\beta, \phi) = \sum_{i=1}^{n} \left\{ \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi) \right\}. \]

右边看不到\(\beta\),它是一个隐函数

\[g(E(y_i)) = {x_i}^{T} \beta\\ E(y_i) = \mu = \frac{\partial b(\theta_i)}{\partial \theta_i} \]

似然估计的常用大样本性质；consistent, efficient,
and asymptotically
小样本：bootstrap,自助法

GLM的统计模拟

课后习题9.6

logistic回归的NLE的表现怎么样
模拟代码：

options(digits=2)
n<-20
#n<-40
#n<-100
b0<-1
b1<-2
y<-rep(0,n)
##一次模拟的结果
x<-rnorm(n,0,1)
#x<-rt(n,3)
#x<-runif(n, min=-1, max=1)
p<-exp(b0+b1*x)/(1+exp(b0+b1*x))
p
for (i in 1:n) y[i]<-rbinom(1,1,p[i])
B<-glm(y~x,family=binomial)$coefficients
B
(bias<-c(B[1]-b0,B[2]-b1))
(mse<-c((B[1]-b0)**2,(B[2]-b1)**2))

##N次模拟的结果
# start_time <- sys.time() #记录运行时间
N<-1000
bias1<-rep(0,N)
bias2<-rep(0,N)
for (j in 1:N){
  x<-rnorm(n,0,1)
  #x<-rt(n,3)
  #x<-runif(n, min=-1, max=1)
  p<-exp(b0+b1*x)/(1+exp(b0+b1*x))
  for (i in 1:n) 
  y[i]<-rbinom(1,1,p[i])
  B<-glm(y~x,family=binomial)$coefficients
  bias1[j]<-B[1]-b0
  bias2[j]<-B[2]-b1
}
(bias1.mean<-mean(bias1))
(bias2.mean<-mean(bias2))
hist(bias1)
hist(bias2)