R語言實現混合模型
普通的線性迴歸只包含兩項影響因素，即固定效應（fixed-effect）和噪聲（noise）。噪聲是我們模型中沒有考慮的隨機因素。而固定效應是那些可預測因素，而且能完整的劃分總體。例如模型中的性別變數，我們清楚只有兩種性別，而且理解這種變數的變化對結果的影響。
那麼為什麼需要 Mixed-effect Model？因為有些現實的複雜資料是普通線性迴歸是處理不了的。例如我們對一些人群進行重複測量，此時存在兩種隨機因素會影響模型，一種是對某個人重複測試而形成的隨機噪聲，另一種是因為人和人不同而形成的隨機效應（random effect）。如果將一個人的測量資料看作一個組，隨機因素就包括了組內隨機因素（noise）和組間隨機因素（random effect）。這種巢狀的隨機因素結構違反了普通線性迴歸的假設條件。
你可能會把人員（組間的隨機效應）看作是一種分類變數放到普通線性迴歸模型中，但這樣作是得不償失的。有可能這個factor的level很多，可能會用去很多自由度。更重要的是，這樣作沒什麼意義。因為人員ID和性別不一樣，我們不清楚它的意義，而且它也不能完整的劃分總體。也就是說樣本資料中的路人甲，路人乙不能完全代表總體的人員ID。因為它是隨機的，我們並不關心它的作用，只是因為它會影響到模型，所以不得不考慮它。因此對於隨機效應我們只估計其方差，不估計其迴歸係數。
混合模型中包括了固定效應和隨機效應，而隨機效應有兩種方式來影響模型，一種是對截距影響，一種是對某個固定效應的斜率影響。前者稱為 Random intercept model，後者稱為 Random Intercept and Slope Model。Random intercept model的函式結構如下
Yij = a0 + a1*Xij + bi + eij
a0: 固定截距 a1: 固定斜率 b: 隨機效應（隻影響截距） X: 固定效應 e: 噪聲
混合線性模型有時又稱為多水平線性模型或層次結構線性模型由兩個部分來決定，固定效應部分+隨機效應部分，
二、R語言中的線性混合模型可用包
1、nlme包
這是一個比較成熟的R包，是R語言安裝時預設的包，它除了可以分析分層的線性混合模型，也可以處理非線性模型。在優勢方面，個人認為它可以處理相對複雜的線性和非線性模型，可以定義方差協方差結構，可以在廣義線性模型中定義幾種分佈函式和連線函式。
它的短板：隨機效應的定義過於呆板；資料量大時速度很慢；預設情況下不能處理系譜資料；不能處理多元資料。
2、lme4包
lme4包是由Douglas Bates開發，他也是nlme包的作者之一，相對於nlme包而言，它的執行速度快一點，對於睡覺效應·隨機效應的結構也可以更復雜一點，但是它的缺點也和nlme一樣：不能處理協方差和相關係數結構；它可以與構建係數的包連線，比如mmpedigree包，但是結合比較脆弱。
3、ASReml-R包
ASReml-R是ASReml的R版本，它的優點：可以處理複雜的隨機因子結構；可以處理多元資料；可以處理系譜資料；可以處理大批量的資料
主要的缺點：它是收費的，當然它對於不發達國家的科研機構是免費的，不過需要申請和被稽核；它的使用者主要是育種公司、科研機構等，它可以在各種平臺上執行，包括Windows、Linux、OS X等。
二、多水平模型案例分析
案例一：
1、首先匯入資料，檢視一下資料的結構
資料來源：一個傳統的裂區資料來說明不同軟體包的用法，這個資料oats是在MASS包中，是研究大麥品種和N肥處理的裂區試驗，其中品種為主區，肥料為裂區。
library(MASS)data(oats)names(oats) = c('block', 'variety', 'nitrogen', 'yield')head(oats) block variety nitrogen yield
1 I Victory 0.0cwt 111
2 I Victory 0.2cwt 130
3 I Victory 0.4cwt 157
4 I Victory 0.6cwt 174
5 I Golden.rain 0.0cwt 117
6 I Golden.rain 0.2cwt 114oats$mainplot=oats$varietyoats$subplot=oats$nitrogen2、nlme包：用這個包很簡單，y-變數寫在左邊，然後是固定因子，然後是隨機因子，注意1|block/mainplot是裂區試驗殘差的寫法，因為裡面有兩個殘差。程式碼如下：
library(lme4)library(nlme)m1.nlme = lme(yield ~ variety*nitrogen,random = ~ 1|block/mainplot,data = oats)summary(m1.nlme)方差分析結果為：anova(m1.nlme) numDF denDF F-value p-value
(Intercept) 1 45 245.14299 <.0001
variety 2 10 1.48534 0.2724
nitrogen 3 45 37.68562 <.0001
variety:nitrogen 6 45 0.30282 0.93223、lme4包lme4包的語法也相似，隨機效應有著和nlme相同的語法，不同的是lme4包它的結果給出了隨機效應的標準差，而不是方差。library(lme4)m1.lme4 = lmer(yield ~ variety*nitrogen + (1|block/mainplot),data = oats)summary(m1.lme4)anova(m1.lme4)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
variety 526.1 263.0 2 10 1.485 0.2724
nitrogen 20020.5 6673.5 3 45 37.686 2.458e-12 ***
variety:nitrogen 321.8 53.6 6 45 0.303 0.9322
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1案例二：利用nlme包下面我們以faraway包中的psid資料來舉例，該資料集是對美國人的收入情況進行調查所得到的，其中包括了年齡、教育、性別、時間和個體ID這幾個變數，我們希望瞭解這些因素對收入的影響。
age educ sex income year person 1 31 12 M 6000 68 1 2 31 12 M 5300 69 1 3 31 12 M 5200 70 1 4 31 12 M 6900 71 1 5 31 12 M 7500 72 1 6 31 12 M 8000 73 1
如果假設認為這些調查物件是同質的，也就是個體間沒有差異性，那麼可以將資料完全彙集（complete pooling）到一起，直接利用lm函式進行迴歸。但這個混合效應模型的同質假設往往不成立，資料彙集導致過度簡化。另一種思路是假設研究的異質性，將不同的個體分別進行迴歸，從而得到針對特定個體的估計值，這稱為不彙集（no pooling）。但這種方法導致每個迴歸所用到的樣本減少，從而難以估計統計量的標準差。
多層迴歸模型的思路是前兩者的折中，所以又稱為部分彙集（partial pooling）。在R語言中我們使用mgcv包中的lmer函式來完成這項工作。首先載入faraway包以便讀取psid資料集，然後載入mgcv包，再將年份資料中心化以方便解釋模型，最後用lmer函式進行建模。
install.packages("faraway")library(faraway)library(nlme) age educ sex income year person cyear
1 31 12 M 6000 68 1 -10
2 31 12 M 5300 69 1 -9
3 31 12 M 5200 70 1 -8
4 31 12 M 6900 71 1 -7
5 31 12 M 7500 72 1 -6
6 31 12 M 8000 73 1 -5library(mgcv)psid$cyear <- psid$year-78head(psid)model1=lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid) lmer函式使用和lm是類似的，一般變量表示固定效應，括號內豎線右側的person表示它是一個隨機效應，它與模型中其它變數相加，而且與年份cyear變數相乘，影響其斜率。這就是一個隨機效應模型。如果認為隨機效應隻影響模型截距，那麼固定效應迴歸模型可以用下面的公式
model2=lmer(log(income) ~ cyear*sex +age+educ+(1|person),psid)
為了判斷哪一個模型更為合適，可以使用anova函式，從結果中觀察到P值很小，判斷應當使用model1 anova(model1,model2)Data: psid
Models:
..1: log(income) ~ cyear * sex + age + educ + (1 | person)
object: log(income) ~ cyear * sex + age + educ + (cyear | person)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
..1 8 3943.8 3987.1 -1963.9 3927.8
object 10 3805.5 3859.6 -1892.7 3785.5 142.3 2 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1得到的模型結果還可以用各種泛型函式如summary、predict、resid進行進一步處理。當響應變數不符合正態分佈假設時，還可以用廣義多層迴歸進行(glmer)建模
案例三：
1、Generate a longitudinal dataset and convert it into long format
library(MASS) dat.tx.a <- mvrnorm(n=250, mu=c(30, 20, 28), Sigma=matrix(c(25.0, 17.5, 12.3, 17.5, 25.0, 17.5, 12.3, 17.5, 25.0), nrow=3, byrow=TRUE)) dat.tx.b <- mvrnorm(n=250, mu=c(30, 20, 22), Sigma=matrix(c(25.0, 17.5, 12.3, 17.5, 25.0, 17.5, 12.3, 17.5, 25.0), nrow=3, byrow=TRUE))
dat <- data.frame(rbind(dat.tx.a, dat.tx.b)) names(dat) <- c('measure.1', 'measure.2', 'measure.3')
根據指定的均值和協方差生成多元正態資料：MASS包中的mvrnorm()函式 mvrnorm(n,mean,sigma) measure.1 measure.2 measure.3 1 25.31761 20.89468 34.65525 2 19.30890 22.82886 33.24505 3 25.53250 16.27554 24.61767 4 27.52445 23.07668 32.62275 5 35.89934 28.24790 36.07344 6 25.71556 13.36878 26.69632
dat <- data.frame(subject.id = factor(1:500), tx = rep(c('A', 'B'), each = 250), dat)
subject.id tx measure.1 measure.2 measure.3
1 1 A 25.31761 20.89468 34.65525
2 2 A 19.30890 22.82886 33.24505
3 3 A 25.53250 16.27554 24.61767
4 4 A 27.52445 23.07668 32.62275
5 5 A 35.89934 28.24790 36.07344
6 6 A 25.71556 13.36878 26.69632rm(dat.tx.a, dat.tx.b) dat <- reshape(dat,varying = c('measure.1','measure.2','measure.3'),idvar = 'subject.id', direction = 'long') subject.id tx time measure 1.1 1 A 1 25.31761 2.1 2 A 1 19.30890 3.1 3 A 1 25.53250 4.1 4 A 1 27.52445 5.1 5 A 1 35.89934 6.1 6 A 1 25.71556 subject.id tx time measure 495.3 495 B 3 29.14236 496.3 496 B 3 17.03742 497.3 497 B 3 18.75601 498.3 498 B 3 20.80353 499.3 499 B 3 31.94328 500.3 500 B 3 25.78684
2、Multilevel growth models
有很多R包可以做多水平分析，其中 lme4對於一般模型（二分類及離散輸出）比較適合，
另外一個nlme包 比較適合連續輸出變數（正態或高斯分佈）
install.packages('lme4') library(Matrix) library(lme4)
m <- lmer(measure ~ time + (1 | subject.id), data = dat)
summary(m)
注：符號 DV ~ IV + (1 | rand.int) ，其中 DV 為輸出變數，IV 為自變數, 1 為自變數的係數或斜率,
rand.int 為隨機截距變數
Likewise, a random slopes model is specified using the syntax DV ~ IV + (rand.slope | rand.int).
結果顯示如下：
Linear mixed model fit by REML ['lmerMod'] Formula: measure ~ time + (1 | subject.id) Data: dat REML criterion at convergence: 9750.6 Scaled residuals: Min 1Q Median 3Q Max -2.41363 -0.69275 0.04628 0.69361 2.47612 Random effects:
Groups Name Variance Std.Dev. subject.id (Intercept) 8.136 2.852 Residual 32.268 5.681 Number of obs: 1500, groups: subject.id, 500 Fixed effects: Estimate Std. Error t value (Intercept) 29.2715 0.4085 71.66 time -2.2832 0.1796 -12.71 Correlation of Fixed Effects: (Intr) time -0.880
3、 Multilevel growth models with approximate p-values
install.packages('lmerTest') library(lmerTest) m <- lmer(measure ~ time + (1 | subject.id), data = dat) summary(m)
結果是一樣的，不過多了P值
4、Calculating 95% CI and PI
dat.new <- data.frame(time = 1:3) dat.new$measure <- predict(m, dat.new, re.form = NA) m.mat <- model.matrix(terms(m), dat.new) dat.new$var <- diag(m.mat %*% vcov(m) %*% t(m.mat)) + VarCorr(m)$subject.id[1] dat.new$pvar <- dat.new$var + sigma(m)^2 dat.new$ci.lb <- with(dat.new, measure - 1.96 * sqrt(var)) dat.new$ci.ub <- with(dat.new, measure + 1.96 * sqrt(var)) dat.new$pi.lb <- with(dat.new, measure - 1.96 * sqrt(pvar)) dat.new$pi.ub <- with(dat.new, measure + 1.96 * sqrt(pvar))
5、Plot the average values
library(ggplot2) ggplot(data = dat.new, aes(x = time, y = measure)) + geom_line(data = dat, alpha = .02, aes(group = subject.id)) + geom_errorbar(width = .02, colour = 'red', aes(x = time - .02, ymax = ci.ub, ymin = ci.lb)) + geom_line(colour = 'red', linetype = 'dashed', aes(x = time-.02)) + geom_point(size = 3.5, colour = 'red', fill = 'white', aes(x = time - .02)) + geom_errorbar(width = .02, colour = 'blue', aes(x = time + .02, ymax = pi.ub, ymin = pi.lb)) + geom_line(colour = 'blue', linetype = 'dashed', aes(x = time + .02)) + geom_point(size = 3.5, colour = 'blue', fill = 'white', aes(x = time + .02)) + theme_bw()