首先先明确几个概念(例子以SLR说明,MLR同理)

1. Y=β₀+β₁X₁+β₂X₂+...β_pX_p+e

此式代表着数据点的值,也就是观测到的值, p代表predictor的个数, X_p在矩阵中实际上为X_i,p, i代表第i个数据
这里β₀,β₁,e等都是未知的,只有观测到的每个数据(Y_i,X₁,X₂,...X_p)是已知的

2. β₀,β₁与 $\hat{Y}$

β₀代表Y轴上的截距,可能并没有含义(具体看定义)
β₁代表X₁每上升1,其它X固定,对应Y上升量为β₁
由于β₀,β₁是未知的,我们采用 $\hat{β_0},\hat{β_1}$ 来模拟β₀,β₁的真实值获得拟合直线.显然当 $\hat{β_0},\hat{β_1}$ 带入1.的公式中,此时的Y也应该用 $\hat{Y}$ 来表示. 这里 $\hat{Y}$ 实际表示 $\hat E(Y|X=x)$ ,也可写成 $\hat μ_{Y|X}$ ,表示X=x时的平均值
例如: 假设我们有10个病人,我们想知道吃完降压药后血压下降和时间是否有线性关系,X₁可以取几个值: 10,20,30,50min.... 通过公式得到的 $\hat E(Y|X=40)$ 实际上代表10个人血压在40min时的平均值.

3. $\hat{Y}_{pred|x_0}$

显然 $\hat{Y}$ 仅能作为预测的样本平均值.对于每一个个体如果我们选取X₀,对应Y用 $\hat{Y}_{pred|x_0}$ 来表示.虽然我们不知道 $\hat{Y}_{pred|x_0}$ 的具体值,但是我们能估计其分布:
(1)均值:均值为 $\hat μ_{Y|X}$ **,且为 $μ_{Y|X}$ 的无偏估计;
(2)95%CI: X=x₀时对应的95%CI 为 $\hat μ_{Y|x_0} ± t_{n-1-p,1-α/2}* S_{\hat Y|x_0}$
(3)方差:variance由两部分组成: 个体的差异 $σ^2$ (如每个人血压差异,同一个人每次量血压也有差异) + 抽样样本对真实均值的偏离 $Var(\hat μ_{Y|X})$
$Var(\hat{Y}_{pred|x_0})=Var(\hat μ_{Y|X})+σ^2$
具体怎么算可以自己推导,这里很容易理解:σ^2不会变化,但抽样的样本数量增加,样本均值会更接近总体均值, $Var(\hat μ_{Y|X})$ 会变小

使用随机生成的5个数据作为示例

x1 <- rnorm(5)
x2 <- rnorm(5) #按标准正态分布生成5个数字
y <- c(1,2,3,4,5)
xy.lm <- lm(y~x1+x2) #选择因变量和自变量

#直接进行分析
summary(xy.lm)
#或
anova(xy.lm)

结果解读:

image.png

先看Summary
Coefficients中β₀,β₁,β₂,即左边对应的系数.estimate就是

\hat β

,然后是

\hat β

的标准差.
t-test和F-test没有本质上的区别.如果T~ T_n-2 那么T²~ F_1,n-2 .这里可以看到t²的值=F,p值也是相同的.
H₀为β_i=0 (Y和X_i没有线性相关关系)

跳过summary先讲Sum Sq和Mean Sq
SSEresidual(残差)的平方和,即 $\sum(y_i-\hat y_i)^2$
SSR来自回归模型(regression) $\sum(\hat y_i-\overline y_i)^2$
SSY/SST 来自Y的总体: $\sum(y_i-\overline y_i)^2$
MSR=SSR/p
MSE=SSE/(n-1-p)
F=MSR/MSE 如果β₁=0,模型会是一条水平直线,导致MSR=0.也就是说越是β₁≠0,F会越大,就可以拒绝H₀,但不能说明拟合程度好
R²: R²=(SSY-SSE)/SSY=SSR/SSY 意思是X回归的方差占Y的总方差的比例,代表在回归直线上的集中程度 (简单理解为误差比例越小,拟合度越好)
Adjusted R²: R²_adj=1-(SSE/(n-p-1))/(SST/(n-1)) 可按预测自变量数目p进行调整. 当X很多时R²会过大导致无法代表拟合程度,但是会限制额外的自变量.
注意anova()使用的是Type I,和变量的顺序有关.更常使用anova(lm, type="III") 或anova(model1, model2) 研究每一个变量对模型的贡献

分步计算过程

找找对应关系
想想怎么计算 $\hat{Y}_{pred|x_0}$ 的均值,95%CI和标准差 (公式见开头)

将XY转化为矩阵,构造design matrix

Y<-matrix(y)
X <- matrix(c(rep(1,5), x1, x2), nrow=5,byrow=F)
Xt <- t(X)

求H matrix;I,J辅助阵

矩阵相乘使用%*%
H <- X%*%solve(Xt%*%X)%*%Xt
J_n<-matrix(rep(1,5*5),nrow=5)
I_n <-diag(5)

Sum square

SStot <- t(Y)%*%(I_n - J_n/5)%*%Y #J/5 对应5个数据
SSreg <- t(Y)%*%(H - J_n/5)%*%Y
SSerr <- t(Y)%*%(I_n - H)%*%Y
print("SST    SSR    SSE");print(paste(SStot,SSreg,SSerr,sep="    "))

print("MSE");print(MSE<- as.numeric(SSerr/(5-3)))
varhate<- MSE*(I_n-H)

image.png

求 $\hat β$ , COV(β), SE of β

betahat <- solve((Xt%*%X))%*%Xt%*%Y 
print("betahat");print(betahat)
print("Cov of beta");print(cov_betas <- MSE*solve(t(X)%*%X))
print("se of beta");print(se_betas <- sqrt((diag(diag(cov_betas)))))

image.png

求 $\hat Y$

yhat<-X%*%betahat
    #or yhat<-H%*%Y
print("Yhat");print(yhat)

F statistic (总体)

print(MSreg <- SSreg/2) #2个X
Fstat <- MSreg/MSE
print(Fstat)
pvalF <- 1-pf(Fstat, 1, 5-3)
print(paste(Fstat,pvalF))

image.png

2.ANOVA 使用R计算详解

2.ANOVA 使用R计算详解

首先先明确几个概念(例子以SLR说明,MLR同理)

1. Y=β₀+β₁X₁+β₂X₂+...β_pX_p+e

2. β₀,β₁与 $\hat{Y}$

3. $\hat{Y}_{pred|x_0}$

使用随机生成的5个数据作为示例

结果解读:

分步计算过程

将XY转化为矩阵,构造design matrix

求H matrix;I,J辅助阵

Sum square

求 $\hat β$ , COV(β), SE of β

求 $\hat Y$

F statistic (总体)

推荐阅读更多精彩内容

2.ANOVA 使用R计算详解

首先先明确几个概念(例子以SLR说明,MLR同理)

1. Y=β0+β1X1+β2X2+...βpXp+e

2. β0,β1与

3.

使用随机生成的5个数据作为示例

结果解读:

分步计算过程

将XY转化为矩阵,构造design matrix

求H matrix;I,J辅助阵

Sum square

求, COV(β), SE of β

求

F statistic (总体)

推荐阅读更多精彩内容

1. Y=β₀+β₁X₁+β₂X₂+...β_pX_p+e

2. β₀,β₁与 $\hat{Y}$

3. $\hat{Y}_{pred|x_0}$

求 $\hat β$ , COV(β), SE of β

求 $\hat Y$