快一個月沒更新文章啦芋簿，今天收到好幾個粉絲的催更私信，好的吧璃饱，實在對不住大家期待的眼神与斤，看樣子不能再拖啦，想想寫啥好呢荚恶，大家咨詢比較多的撩穿，混合模型算一個，今天就繼續(xù)給大家寫寫混合模型如何做吧裆甩。

混合模型一般都可以用lme4這個包解決冗锁，lme4既可以做線性混合模型齐唆，也可以做廣義線性混合模型還可以做非線性混合模型嗤栓，大家有需要可以只研究這一個包就行。

所謂混合模型就是既有固定效應又有隨機效應的模型：

“mixedeffects”, denotes a model that incorporates both fixed- and random-effects terms in a linear predictor expression from which the conditional mean of the response can be evaluated

第一部分 線性混合模型

直接上例子箍邮，數(shù)據(jù)是來自一篇研究睡眠剝奪的文獻茉帅，整個數(shù)據(jù)大概長下圖這樣，其中我們的受試者在day0的時候可以睡到自然醒锭弊，在之后的日子里所有的受試者就只能睡3個小時了堪澎，我們的響應變量是Reaction，就是對受試者做的測驗的響應時間味滞，我現(xiàn)在關心睡眠剝奪后樱蛤，響應時間的變化情況：

image

對于這么一個縱向數(shù)據(jù)钮呀，

我們來捋一捋：我們只有18個人受試者，每個受試者隨訪10次昨凡，我們需要明白的是爽醋，此時我們的每一次測量是嵌套在人的水平上的，我們可以認為便脊，不同人自己的10次測量是有強烈的相關性的蚂四，而不同人之間的這種關系又不一定是相同的。

直觀一點哪痰，我們可以畫出來每一天所有人響應時間和睡眠剝奪的變化遂赠，畫出來就是下圖：

image

可以看到我們上面的這個大圖是由很多個小圖組成的，每一個小圖中橫軸就是睡眠剝奪的時間晌杰，縱軸是反應時間跷睦。每個小圖就代表著我們要研究的睡眠剝奪和反應時間的關系（具體到人），但是我們也應該注意到這種關系在不同的人上是不同的肋演，體現(xiàn)在：關系的斜率不同和截距不同送讲。（這個關系的不同可以很明顯的在圖中看出來）

所以我們就可以擬合一個帶有隨機效應的混合模型：

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

運行代碼后得到下面的結(jié)果：

image

結(jié)果中有隨機效應的標準差和固定效應的β估計，我得到的截距是251.4惋啃，斜率是10.5哼鬓，這兩個系數(shù)就是我們研究的總體關系的表示，通常需要在文獻中匯報边灭，就意味著睡眠不剝奪的時候人的反應時間是251.4异希，而睡眠每剝奪多一天反應時間增長10.5。

上面這個是最簡單的混合模型绒瘦。我們繼續(xù)看：

lme4包高水平設置介紹

• 混合模型公式

對于一個常見的混合模型称簿，我們可以在lme4包中寫出來如下差不多的混合模型公式：

resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...

這個公式中FEexpr就是固定效應，(REexpr1 | factor1)and(REexpr2 | factor2)都是隨機效應惰帽，理論上你可以弄很多個隨機效應但是實際操作中我們不會關心那么多憨降。

• 理解混合模型公式

我們看到每一個隨機效應在公式中的表達都是(expr | factor)這樣的。豎杠前面的expr就是一個常規(guī)的回歸公式该酗，豎杠后面的factor就是一個常規(guī)的因子授药，你可以把豎杠想象成回歸公式和因子的交互：

One way to think about the vertical bar operator is as a special kind of interaction betweenthe model matrix and the grouping factor。This interaction ensures that the columns of themodel matrix have different effects for each level of the grouping factor.

這種交互的意思就是在因子的不同水平呜魄，我們的回歸是不一樣的悔叽，這也正好和我們前面的解釋相對應，就是在不同的人的水平睡眠剝奪和響應時間的關系不一樣爵嗅。

寫到這娇澎，希望大家能記住下面這張表：

image

這個表就給我們展示了常見的隨機效應的設置，比如(1 | g)睹晒，就是說在因子g的不同水平趟庄，我們響應變量的截距都不一樣括细。表中的第二行有個offset，表示沒有固定效應戚啥。如果我們的數(shù)據(jù)是一個三層嵌套數(shù)據(jù)勒极，我們可以用第三行的設定來表示隨機截距；如果你的數(shù)據(jù)沒有直接嵌套但是在g1和g2的不同水平上存在相關虑鼎，那么可以用第四行的設定辱匿，這個在項目反應理論中比較常見。

在lme4中炫彩，默認認為同一個模型的截距和斜率是存在相關的匾七，如果你確定截距和斜率無關那么設定隨機效應的時候就可以用兩個豎杠，或者把截距和斜率分開來寫江兢，就是說(x || g)和x +(1 | g) + (0 + x | g)表達的隨機效應都是一樣的昨忆。

比如如果我認為睡眠剝奪和反應時間隨機效應的截距和斜率無關，我便可以做如下設定：

fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)#截距和斜率無關的設定

有時候我們擬合一個后又想嘗試對模型進行改變杉允，但又不想重寫邑贴，此時就可以直接對相似的模型基礎上進行更新：

• 模型的更新

比如我想在fm1的基礎上去掉隨機斜率只留隨機截距，我就可以用updata寫出如下代碼：

fm3 <- update(fm1, . ~ . - (Days | Subject) + (1 | Subject))#模型的更新

到底哪一個模型更好呢叔磷？

可以用anova方法進行模型間的比較：

anova(fm1, fm2, fm3)

運行代碼會輸出比較的結(jié)果：

image

其中拢驾，從模型比較的結(jié)果可以看出，給模型增加一個截距和斜率無關的隨機效應相比會使得模型的deviance變小改基，進一步將隨機效應設定為相關繁疤，并不能夠顯著地減小deviance，從而我們就可以知道fm2才是對數(shù)據(jù)擬合最好的模型秕狰。

第二部分 非線性混合模型

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

The fixed-effects parameters describe the general patterns of the data and random-effects parameters describe specific clusters. If the model is nonlinear in the parameters,it is called a nonlinear mixed-effects model (Davidian &Giltinan, 2003)

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

看下面這個圖叹放，這個圖描繪了不同人用了茶堿過后的反應，時間是橫軸低飒，殘留是縱軸许昨，和開篇線性模型中睡眠剝奪和反應時間的例子一樣懂盐，我們把每個人的關系都做了圖出來褥赊，不過從圖中可以明顯看出這種關系并不是簡單線性的。

image

其實這種不是線性的關系存在的情況很多莉恼。

比如漸進回歸：

image

再比如邏輯增長：

image

此時我們要注意到像這兩非線性關系模型的參數(shù)都不是簡單的一個斜率加個截距了拌喉。都有φ1速那，φ2，φ3三個額外參數(shù)尿背。

這兒先給大家寫一個邏輯增長的實際例子：我現(xiàn)在有一個關于樹木周徑的數(shù)據(jù)集端仰，每棵樹隨訪了7次，每次隨訪記錄數(shù)的年齡age田藐，和周徑荔烧，我現(xiàn)在想研究在所有樹木中時間和周徑的關系。

image

很自然汽久，我們可以想到不同的樹這個關系應該是不一樣的鹤竭，我們想探求的一定是考慮了樹水平的變異之后的總體關系，所以不妨先畫出來每個樹的關系：

image

從圖中可以看到我們總共有5棵樹景醇，基本關系是一致的臀稚，但存在些許變異相關（所以考慮混合模型），而且這個關系并不是線性的（時間越大周徑基本不改變了）三痰，所以我們應該考慮非線性的混合模型吧寺。

具體地，我們可以用nlmer方法來擬合非線性混合模型散劫，方法參數(shù)包括3部分：首先是響應變量稚机，然后是非線性函數(shù)，然后是混合效應公式：

The formula argument fornlmeris in three parts: the response, the nonlinear model function depending on covariates and a set of nonlinear model (nm) parameters, and the mixed-effects formula.

比如對我們的數(shù)據(jù)我就可以寫出如下SSlogis方法的代碼：

print(nm1 <- nlmer(circumference ~ SSlogis(age,    Asym, xmid, scal) ~ Asym | Tree, 
                   Orange, start = c(Asym = 200,    xmid = 770, scal = 120)), corr = FALSE)

此時我們選擇的非線性函數(shù)是邏輯增長函數(shù)SSlogis获搏，剛剛給大家解釋了這個函數(shù)是有3個參數(shù)的抒钱，在上面的代碼中，age是我們的預測變量颜凯，Asym, xmid, scal分別是額外的三個參數(shù)（之前的邏輯增長的式子和Asym/(1+exp((xmid-input)/scal))等同）：

image

進一步谋币，擬合邏輯增長是要我們給出這些參數(shù)的初始值的，然后從初始值通過梯度下降尋找各個參數(shù)的最優(yōu)解：

SSlogis has an attribute called "initial", which is a function that nls can call to compute reasonable starting values for fitting a logistic function to the input data.

所以我們看到代碼中都給出了響應參數(shù)的初始值症概。

運行上面代碼后輸出如下結(jié)果：

image

我們可以看到結(jié)果中的固定效應里面有Asym蕾额，Xmid，scal參數(shù)的估計結(jié)果彼城。

那么這些參數(shù)如何解釋呢诅蝶？

留個懸念我們下期再更。關注關注關注募壕，嘿嘿调炬。

小結(jié)

很久沒更新了，今天給大家寫了如何用lme4做混合模型舱馅，包括線性和非線性的例子缰泡，感謝大家耐心看完，自己的文章都寫的很細代嗤，代碼都在原文中棘钞，希望大家都可以自己做一做缠借，請關注后私信回復“數(shù)據(jù)鏈接”獲取所有數(shù)據(jù)和本人收集的學習資料。如果對您有用請先收藏宜猜，再點贊轉(zhuǎn)發(fā)泼返。

也歡迎大家的意見和建議，大家想了解什么統(tǒng)計方法都可以在文章下留言姨拥，說不定我看見了就會給你寫教程哦绅喉，另咨詢代做請私信。