美文网首页
PH525x series - Expressing desig

PH525x series - Expressing desig

作者: 3between7 | 来源:发表于2019-11-27 16:17 被阅读0次

    设计矩阵

    这个章节作者主要是在讲如何利用R中的formula()和model.matrix()函数为各种线性模型生成设计矩阵（也被称为模型矩阵）。比如说，在小鼠饲料与体重的例子中，作者将模型写成：
    Y_i = \beta_0 + \beta_1x_i + ε,i = 1,\cdots,N

    其中，Y_i是小鼠体重，只有当小鼠吃的饲料是高脂饲料时，x_i才等于1。我们将在一次测量中获得的N个不同的实体（different entities）称为试验单元（experimental units）。在本例中，小鼠就是试验单元。

    试验单元这种类型的变量被称作指示变量（indicator variables）,原因是它们可以指示试验单元是否有某种特征。如前所述，我们利用线性代数去表示这个模型：

    假设：
    Y = \left\{\begin{matrix}Y_1 \\Y_2 \\ \vdots \\ Y_N \end{matrix}\right\},X = \left\{\begin{matrix}1 &X_1 \\1&X_2 \\ \vdots& \vdots \\ 1&X_N \end{matrix}\right\},\beta = \left\{\begin{matrix}\beta_0 \\\beta_1 \end{matrix}\right\} ,ε = \left\{\begin{matrix}ε_1 \\ε_2 \\ \vdots \\ ε_N \end{matrix}\right\}
    那么模型Y_i = \beta_0 + \beta_1x_i + ε,i = 1,\cdots,N便可以被写成：
    \left\{\begin{matrix}Y_1 \\Y_2 \\ \vdots \\ Y_N \end{matrix}\right\} = \left\{\begin{matrix}1 &X_1 \\1&X_2 \\ \vdots& \vdots \\ 1&X_N \end{matrix}\right\} \left\{\begin{matrix}\beta_0 \\\beta_1 \end{matrix}\right\} + \left\{\begin{matrix}ε_1 \\ε_2 \\ \vdots \\ ε_N \end{matrix}\right\}

    更简单的一种写法是：
    \mathbf{Y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}

    其中，矩阵\mathbf{X}就是所谓的设计矩阵。一旦我们定义了一个设计矩阵，我们便可以着手去计算最小二乘估计了，这一过程被称作拟合模型（fitting the model）。在R语言中，直接给lm()函数传递一个公式就可以进行拟合了，而lm()内部又会调用函数model.matrix()从而将你传入的公式与矩阵\mathbf{X}联系起来。

    设计矩阵的选择

    设计矩阵的选择是构建线性模型中至关重要的一步，因为它决定了那个系数会被拟合进模型当中，它还决定了样本之间的内部关系。一个常见的误解是，选择设计矩阵是直接看试验中样本如何描述而开始的，也就是在设计矩阵时，仅考虑一些诸如组别、试验批次等等基本信息是不够的，还需要考虑\mathbf{X}中的变量是如何解释\mathbf{Y}中的观测值这种类似的问题。

    在上述例子中，我们是使用线性模型在不同的组别中进行比较。因此，最终我们得到的设计矩阵最少应该有两列：第一列是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

    如果我们假设饲料与性别的影响是相同的（也就是互相独立），那么我们的线性模型就会变成：
    Y_i = \beta_0 + \beta_1x_{i,1} + \beta_2x_{i,2} + ε_i

    那么在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()函数的作用其实就是保持参数自身类型不变。

    但是需要注意的是，除非数据支持所使用的模型，否则最好不要使用连续型变量。比如说，正是因为我们知道自由落体问题有引力理论做支持，所以才去建立的模型，若提前并不知道的话，这个过程就会比较困难。

    原文链接

    相关文章

      网友评论

          本文标题：PH525x series - Expressing desig

          本文链接：https://www.haomeiwen.com/subject/gctpwctx.html

          延伸阅读

          深度阅读

          • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

          栏目导航

          热点阅读

          关于我们|服务条款|联系我们|PH525x series - Expressing desig|投稿指南|网站地图|RSS订阅|排版工具|手机版

          提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

          Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

          备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

          本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

          所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！