设计矩阵
这个章节作者主要是在讲如何利用R中的formula()和model.matrix()函数为各种线性模型生成设计矩阵（也被称为模型矩阵）。比如说，在小鼠饲料与体重的例子中，作者将模型写成：
其中，是小鼠体重，只有当小鼠吃的饲料是高脂饲料时，
才等于1。我们将在一次测量中获得的N个不同的实体（different entities）称为试验单元（experimental units）。在本例中，小鼠就是试验单元。
试验单元这种类型的变量被称作指示变量（indicator variables）,原因是它们可以指示试验单元是否有某种特征。如前所述，我们利用线性代数去表示这个模型：
假设：
那么模型便可以被写成：
更简单的一种写法是：
其中，矩阵就是所谓的设计矩阵。一旦我们定义了一个设计矩阵，我们便可以着手去计算最小二乘估计了，这一过程被称作拟合模型（fitting the model）。在R语言中，直接给
lm()函数传递一个公式就可以进行拟合了，而lm()内部又会调用函数model.matrix()从而将你传入的公式与矩阵联系起来。
设计矩阵的选择
设计矩阵的选择是构建线性模型中至关重要的一步，因为它决定了那个系数会被拟合进模型当中，它还决定了样本之间的内部关系。一个常见的误解是，选择设计矩阵是直接看试验中样本如何描述而开始的，也就是在设计矩阵时，仅考虑一些诸如组别、试验批次等等基本信息是不够的，还需要考虑中的变量是如何解释
中的观测值这种类似的问题。
在上述例子中，我们是使用线性模型在不同的组别中进行比较。因此，最终我们得到的设计矩阵最少应该有两列：第一列是intercept列`，该列全部是1；第二列则指定哪些样本在第二组。
在小鼠体重这个例子里，有两个系数需要在线性模型中拟合，一个是intercept，代表第一个组的群体均值；另一个系数则代表两组的群体均值之差。
设计矩阵的代码是由波浪符~以及变量名组成，比如~ group，其中的~意思是我们要用~右边的变量去建模，它右边的变量的作用是告诉我们哪些样本在哪个组。
举例说明，假设我们有对照组和高脂饲料组，分别用1和2表示组别，且每组两个样本。我们首先要做的是告诉R这个值是因子水平的变量，然后在用~ group告诉R我们要用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"
也可以加上formula()函数，它的作用是告诉R传入它的表达式是公式：
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)
## (Intercept) group
## 1 1 1
## 2 1 1
## 3 1 2
## 4 1 2
## attr(,"assign")
## [1] 0 1
可以看到，得到的矩阵并不是我们想要的。这是因为传入函数中的group现在是数值型变量，但实际上我们是希望它能够指代组别，也就是第二列只有0或者1。
另外需要注意的是，传入model.matrix()和lm()函数的因子的具体水平的名字是啥都行，最重要的是哪些水平的顺序，举例来说：
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"
得到的设计矩阵与之前得到的相同。
多组别的设计矩阵
利用相同的公式，我们可以设计多组别的设计矩阵，现假设我们有第三种饲料：
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"
与之前得到的矩阵相比，上面的矩阵多了第三列，它的作用是指定哪些样本是第三组。在公式里我们还可以再+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"
这种矩阵给每一个组都拟合了一个系数，后续我们会对此进行更加深入的探索。
多变量设计矩阵
在生命科学领域，有很多实验是有多种变量的。比如，我们在小鼠体重的问题中再加上一个性别变量：
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)
## 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)
## (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"
之前我们提到过，上面的这个设计矩阵有个前提假设：两个变量的效应相同，但是往往在实验中常常会有累加效应（additive effect），如果在设计矩阵时，你觉得两个变量之间可能会有潜在的内在关系，你可以使用:：
model.matrix(~ diet + sex + diet:sex)
或者：
model.matrix(~ 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"
参考水平
在上个设计矩阵中，像diet2、sexm是参考水平(reference level)，它就是做对比用的。默认情况下，R会将按照字母顺序排序的第一个水平当做参考水平。如果想更换reference level，使用relevel()函数即可：
group <- factor(c(1,1,2,2))
group <- relevel(group, "2")
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"
model.matrix从哪里寻找数据？
除非在model.matrix()函数中传入data参数，否则的话该函数会从R的全局环境中寻找变量：
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
连续变量的设计矩阵
在一些特定的情况，我们会使用数值型变量去设计公式，比如说自由落体问题，在那个例子中，时间和时间的平方在模型中都是连续变量：
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()函数的作用其实就是保持参数自身类型不变。
但是需要注意的是，除非数据支持所使用的模型，否则最好不要使用连续型变量。比如说，正是因为我们知道自由落体问题有引力理论做支持，所以才去建立的模型，若提前并不知道的话，这个过程就会比较困难。
网友评论