原文鏈接：http://tecdat.cn/?p=23426?

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

序言

混合線性模型，又名多層線性模型(Hierarchical linear model)。它比較適合處理嵌套設(shè)計(jì)(nested)的實(shí)驗(yàn)和調(diào)查研究數(shù)據(jù)。此外，它還特別適合處理帶有被試內(nèi)變量的實(shí)驗(yàn)和調(diào)查數(shù)據(jù)，因?yàn)樵撃Ｐ筒恍枰僭O(shè)樣本之間測(cè)量獨(dú)立，且通過設(shè)置斜率和截距為隨機(jī)變量，可以分離自變量在不同情境中(被試內(nèi)設(shè)計(jì)中常為不同被試)對(duì)因變量的作用。

簡(jiǎn)單的說，混合模型中把研究者感興趣的自變量對(duì)因變量的影響稱為固定效應(yīng)，把其他控制的情景變量稱為隨機(jī)效應(yīng)。由于模型中包括固定和隨機(jī)效應(yīng)，故稱為混合線性模型。無論是用方差分析進(jìn)行差異比較，還是回歸分析研究自變量對(duì)因變量的影響趨勢(shì)，混合線性模型比起傳統(tǒng)的線性模型都有更靈活的表現(xiàn)。
?

非線性混合模型就是通過一個(gè)連接函數(shù)將線性模型進(jìn)行拓展，并且同時(shí)再考慮隨機(jī)效應(yīng)的模型。

非線性混合模型常常在生物制藥領(lǐng)域的分析中會(huì)用到，因?yàn)楹芏鄤┝糠磻?yīng)并不是線性的，如果這個(gè)時(shí)候數(shù)據(jù)再有嵌套結(jié)構(gòu)，那么就需要考慮非線性混合模型了。

?本文中我們用（非）線性混合模型分析藻類數(shù)據(jù)。這個(gè)問題的參數(shù)是：已知截距（0日值）在各組和樣本之間是相同的。

數(shù)據(jù)

用lattice和ggplot2繪制數(shù)據(jù)。

xyplot(jitter(X)~Day, groups=Group)

ggplot版本有兩個(gè)小優(yōu)勢(shì)。1. 按個(gè)體和群體平均數(shù)添加線條[用stat_summary應(yīng)該和用xyplot的type="a "一樣容易])；2.調(diào)整點(diǎn)的大小，使重疊的點(diǎn)可視化。(這兩點(diǎn)當(dāng)然可以用自定義的 panel.xyplot 來實(shí)現(xiàn) ...)

1. ## 必須用手進(jìn)行匯總

2. ggplot(d,aes(x=Day,y=X,colour=Group))

從這些圖片中得出的主要結(jié)論是：（1）我們可能應(yīng)該使用非線性模型，而不是線性模型；（2）可能存在一些異方差（在較低的平均值上有較大的方差，好像在 X=0.7的數(shù)據(jù)有一個(gè) "天花板"）；看起來可能存在個(gè)體間的變化（特別是基于t2的數(shù)據(jù)，其中個(gè)體曲線近乎平行）。然而，我們也將嘗試線性擬合來說明問題。

使用nlme

用lme的線性擬合失敗。

LME <- lme(X ~ 1, random = ~Day|Individual, data=d)

如果我們用control=lmeControl(msVerbose=TRUE))運(yùn)行這個(gè)程序，就會(huì)得到輸出，最后是。?

可以看到考慮到組*日效應(yīng)的模型也失敗了。?

LME1 <- lme(X ~ Group*Day, random = ~Day|Individual, data=d)

我試著用SSfpl擬合一個(gè)非線性模型，一個(gè)自啟動(dòng)的四參數(shù)Logistic模型（參數(shù)為左漸近線、右漸近線、中點(diǎn)、尺度參數(shù)）。這對(duì)于nls擬合來說效果不錯(cuò)，給出了合理的結(jié)果。?

1. nlsfit1 <- nls(X ~ SSfp)

2. coef(nlsfit1)

可以用gnls來擬合組間差異（我需要指定起始值

我的第一次嘗試不太成功。?

1. gnls(

2. X ~ SSfpl)

?但如果我只允許asymp.R在各組之間變化，就能運(yùn)行成功。?

params=symp.R~Group

繪制預(yù)測(cè)值。

1. 


2. g1 + geom_line()

這些看起來很不錯(cuò)（如果能得到置信區(qū)間就更好了--需要使用delta法或bootstrapping）。?

1. dp <- data.frame(d,res=resid(gnlsfit2),fitted=fitted(gnlsfit2))

2. (diagplot1 <- ggplot(dp,aes(x=factor(Individual),

3. y=res,colour=Group))+

4. geom_boxplot(outlier.colour=NULL)+

5. scale_colour_brewer(palette="Dark2"))

6. 


除了7號(hào)樣本外，沒有很多證據(jù)表明個(gè)體間的變異......如果我們想忽略個(gè)體間的變異，可以用?

anova(lm(res~Individual))

大的（p\）值可以接受個(gè)體間不存在變異的無效假設(shè)...

更一般的診斷圖--殘差與擬合，同一個(gè)體的點(diǎn)用線連接?？梢园l(fā)現(xiàn)，隨著平均數(shù)的增加，方差會(huì)逐漸減小。

plot(dp,(x=fitted,y=res,colour=Group))

我不能用nlme來處理三個(gè)參數(shù)因組而異模型，但如果我只允許asymp變化，就可以運(yùn)行。?

nlme(model=list(fixed=with(c(asymp.R,xmid,scale,asymp.L),...)

右側(cè)漸近線中的方差估計(jì)值是非零的。

加入隨機(jī)效應(yīng)后，參數(shù)根本就沒有什么變化。?

最大的比例差異是3.1%（在比例參數(shù)中）。?

1. nlmefit2 <- update(list(asyR+xmd+scal+asp ~1),

2. start )

我們可以通過AIC或似然比檢驗(yàn)來比較模型

AICtab(nlmefit1,nlmefit2,weights=TRUE)

anova(nlmefit1,nlmefit2)

可以做一個(gè)F測(cè)試而不是 LRT（即考慮到有限大小的修正）。
?

pchisq(iff,df=2,lower.tail=FALSE)

?

1. ##分母非常大的F檢驗(yàn)。

2. pf(diff/2,df1=2,df2=1000000,lower.tail=FALSE)

我們不知道真正相關(guān)的df，但上面的總結(jié)表明df是40。?

nlmer

我想現(xiàn)在可以為nlmer得到正確的模型規(guī)范，但我找不到一個(gè)方便的語法來進(jìn)行固定效應(yīng)建模（即在這種情況下允許一些參數(shù)因組而異）--當(dāng)我構(gòu)建了正確的語法，nlmer無法得到答案。

基本的RE模型（沒有群體效應(yīng)）運(yùn)行良好。

1. nlmer(

2. X ~ SSfpl(Day, asy, as, x, s) ~

3. asy|Indi,)

根據(jù)我的理解，人們只需要構(gòu)建自己的函數(shù)來封裝固定效應(yīng)結(jié)構(gòu)；為了與nlmer一起使用，該函數(shù)還需要計(jì)算相對(duì)于固定效應(yīng)參數(shù)的梯度。這有點(diǎn)麻煩，但可以通過修改派生函數(shù)生成的函數(shù)，使之稍微自動(dòng)化。

1. 構(gòu)建虛擬變量:

1. mm <- model.matrix(~Group,data=d)

2. grp2 <- mm[,2]

1. 構(gòu)建一個(gè)函數(shù)來評(píng)估預(yù)測(cè)值及其梯度；分組結(jié)構(gòu)是硬編碼的。

deriv(~A+((B0+B1*grp2+B2*grp3-A)/(1+exp((x-xmid)/scale)

1. 通過插入與傳遞給函數(shù)的參數(shù)名稱相匹配的行來查看所產(chǎn)生的函數(shù)，并將這些參數(shù)名稱分配給梯度矩陣。

1. L1 <- grep("^ +\\.value +<-")

2. L2 <- grep("^ +attr\\(\\.value",)

3. eval(parse(text))

嘗試一下擬合:

1. nlmer(

2. X ~ fpl(Day, asym, as, asymp, asR3, xmi, sca) ~

3. as|Indi,

4. start = ?list(nlpars)),data=d)

失敗了（但我認(rèn)為這是由于nlmer本身造成的，而不是設(shè)置有什么根本性的問題）。為了確定，我應(yīng)該按照同樣的思路生成一個(gè)更大的人工數(shù)據(jù)集，看看我是否能讓它工作起來。

現(xiàn)在我們可以用穩(wěn)定版（lme4.0）得到一個(gè)答案。

結(jié)果不理想

fixef(nlmerfit2)

range(predict(nlmerfit2))

我不能確定，在nlmer中是否有更簡(jiǎn)單的方法來做固定效果。

AD模型生成器

我們還可以使用AD模型生成器來解決這個(gè)問題。它可以處理更復(fù)雜的模型，比如擬合更多參數(shù)的群體效應(yīng)。

部分原因是我對(duì)ADMB的熟悉程度較低，這有點(diǎn)費(fèi)勁，最后我通過循序漸進(jìn)的步驟才成功。

最小的例子

首先嘗試沒有隨機(jī)效應(yīng)、分組變量等。(即等同于上面的nls擬合）。)

1. ##設(shè)置數(shù)據(jù)：調(diào)整名稱，等等

2. d0 <- c(list(nobs=nrow(d)),as.list(d0))

3. ##起始值：調(diào)整名稱，增加數(shù)值

4. names(svec3) <- gsub("\\.","",names(svec3)) ?## 移除點(diǎn)

5. svec3$asympR <- 0.6 ## 單一值

6. ## 運(yùn)行

7. do_admb("algae0",

8. data,

9. params,

10. run.opts)

結(jié)果不錯(cuò)

固定效應(yīng)模型

現(xiàn)在嘗試用固定效應(yīng)分組，使用上面構(gòu)建的虛擬變量（也可以使用if語句，或者用R[Group[i]]的for循環(huán)中的R值向量，或者（最佳選擇）為R傳遞一個(gè)模型矩陣...）。我們必須使用elem_div而不是/來對(duì)兩個(gè)向量進(jìn)行元素除法。

1. 


2. model1 <- "

3. 參數(shù)部分

4. 向量 pred(1,nobs) // 預(yù)測(cè)值

5. 向量Rval(1,nobs) //預(yù)測(cè)值

6. 


7. 過程部分

8. pred = as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal)

9. "

10. 


試著用模型矩陣來代替它。

1. 


2. model1B <- "

3. 參數(shù)部分

4. 向量 pred(1,nobs) // 預(yù)測(cè)值

5. 向量Rval(1,nobs) //預(yù)測(cè)值

6. 


7. 過程部分

8. pred = asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc)) 。

9. "

10. 


當(dāng)然，在參數(shù)相同的情況下，也可以工作。

隨機(jī)效應(yīng)

現(xiàn)在添加隨機(jī)效應(yīng)?；貧w函數(shù)并沒有完全實(shí)現(xiàn)隨機(jī)效應(yīng)模型（盡管這應(yīng)該在即將到來的版本中被修復(fù)），所以我們用公式減去（n/2 log({RSS}/n)），其中RSS是殘差平方和。

1. model2 <- "

2. 參數(shù)部分

3. 向量 pred(1,nobs) // 預(yù)測(cè)值

4. 向量Rval(1,nobs) //預(yù)測(cè)值

5. 


6. 過程部分

7. pred = asym+elem

8. f = 0.5*no*log(norm2(X-pr)/n)+norm2(R)。

9. "

由于ADMB不處理稀疏矩陣，也不懲罰循環(huán)，如果將隨機(jī)效應(yīng)實(shí)現(xiàn)為(i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]]，效率會(huì)略高，但我是懶人/我喜歡矩陣表示的緊湊性和可擴(kuò)展性.

現(xiàn)在我們終于可以測(cè)試R以外的參數(shù)的固定效應(yīng)差異了。

1. 


2. model3 <- "

3. 參數(shù)部分

4. 向量 prd(1,nobs) // 預(yù)測(cè)值

5. 向量Rl(1,nobs) // 預(yù)測(cè)值

6. 向量 scalal(1,nobs)

7. 向量xmal(1,nobs)

8. sdror opr(1,nobs) //輸出預(yù)測(cè)值

9. 


10. 程序部分

11. Rval = XR*Rve+Rsma*（Z*Ru）。

12. xmval = Xd*xdvec;....

13. f = 0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru)

14. "

結(jié)果:

summary(admbfit3)

有一個(gè)非常大的AIC差異。如上文所示，對(duì)nlme擬合的似然比F測(cè)試是作為一種練習(xí)......

對(duì)于該圖，最好是按組指定參數(shù)重新進(jìn)行擬合，而不是按基線+對(duì)比度進(jìn)行擬合。

1. 


2. fit3B <- do_admb(,

3. data,

4. params,

5. re,

6. run.opts=run.control)

1. plot2(list(cc),intercept=TRUE)

2. 
   

現(xiàn)在我們對(duì)標(biāo)準(zhǔn)化的問題很困擾，所以（經(jīng)過一番折騰）我們可以在不同的面板上重新畫出群體變化的參數(shù)。

診斷圖

1. ##放棄條件模式/樣本-R估計(jì)值

2. 


3. diagplot1 %+% dp2

也許這暗示了兩個(gè)實(shí)驗(yàn)組中更大的差異？

擬合與殘差

diagplot2 %+% dp2

疊加預(yù)測(cè)（虛線）:

g1 + geom_line

如果能生成平滑的預(yù)測(cè)曲線（即對(duì)中間的日值），那就更好了，但也更繁瑣。

結(jié)論

• 從參數(shù)估計(jì)中得出的主要結(jié)論是，第三組下降得更早一些（xmidvec更?。?，同時(shí)下降得更遠(yuǎn)（Rvec更低）。

似然分析

計(jì)算一個(gè)( sigma^2_R ) 似然函數(shù)的代碼并不難，但運(yùn)行起來有點(diǎn)麻煩：它很慢，而且計(jì)算在置信度下限附近的幾個(gè)點(diǎn)上出現(xiàn)了非正-無限矩陣；我運(yùn)行了另一組值，試圖充分覆蓋這個(gè)區(qū)域。

1. lapply(Rsigmavec,fitfun)

2. ## 嘗試填補(bǔ)漏洞

3. lapply(Rsigmavec2,fitfun)

帶有插值樣條的剖面圖和似然比檢驗(yàn)分界線。?

在sigma^2_R 上的95%剖面置信區(qū)間是{0.0386,0.2169}。

我沒有計(jì)算過，但轉(zhuǎn)換后的剖面圖（在對(duì)應(yīng)于偏離度與最小偏離度的平方根偏差的 y ）上，所以二次剖面將是一個(gè)對(duì)稱的V）顯示，二次近似對(duì)這種情況相當(dāng)糟糕 ...

1. ggplot(sigma,sqrt(2*(NLL-min(NLL))+

2. geom_point()

擴(kuò)展

• 更多地討論分母df問題。參數(shù)引導(dǎo)法/MCMC？

• 我們可以嘗試在xmid和scale參數(shù)中加入隨機(jī)效應(yīng)。

• 在組間或作為X的函數(shù)的方差（無論是殘差還是個(gè)體間的方差）中可能有額外的模式。

最受歡迎的見解

1.基于R語言的lmer混合線性回歸模型

2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

3.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

4.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例2

5.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

6.線性混合效應(yīng)模型Linear Mixed-Effects Models的部分折疊Gibbs采樣

7.R語言LME4混合效應(yīng)模型研究教師的受歡迎程度

8.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長(zhǎng)

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM