1. 數(shù)據(jù)這里使用sleepstudy數(shù)據(jù)集，看一下免費(fèi)的R包lme4和付費(fèi)包asreml如何處理不同的混合線性模型，以加深對(duì)混合線性模型的理解。 「數(shù)據(jù)描述：」 ?睡眠剝奪研究中受試者每天的平均反應(yīng)時(shí)間。第0天，受試者有正常的睡眠時(shí)間。從那天晚上開始，他們每晚只能睡3個(gè)小時(shí)，依次進(jìn)行0~9天。觀察結(jié)果（y變量）代表了每天對(duì)每個(gè)受試者進(jìn)行的一系列測(cè)試的平均反應(yīng)時(shí)間。 ? 「數(shù)據(jù)預(yù)覽：」 > head(dat)
Reaction Days Subject
1 249.5600 0 308
2 258.7047 1 308
3 250.8006 2 308
4 321.4398 3 308
5 356.8519 4 308
6 414.6901 5 308
- Reaction:為觀測(cè)值，遇到刺激的反應(yīng)時(shí)間
- Subject：實(shí)驗(yàn)對(duì)象（ID）
「原數(shù)據(jù)可視化：」 這里，Subject為每個(gè)實(shí)驗(yàn)對(duì)象（人），做一下折線圖，看一下不同人在不同天數(shù)的反應(yīng)時(shí)間。 library(lme4)
data("sleepstudy")
dat = sleepstudy
ggplot(dat,aes(x = Days, y = Reaction, group = Subject, color = Subject)) +
geom_point() + geom_line()
 可以看到，不同的人差別比較大，不同的處理天數(shù)差別比較大，但是具體到每個(gè)人變化是不同的。
- 有些人在0天時(shí)，反應(yīng)時(shí)間比如高（截距），后面隨著天數(shù)的增加，增加得快，或者增加的慢（斜率）
- 有些人在0天時(shí)，反應(yīng)時(shí)間比較短，后面隨著天數(shù)的增加，增加得快，或者增加得慢
- 截距：intercept，為0天的反應(yīng)
lmer常用模型公式如下： mod= lmer(data = , formula = y ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor))
- Random_intercept，為隨機(jī)截距，即認(rèn)為不同群體因變量的分布不同（通俗的解釋：有的人生下來(lái)起點(diǎn)高，是富二代，有的人是一般群眾，起點(diǎn)低）
- Random_Slope，為隨機(jī)斜率，即認(rèn)為不同群體受固定因子的影響不同（通俗解釋：有的人是學(xué)霸，學(xué)習(xí)能力強(qiáng)，2個(gè)小時(shí)學(xué)會(huì)，斜率高；有的人是學(xué)渣，2天才能學(xué)會(huì)，斜率低）
?參考: https://zhuanlan.zhihu.com/p/63092231 ? 2. 模型1：截距隨機(jī)，斜率固定這種模型，認(rèn)為不同人起點(diǎn)有差異，但是隨著剝奪睡眠，他們的變化趨勢(shì)沒(méi)有差異（平行的） 這里，隨機(jī)因子為(1 | Subject)，可以擴(kuò)展為（1 + 1|Subject）認(rèn)為斜率是固定的，截距是隨機(jī)的。 library(lme4)
mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)
summary(mod1a)
library(asreml)
mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)
兩者結(jié)果一致： 使用asreml軟件的predict函數(shù)進(jìn)行預(yù)測(cè)，查看預(yù)測(cè)值的分布： pre = predict(mod1b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)
ggplot(pre,aes(x = Days, y = pre_value,
group = Subject, color = Subject)) +
geom_point() + geom_line()
 可以看出，每個(gè)人的斜率是一樣的，截距不一樣。
由結(jié)果可以寫出擬合的函數(shù)，比如fixed的值： - Intercept：251.4051，為整體均值，加上各個(gè)Subject的效應(yīng)值，即為各個(gè)Subject的斜率
 比如 Subject308，它的截距為：251.4051 + 40.7786 = 292.1837，那么它的公式為： 如框內(nèi)所示： 3. 模型2：截距隨機(jī)，斜率隨機(jī)，沒(méi)有相關(guān)性這里，不考慮相關(guān)性的寫法是|| mod2a = lmer(Reaction ~ Days + (Days || Subject), data=dat)
summary(mod2a)
mod2b = asreml(Reaction ~ Days, random = ~ Subject/Days, data=dat)
summary(mod2b)$varcomp
summary(mod2b,coef=T)$coef.fixed
兩者結(jié)果一致： 使用asreml軟件的predict函數(shù)進(jìn)行預(yù)測(cè)，查看預(yù)測(cè)值的分布： pre = predict(mod2b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)
ggplot(pre,aes(x = Days, y = pre_value,
group = Subject, color = Subject)) +
geom_point() + geom_line()
 可以看出，不同Subject個(gè)體，截距不一樣，斜率也不一樣。
4. 模型3：截距隨機(jī)，斜率隨機(jī)，考慮相關(guān)性mod3a = lmer(Reaction ~ Days + (Days | Subject), data=dat)
summary(mod3a) # 它沒(méi)有估計(jì)出協(xié)方差，asreml可以估計(jì)出協(xié)方差
mod3b = asreml(Reaction ~ Days,
random = ~ str(~Subject + Subject:Days, ~us(2):id(Subject)),
data=dat)
summary(mod3b)$varcomp
summary(mod3b,coef=T)$coef.fixed
summary(mod3b,coef=T)$coef.random
結(jié)果一致，asreml結(jié)果更完整： 使用asreml軟件的predict函數(shù)進(jìn)行預(yù)測(cè)，查看預(yù)測(cè)值的分布： pre = predict(mod3b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)
ggplot(pre,aes(x = Days, y = pre_value,
group = Subject, color = Subject)) +
geom_point() + geom_line()
 4. 模型3：斜率隨機(jī)，截距固定mod4a = lmer(Reaction ~ Days + (0 + Days | Subject), data=dat)
summary(mod4a)
mod4b = asreml(Reaction ~ Days,
random = ~ Subject:Days,
data=dat)
summary(mod4b)$varcomp
summary(mod4b,coef=T)$coef.fixed
ranef(mod4a)
summary(mod4b,coef=T)$coef.random
結(jié)果一致： 使用asreml軟件的predict函數(shù)進(jìn)行預(yù)測(cè)，查看預(yù)測(cè)值的分布： pre = predict(mod4b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)
ggplot(pre,aes(x = Days, y = pre_value,
group = Subject, color = Subject)) +
geom_point() + geom_line()
 可以看到，該模型截距都一樣，截距固定。
5. lme4包模型比較「模型1：」 mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)
「模型2：」 mod2a = lmer(Reaction ~ Days + (Days || Subject), data=dat)
「模型3：」 mod3a = lmer(Reaction ~ Days + (Days | Subject), data=dat)
「模型4：」 mod4a = lmer(Reaction ~ Days + (0 + Days | Subject), data=dat)
「模型比較：」 > anova(mod1a,mod2a,mod3a,mod4a)
refitting model(s) with ML (instead of REML)
Data: dat
Models:
mod1a: Reaction ~ Days + (1 | Subject)
mod4a: Reaction ~ Days + (0 + Days | Subject)
mod2a: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
mod3a: Reaction ~ Days + (Days | Subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mod1a 4 1802.1 1814.8 -897.04 1794.1
mod4a 4 1782.1 1794.8 -887.04 1774.1 19.9983 0
mod2a 5 1762.0 1778.0 -876.00 1752.0 22.0771 1 2.619e-06 ***
mod3a 6 1763.9 1783.1 -875.97 1751.9 0.0639 1 0.8004
---
Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1
可以看出，mod2a為最優(yōu)模型，它與mod3a不顯著，與其他模型達(dá)到極顯著水平。它的AIC和BIC也最低，是最優(yōu)模型。即模型中截距隨機(jī)，斜率隨機(jī)的模型最優(yōu)： 6. asreml包模型比較「模型1：」 mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)
「模型2：」 mod2b = asreml(Reaction ~ Days, random = ~ Subject/Days, data=dat)
「模型3：」 mod3b = asreml(Reaction ~ Days,
random = ~ str(~Subject + Subject:Days, ~us(2):id(Subject)),
data=dat)
「模型4：」 mod4b = asreml(Reaction ~ Days,
random = ~ Subject:Days,
data=dat)
「模型比較：」 aic1 = summary(mod1b)$bic
aic2 = summary(mod2b)$bic
aic3 = summary(mod3b)$bic
aic4 = summary(mod4b)$bic
bic1 = summary(mod1b)$bic
bic2 = summary(mod2b)$bic
bic3 = summary(mod3b)$bic
bic4 = summary(mod4b)$bic
re = data.frame(Model = paste0("Model: ",1:4), AIC = c(aic1,aic2,aic3,aic4),
BIC = c(bic1,bic2,bic3,bic4))
re
結(jié)果： > re
Model AIC BIC
1 Model: 1 1469.687 1469.687
2 Model: 2 1432.073 1432.073
3 Model: 3 1437.213 1437.213
4 Model: 4 1449.746 1449.746
可以看出，模型2最優(yōu)，它的AIC和BIC結(jié)果最低。 ?注意，asreml和lme4計(jì)算AIC和BIC的方法不一樣，所以結(jié)果有所差異。 ? 使用LRT似然比檢驗(yàn)?zāi)Ｐ停?/p>lrt.asreml(mod1b,mod2b)
lrt.asreml(mod1b,mod3b)
lrt.asreml(mod2b,mod3b)
結(jié)果： > lrt.asreml(mod1b,mod3b)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
mod3b/mod1b 2 42.837 1.545e-10 ***
---
Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1
> lrt.asreml(mod2b,mod3b)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
mod3b/mod2b 1 0.041056 0.4197
模型2最優(yōu)，模型2和模型3不顯著，模型2和模型1達(dá)到極顯著。 7. 模型2的散點(diǎn)圖和擬合圖合并這個(gè)模型最優(yōu)！
|
