DALS012-Linear Models线性模型02

title: DALS012-Linear Models线性模型02
date: 2019-08-17 12:0:00
type: "tags"
tags:

线性代数
categories:
生物统计

前言

这一部分是《Data Analysis for the life sciences》的第5章线性模型的第2小节。这一部分的主要内容涉及矩阵的数学计算原理。

先来看一下前面的一个案例。

小鼠饲料案例

在前面的内容中，我们使用了t检验来研究高脂饲料(hf)饲喂的小鼠与普通饲料(chow)饲喂的小鼠体重的差异，现在我们使用线性模型来研究一下这两种小鼠的体重是否有差异，最终我们会发现，这两种统计方法在本质上是一样的，如下所示：

dir <- system.file(package = "dagdata")
list.files(dir)
dir
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
# input raw data
head(dat)
str(dat)
stripchart(dat$Bodyweight ~ dat$Diet, vertical=TRUE, method="jitter",
main="Bodyweight over Diet")

结果如下所示：

> head(dat)
  Diet Bodyweight
1 chow      21.51
2 chow      28.14
3 chow      24.04
4 chow      23.45
5 chow      23.68
6 chow      19.79
> str(dat)
'data.frame':   24 obs. of  2 variables:
 $ Diet      : Factor w/ 2 levels "chow","hf": 1 1 1 1 1 1 1 1 1 1 ...
 $ Bodyweight: num  21.5 28.1 24 23.4 23.7 ...

image

从这张图上我们大致可以看出来，hf组的小鼠体重比chow组的小鼠体重高一些。现在我们使用设计矩阵 $\mathbf{X}$ （公式为~Diet）来计算一下这个结果，在设计矩阵中，第2列是分组信息，也就是饲料的类型，如下所示：

levels(dat$Diet)
X <- model.matrix(~ Diet, data=dat)
X

结果如下所示：

> X
   (Intercept) Diethf
1            1      0
2            1      0
3            1      0
4            1      0
5            1      0
6            1      0
7            1      0
8            1      0
9            1      0
10           1      0
11           1      0
12           1      0
13           1      1
14           1      1
15           1      1
16           1      1
17           1      1
18           1      1
19           1      1
20           1      1
21           1      1
22           1      1
23           1      1
24           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Diet
[1] "contr.treatment"

`lm()`函数背后的数学原理

现在我们先不用lm()来进行计算，先用设计矩阵%\mathbf{X} $来计算一下$ \beta$，这个值能够使线性模型的平方和最小，公式如下所示：
$\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y}$
在R中，可以使用矩阵相乘的符号%*%，以及solve()函数来求解上述方程，如下所示：

Y <- dat$Bodyweight
X <- model.matrix(~ Diet, data=dat)
solve(t(X) %*% X) %*% t(X) %*% Y

结果如下所示：

> Y <- dat$Bodyweight
> X <- model.matrix(~ Diet, data=dat)
> solve(t(X) %*% X) %*% t(X) %*% Y
                 [,1]
(Intercept) 23.813333
Diethf       3.020833

这些系数是对照组的平均值（23.813333），以及两组的差值（3.020833），计算一下，如下所示：

s <- split(dat$Bodyweight, dat$Diet)
mean(s[["chow"]])
mean(s[["hf"]]) - mean(s[["chow"]])

结果如下所示：

> s <- split(dat$Bodyweight, dat$Diet)
> mean(s[["chow"]])
[1] 23.81333
> mean(s[["hf"]]) - mean(s[["chow"]])
[1] 3.020833

现在我们使用lm()函数两周来计算一下，如下所示：

fit <- lm(Bodyweight ~ Diet, data=dat)
summary(fit)
(coefs <- coef(fit))

结果如下所示：

> summary(fit)

Call:
lm(formula = Bodyweight ~ Diet, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1042 -2.4358 -0.4138  2.8335  7.1858 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   23.813      1.039  22.912   <2e-16 ***
Diethf         3.021      1.470   2.055   0.0519 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.6 on 22 degrees of freedom
Multiple R-squared:  0.1611,    Adjusted R-squared:  0.1229 
F-statistic: 4.224 on 1 and 22 DF,  p-value: 0.05192
> (coefs <- coef(fit))
(Intercept)      Diethf 
  23.813333    3.020833

线性回归系数

下图说明了前面我们计算出来的2个系数，即对照组小鼠的平均体重，以及hf组的小鼠体重与对照组的差异，如下所示：

stripchart(dat$Bodyweight ~ dat$Diet, vertical=TRUE, method="jitter",
           main="Bodyweight over Diet", ylim=c(0,40), xlim=c(0,3))
a <- -0.25
lgth <- .1
library(RColorBrewer)
cols <- brewer.pal(3,"Dark2")
abline(h=0)
arrows(1+a,0,1+a,coefs[1],lwd=3,col=cols[1],length=lgth)
abline(h=coefs[1],col=cols[1])
arrows(2+a,coefs[1],2+a,coefs[1]+coefs[2],lwd=3,col=cols[2],length=lgth)
abline(h=coefs[1]+coefs[2],col=cols[2])
legend("right",names(coefs),fill=cols,cex=.75,bg="white")

image

在这里使用前面的小鼠案例主要是为了说明了回归与t检验在本质上是一样的，这个简单的线性回归模型给出了与特定类型t检验相同的结果（t统计量与p值）。这是两组之间的t检验，前提是两组的总体标准差相同。当我们假设误差 $\boldsymbol{\varepsilon}$ 都是均匀分布时，这些误差被编码到线性模型中。虽然在这种情况下的线性模相当于t检验，但我们可以对线性模型进行拓展，将其扩展到复杂的形式。在下面的内容中，我们将会说明线性模型能够得到几乎相同的结果：

> summary(fit)$coefficients
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 23.813333   1.039353 22.911684 7.642256e-17
Diethf       3.020833   1.469867  2.055174 5.192480e-02

t检验的结果与之相同：

ttest <- t.test(s[["hf"]], s[["chow"]], var.equal=TRUE)
summary(fit)$coefficients[2,3]
ttest$statistic

结果如下所示：

> summary(fit)$coefficients[2,3]
[1] 2.055174
> ttest$statistic
       t 
2.055174

练习

P189

标准误

使用矩阵代数来计算最小二乘估计时，所得到的估计值是取机变量。我们还能计算这些值的标准误。线性代数也能满足这些任务，先看一些案例。

自由落体

在学习统计学的过程中，了解随机源于何处是有必要的。在自由落体这个案例中，当我们对落下来的球进行测量时，就会出现检测误差，这就是自由落体运行中随机的来源。假设我们做了好几次实验，每次实验我们都会得到一组测量误差。这就意味着，我们得到的数据基本是随机变化的，这反过来，也意味着，我们计算得到的估计值也是随机变化的。在这个案例中，每次我们进行实验时，引力常数的估计值也都在发生变化。引力常数是一个确定的值，但是我们对它的估计不是。为了说明这一步，我们可以使用Monte Carlo模拟一下。具体来说，我们将会重复地生成数据，并对每次的二次项(quadratic term)进行估计，如下所示：

set.seed(1)
B <- 10000
h0 <- 56.67
v0 <- 0
g <- 9.8 ##meters per second
n <- 25
tt <- seq(0,3.4,len=n) ##time in secs, t is a base function
X <-cbind(1,tt,tt^2)
##create X'X^-1 X'
A <- solve(crossprod(X)) %*% t(X)
betahat<-replicate(B,{
y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
betahats <- A%*%y
return(betahats[3])
})
head(betahat)

结果如下所示：

> head(betahat)
[1] -5.038646 -4.894362 -5.143756 -5.220960 -5.063322 -4.777521

从上面结果我们可以发现，每次的估计值都不一样，这是因为估计值 $\hat{\beta}$ 是一个随机变量，它其实是服从正态分布的，如下所示：

library(rafalib)
mypar(1,2)
hist(betahat)
qqnorm(betahat)
qqline(betahat)

image

上图显示的是，通过Monte Carlo模拟生成的自由落体运行数据对回归系数的估计值的分布。左侧是直方图，右侧是qq图。

由于 $\hat{\beta}$ 是我们生成的那些服从正态分布的数据的线性组合，因此从QQ图上我们也可以看出， $\hat{\beta}$ 也服从正态分布。此外，这个分布的均值的真实参数0.5g，可以通过Monte Carlo模拟来所证实，如下所示：

round(mean(betahat),1)

结果如下所示：

> round(mean(betahat),1)
[1] -4.9

但是，当我们对估计值进行计算时，无法得到精确的数值，这是因为我们估计值的标准误大约是：

> sd(betahat)
[1] 0.2129976

在接下来的案例中，我们将会演示一下，不通过Monte Carlo模拟来计算标准误的方法，因为在实际的分析中，我们无法精确地知道误差是如何产生的。

父子身高

在父母身高的案例中，我们也遇到了随机问题，因为我们是对父子身高的总体进行随机抽样。为了说明这个问题，现在假设我们已经有了这个总体：

library(UsingR)
x <- father.son$fheight
y <- father.son$sheight
n <- length(y)

现在我们使用Monte Carlo模拟来生成一个含有50个数据的样本，如下所示：

N <- 50
B <-1000
betahat <- replicate(B,{
index <- sample(n,N)
sampledat <- father.son[index,]
x <- sampledat$fheight
y <- sampledat$sheight
lm(y~x)$coef
})
betahat <- t(betahat) #have estimates in two columns

现在绘制一个QQ图，我们看到，我们的估计值大致是服从标准正态分布的，如下所示：

mypar(1,2)
qqnorm(betahat[,1])
qqline(betahat[,1])
qqnorm(betahat[,2])
qqline(betahat[,2])

image

现在看一下估计值之间的相关性，如下所示：

> cor(betahat[,1],betahat[,2])
[1] -0.9992293

当我们计算我们估计值的线性组合时，我们需要知道这个信息能够正确地计算出这些线性组合的标准误差。在接下来的部分中，我们会提到方差-协方差(variance-covariance)矩阵。两个随机变量的协方差定义如下：

mean( (betahat[,1]-mean(betahat[,1] ))* (betahat[,2]-mean(betahat[,2])))

结果为：

> mean( (betahat[,1]-mean(betahat[,1] ))* (betahat[,2]-mean(betahat[,2])))
[1] -1.035291

协方差是相关系数乘以每个随机变量的标准差，如下所示：

$\mbox{Corr}(X,Y) = \frac{\mbox{Cov}(X,Y)}{\sigma_X \sigma_Y}$

除此之外，这个量在实际分析中没有一个现实意义上解释。但是，它对于数学推导的过程来说，有着重要意义。在接下来的部分中，我们会描述用于计算线性模型估计的标准差的矩阵代数计算方法。

方差-协方差矩阵

我们先来定义一什么是方差-协方差矩阵(variance-covariance matrix) $\Sigma$ 。地于一个元素是随机变量的向量 $\mathbf{Y}$ 来说，我们定义 $\Sigma$ 是含有 $i,j$ 维的矩阵，如下所示：
$\Sigma_{i,j} \equiv \mbox{Cov}(Y_i, Y_j)$

如果 $i=j$ ，那么协方差就等于方差，如果变量都是独立的，则协方差都为0。在我们现在所学到的各种向量里， $\mathbf{Y}$ 向量的每个观测值 $Y_{i}$ 是从总体中抽样得到的，我们可以假设它的每个元素都是独立的，并且 $Y_{i}$ 有着相同的方差 $\sigma^2$ ，因此方差-协方差矩阵只有两类元素，如下所示：
$\mbox{Cov}(Y_i, Y_i) = \mbox{var}(Y_i) = \sigma^2$
$\mbox{Cov}(Y_i, Y_j) = 0, \mbox{ for } i \neq j$

这就是暗示了 $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}$ ，其中 $\mathbf{I}$ 是单位矩阵(Identity matrix)。

线性组合的方差

线性代数提供的一个有用结果就是 $\mathbf{Y}$ 的 $\mathbf{AY}$ 线性结合的方差-协方差矩阵，它的计算如下所示：
$\mbox{var}(\mathbf{AY}) = \mathbf{A}\mbox{var}(\mathbf{Y}) \mathbf{A}^\top$
例如，如果 $Y_{1}$ 和 $Y_{2}$ 是独立的，并且其方差为 $\sigma^2$ ，那么：

$\mbox{var}\{Y_1+Y_2\} = \mbox{var}\left\{\begin{pmatrix}1&1\end{pmatrix}\begin{pmatrix} Y_1\\Y_2\end{pmatrix}\right\}=\begin{pmatrix}1&1\end{pmatrix} \sigma^2 \mathbf{I}\begin{pmatrix} 1\\1\\ \end{pmatrix}=2\sigma^2$

我们会使用上述的结果来获取LSE的标准误。

LES标准误

$\hat{\beta}$ 是一个线性组合，即 $\mathbf{Y}$ : $\mathbf{AY}$ with $\mathbf{A}=\mathbf{(X^\top X)^{-1}X}^\top$ 的线性组合，因此我们可以使用上面的公式来计算一个我们估计值的方差：

$$
\mbox{var}(\boldsymbol{\hat{\beta}}) = \mbox{var}( \mathbf{(X^\top X)^{-1}X\top Y} ) =

\\mathbf{(X^\top X)^{-1} X^\top} \mbox{var}(Y) (\mathbf{(X^\top X)^{-1} X^\top})\top = \

\mathbf{(X^\top X)^{-1} X^\top} \sigma^2 \mathbf{I} (\mathbf{(X^\top X)^{-1} X^\top})\top = \

\sigma^2 \mathbf{(X^\top X)^{-1} X^\top}\mathbf{X} \mathbf{(X^\top X)^{-1}} = \
\sigma^2\mathbf{(X\top X)^{-1}}
$$

其中，这个矩阵的平方根的对角线含有我们估计值的标准误。

估计 $\sigma^2$

为了从上述公式中获得一个精确的估计值，我们需要估计 $\sigma^2$ 。在前面我们通过抽样估计了标准误。但是 $Y$ 的抽样标准误不是 $\sigma$ ，因为 $Y$ 还包含了变异，这些变量是由模型 $\mathbf{X\beta}$ 的确定部分引入的。此时我们采用的方法是利用残差。

我们可以按下面的公式构建残差：
$\mathbf{r}\equiv\boldsymbol{\hat{\varepsilon}} = \mathbf{Y}-\mathbf{X}\boldsymbol{\hat{\beta}}$
其中， $r$ 和 $\hat{\epsilon}$ 都可以用来表示残差(residuals)。然后我们使用上述的公式来估计残差，这种方式与单变量情况下类似：
$s^2 \equiv \hat{\sigma}^2 = \frac{1}{N-p}\mathbf{r}^\top\mathbf{r} = \frac{1}{N-p}\sum_{i=1}^N r_i^2$

其中， $N$ 是样本数目， $p$ 是 $\mathbf{X}$ 的列数，或者是参数的数目（包括截矩 $\beta_{0}$ ）。公式中还要除以 $N-p$ ，这是因为根据数学理论（不用深究，这个跟自由度能关），除以 $N-p$ ，表示无偏估计。下面在R中计算一下，观察一下我们是否能够获得与Monte Carlo模拟一样的数值：

n <- nrow(father.son)
N <- 50
index <- sample(n,N)
sampledat <- father.son[index,]
x <- sampledat$fheight
y <- sampledat$sheight
X <- model.matrix(~x)
N <- nrow(X)
p <- ncol(X)
XtXinv <- solve(crossprod(X))
resid <- y - X %*% XtXinv %*% crossprod(X,y)
s <- sqrt( sum(resid^2)/(N-p))
ses <- sqrt(diag(XtXinv))*s
summary(lm(y~x))$coef[,2]
ses
apply(betahat,2,sd)

结果如下所示：

> summary(lm(y~x))$coef[,2]
(Intercept)           x 
  8.3899781   0.1240767 
> ses
(Intercept)           x 
  8.3899781   0.1240767 
> apply(betahat,2,sd)
(Intercept)           x 
  8.3817556   0.1237362

从上面我们可以看出来，我们的计算结果与Monte Carlo模拟的结果几乎一样。

估计值的线性组合

多数情况下，我们会计算估计值的一个线性组合，例如 $\hat{\beta_{2}}-\hat{\beta_{1}}$ 。这就是 $\hat{\beta}$ 的一个线性组合，如下所示：
$\hat{\beta}_2 - \hat{\beta}_1 = \begin{pmatrix}0&-1&1&0&\dots&0\end{pmatrix} \begin{pmatrix} \hat{\beta}_0\\ \hat{\beta}_1 \\ \hat{\beta}_2 \\ \vdots\\ \hat{\beta}_p \end{pmatrix}$
通过上述的公式，我们就知道中如何计算 $\hat{\beta}$ 的方差-协方差矩阵。

CLT与t分布

我们已经了解了如何过计算估计值的标准误。但是，我们在第1章中了解到，在计算置信区间时，我们需要知道随机变量的分布。我们努力计算标准误的原因是因为，CLT适用于线性模型。如果 $N$ 足够大，那么LSE就会服从正态分布，其中这个正态分布的均值为 $\beta$ ，标准误就是前面计算的标准误。对于小样本来说，如果 $\varepsilon$ 服从正态分布，那么 $\hat{\beta}-\hat{\beta}$ 服从t分布。在这里，我们不需要知道空上计算过程，但是这个结果很有，因为这是我们在线性模型中计算p值，置信区间的基础。

代码与数学

构建线性模型的标准做法要么是假设 $\mathbf{X}$ 是固定的，要么我们给它们设置条件。因此 $\mathbf{X\beta}$ 没有像 $\mathbf{X}$ 那样固定的方差。这就浊为什么我们要写 $\mbox{var}(Y_i) = \mbox{var}(\varepsilon_i)=\sigma^2$ 。在实际计算中，这有可能造成误解，如下所示：

x = father.son$fheight
beta = c(34,0.5)
var(beta[1]+beta[2]*x)

结果如下所示：

> var(beta[1]+beta[2]*x)
[1] 1.883576

这个值并不等于0，也不接近于0。这就是为什么我们要谨慎区分代码与数学的一个案例。var()函数仅仅是用于计算我们输入数值的方差，而数学对方差的定义则是要考虑到这些数字是否是随机变量。在上述的R代码中，x并没有固定：我们只是让它发生变化，但是，当我们写下 $\mbox{var}(Y_i) = \sigma^2$ 时，我们是把数学上的x强制固定下来。同样，如果我们使用R来计算自由落体运行案例中 $Y$ 的方差，我们也会获得一个与 $\sigma^2=1$ （已知方差）明显不同的数值，如下所示：

n <- length(tt)
y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
var(y)

结果如下所示：

> n <- length(tt)
> y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
> var(y)
[1] 329.5136

练习

P199

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,793评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,567评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,342评论 0赞 338
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,825评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,814评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,680评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,033评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,687评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 42,175评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,668评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,775评论 1赞 332
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,419评论 4赞 321
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,020评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,978评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,206评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,092评论 2赞 351
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,510评论 2赞 343

DALS012-Linear Models线性模型02

前言

小鼠饲料案例

lm()函数背后的数学原理

线性回归系数

练习

标准误

自由落体

父子身高

方差-协方差矩阵

线性组合的方差

LES标准误

估计

估计值的线性组合

CLT与t分布

代码与数学

练习

推荐阅读更多精彩内容

`lm()`函数背后的数学原理

估计 $\sigma^2$