title: DALS011-Linear Models线性模型01MatrixDescription
date: 2019-08-11 12:0:00
type: "tags"
tags:

• 线性代数
  categories:
• Data Analysis for the life sciences

前言

这一部分是《Data Analysis for the life sciences》的第5章线性模型的第1小节。这一部分的主要内容涉及矩阵的表示方法。

背景知识

数据分析中的很多模型都可以用矩阵代数来表示。我们指的这些模型都是线性模型(linear models)。线性(Linear)并不是说它们是直线，而是说是线性组合。使用矩阵代数来描述线性模型非常方便，因为我们可以更好地构建这些模型，并使用数学工具来方便计算。在这一部分中，我们将详细地描述我们如何使用矩阵代数来构建并拟合线性模型。

在本书中，我们关注的线性模型是基于二分法的线性模型：例如治疗组与对照组。在前面提到的小鼠饲料的案例中就是这种线性模型的一个典型案例。在这一部分中，我们会描述一些比较复杂的案例，不过还会继续关注二分法变量的线性模型。

当我们学习线性模型时，需要记住，我们使用的变量仍然是随机变量。这就意味着，我们使用线性模型获得的估计值也是随机变量。虽然数学上这理解起来比较复杂，但是我们前面学到的一些概念在这里仍然适用。现在可以先做一些练习，在线性模型的练习中复习一些随机变量的概念。

练习

P172

设计矩阵(The Design Matrix)

R中的formual()和model.matrix()可以用于生成各种线性模型的设计矩阵(design matrices)（也被称为模型矩阵，即model matrices)。例如，在小鼠的饲料实验中，我们可以这样构建模型，如下所示：

其中，表示体重，而等于1的时候，表示小鼠吃的是高脂饲料(hf)。我们可以使用实验单位(Experimental unit)来表示N个不同的实验动物，这里的实验单位是小鼠。我们称这类变量为指示变量(indicator variables)，因此这类变量仅用于标明实验单位，表示某种干预的有无。现在用矩阵代数来描述就是以下形式：

或为：

简化形式为：

其中是设计矩阵，一旦我们定义了设计矩阵，那么我们就能计算出最小二乘估计。这个过程称为拟合模型(fitting the model)。在R中，我们可以直接在lm()函数中输入一个公式(formula)直接进行拟合。在后面的一个脚本中，我们会使用model.matrix()函数，这个函数是在lm()函数内部使用的。这种使用方式哦可以方便地让我们将R的formula与矩阵连接起来，有助于我们对lm()的结果进行解释。

设计方案

设计矩阵的选择是线性模型中的关键步骤，因为它编码了哪些系数将适合用于建模，以及样本之间的相互关系。一个常见的误解就是，选择实验设计要遵循简单明了的原则，即要在描述中直接标明样本。事实并非如此。关于每个样本（无论是对照组还是处理组，实验批次等）的基础信息里并不意味着这是一个“正确”的设计矩阵。设计矩阵还应该对中的变量是如何解释的值这些信息进行编码，这是研究者们必须要考虑的。

在这一部分的案例中，我们使用线性模型来比较不同组之间的差异。因此，最终进行计算的设计矩阵应该至少包含两列，即截矩(intercept)列，截矩列是第1列，另外就是第2列，第2列指定样本。在这种情况下，线性模型中的两个系数(coefficient)分别为：第1个系数是截矩(interccept)，它代表了第1组的总体均数，第2个系数是差值，这个差值表示的是第2组与第1组总体的差值。当我们进行统计学分析时，往往第2个系数是我们最感兴趣的，因为我们想知道这两组数据是否有差异。

我们可以在R中通过2行代码来编制实验设计。我们先是使用波浪线~来构建一个公式。这种形式的公式意味着，我们想使用波浪线右侧的变量的观测值来构建统计模型。然后我们再放进一个变量的名称，就能知道样本来源于哪些组。

看一个案例。
现在我们有两组数据，分别是对照组和高脂饲料组(hf)小鼠的体重数据。现在我们用1来表示对照组小鼠，2来表示hf组。我们首先要告诉R，这些值不能被当成数值，而是应该当成因子。然后使用公式~group来构建模型，如下所示：

group <- factor( c(1,1,2,2) )
model.matrix(~ group)

结果如下所示：

> group <- factor( c(1,1,2,2) )
> model.matrix(~ group)
  (Intercept) group2
1           1      0
2           1      0
3           1      1
4           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

这里面有很多信息，attr后面的信息不用管。

其实上面的计算过程也可以使用model.matrix()来构建，如下所示：

group <- factor( c(1,1,2,2) )
model.matrix(formula(~ group))

如下所示：

> group <- factor( c(1,1,2,2) )
> model.matrix(formula(~ group))
  (Intercept) group2
1           1      0
2           1      0
3           1      1
4           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

如果我们不对组别(group)进行因子转换，那么就会出现下面的情况：

group <- c(1,1,2,2)
model.matrix(~ group)

即：

> model.matrix(~ group)
  (Intercept) group
1           1     1
2           1     1
3           1     2
4           1     2
attr(,"assign")
[1] 0 1

这种形式不是设计矩阵。因为没有经过factor转换的数值就是单纯的数值，而非指示变量(indicator variable)，因子类型变量的数值与具体的数字无关，字符串也能生成这些变量，如下所示：

group <- factor(c("control","control","highfat","highfat"))
model.matrix(~ group)

结果如下所示：

> group <- factor(c("control","control","highfat","highfat"))
> model.matrix(~ group)
  (Intercept) grouphighfat
1           1            0
2           1            0
3           1            1
4           1            1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

多组设计

现在我们设计了3组数据，如下所示：

group <- factor(c(1,1,2,2,3,3))
model.matrix(~ group)

结果如下所示：

> group <- factor(c(1,1,2,2,3,3))
> model.matrix(~ group)
  (Intercept) group2 group3
1           1      0      0
2           1      0      0
3           1      1      0
4           1      1      0
5           1      0      1
6           1      0      1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

现在我们的设计矩阵中有了第3列，这个第3列属于第3组。另外一种方法就是通过在公式中指定+0来添加第3组，如下所示：

group <- factor(c(1,1,2,2,3,3))
model.matrix(~ group + 0)

结果如下所示：

> group <- factor(c(1,1,2,2,3,3))
> model.matrix(~ group + 0)
  group1 group2 group3
1      1      0      0
2      1      0      0
3      0      1      0
4      0      1      0
5      0      0      1
6      0      0      1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

注：在这一处我不太理解，原文是This group now fits a separate coefficient for each group.翻译为汉语则是这一组现在适用于每组的单独系数，等后文中有介绍的时候再补上，以下是其它的参考书（数学建模：基于R.薛毅 著）对这一内容的理解：

lm()中的参数formula是模型公式：

1. y ~ 1+x1 + x2这种形式，表示常数项，，的系数；
2. y ~ x1 + x2这种形式，去掉了1，其意义与原来加1的一样。
3. y ~ 0 + x1 + x2这种形式，添加了0，表示拟合成就齐次线性模型（或者是改为y ~ -1 + x1 + x2）。

多变量设计

前面我们构建的线性模型只有1个变量，即饲料。而在生命科学中，在一次实验中，往往会涉及多个变量，例如我们也许会研究饲料与性别这两个因素对体重的影响，在这种情况下，我们也许就会设计4组，如下所示：

diet <- factor(c(1,1,1,1,2,2,2,2))
sex <- factor(c("f","f","m","m","f","f","m","m"))
table(diet,sex)

组别如下所示：

> table(diet,sex)
    sex
diet f m
   1 2 2
   2 2 2

假如我们认为饲料对体重的影响和性别对体重的影响（这只是假设）是相同的，那么我们就可以构建一个线性模型，如下所示：

现在我们在R中对这个模型进行拟合，此时我们需要添加另外一个变量来构建一个设计矩阵，如下所示：

diet <- factor(c(1,1,1,1,2,2,2,2))
sex <- factor(c("f","f","m","m","f","f","m","m"))
model.matrix(~ diet + sex)

如下所示：

> diet <- factor(c(1,1,1,1,2,2,2,2))
> sex <- factor(c("f","f","m","m","f","f","m","m"))
> model.matrix(~ diet + sex)
  (Intercept) diet2 sexm
1           1     0    0
2           1     0    0
3           1     0    1
4           1     0    1
5           1     1    0
6           1     1    0
7           1     1    1
8           1     1    1
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$diet
[1] "contr.treatment"

attr(,"contrasts")$sex
[1] "contr.treatment"

这个设计矩阵包含一个截矩(intercept)，一个饲料变量(diet)，一个性别变量(sex)。我们可以说这个线性模型解释了组间差异与条件变量的差异。但是，如上所述，这个模型成立的前提是，我们假设了饮食和性别的差异会对小鼠的体重产生相同的影响。我们称这些情况为累加效应(additive effect)。对于每一个变量，我们添加一个效应时，不用考虑其它的变量是什么。针对这种情况，还有一种线性模型，这种线性模型里还会考虑两个变量之间的潜在交互作用。

如果我们要考虑交互作用的话，模型的构建如下所示：

model.matrix(~ diet + sex + diet:sex)
# OR
# model.matrix(~ diet*sex)

结果如下所示：

> 
> model.matrix(~ diet + sex + diet:sex)
  (Intercept) diet2 sexm diet2:sexm
1           1     0    0          0
2           1     0    0          0
3           1     0    1          0
4           1     0    1          0
5           1     1    0          0
6           1     1    0          0
7           1     1    1          1
8           1     1    1          1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$diet
[1] "contr.treatment"

attr(,"contrasts")$sex
[1] "contr.treatment"

指定比较水平relevel()

在构建设计矩阵时，我们所选水平的参考水平用于比较。默认情况下，我们都是按照第1个水平因素进行比较，但是我们可以使用relevel()函数来指定哪个水平因素作用为对照，例如我们选择第2组作为参考水平，如下所示：

group <- factor(c(1,1,2,2))
group <- relevel(group, "2")
model.matrix(~ group)

结果如下所示：

> model.matrix(~ group)
  (Intercept) group1
1           1      1
2           1      1
3           1      0
4           1      0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

这里涉及了relevel函数，此函数信息如下所示：

relevel()函数

全称reorder levels of factor，函数的功能是：对一个因子的水平重新排序，以通过指定ref来指定第1个水平作为参考水平（例如在回归中使用二元解释变量，可以指定哪个水平作为参考水平），其余的往后移。在contr.treatment对照中非常有用，contr.treatment通常会将第1个水平当作参考水平。
用法：
relevel(x,ref,...)，其中x是一个未排序的因子，ref：指定的参考水平，通常是一个字符串。

或者是在factor()中直接明确地指出，如下所示：

group <- factor(c(1,1,2,2))
group <- factor(group, levels=c("1","2"))
model.matrix(~ group)

结果如下所示：

> group <- factor(c(1,1,2,2))
> group <- factor(group, levels=c("1","2"))
> model.matrix(~ group)
  (Intercept) group2
1           1      0
2           1      0
3           1      1
4           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

model.matrix的计算

model.matrix()函数会捕获R全局环境中的变量，除非这些数据已经作为数据框的形式提供给到model.matrix()函数，如下所示：

group <- 1:4
model.matrix(~ group, data=data.frame(group=5:8))

结果如下所示：

> group <- 1:4
> model.matrix(~ group, data=data.frame(group=5:8))
  (Intercept) group
1           1     5
2           1     6
3           1     7
4           1     8
attr(,"assign")
[1] 0 1

连续型变量

前面内容中的变量都是指示数据（有的统计学书中指的是分类变量，哑变量），而在一些实验设计中，我们研究的是数值型变量（连续型变量），不需要将它们转换到第1列。例如在自由落体实验中，时间是一个连续型变量，如下所示：

tt <- seq(0,3.4,len=4)
model.matrix(~ tt + I(tt^2))

结果如下所示：

> tt <- seq(0,3.4,len=4)
> model.matrix(~ tt + I(tt^2))
  (Intercept)       tt   I(tt^2)
1           1 0.000000  0.000000
2           1 1.133333  1.284444
3           1 2.266667  5.137778
4           1 3.400000 11.560000
attr(,"assign")
[1] 0 1 2

在这个模型里，我们用到了I()函数，I()的功能是，对一个变量进行数学变换，在上面的公式里，I(tt^2)的意思是就是，通过I()，将tt^2这个变量当作是一个变量，如果不加I()，则是单纯地计算tt^2的值。查看I()的文档就可以发现，这个函数的全称是Inhibit Interpretation/Conversion of Objects。

在生命科学中，我们还有可能会涉及这样一种情况，例如我们想研究一个药物的量效关系，例如研究这个药物在0mg，10mg和20mg时的效果。
将连续型数据强制作为变量很难进行分析。这就是为什么指示变量只是假设2组的均值不同，而连续型变量会假设预测测变量和结果之间存在着某种特定的关系。

在自由落体实验中，我们有重力加速度这个理论的支持。而在父子身高案例中，由于这些数据是二元正态分布变量，如果我们加上一些条件的话，父子之间的身高存在着线性关系。但是，我们发现，如果不对年龄等这类变量进行“调整(adjust)”的话，并将它们包含在线性模型中，其统计结果就值得商榷。我们不支持将这样的变量加到线性模型中，除非已经有支持这种数据的模型出现。

练习

P183