StatQuest学习笔记05——线性模型 前言——主要内容

StatQuest学习笔记05——线性模型

前言——主要内容

篇笔记是StatQuest视频教程的第16到22。其中第16内容是线性回归（Linear Regression，它也叫一般线性模型 即General Linear Models）的基本原理；第17内容是在R语言中进行线性回归；第18内容是多重线性回归；第19内容是R主语言与多重线性回归；第20内容是一般线性模型中的t检验与ANOVA；第21内容是设计矩阵；第22内容是在R语言中设计矩阵。

线性回归主要思想

线性回归的主要思想就是利用最小二乘法找到最佳拟合曲线，计算，计算p值，如下所示：

在第4篇笔记中，我们利用了最小二乘法计算了小鼠的体重与大小的线性关系，具体的计算过程可以参见上一篇笔记，最终的线性回归方程如下所示：

由于这个方程的斜率不为0，因此我们知道小鼠的体重的话，就能够预测小鼠的大小。此时我们要看一下这个回归方程的好坏程度，也就是看一下的意义，第一步，计算小鼠的平均大小，如下所示：

将这些小鼠大小的点投射到y轴上，如下所示：

下图中的黑色横线就是小鼠大小的平均值，如下所示：

接着计算离均差平方和（SS），这里需要与残差平方和进行区分。

离均差平方和是指一组数据中，每个数据与均值差值的平方，英文就是SS，全称为sum of squares of deviation from mean ），而残差平方和（SSE）是指观察值与拟合值之间的差异的平方，图计算的是SS：

SS的公式，而这组数据的方差则是指SS除以样本数，如下所示：

此时再回到原始数据，我们将通过最小二乘得到的曲线与每个原始值差值的平方加起来，这个是残差平方和SS(fit)（这个其实就是上篇笔记中提到的SSE），残差平方和除以样本数，就是残差的平均方差（即Var(fit)=SS(fit)/n，这个值没有找到相应的术语），它的公式如下所示：

通常来说，变异等于总的离均差平方和除以样本数，也就是均方差，如下所示：

下图的左图是小鼠的大小的变异（variance），右图是原始数据与回归曲线之间的变异：

如果将上述左图的点与自身均值差值的平方和与右图中原始数据的点和拟合的曲线的距离相比，右边这种通过最小二乘法得到的曲线的变异更小。因此我们在右图中可以知道，如果考虑了小鼠的体重，那么这条曲线可以解释小鼠大小的变异，也就是说，小鼠体重越大，它的大小就越大，反之，也是如此，如下所示：

指的就是小鼠的体重可以解释小鼠大小的程度，它的公式如下所示：

经计算，发现原始数据的均言状为11.1，拟合曲线的均方差为4.4，那么的结果如下所示：

从上图我们可以说，当考虑了小鼠的体重时，小鼠的大小的变异程度会降低60%（这是与原始数据只考虑小鼠的大小，不考虑小鼠的体重相比），也可以说，小鼠大小的变异可以解释小鼠体重变异的60%。

我们也可以采用平方和来计算，如下所示：

我们再看一个案例，如下所示：

在这个案例中，拟合曲线的变异是0，因此小鼠的体重可以解释小鼠大小100%的变异。

再看一个案例，在这个案例中，我们知道小鼠的体重，但是无法根据小鼠的体重来预测小鼠的大小，因为原始数据的变异为11.1，拟合曲线的变异也是11.1，那么就等于0，如下所示：

当我们计算离均差平方和（SS(mean)）时，我们会将忽略小鼠的体重，将每个数据点投射到y轴（也就是小鼠的大小上），如下所示：

此时我们绘制一条小鼠大小的平均值曲线，如下图所示：

此时计算这些数据点距离这个平均值曲线的距离，也就是SS(mean）)，最终我们会得到这个回归曲线，如下所示：

以上就是最小二乘法求线性回归的基本思想，它可以应用到更复杂的公式中去，它的思想就是，先计算原始数据的离均差平方和SS(mean)，再计算这些数据与拟合曲线的的差值平方和，如下所示：

复杂案例

此时我们再看一个稍微复杂的案例。

还是以小鼠为例说明，如果我们知道小鼠的体重（weight），尾巴长度（tail length）和小鼠的身体长度（body length），我们就可以根据前两个变量来预测第三个变量，如下所示：

此时，我们将这3个变量投射到一个三维坐标系中，如下所示：

例如我们有一个数据，它在坐标中的位置就如下所示：

现在把所有的数据点都放到这个坐标中，则如下所示：

此时，我们利用最小二乘法来计算拟合曲线，由于这是三维坐标，这个曲线其实是不是曲线（在二维坐标中是曲线，在三维坐标中就是一个平面了），而是一个平面，如下所示：

这个平面的方程就如下所示：

其中，y是身体长度，x是体重，z是尾巴长度，而0.1则是在y轴上的截矩。如果我们知道一个小鼠的体重和尾巴长度，那么我们就能计算出小鼠的身体长度，如下所示：

此时，我们按照二维坐标的计算方式来计算，如下所示：

此时，这个z轴数据是没用的，它并不会用于计算SS(fit)，此时可以将z值视为0并不会影响预测小鼠的身体长度，然后使用最小二乘法，如下所示（这一步的推导过程我不理解）：

在上图中，把z视为0，意味着，与更少的参数相比，更多的参数并不会让SS(fit)更差（字面翻译，具体原理还没弄懂）。换句话说，像下面的这个公式那样：

虽然Mouse size这个等号后面有很多参数，但是这个公式与下面Mouse size = 0.3 + mouse weight是差不多的，剔除的那几个参数并不会影响Mouse size，因此可以把那些不必要的参数系数设为0，如下所示：

由于存在一些随机因素，小老鼠在这个数据集中得到硬币朝上（get heads）概率并不比大老鼠的高，如下所示：

因此在这个公式中，如果有更多的参数，那么我们就有了更多的随机事情来降低SS(fit)，导致产生一个更好的，因此人们会使用一个调整的（adjusted ，从本质上来讲，就是参数的数据来对进行校正。

虽然是一个很重要的指标，但是它也有一些局限，看下面的案例：

在上图中，只有2个点，计算的SS(mean)=10，SS(fit)=0，后者等于0是因为你总能画出经过两个点的直线，并求出其方程，因此求出的R的平方就是100%，如下所示：

虽然R的平方是100%，这是一个很好的数字，但是任意两个数字都能得到同样的R的平方（其实在利用Qubit测文库的时候，就只是使用了2个标准品，标准品虽然含量非常准，但也是一个非常粗糙的估计，因此在对文库进行定量时，还要配合使用qPCR），下面的这两个数字得到的R的平方也是100%，如下所示：

为了解决这个问题，我们需要研究一下R的平方是否有统计学意义，如下所示：

在计算p值之前，我们先看一下是如何定义的，如下所示：

在小鼠的体重与大小的数据中，我们得到的是0.6，如下所示：

在计算方面，我们会使用一个统计量，即F，F的公式如下所示：

的计算公式和F统计量的计算公式的分子是一样的，如下所示：

这个分子的意义就是小鼠体重的变异，如下所示：

F的分母有点不太一样，它的意义如下所示：

在上图中，下图中的蓝色虚线就是拟合曲线与实际值的差值，它的平方和就变异程度，也就是无法被体重解释的小鼠大小的程度，现在再看一下和F统计量的公式，如下所示：

这两个公式的有一部分是一样的，也就是那些平方和，如下所示：

转换到图形上，则是这个样子：

F统计量分母中的一部分是拟合后的平方和，如下所示：

F后面的那一坨东西是自由度，如下所示：

我们此时看一下什么是自由度，就先从红色方框中的部分开始，如下所示：

拟合曲线中有2个参数，即y=y-intercept + slope x中有intercept和slope这两个参数，它就有2个自由度，而均值的公式是y=y-intercept，它只有一个参数，因此它就是1个自由度。换话句话说，与均值的公式相比，拟合曲线多了一个自由度，即斜率，斜率代表了体重与大小的关系，如下所示：

那么原来红色方框中的那一部分的值就是1，如下所示：

那么F统计量的分子就是由这个多出来的参数解释的变异程度，在这个案例中，就是由小鼠的体重解释小鼠大小的变异程度，如下所示：

如果我们使用小鼠的体重，尾巴长度来解释小鼠大小的变异，就是这个样子：

它的拟合曲线的自由度就是3，如下所示：

此时我们再看一下F统计量的分子，这个分母就是由拟合曲线无法解释的程度，这个分母中的有一个参数n，随后我们会解释，如下所示：

直观上讲，如果我们的公式中有更多的的参数，那么我们就需要多的数据来估计它们，例如我们只需要2个点就能绘制出一条直线，但是我们如果要绘制一个平面，那么至少需要3个点。再接着看F统计量，如果拟合曲线比较好，那么F的值就是一个比较大的数字（从公式上也能看出来，拟合比较好的话，肯定是我们解释的程度要远远大于不能解释的程度），如下图所示：

此时，我们需要把F值转换为p值，先看一个小案例，下面的是一组随机数据，它们计算出来的F值是2，把F值绘制到一个直方图上，如下所示：

再来看一组数据，计算其SS(mean)和SS(fit)，计算F值，再添加到直方图上，如下所示：

再来计算一组随机数据，方法同上，得到F值为1，如下所示：

不断地重复这个过程，我们会得到很多F值，然后把它们都绘制到直方图上，例如其中一个F值是6，如下所示：

从上面的这个图中我们可以知道这些信息：

1. x轴上对应的是F值，从左到右依次增大；

2. F值越大，这个F值的点占到整个直方图中的比例也就越低，换句话讲，就是F值更大，它就是一个越不容易得到的值，如果把得到F值的这个概率定为p，那么，F值越大，p值就越小（也就是说F值越大，它越是一个极端的值，不太容易出现）。

此时，把上面的直方图用一条曲线拟合起来，事实上，我们在计算F值对应的p值时，通常是使用曲线，而不是直方图，如下所示：

此时我们看一个F值的分布，F值的分布曲线是由其自由度决定的，如下所示：

再看由两个F分布叠加的图形，如下所示：

从这张图形上我们可以看出，自由度越大，F分布的曲线下降得越厉害，也就是说，拟合曲线中的参数越多（自由度高），样本数目赵多，p值越小。其实这也是F分布的原理。

R语言与线性回归

此时我们看一下原始数据，如下所示：

weight = 0.9, 1.8, 2.4, 3.5, 3.9, 4.4, 5.1, 5.6, 6.3size = 1.4, 2.6, 1.0, 3.7, 5.5, 3.2, 3.0, 4.9, 6.3

此时绘制出小鼠的体重与其大小的散点图，如下所示：

weight <- c(0.9, 1.8, 2.4, 3.5, 3.9, 4.4, 5.1, 5.6, 6.3) size <- c(1.4, 2.6, 1.0, 3.7, 5.5, 3.2, 3.0, 4.9, 6.3) mouse.data <- data.frame(weight, size) head(mouse.data) plot(mouse.data$weight,mouse.data$size) mouse.regression <- lm(size~weight, data=mouse.data) summary(mouse.regression)

图形如下所示：

线性回归统计结果如下所示：

> summary(mouse.regression) Call: lm(formula = size ~ weight, data = mouse.data) Residuals:    Min      1Q  Median      3Q     Max -1.5482 -0.8037  0.1186  0.6186  1.8852 Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)   0.5813     0.9647   0.603   0.5658  weight        0.7778     0.2334   3.332   0.0126 * --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.19 on 7 degrees of freedom Multiple R-squared:  0.6133,    Adjusted R-squared:  0.558 F-statistic:  11.1 on 1 and 7 DF,  p-value: 0.01256

现在解释一下这个结果：

第一部分：

Call: lm(formula = size ~ weight, data = mouse.data)

这一行表示的是使用lm()函数进行线性回归的公式，这里的size是自变量，weight是因变量，数据集来源于mouse.data。

第二部分：

Residuals:    Min      1Q  Median      3Q     Max -1.5482 -0.8037  0.1186  0.6186  1.8852

这是残差的结果，它表示原始数据距离拟合曲线的距离的 一些统计情况，如果是理想情况，那么正值与负值的频率是一样的，也就是说原始值在拟合曲线上方与下方的点是一样的。

第三部分：

Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)   0.5813     0.9647   0.603   0.5658  weight        0.7778     0.2334   3.332   0.0126 * Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.19 on 7 degrees of freedom

这一部分是经过最小二乘法计算得到出的拟合曲线的系数，其中第1列第1个数字0.5813，就是截矩，截矩下方的0.7778是斜率，那么得到的拟合曲线公式就是：

size = 0.5813 + 0.7778 x weight

第2列与第3列是回归曲线的截矩与斜率的标准差与t值。第4列，也就是最后一列，是这些估计值的p值，通常来说，我们并不关注截矩的这个值。也就是说，它的p值大小我们不乎。但是，斜率，也就是说weight的这个p值必须要小于0.05，也就是说，它必须要有统计学意义，只有这样，我们才能对小鼠的大小进行一个可靠的估计，如下所示：

这个案例的线性回归结果中，weight的p值是小于0.05的，也就是说它有意义，它后面有一个*号，这个*表示p值大于0.01，小于0.5，我们可以在统计结果中看到*的数量对应的p值大小。

最后一行是下面的内容：

Residual standard error: 1.19 on 7 degrees of freedom

它表示，拟合曲线有7个自由度，残差的标准差为1.19。

第四部分：

Multiple R-squared:  0.6133,    Adjusted R-squared:  0.558 F-statistic:  11.1 on 1 and 7 DF,  p-value: 0.01256

在这一部分的结果中，我们可以得到这些信息：

1. Multiple R-squared:  0.6133表示的是，它表示体重可以解释61%的小鼠大小的变异，这已经是很好的结果了；

2. Adjusted R-squared:  0.558这个是由参数的多少调控后的的值；

3. F-statistic:  11.1 on 1 and 7 DF,这表示F的统计量为11.1，1和7是两个自由度；

4. p-value: 0.01256这是F统计量对应的结果。

最后我们可以在这个散点图上添加拟合曲线，输入命令abline(mouse.regression, col="blue")即可，如下所示：

多元线性回归

前面我们讲的线性回归其实只有2个变量（小鼠的体重与小鼠的大小），这种线性回归可以视为简单线性回归，其实线性回归还可以涉及到多个变量，这就是多元线性回归。而很多不太理解线性回归的人通常认为简单线性回归与多元线性回归差异很大，其实它们在本质上没多少区别。下图的左图是简单线性回归，右图是多元线性回归，在计算方面，它们没有区别，如下所示：

在多元线性回归中，需要对进行调整，从而与多元线性回归方程中的多个参数对应，如下所示：

在计算F值时，这两种回归的原理是一样的，只是简单线性回归的p(fit)=2，如下所示：

多元线性回归的p(fit)=3（只是在这个案例中，多元线性回归的p(fit)=3，如果是其他的案例，还要具体情况具体分析），如下所示：

它们的截矩都是一样的，如下所示：

此时我们分别将这两种回归方程中各个数据点与截矩进行对比，来找出这两个回归的优劣：

比较均值可能不太容易看出来，我们可以将这两个回归直接相互进行对比，我们就能知道，我们是否值得花时间和精力收集额外的一个数据，并将其添加到回归方程中，即一个小鼠尾巴（tail length）这个数据，如下所示：

此时，我们按照之前计算F的方法来计算一下这两个线性回归方和的F值，跟先前计算F值不太相同的地方在于，我们使用了简单线性回归的来替换原来的均值，如下所示：

SS(simple)如下所示：

此时，p(simple)=2，如下所示：

SS(multiple)如下所示：

而p(multiple)=3，如下所示：

如果简单线性回归和多元线性回归的差值很大，并且p值很小，那么向简单线性回归中添加一个小鼠尾巴变量（tail length）就是有必要的，如下所示：

R语言与多重线性回归

此时，我们添加上小鼠尾巴的数据，如下所示：

tail = 0.7, 1.3, 0.7, 2.0, 3.6, 3.0, 2.9, 3.9, 4.0

构建小鼠的三个变量的数据框，如下所示：

weight <- c(0.9, 1.8, 2.4, 3.5, 3.9, 4.4, 5.1, 5.6, 6.3) size <- c(1.4, 2.6, 1.0, 3.7, 5.5, 3.2, 3.0, 4.9, 6.3) tail <- c(0.7, 1.3, 0.7, 2.0, 3.6, 3.0, 2.9, 3.9, 4.0) mouse.data <- data.frame(size,weight,tail) head(mouse.data)

构建好的数据框如下所示：

> head(mouse.data)  size weight tail1  1.4    0.9  0.72  2.6    1.8  1.33  1.0    2.4  0.74  3.7    3.5  2.05  5.5    3.9  3.66  3.2    4.4  3.0

简单线性回归

我们先回顾一下简单线性回归，先绘制小鼠体重与其大小的散点图，如下所示：

plot(mouse.data$weight,mouse.data$size)

然后利用lm()函数来对小鼠的这体重和大小这两个变量做回归分析，再用aline()函数添加回归曲线，代码如下：

simple.regression <- lm(size ~ weight, data=mouse.data) summary(simple.regression) abline(simple.regressin, col="red",lwd = 2)

简单线性回归结果如下所示：

> summary(simple.regression) Call: lm(formula = size ~ weight, data = mouse.data) Residuals:    Min      1Q  Median      3Q     Max -1.5482 -0.8037  0.1186  0.6186  1.8852 Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)   0.5813     0.9647   0.603   0.5658  weight        0.7778     0.2334   3.332   0.0126 * --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.19 on 7 degrees of freedom Multiple R-squared:  0.6133,    Adjusted R-squared:  0.558 F-statistic:  11.1 on 1 and 7 DF,  p-value: 0.01256

回归曲线如下所示：

多元线性回归

绘制多元线性回归的散点图，如下所示：

plot(mouse.data)

图片如下所示：

在绘制多元线性回归的图形时，不用指定x轴与y轴，R会把所有的变量进行排列组合，绘制出一多个图形，现在解释一下这个图形；

第一：size与weight，x轴是weight，y轴是size

第二，看下面的图形，x轴是size，y轴是weigth：

上述的这两个图形是非常相似的，只是它们的坐标轴颠倒了一下：

第三，size与tail，这两个变量位于整张图形的对角上，它们的图形也非常相似，也是xy轴颠倒了一下：

同样类似的还有weight和tail变量：

从图形上我们可以看出来，weight和tail与size相关，如下所示：

weigth和tail也有一定的相关关系，这表明，这两个变量也能够为我们提供一些类似的信息，但我们的模型有可能不太需要这些信息：

此时，我们使用lm()函数进行多元线性回归分析，如下所示：

multiple.regression <- lm(size ~weight + tail, data=mouse.data)

在R中，我们指定各种变量的位置的依据如下所示：

默认情况下，y轴截矩和两个斜率如下所示：

查看一下多元线性回归的结果，如下所示：

> summary(multiple.regression) Call: lm(formula = size ~ weight + tail, data = mouse.data) Residuals:     Min       1Q   Median       3Q      Max -0.99928 -0.38648 -0.06967  0.34454  1.07932 Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)   0.7070     0.6510   1.086   0.3192  weight       -0.3293     0.3933  -0.837   0.4345  tail          1.6470     0.5363   3.071   0.0219 * --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.8017 on 6 degrees of freedom Multiple R-squared:  0.8496,    Adjusted R-squared:  0.7995 F-statistic: 16.95 on 2 and 6 DF,  p-value: 0.003399

上述结果与简单线性回归类似，我们看一下最后一部分：

Multiple R-squared:  0.8496,    Adjusted R-squared:  0.7995 F-statistic: 16.95 on 2 and 6 DF,  p-value: 0.003399

的值为0.8496，调整后的值为0.7995，F值为16.95，p值为0.003399，这个拟合效果是非常好的。

我们再看一下各个自变量的系数，如下所示：

Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)   0.7070     0.6510   1.086   0.3192  weight       -0.3293     0.3933  -0.837   0.4345  tail          1.6470     0.5363   3.071   0.0219 * --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我们看一下weight这个自变量的系数，并将它与简单线性回归中的weight进行对比，也就是如下所示：

在多元线性回归中，weight的系数为-0.3293，Pr(>|t|)为0.4345；而在简单线性回归中，weight的系数为0.7778，Pr(>|t|)为0.0126。

这就说明，与只使用tail这一个变量来预测小鼠的大小相比，我们使用weight和tail这两个变量并不是非常显著。

我们再看一下tail这个变量，如下所示：

tail这个变量的p值为0.0219。这就说明，我们使用weight和tail这两个变量比只使用weight这个变量来预测小鼠的大小（size），有显著的区别。

结论就是，我们使用weight和tail这两个变量来预测小鼠的size是比较好的：

但是，如果我们精力与时间有限，那么我们可以中使用小鼠的tail来预测其size，不用考虑weight变量：

t检验与ANOVA

t检验与一般线性模型

统计学初学者通常会首选接触到t检验、方差分析、线性回归等方法，不少人的感觉就是，t检验用于两组均值比较，方差分析用于多组均值比较，而线性回归则用于自变量对因变量的影响分析。看起来似乎没有什么关系，但它们却统一在一个模型下，这就是一般线性模型（General Linear Model）。一般线性模型并不是一个具体的模型，而是多种方法的统称，像t检验、方差分析、线性回归等都从属于一般线性模型的范畴。

——《白话统计》（冯国双.P22.2018）

回顾线性回归

线性回归就是对两个变量采用最小二乘法，求出其线性回归方程，从而使用一个变量预测另外一个变量，在前面的案例中，我们就是检测了小鼠的体重和小鼠的大小，从而求出线性回归方向，从线性回归方程中，我们知道这两个特点：

1. 小鼠的体重对于小鼠的大小预测的效果怎么样（这个要看，这个值越大，表示预测效果越好）；

2. 我们得到的这个线性回归方程是否有统计学意义（这个要看p值，p值小于某个阈值（通常是0.05）就表示有意义）。

如下所示：

将线性回归方程应用于t检验

此时，我们研究一下能否将线性回归方程的思想应用到t检验上（其实是能的），看下面的案例，在这个案例中，我们有一个对照组，有一个突变组，研究了它们某个基因的表达情况，如下所示：

t检验的目标就在于比较两且数据的均值，并观察它们是否有统计学上的差异，如下所示：

如果有一个相同的方法可以计算线性回归和t检验的p值，那么我们就可以在更加复杂的情况下，使用这种方法：

下图分别是线性回归和t检验的散点图：

此时， 我们开始计算F值（由F值就能计算出p值）。

第一步：先忽略x轴坐标的意义，分别计算线性回归图形中所有点的均值和t检验中所有点的均值（即把control和mutant这两组数据混在一块，计算出一个均值），如下所示：

第二步：计算残差平方和，即原始数据点的对应的y轴坐标减去均值，平方后的和，如下所示：

这个计算过程比较简单。

第三步：找出拟合曲线，此时就需要x轴的坐标了：

其中线性回归这个图形的拟合曲线很好找，但是t检验这一组的拟合曲线不太好找，如下所示：

此时，我们先针对control组找到一个拟合曲线，其实就是它的均值所在位置：

得到的曲线为y=2.2，如下所示：

对于mutant组同样如此，它的回归曲线方程为y=3.6，如下所示：

此时，我们就得到了线性回归这一组数据的拟合曲线，以及t检验这一组中的两条回归曲线，如下所示：

此时，我们将t检验这一组图形中的两条曲线合并为一个方程，接着计算F值：

虽然t检验组中的这两条曲线看起来比较诡异，但我们还是能接受的：

将这两条曲线的方程合并为一条，那么这个方程就如下所示：

其中，红色部分表示的是control组的均值，如下所示：

绿色的3.6是mutant组的均值：

残差为蓝色的部分，但这条曲线也是比较诡异的，我们计算control组时，它的mutant这组的系数为0，其实就相当于y=1 x 2.2 + residual了，如下所示：

那么control组中的这几个点的计算过程如下所示：

对于mutant这一组，当我们计算它时，它的control的系数为0，那么这个公式就变成了y = 1 x 3.6 + residual了，如下所示：

计算mutant组的数据时，就是0乘以control，加上1乘以mutant，如下所示：

还有：

现在得到了所有的数据点：

现在我们把这个计算过程中的系数提取出来，看一下：

这个系数看起来给人的印象就是这样的：在分别计算两组的数据时，它们相应的系数（一个是2.2，一个是3.6）就是像是有一对开关在控制，在control组中，2.2的系数为1，3.6的系数为0，mutant正好相反。因此，提取它们，构成一个设计矩阵，如下所示：

这个矩阵我们称为设计矩阵（design matrix），它代表了原始的数据，如下所示：

其中column1表示control组的开或关（on or off）：

column2表示了mutant组的开或关：

在实际运用中，每个column都已经是默认存在了，那么这个公式就会写为如下形式：

此时，我们就拟合了control和mutant数据，把它们放到了一个公式里（其实是一个矩阵），现在计算F值和p值了，如下所示：

第四步：计算SS(fit)，即实际值与拟合曲线差值的平方和：

计算结果如下所示：

F值的计算结果如下所示：

在左图的线性回归图形中，ｙ轴表示的是小鼠的平均大小，ｐ(mean)=1，如下所示：

在右图的t检验中，y轴表示的是平均基因的表达水平，p(mean)=1，如下所示：

同理，线性回归中的p(fit)=2，如下所示：

t检验中的p(fit)=2，如下所示：

现在计算p值，先看一下过程，先计算一SS(mean)，然后计算SS(fit)，如下所示：

t检验的p值和F值的计算至此结束，需要注意的是，我们在前面构建的矩阵并不是t检验的一个标准矩阵（这里先不管这个），接着看ANOVA。

将线性回归方程应用于ANOVA

为了方便理解，我们可以将ANOVA看成升级版的t检验，t检验用于比较两组的数据，ANOVA用于比较多组之间的数据，如下所示：

第一步：计算SS(mean)：

第二步：计算SS(fit)：

第三步：设计矩阵：

第四步：计算F值

设计矩阵

矩阵的查看

在后文，一般线性模型（General Linear Models）都缩写为GLM。在前面的部分中，我们提到了，我们设计的那个矩阵并不是一个标准的矩阵（下图左图），右图才是一个标准的t检验矩阵，如下所示：

这个标准的矩阵对应的这个方程有点不太一样，它的后面是difference(mutant - control)，我们看一下这个新的矩阵：

在这个矩阵中，第1列所有的数据（control和mutant）都变成了（mean(control)），如下所示：

不过，第2列中，只有mutant的数据变成了difference(mutant - control)，如下所示：

例如，下图的数字1对应的数据点（我暂把它叫做A点，它是control组中的点），就是mean(control)的“开”，如下所示：

同理，数字0表示A点对应的就是difference(mutant - control)的“关”，如下所示：

这样我们就很好理解了，control组中的A这个点的对应的mean(control)这个系数为1，A点对应的difference(mutant-control)这个系数为0，加起来，就是A点的这个数字。同时，在mutant中的某个点（我把它称为B点，它位于mutant组），它对应的mean(control)这个系数为1（这个好理解，因为它后面还有mutant - control），对应的difference(mutant-control)这个系数也为1，如下所示：

以下两个方程（设计矩阵和标准矩阵）的残差是相同的，它们也有相同的参数，因此p(fit)也是相同的：

因此，它们的F值也是相同的：

现在我们面临一个疑问：即使这两个方程的结果都一样，为什么右边的这个方程更普遍？原作者认为这可能与回归有关：

我们先看一下这个设计矩阵是如何工作的，如下所示：

设计矩阵的第1列的这个系数是与mean(control)相乘的，如下所示：

第2列的数字是与difference(mutant - control)相乘的，如下所示：

将mean(control)乘以1，那就是它的原始数字，如下所示：

将difference(mutant - control)乘以0，就是没有，如下所示：

这个设计矩阵中的数字是1和1，因此它非常适合做t检验或ANOVA，因为这些数字可以代表不同的分类（其实我们也可以使用其它的数字）：

例如，下面是一个线性回归的设计矩阵，这个矩阵就与左边的这个线性回归方程是等价的：

在这个矩阵的第1列，都是数字1，如下所示：

第2列是这些点在对应的x轴上的坐标，如下所示：

我们现在看一下第1行，它对应的是第1个点，如下所示：

此时用1乘以这个这个公式中的系数（y-intercept），这就表示，这个系数（y-intercept）处于打开（on）的状态，如下所示：

再用第1行的第2个数字（也就是0.9）乘以斜率（slope），如下所示：

为了这让个过程更加具体，我们把这个方程中的矩阵和斜率用具体的数字替换一下，计算出拟合曲线上的真实数据，如下所示：

拟合的曲线斜率非常小，它为0.01（注意：这条拟合曲线与前面我们做的简单线性回归的曲线不同，因为这里面的数据跟前面的貌似不是一批的），如下所示：

它的斜率为0.8，与第1个数据点相乘，如下所示：

经计算，得到0.73，它对应的是就是拟合曲线上的第一个点，就是绿叉的那个点，如下所示：

现在我们看第2行，它对应的是第2个点，如下所示：

现在我们计算拟合曲线上对应的第2个点，下图绿叉所在位置，如下所示：

现在把这个矩阵中所有的点都在拟合曲线上计算出来，如下所示：

一旦计算出这些拟合曲线上的点，那么我们就可以计算其残差，并计算p值，如下所示：

用这个案例就是为了说明，设计矩阵不一定只含有１和０，也可以含有一组其它的数字，将这数字投射到方程中，一次投射一行：

我们再看t检验的案例，t检验的这种案例更加常见，因此我们后文中的案例都按照这种格式进行，如下的所示：

现在我们知道了，我们可以把任何数字都放到一个设计矩阵里。

第1个案例

现在我们再看一些案例，这些案例混合了t检验与线性回归，下图是两种不同种类的小鼠的体重与大小，不同颜色用于区分不同种类的小鼠，如下所示：

其中，红色是对照组小鼠，如下所示：

绿色的是突变组小鼠（mutant）如下所示：

从直观上我们就能够看出，突变组小鼠的比较大，如下所示：

换句话说，突变组与对照组小鼠趋势是这个样子的：

我们现在面临的问题就是：如何使用统计学的方法来判断这两组小鼠是否有统计学意义上的差别？

如果我们对整体小鼠（包括mutant和control）做一个回归，那么我们也无法找出这两组小鼠是否有差异，如下所示：

如果我们忽略小鼠的体重与大小之间的关系，只对小鼠的大小做一个统计，我们就发现，小鼠的大小（size）并没有统计学意义上的差别，如下所示：

由于小鼠的体重（weight）和大小（size）之间的这个关系与小鼠的类型有关，因此我们可以把这两种因素归到一种统计学方法中来，如下所示：

换句话讲，我们并不直接比较小鼠的大小（这个比较是用t检验来进行的），而是用另外的一种方法来比较，如下所示：

这种方法就是，比较这两种小鼠的回归曲线，如下所示：

为了比较这两个曲线，我们需要一个含有对照组小鼠方程截矩的方程，如下所示：

这个方程还要含有突变组小鼠与对照组小鼠的差值，如下所示：

这个方程还要含有斜率（在这个案例中，这两条曲线的斜率是相同的，但在实际运用中，可能不同），如下所示：

这就意味着，我们需要一个设计矩阵，这个矩阵的第1列是数字1，它们是y轴上的某些点，如下所示：

设计矩阵的第2列是0和1，它们代表mutant off（也就是突变组与对照组的差值），如下所示：

其中0代表off，它代表control的off，如下所示：

1代表mutant小鼠，如下所示：

第3列是小鼠的体重，如下所示：

绿色的是突变组，如下所示：

现在看第1行，1，0，2.4对应的方程的数字如下所示：

它投射到拟合曲线上的点如下所示：

第2行的数字对应到方程上如下所示：

第2个点投射到拟合曲线上就是这个样子：

把所有的点投射到相应的曲线上，就如下所示：

一旦我们确定了拟合曲线上的点的位置，就可以计算其残差，如下所示：

第1个simpler model

此时我们把这个复杂的模型（用fancy model表示）与简单的模型（simpler model）比较，如下所示：

此时，把这个fancy model的残差平方和代入到下面的公式中，如下所示：

其中fancy model有3个驎肆，也就是p(fancy)为3，如下所示：

把simple的SS代入到下面方程，如下所示：

simple方程中的参数只有1个，也就是说p(simple)为1，代入到下面的方程，如下所示：

最终得到的F值为21.88，p值为0.003，如下所示：

这里需要提一下的是，simpler model可以是任何的模型，下图中这个simpler model就是一个简单线性回归模型，如下所示：

第2个simpler model

我们用这个simple model来计算一下，如下所示：

得到的结果为F值是32.6，p值为0.0023，如下所示：

从上图的结果来看，使用小鼠的体重(weight)和种类(type)能够比只使用小鼠的体重(weigth)更好地预测小鼠的大小(size)。再看下面的一个simpler model，在这个模型中，我们忽略了小鼠的体重，如下所示：

第3个simpler model

现在把这个simpler model的平方和加入到F值的公式中，如下所示：

这次得到的F值与p值如下所示：

这个结果的p值与第一个simpler model（即考虑全部的小鼠平均值的那个模型）的p值相比，它更小，这就说明，使用小鼠的体重(weight)和类型(type)来预测小鼠的大小，比只使用类型(type)更好。

第2个案例：批次效应

现在A实验室（Lab A）做了一个实验，如下图左图所示，而B实验室（Lab B）实验重复了A实验室的实验。现在我们要综合这两次的实验结果，观察一下，mutant组与control相比，有没有差异，但是，这两次实验不是一次做的，因此存在批次效果（batch effect），我们现在要解决的问题是，如何消除这个批次效应，如下所示：

首选，我们先计算一下Ａ实验的control组的均值，如下所示：

第二，我们计算一下B实验的control这个结果与A实验的control这个均值的差值，这用于研究批次效应，如下所示：

第三，我们再加入mutant和control数据的差值，如下所示：

此时，设计一个矩阵，如下所示：

从本质上来讲， 我们想知道上述方程最后一项重要不重要。也就是说矩阵最后一项重要与否，如下所示：

现在我们比较这个复杂方程（simple equation）与简单方程（simpler equation，简单方程忽略了control与mutant的分组）的区别，如下所示：

如果得到一个很小的p值，那么这就告诉我们要研究control与mutant的差异，这就意味着我control和mutant之间存在着显著性差异。

R语言与设计矩阵

第1案例：小鼠体重

还看前文的案例，我们想知道mutant组的小鼠和control的小鼠在size上是否有统计学意义上的差别。我们分别做了这两种小鼠的回归曲线，如下所示：

此时我们使用R语言来进行分析，输入以下代码：

Type <- factor(c(  rep("Control", times=4),  rep("Mutant", times=4))) Weight <- c(2.4, 3.5, 4.4, 4.9, 1.7, 2.8, 3.2, 3.9) Size <- c(1.9, 3, 2.9, 3.7, 2.8, 3.3, 3.9, 4.8) model.matrix(~Type+Weight)

结果如下所示：

> Type <- factor(c( +   rep("Control", times=4), +   rep("Mutant", times=4))) > > Weight <- c(2.4, 3.5, 4.4, 4.9, 1.7, 2.8, 3.2, 3.9) > Size <- c(1.9, 3, 2.9, 3.7, 2.8, 3.3, 3.9, 4.8) > model.matrix(~Type+Weight)  (Intercept) TypeMutant Weight1           1          0    2.42           1          0    3.53           1          0    4.44           1          0    4.95           1          1    1.76           1          1    2.87           1          1    3.28           1          1    3.9attr(,"assign") [1] 0 1 2attr(,"contrasts") attr(,"contrasts")$Type [1] "contr.treatment"

其中model.matrix()这个函数是用于构建设计矩阵的，它能自动将因子型数据转换为哑变量（dumpy variables）。构建好的矩阵的所代表的第1列是对照组的截矩，如下所示：

第2列代表的数据是mutant与control的差值，如下所示：

最后1列是如下所示：

现在构建线性回归模型，结果如下所示：

> model <- lm(Size~Type + Weight) > summary(model) Call: lm(formula = Size ~ Type + Weight) Residuals:       1        2        3        4        5        6        7        8 0.05455  0.34562 -0.41623  0.01607 -0.01753 -0.32646 -0.02062  0.36461 Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)  0.08052    0.52744   0.153  0.88463   TypeMutant   1.48685    0.26023   5.714  0.00230 ** Weight       0.73539    0.13194   5.574  0.00256 ** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.3275 on 5 degrees of freedom Multiple R-squared:  0.8975,    Adjusted R-squared:  0.8564 F-statistic: 21.88 on 2 and 5 DF,  p-value: 0.003367

现在看这个结果，第1部分如下所示：

Call: lm(formula = Size ~ Type + Weight) Residuals:       1        2        3        4        5        6        7        8 0.05455  0.34562 -0.41623  0.01607 -0.01753 -0.32646 -0.02062  0.36461

这是每个数据与拟合曲线的残差。

第2部分，我们看最后的信息：

Multiple R-squared:  0.8975,    Adjusted R-squared:  0.8564

这些是多重，因为我们的拟合曲线是一个比复杂的方程（它有2个参数），后面的是根据参数的数量调整后的。最后就是F值与p值了，p值为0.003367，这个p值是线性回归的p值整体数据（其实就是simpler model）对比后的p值，如下所示：

接着看各项的系数，如下所示：

Coefficients:            Estimate Std. Error t value Pr(>|t|)   (Intercept)  0.08052    0.52744   0.153  0.88463   TypeMutant   1.48685    0.26023   5.714  0.00230 ** Weight       0.73539    0.13194   5.574  0.00256 ** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结果中Weight这个变量的p值为0.00256，它表示的是，如果在这个复杂模型（fancy model）中剔除Weight这个变量，构建成的简单模型（simpler model）与原始数据的拟合效果。它反映的是一个复杂模型的统计效果和只用常规t检验的统计效果的比较（p值小于0.05，说明前者的统计效果更好）。换句话讲，它反映了复杂模型中剔除了Weight这个变量后，使用最小二乘法进行的检验与原复杂模型检验的比较，由于p值很小，这就说明了Weight这个变量的重要性，如下所示：

同理，TypeMutant这个变量对应的p值为0.00230，它反映的是原复杂方程（fancy model）与该方程剔除小鼠类型这个变量后的比较（剔除了小鼠类型的这个变量，就变成了普通线性回归方程），如下所示：

在图形上反映的效果就是如下所示：

第2案例：批次效果

这个案例还是前面的那个案例，在这个案例中，我们会看到两个实验室做的同一个实验，它们存在着批次效应，如下所示：

现在使用R语言进行分析，如下所示：

Lab <- factor(c(rep("A",times=6),rep("B",times=6))) Type <- factor(c(rep("Control", times = 3),rep("Mutant", times = 3),rep("Control", times = 3),rep("Mutant", times = 3))) Expression <- c(1.7, 2, 2.2,3.1, 3.6, 3.9, 0.9, 1.2, 1.9, 1.8, 2.2, 2.9) model.matrix(~Lab+Type)

构建好的设计矩阵如下所示：

> model.matrix(~Lab+Type)   (Intercept) LabB TypeMutant1            1    0          02            1    0          03            1    0          04            1    0          15            1    0          16            1    0          17            1    1          08            1    1          09            1    1          010           1    1          111           1    1          112           1    1          1attr(,"assign") [1] 0 1 2attr(,"contrasts") attr(,"contrasts")$Lab [1] "contr.treatment"attr(,"contrasts")$Type [1] "contr.treatment"

第1列代表的是Lab A control mean，如下所示：

第2列是lab B offset，如下所示：

第3列是difference(mutant - control)，如下所示：

现在计算一下批次效应，如下所示：

batch.lm <- lm(Expression ~ Lab + Type) summary(batch.lm)

结果如下所示：

> summary(batch.lm) Call: lm(formula = Expression ~ Lab + Type) Residuals:    Min      1Q  Median      3Q     Max -0.6500 -0.2833 -0.0500  0.2750  0.7167 Coefficients:            Estimate Std. Error t value Pr(>|t|)     (Intercept)   2.1167     0.2279   9.287  6.6e-06 *** LabB         -0.9333     0.2632  -3.546 0.006250 ** TypeMutant    1.2667     0.2632   4.813 0.000956 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.4558 on 9 degrees of freedom Multiple R-squared:  0.7989,    Adjusted R-squared:  0.7542 F-statistic: 17.87 on 2 and 9 DF,  p-value: 0.0007342

这个结果跟先前多元回归的结果类似，最后一部分就是以及调整后的，它的p值为0.0007342，但是这个p值不是我们所需要的，这个p值代表的是复杂方程与一个简单方程（也就是下图中y = overall mean）的统计学差异，这个p值是不能用于判断mutant和control差异的，如下所示：

我们再看一下TypeMutant这个变量对应的p值，这个p值是0.000956，这个才是我们所需要的p值，它表示的是，我们去掉difference(mutant - control)这个变量后，我们是否还能拟合出一条好的曲线（我个人理解就是：原假设就是去掉`TypeMutant后，剩下的数据能够拟合出一条好的曲线方程，但是结果却是p值小于0.05，是一个非常小的值，因此原假设不成立，也就是说这个变量不能去掉），如下所示：

总之，这个变量（TypeMutant）是我们所需要的，如下所示：

这个案例说明的是，我们在做线性回归时，不一定总要看最终的p值（也就是统计结果右下角的p值），而是针对自己具体的容易来查找相应的p值。