R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數據實例

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

序言

混合線性模型，又名多層線性模型(Hierarchical linear model)。它比較適合處理嵌套設計(nested)的實驗和調查研究數據。此外，它還特別適合處理帶有被試內變量的實驗和調查數據，由於該模型不須要假設樣本之間測量獨立，且經過設置斜率和截距爲隨機變量，能夠分離自變量在不一樣情境中(被試內設計中常爲不一樣被試)對因變量的做用。

簡單的說，混合模型中把研究者感興趣的自變量對因變量的影響稱爲固定效應，把其餘控制的情景變量稱爲隨機效應。因爲模型中包括固定和隨機效應，故稱爲混合線性模型。不管是用方差分析進行差別比較，仍是迴歸分析研究自變量對因變量的影響趨勢，混合線性模型比起傳統的線性模型都有更靈活的表現。
 

非線性混合模型就是經過一個鏈接函數將線性模型進行拓展，而且同時再考慮隨機效應的模型。

非線性混合模型經常在生物製藥領域的分析中會用到，由於不少劑量反應並非線性的，若是這個時候數據再有嵌套結構，那麼就須要考慮非線性混合模型了。

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

數據

用lattice和ggplot2繪製數據。

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

ggplot版本有兩個小優點。1. 按個體和羣體平均數添加線條[用stat_summary應該和用xyplot的type="a "同樣容易])；2.調整點的大小，使重疊的點可視化。(這兩點固然能夠用自定義的 panel.xyplot 來實現 ...)

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

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

使用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擬合一個非線性模型，一個自啓動的四參數Logistic模型（參數爲左漸近線、右漸近線、中點、尺度參數）。這對於nls擬合來講效果不錯，給出了合理的結果。 

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

能夠用gnls來擬合組間差別（我須要指定起始值

個人第一次嘗試不太成功。 

gnls(
    X ~ SSfpl)

 但若是我只容許asymp.R在各組之間變化，就能運行成功。 

params=symp.R~Group

繪製預測值。

g1 + geom_line()

這些看起來很不錯（若是能獲得置信區間就更好了--須要使用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號樣本外，沒有不少證據代表個體間的變異......若是咱們想忽略個體間的變異，能夠用 

anova(lm(res~Individual))

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

更通常的診斷圖--殘差與擬合，同一個體的點用線鏈接。能夠發現，隨着平均數的增長，方差會逐漸減少。

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

我不能用nlme來處理三個參數因組而異模型，但若是我只容許asymp變化，就能夠運行。 

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

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

加入隨機效應後，參數根本就沒有什麼變化。 

最大的比例差別是3.1%（在比例參數中）。 

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

我想如今能夠爲nlmer獲得正確的模型規範，但我找不到一個方便的語法來進行固定效應建模（即在這種狀況下容許一些參數因組而異）--當我構建了正確的語法，nlmer沒法獲得答案。

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

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

根據個人理解，人們只須要構建本身的函數來封裝固定效應結構；爲了與nlmer一塊兒使用，該函數還須要計算相對於固定效應參數的梯度。這有點麻煩，但能夠經過修改派生函數生成的函數，使之稍微自動化。

2. 構建虛擬變量:

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

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

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

2. 經過插入與傳遞給函數的參數名稱相匹配的行來查看所產生的函數，並將這些參數名稱分配給梯度矩陣。

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自己形成的，而不是設置有什麼根本性的問題）。爲了肯定，我應該按照一樣的思路生成一個更大的人工數據集，看看我是否能讓它工做起來。

如今咱們能夠用穩定版（lme4.0）獲得一個答案。

結果不理想

fixef(nlmerfit2)

range(predict(nlmerfit2))

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

AD模型生成器

咱們還可使用AD模型生成器來解決這個問題。它能夠處理更復雜的模型，好比擬合更多參數的羣體效應。

部分緣由是我對ADMB的熟悉程度較低，這有點費勁，最後我經過按部就班的步驟才成功。

最小的例子

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

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

結果不錯

固定效應模型

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

model1 <- "
參數部分
   向量 pred(1,nobs) // 預測值
   向量Rval(1,nobs) //預測值

過程部分
   pred = as+elem(Rval-asy,1.0+exp(-(Day-xmid)/scal) 
"

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

model1B <- "
參數部分
   向量 pred(1,nobs) // 預測值
   向量Rval(1,nobs) //預測值

過程部分
   pred = asym+ele(Rv-asy,1.0+exp(-(Da-xmi)/sc)) 。
"

固然，在參數相同的狀況下，也能夠工做。

隨機效應

如今添加隨機效應。迴歸函數並無徹底實現隨機效應模型（儘管這應該在即將到來的版本中被修復），因此咱們用公式減去（n/2 log({RSS}/n)），其中RSS是殘差平方和。

model2 <- "
參數部分
   向量 pred(1,nobs) // 預測值
   向量Rval(1,nobs) //預測值

過程部分
   pred = asym+elem
   f = 0.5\*no\*log(norm2(X-pr)/n)+norm2(R)。
"

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

如今咱們終於能夠測試R之外的參數的固定效應差別了。

model3 <- "
參數部分
   向量 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測試是做爲一種練習......

對於該圖，最好是按組指定參數從新進行擬合，而不是按基線+對比度進行擬合。

fit3B <- do_admb(,
        data,
        params,
        re,
        run.opts=run.control)

plot2(list(cc),intercept=TRUE)

如今咱們對標準化的問題很困擾，因此（通過一番折騰）咱們能夠在不一樣的面板上從新畫出羣體變化的參數。

診斷圖

##放棄條件模式/樣本-R估計值

diagplot1 %+% dp2

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

擬合與殘差

diagplot2 %+% dp2

疊加預測（虛線）:

g1 + geom_line

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

結論

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

似然分析

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

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

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

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

我沒有計算過，但轉換後的剖面圖（在對應於偏離度與最小偏離度的平方根誤差的 y ）上，因此二次剖面將是一個對稱的V）顯示，二次近似對這種狀況至關糟糕 ...

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

擴展

• 更多地討論分母df問題。參數引導法/MCMC？
• 咱們能夠嘗試在xmid和scale參數中加入隨機效應。
• 在組間或做爲X的函數的方差（不管是殘差仍是個體間的方差）中可能有額外的模式。

---

 

最受歡迎的看法

1.基於R語言的lmer混合線性迴歸模型

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

3.R語言線性混合效應模型實戰案例

4.R語言線性混合效應模型實戰案例2

5.R語言線性混合效應模型實戰案例

6.線性混合效應模型Linear Mixed-Effects Models的部分摺疊Gibbs採樣

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

8.R語言中基於混合數據抽樣(MIDAS)迴歸的HAR-RV模型預測GDP增加

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