R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例
原文鏈接：http://tecdat.cn/?p=23426
最近我們被客戶要求撰寫關(guān)于線性混合模型的研究報告，包括一些圖形和統(tǒng)計輸出。
混合線性模型，又名多層線性模型(Hierarchical linear model)。它比較適合處理嵌套設(shè)計(nested)的實驗和調(diào)查研究數(shù)據(jù)
序言
此外，它還特別適合處理帶有被試內(nèi)變量的實驗和調(diào)查數(shù)據(jù)，因為該模型不需要假設(shè)樣本之間測量獨立，且通過設(shè)置斜率和截距為隨機變量，可以分離自變量在不同情境中(被試內(nèi)設(shè)計中常為不同被試)對因變量的作用。
簡單的說，混合模型中把研究者感興趣的自變量對因變量的影響稱為固定效應(yīng)，把其他控制的情景變量稱為隨機效應(yīng)。由于模型中包括固定和隨機效應(yīng)，故稱為混合線性模型。無論是用方差分析進(jìn)行差異比較，還是回歸分析研究自變量對因變量的影響趨勢，混合線性模型比起傳統(tǒng)的線性模型都有更靈活的表現(xiàn)。
非線性混合模型就是通過一個連接函數(shù)將線性模型進(jìn)行拓展，并且同時再考慮隨機效應(yīng)的模型。
非線性混合模型常常在生物制藥領(lǐng)域的分析中會用到，因為很多劑量反應(yīng)并不是線性的，如果這個時候數(shù)據(jù)再有嵌套結(jié)構(gòu)，那么就需要考慮非線性混合模型了。
本文中我們用（非）線性混合模型分析藻類數(shù)據(jù)。這個問題的參數(shù)是：已知截距（0日值）在各組和樣本之間是相同的。
數(shù)據(jù)
用lattice和ggplot2繪制數(shù)據(jù)。
xyplot(jitter(X)~Day,?groups=Group)
ggplot版本有兩個小優(yōu)勢。1. 按個體和群體平均數(shù)添加線條[用stat_summary應(yīng)該和用xyplot的type="a "一樣容易])；2.調(diào)整點的大小，使重疊的點可視化。(這兩點當(dāng)然可以用自定義的 panel.xyplot 來實現(xiàn) ...)
##?必須用手進(jìn)行匯總ggplot(d,aes(x=Day,y=X,colour=Group))
從這些圖片中得出的主要結(jié)論是：（1）我們可能應(yīng)該使用非線性模型，而不是線性模型；（2）可能存在一些異方差（在較低的平均值上有較大的方差，好像在 X=0.7的數(shù)據(jù)有一個 "天花板"）；看起來可能存在個體間的變化（特別是基于t2的數(shù)據(jù)，其中個體曲線近乎平行）。然而，我們也將嘗試線性擬合來說明問題。
使用nlme
用lme的線性擬合失敗。
LME?<-?lme(X?~?1,?random?=?~Day|Individual,?data=d)
如果我們用control=lmeControl(msVerbose=TRUE))運行這個程序，就會得到輸出，最后是。?
可以看到考慮到組*日效應(yīng)的模型也失敗了。
LME1?<-?lme(X?~?Group*Day,?random?=?~Day|Individual,?data=d)
我試著用SSfpl擬合一個非線性模型，一個自啟動的四參數(shù)Logistic模型（參數(shù)為左漸近線、右漸近線、中點、尺度參數(shù)）。這對于nls擬合來說效果不錯，給出了合理的結(jié)果。
nlsfit1?<-?nls(X?~?SSfp)coef(nlsfit1)
可以用gnls來擬合組間差異（我需要指定起始值
我的第一次嘗試不太成功。
gnls(????X?~?SSfpl)
但如果我只允許asymp.R在各組之間變化，就能運行成功。
params=symp.R~Group
繪制預(yù)測值。
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\）值可以接受個體間不存在變異的無效假設(shè)...
更一般的診斷圖--殘差與擬合，同一個體的點用線連接?？梢园l(fā)現(xiàn)，隨著平均數(shù)的增加，方差會逐漸減小。
plot(dp,(x=fitted,y=res,colour=Group))
點擊標(biāo)題查閱往期內(nèi)容
非線性混合效應(yīng) NLME模型對抗哮喘藥物茶堿動力學(xué)研究
左右滑動查看更多
01
02
03
04
我不能用nlme來處理三個參數(shù)因組而異模型，但如果我只允許asymp變化，就可以運行。
nlme(model=list(fixed=with(c(asymp.R,xmid,scale,asymp.L),...)
右側(cè)漸近線中的方差估計值是非零的。
加入隨機效應(yīng)后，參數(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)
我們不知道真正相關(guān)的df，但上面的總結(jié)表明df是40。?
nlmer
我想現(xiàn)在可以為nlmer得到正確的模型規(guī)范，但我找不到一個方便的語法來進(jìn)行固定效應(yīng)建模（即在這種情況下允許一些參數(shù)因組而異）--當(dāng)我構(gòu)建了正確的語法，nlmer無法得到答案。
基本的RE模型（沒有群體效應(yīng)）運行良好。
?nlmer(??X?~?SSfpl(Day,?asy,?as,?x,?s)?~?????????asy|Indi,)
根據(jù)我的理解，人們只需要構(gòu)建自己的函數(shù)來封裝固定效應(yīng)結(jié)構(gòu)；為了與nlmer一起使用，該函數(shù)還需要計算相對于固定效應(yīng)參數(shù)的梯度。這有點麻煩，但可以通過修改派生函數(shù)生成的函數(shù)，使之稍微自動化。
構(gòu)建虛擬變量:
mm?<-?model.matrix(~Group,data=d)grp2?<-?mm[,2]
構(gòu)建一個函數(shù)來評估預(yù)測值及其梯度；分組結(jié)構(gòu)是硬編碼的。
deriv(~A+((B0+B1*grp2+B2*grp3-A)/(1+exp((x-xmid)/scale)
通過插入與傳遞給函數(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)
失敗了（但我認(rèn)為這是由于nlmer本身造成的，而不是設(shè)置有什么根本性的問題）。為了確定，我應(yīng)該按照同樣的思路生成一個更大的人工數(shù)據(jù)集，看看我是否能讓它工作起來。
現(xiàn)在我們可以用穩(wěn)定版（lme4.0）得到一個答案。
結(jié)果不理想
fixef(nlmerfit2)
range(predict(nlmerfit2))
我不能確定，在nlmer中是否有更簡單的方法來做固定效果。
AD模型生成器
我們還可以使用AD模型生成器來解決這個問題。它可以處理更復(fù)雜的模型，比如擬合更多參數(shù)的群體效應(yīng)。
部分原因是我對ADMB的熟悉程度較低，這有點費勁，最后我通過循序漸進(jìn)的步驟才成功。
最小的例子
首先嘗試沒有隨機效應(yīng)、分組變量等。(即等同于上面的nls擬合）。)
##設(shè)置數(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)
結(jié)果不錯
固定效應(yīng)模型
現(xiàn)在嘗試用固定效應(yīng)分組，使用上面構(gòu)建的虛擬變量（也可以使用if語句，或者用R[Group[i]]的for循環(huán)中的R值向量，或者（最佳選擇）為R傳遞一個模型矩陣...）。我們必須使用elem_div而不是/來對兩個向量進(jìn)行元素除法。
model1?<-?"參數(shù)部分???向量?pred(1,nobs)?//?預(yù)測值???向量Rval(1,nobs)?//預(yù)測值過程部分???pred?=?as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal)?"
試著用模型矩陣來代替它。
model1B?<-?"參數(shù)部分???向量?pred(1,nobs)?//?預(yù)測值???向量Rval(1,nobs)?//預(yù)測值過程部分???pred?=?asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc))?。"
當(dāng)然，在參數(shù)相同的情況下，也可以工作。
隨機效應(yīng)
現(xiàn)在添加隨機效應(yīng)?；貧w函數(shù)并沒有完全實現(xiàn)隨機效應(yīng)模型（盡管這應(yīng)該在即將到來的版本中被修復(fù)），所以我們用公式減去（n/2 log({RSS}/n)），其中RSS是殘差平方和。
model2?<-?"參數(shù)部分???向量?pred(1,nobs)?//?預(yù)測值???向量Rval(1,nobs)?//預(yù)測值過程部分???pred?=?asym+elem???f?=?0.5*no*log(norm2(X-pr)/n)+norm2(R)。"
由于ADMB不處理稀疏矩陣，也不懲罰循環(huán)，如果將隨機效應(yīng)實現(xiàn)為(i=1; i<=nobs; i++) Rval[i] += Rsigma*Ru[Group[i]]，效率會略高，但我是懶人/我喜歡矩陣表示的緊湊性和可擴展性.
現(xiàn)在我們終于可以測試R以外的參數(shù)的固定效應(yīng)差異了。
model3?<-?"參數(shù)部分???向量?prd(1,nobs)?//?預(yù)測值???向量Rl(1,nobs)?//?預(yù)測值???向量?scalal(1,nobs)???向量xmal(1,nobs)???sdror?opr(1,nobs)?//輸出預(yù)測值程序部分???Rval?=?XR*Rve+Rsma*（Z*Ru）。???xmval?=?Xd*xdvec;....???f?=?0.5*nobs*log(norm2(X-pred)/nobs)+norm2(Ru)"
結(jié)果:
summary(admbfit3)
有一個非常大的AIC差異。如上文所示，對nlme擬合的似然比F測試是作為一種練習(xí)......
對于該圖，最好是按組指定參數(shù)重新進(jìn)行擬合，而不是按基線+對比度進(jìn)行擬合。
fit3B?<-?do_admb(,????????data,????????params,????????re,????????run.opts=run.control)
plot2(list(cc),intercept=TRUE)
現(xiàn)在我們對標(biāo)準(zhǔn)化的問題很困擾，所以（經(jīng)過一番折騰）我們可以在不同的面板上重新畫出群體變化的參數(shù)。
診斷圖
##放棄條件模式/樣本-R估計值diagplot1?%+%?dp2
也許這暗示了兩個實驗組中更大的差異？
擬合與殘差
diagplot2?%+%?dp2
疊加預(yù)測（虛線）:
g1?+?geom_line
如果能生成平滑的預(yù)測曲線（即對中間的日值），那就更好了，但也更繁瑣。
結(jié)論
從參數(shù)估計中得出的主要結(jié)論是，第三組下降得更早一些（xmidvec更?。瑫r下降得更遠(yuǎn)（Rvec更低）。
似然分析
計算一個( sigma^2_R ) 似然函數(shù)的代碼并不難，但運行起來有點麻煩：它很慢，而且計算在置信度下限附近的幾個點上出現(xiàn)了非正-無限矩陣；我運行了另一組值，試圖充分覆蓋這個區(qū)域。
lapply(Rsigmavec,fitfun)##?嘗試填補漏洞lapply(Rsigmavec2,fitfun)
帶有插值樣條的剖面圖和似然比檢驗分界線。?
在sigma^2_R 上的95%剖面置信區(qū)間是{0.0386,0.2169}。
我沒有計算過，但轉(zhuǎn)換后的剖面圖（在對應(yīng)于偏離度與最小偏離度的平方根偏差的 y ）上，所以二次剖面將是一個對稱的V）顯示，二次近似對這種情況相當(dāng)糟糕 ...
ggplot(sigma,sqrt(2*(NLL-min(NLL))+??geom_point()
擴展
更多地討論分母df問題。參數(shù)引導(dǎo)法/MCMC？
我們可以嘗試在xmid和scale參數(shù)中加入隨機效應(yīng)。
在組間或作為X的函數(shù)的方差（無論是殘差還是個體間的方差）中可能有額外的模式。
點擊文末?“閱讀原文”
獲取全文完整代碼數(shù)據(jù)資料。
本文選自《R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例》。
點擊標(biāo)題查閱往期內(nèi)容
R語言混合線性模型、多層次模型、回歸模型分析學(xué)生平均成績GPA和可視化
R語言線性混合效應(yīng)模型（固定效應(yīng)&隨機效應(yīng)）和交互可視化3案例
R語言用lme4多層次（混合效應(yīng)）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)R語言 線性混合效應(yīng)模型實戰(zhàn)案例
R語言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應(yīng)模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言建立和可視化混合效應(yīng)模型mixed effect model
R語言LME4混合效應(yīng)模型研究教師的受歡迎程度
R語言 線性混合效應(yīng)模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
SPSS中的多層（等級）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)
用SPSS估計HLM多層（層次）線性模型模型
