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

• 線性代數(shù)
  categories:
• Data Analysis for the life sciences

前言

這一部分是《Data Analysis for the life sciences》的第5章線性模型的第1小節(jié)。這一部分的主要內(nèi)容涉及矩陣的表示方法。

背景知識

數(shù)據(jù)分析中的很多模型都可以用矩陣代數(shù)來表示恩敌。我們指的這些模型都是線性模型(linear models)窗慎。線性(Linear)并不是說它們是直線，而是說是線性組合酥郭。使用矩陣代數(shù)來描述線性模型非常方便，因為我們可以更好地構(gòu)建這些模型，并使用數(shù)學工具來方便計算办悟。在這一部分中，我們將詳細地描述我們?nèi)绾问褂镁仃嚧鷶?shù)來構(gòu)建并擬合線性模型滩褥。

在本書中病蛉，我們關(guān)注的線性模型是基于二分法的線性模型：例如治療組與對照組。在前面提到的小鼠飼料的案例中就是這種線性模型的一個典型案例瑰煎。在這一部分中铺然，我們會描述一些比較復雜的案例，不過還會繼續(xù)關(guān)注二分法變量的線性模型酒甸。

當我們學習線性模型時魄健，需要記住，我們使用的變量仍然是隨機變量插勤。這就意味著沽瘦，我們使用線性模型獲得的估計值也是隨機變量革骨。雖然數(shù)學上這理解起來比較復雜，但是我們前面學到的一些概念在這里仍然適用∥隽担現(xiàn)在可以先做一些練習良哲，在線性模型的練習中復習一些隨機變量的概念。

練習

P172

設計矩陣(The Design Matrix)

R中的formual()和model.matrix()可以用于生成各種線性模型的設計矩陣(design matrices)（也被稱為模型矩陣助隧，即model matrices)臂外。例如，在小鼠的飼料實驗中喇颁，我們可以這樣構(gòu)建模型漏健，如下所示：

其中，表示體重橘霎，而等于1的時候蔫浆，表示小鼠吃的是高脂飼料(hf)。我們可以使用實驗單位(Experimental unit)來表示N個不同的實驗動物姐叁，這里的實驗單位是小鼠瓦盛。我們稱這類變量為指示變量(indicator variables)，因此這類變量僅用于標明實驗單位外潜，表示某種干預的有無≡罚現(xiàn)在用矩陣代數(shù)來描述就是以下形式：

或為：

簡化形式為：

其中是設計矩陣，一旦我們定義了設計矩陣处窥，那么我們就能計算出最小二乘估計嘱吗。這個過程稱為擬合模型(fitting the model)。在R中滔驾，我們可以直接在lm()函數(shù)中輸入一個公式(formula)直接進行擬合谒麦。在后面的一個腳本中，我們會使用model.matrix()函數(shù)哆致，這個函數(shù)是在lm()函數(shù)內(nèi)部使用的绕德。這種使用方式哦可以方便地讓我們將R的formula與矩陣連接起來，有助于我們對lm()的結(jié)果進行解釋摊阀。

設計方案

設計矩陣的選擇是線性模型中的關(guān)鍵步驟耻蛇，因為它編碼了哪些系數(shù)將適合用于建模，以及樣本之間的相互關(guān)系胞此。一個常見的誤解就是臣咖，選擇實驗設計要遵循簡單明了的原則，即要在描述中直接標明樣本豌鹤。事實并非如此亡哄。關(guān)于每個樣本（無論是對照組還是處理組枝缔，實驗批次等）的基礎信息里并不意味著這是一個“正確”的設計矩陣布疙。設計矩陣還應該對中的變量是如何解釋的值這些信息進行編碼蚊惯，這是研究者們必須要考慮的。

在這一部分的案例中灵临，我們使用線性模型來比較不同組之間的差異截型。因此，最終進行計算的設計矩陣應該至少包含兩列儒溉，即截矩(intercept)列宦焦，截矩列是第1列，另外就是第2列顿涣，第2列指定樣本波闹。在這種情況下，線性模型中的兩個系數(shù)(coefficient)分別為：第1個系數(shù)是截矩(interccept)涛碑，它代表了第1組的總體均數(shù)精堕，第2個系數(shù)是差值，這個差值表示的是第2組與第1組總體的差值蒲障。當我們進行統(tǒng)計學分析時歹篓，往往第2個系數(shù)是我們最感興趣的，因為我們想知道這兩組數(shù)據(jù)是否有差異揉阎。

我們可以在R中通過2行代碼來編制實驗設計庄撮。我們先是使用波浪線~來構(gòu)建一個公式。這種形式的公式意味著毙籽，我們想使用波浪線右側(cè)的變量的觀測值來構(gòu)建統(tǒng)計模型洞斯。然后我們再放進一個變量的名稱，就能知道樣本來源于哪些組坑赡。

看一個案例巡扇。
現(xiàn)在我們有兩組數(shù)據(jù)，分別是對照組和高脂飼料組(hf)小鼠的體重數(shù)據(jù)】逯裕現(xiàn)在我們用1來表示對照組小鼠厅翔，2來表示hf組。我們首先要告訴R搀突，這些值不能被當成數(shù)值刀闷，而是應該當成因子。然后使用公式~group來構(gòu)建模型仰迁，如下所示：

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

結(jié)果如下所示：

> 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()來構(gòu)建徐许，如下所示：

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)進行因子轉(zhuǎn)換施蜜，那么就會出現(xiàn)下面的情況：

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

這種形式不是設計矩陣。因為沒有經(jīng)過factor轉(zhuǎn)換的數(shù)值就是單純的數(shù)值雌隅，而非指示變量(indicator variable)翻默，因子類型變量的數(shù)值與具體的數(shù)字無關(guān)缸沃，字符串也能生成這些變量，如下所示：

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

結(jié)果如下所示：

> 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"

多組設計

現(xiàn)在我們設計了3組數(shù)據(jù)修械，如下所示：

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

結(jié)果如下所示：

> 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"

現(xiàn)在我們的設計矩陣中有了第3列趾牧，這個第3列屬于第3組。另外一種方法就是通過在公式中指定+0來添加第3組肯污，如下所示：

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

結(jié)果如下所示：

> 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.翻譯為漢語則是這一組現(xiàn)在適用于每組的單獨系數(shù)，等后文中有介紹的時候再補上蹦渣，以下是其它的參考書（數(shù)學建模：基于R.薛毅 著）對這一內(nèi)容的理解：

lm()中的參數(shù)formula是模型公式：

1. y ~ 1+x1 + x2這種形式哄芜，表示常數(shù)項，柬唯，的系數(shù)忠烛；
2. y ~ x1 + x2這種形式，去掉了1权逗，其意義與原來加1的一樣美尸。
3. y ~ 0 + x1 + x2這種形式，添加了0斟薇，表示擬合成就齊次線性模型（或者是改為y ~ -1 + x1 + x2）师坎。

多變量設計

前面我們構(gòu)建的線性模型只有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

假如我們認為飼料對體重的影響和性別對體重的影響（這只是假設）是相同的盟萨，那么我們就可以構(gòu)建一個線性模型，如下所示：

現(xiàn)在我們在R中對這個模型進行擬合了讨，此時我們需要添加另外一個變量來構(gòu)建一個設計矩陣捻激，如下所示：

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)胞谭。我們可以說這個線性模型解釋了組間差異與條件變量的差異。但是男杈，如上所述丈屹，這個模型成立的前提是，我們假設了飲食和性別的差異會對小鼠的體重產(chǎn)生相同的影響伶棒。我們稱這些情況為累加效應(additive effect)旺垒。對于每一個變量彩库，我們添加一個效應時，不用考慮其它的變量是什么袖牙。針對這種情況，還有一種線性模型舅锄，這種線性模型里還會考慮兩個變量之間的潛在交互作用鞭达。

如果我們要考慮交互作用的話，模型的構(gòu)建如下所示：

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

結(jié)果如下所示：

> 
> 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()

在構(gòu)建設計矩陣時皇忿，我們所選水平的參考水平用于比較畴蹭。默認情況下，我們都是按照第1個水平因素進行比較鳍烁，但是我們可以使用relevel()函數(shù)來指定哪個水平因素作用為對照叨襟，例如我們選擇第2組作為參考水平，如下所示：

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

結(jié)果如下所示：

> 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函數(shù)幔荒，此函數(shù)信息如下所示：

relevel()函數(shù)

全稱reorder levels of factor糊闽，函數(shù)的功能是：對一個因子的水平重新排序，以通過指定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)

結(jié)果如下所示：

> 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()函數(shù)會捕獲R全局環(huán)境中的變量君编，除非這些數(shù)據(jù)已經(jīng)作為數(shù)據(jù)框的形式提供給到model.matrix()函數(shù)，如下所示：

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

結(jié)果如下所示：

> 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

連續(xù)型變量

前面內(nèi)容中的變量都是指示數(shù)據(jù)（有的統(tǒng)計學書中指的是分類變量川慌，啞變量）啦粹，而在一些實驗設計中，我們研究的是數(shù)值型變量（連續(xù)型變量）窘游，不需要將它們轉(zhuǎn)換到第1列唠椭。例如在自由落體實驗中，時間是一個連續(xù)型變量忍饰，如下所示：

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

結(jié)果如下所示：

> 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()函數(shù)，I()的功能是艾蓝，對一個變量進行數(shù)學變換力崇，在上面的公式里斗塘，I(tt^2)的意思是就是，通過I()亮靴，將tt^2這個變量當作是一個變量馍盟，如果不加I()，則是單純地計算tt^2的值茧吊。查看I()的文檔就可以發(fā)現(xiàn)贞岭，這個函數(shù)的全稱是Inhibit Interpretation/Conversion of Objects。

在生命科學中搓侄，我們還有可能會涉及這樣一種情況瞄桨，例如我們想研究一個藥物的量效關(guān)系，例如研究這個藥物在0mg讶踪，10mg和20mg時的效果芯侥。
將連續(xù)型數(shù)據(jù)強制作為變量很難進行分析。這就是為什么指示變量只是假設2組的均值不同乳讥，而連續(xù)型變量會假設預測測變量和結(jié)果之間存在著某種特定的關(guān)系柱查。

在自由落體實驗中，我們有重力加速度這個理論的支持云石。而在父子身高案例中物赶，由于這些數(shù)據(jù)是二元正態(tài)分布變量，如果我們加上一些條件的話留晚，父子之間的身高存在著線性關(guān)系酵紫。但是，我們發(fā)現(xiàn)错维，如果不對年齡等這類變量進行“調(diào)整(adjust)”的話奖地，并將它們包含在線性模型中，其統(tǒng)計結(jié)果就值得商榷赋焕。我們不支持將這樣的變量加到線性模型中参歹，除非已經(jīng)有支持這種數(shù)據(jù)的模型出現(xiàn)。

練習

P183