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

最近我們被客戶要求撰寫關于線性混合模型的研究報告，包括一些圖形和統(tǒng)計輸出。

混合線性模型，又名多層線性模型(Hierarchical linear model)。它比較適合處理嵌套設計(nested)的實驗和調(diào)查研究數(shù)據(jù)

序言

此外，它還特別適合處理帶有被試內(nèi)變量的實驗和調(diào)查數(shù)據(jù)，因為該模型不需要假設樣本之間測量獨立，且通過設置斜率和截距為隨機變量，可以分離自變量在不同情境中(被試內(nèi)設計中常為不同被試)對因變量的作用。

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

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

非線性混合模型常常在生物制藥領域的分析中會用到，因為很多劑量反應并不是線性的，如果這個時候數(shù)據(jù)再有嵌套結構，那么就需要考慮非線性混合模型了。

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

數(shù)據(jù)

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

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

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

##?必須用手進行匯總ggplot(d,aes(x=Day,y=X,colour=Group))

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

使用nlme

用lme的線性擬合失敗。

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

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

可以看到考慮到組*日效應的模型也失敗了。

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

我試著用SSfpl擬合一個非線性模型，一個自啟動的四參數(shù)Logistic模型（參數(shù)為左漸近線、右漸近線、中點、尺度參數(shù)）。這對于nls擬合來說效果不錯，給出了合理的結果。

nlsfit1?<-?nls(X?~?SSfp)coef(nlsfit1)

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

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

gnls( ????X?~?SSfpl)

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

params=symp.R~Group

繪制預測值。

g1?+?geom_line()

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

dp?<-?data.frame(d,res=resid(gnlsfit2),fitted=fitted(gnlsfit2)) (diagplot1?<-?ggplot(dp,aes(x=factor(Individual), ??????????????y=res,colour=Group))+ ??????geom_boxplot(outlier.colour=NULL)+ ??scale_colour_brewer(palette="Dark2"))

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

anova(lm(res~Individual))

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

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

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

點擊標題查閱往期內(nèi)容

非線性混合效應 NLME模型對抗哮喘藥物茶堿動力學研究

左右滑動查看更多

01

02

03

04

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

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

右側漸近線中的方差估計值是非零的。

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

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

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

我們可以通過AIC或似然比檢驗來比較模型

AICtab(nlmefit1,nlmefit2,weights=TRUE)

anova(nlmefit1,nlmefit2)

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

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

##分母非常大的F檢驗。pf(diff/2,df1=2,df2=1000000,lower.tail=FALSE)

我們不知道真正相關的df，但上面的總結表明df是40。?

nlmer

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

基本的RE模型（沒有群體效應）運行良好。

?nlmer( ??X?~?SSfpl(Day,?asy,?as,?x,?s)?~ ?????????asy|Indi,)

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

1. 構建虛擬變量:

mm?<-?model.matrix(~Group,data=d) grp2?<-?mm[,2]

1. 構建一個函數(shù)來評估預測值及其梯度；分組結構是硬編碼的。

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

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

L1?<-?grep("^?+\.value?+<-") L2?<-?grep("^?+attr\(\.value",)eval(parse(text))

嘗試一下擬合:

nlmer( ??X?~?fpl(Day,?asym,?as,?asymp,?asR3,?xmi,?sca)?~ ?????????as|Indi, ?????start?=??list(nlpars)),data=d)

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

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

結果不理想

fixef(nlmerfit2)

range(predict(nlmerfit2))

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

AD模型生成器

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

部分原因是我對ADMB的熟悉程度較低，這有點費勁，最后我通過循序漸進的步驟才成功。

最小的例子

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

##設置數(shù)據(jù)：調(diào)整名稱，等等d0?<-?c(list(nobs=nrow(d)),as.list(d0))##起始值：調(diào)整名稱，增加數(shù)值names(svec3)?<-?gsub("\.","",names(svec3))??##?移除點svec3$asympR?<-?0.6?##?單一值##?運行?do_admb("algae0", ????????data, ????????params, ????????run.opts)

結果不錯

固定效應模型

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

model1?<-?" 參數(shù)部分 ???向量?pred(1,nobs)?//?預測值???向量Rval(1,nobs)?//預測值過程部分 ???pred?=?as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal)? "

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

model1B?<-?" 參數(shù)部分 ???向量?pred(1,nobs)?//?預測值???向量Rval(1,nobs)?//預測值過程部分 ???pred?=?asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc))?。 "

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

隨機效應

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

model2?<-?" 參數(shù)部分 ???向量?pred(1,nobs)?//?預測值???向量Rval(1,nobs)?//預測值過程部分 ???pred?=?asym+elem ???f?=?0.5*no*log(norm2(X-pr)/n)+norm2(R)。 "

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

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

model3?<-?" 參數(shù)部分 ???向量?prd(1,nobs)?//?預測值???向量Rl(1,nobs)?//?預測值???向量?scalal(1,nobs) ???向量xmal(1,nobs) ???sdror?opr(1,nobs)?//輸出預測值程序部分 ???Rval?=?XR*Rve+Rsma*（Z*Ru）。 ???xmval?=?Xd*xdvec;.... ???f?=?0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru) "

結果:

summary(admbfit3)

有一個非常大的AIC差異。如上文所示，對nlme擬合的似然比F測試是作為一種練習......

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

fit3B?<-?do_admb(, ????????data, ????????params, ????????re, ????????run.opts=run.control)plot2(list(cc),intercept=TRUE)

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

診斷圖

##放棄條件模式/樣本-R估計值diagplot1?%+%?dp2

也許這暗示了兩個實驗組中更大的差異？

擬合與殘差

diagplot2?%+%?dp2

疊加預測（虛線）:

g1?+?geom_line

如果能生成平滑的預測曲線（即對中間的日值），那就更好了，但也更繁瑣。

結論

• 從參數(shù)估計中得出的主要結論是，第三組下降得更早一些（xmidvec更小），同時下降得更遠（Rvec更低）。

似然分析

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

lapply(Rsigmavec,fitfun) ##?嘗試填補漏洞lapply(Rsigmavec2,fitfun)

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

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

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

ggplot(sigma,sqrt(2*(NLL-min(NLL))+ ??geom_point()

擴展

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

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

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

點擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例》。

點擊標題查閱往期內(nèi)容

R語言混合線性模型、多層次模型、回歸模型分析學生平均成績GPA和可視化
R語言線性混合效應模型（固定效應&隨機效應）和交互可視化3案例
R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)R語言 線性混合效應模型實戰(zhàn)案例
R語言混合效應邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言建立和可視化混合效應模型mixed effect model
R語言LME4混合效應模型研究教師的受歡迎程度
R語言 線性混合效應模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對學術能力測驗（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
SPSS中的多層（等級）線性模型Multilevel linear models研究整容手術數(shù)據(jù)
用SPSS估計HLM多層（層次）線性模型模型