更多統(tǒng)計(jì)方法請點(diǎn)擊：醫(yī)學(xué)統(tǒng)計(jì)系列教程，可能是中文互聯(lián)網(wǎng)領(lǐng)域最全面的R語言醫(yī)學(xué)統(tǒng)計(jì)教程！持續(xù)更新中。所有數(shù)據(jù)和代碼都可以免費(fèi)獲取，數(shù)據(jù)文件在QQ群文件，代碼可直接復(fù)制粘貼。為了方便大家，還制作了以下合集： 本文目錄： 理論知識 數(shù)據(jù)探索 空模型 添加1級水平的固定效應(yīng) 添加2級水平的固定效應(yīng) 具有隨機(jī)斜率的MLM 具有交互效應(yīng)的MLM 重復(fù)測量數(shù)據(jù)的MLM 廣義混合效應(yīng)模型 參考資料
理論知識理論知識我就不獻(xiàn)丑了，直接給大家推薦高手的解讀！ 關(guān)于多水平模型（multi-level models，MLM）的概念和理論知識，強(qiáng)烈推薦閱讀馮國雙老師的幾篇文章，這是我目前見過的寫的最通俗易懂的。多水平模型、混合模型、隨機(jī)效應(yīng)模型、固定效應(yīng)模型、隨機(jī)系數(shù)模型、方差成分模型等等，全都詳細(xì)介紹到了： 我之前接觸多水平模型很少，對一些概念很模糊，但是讀完這些推文后，我感到豁然開朗、恍然大悟，強(qiáng)烈推薦大家仔細(xì)讀一讀，對于初學(xué)者幫助很大。 上面是理論知識，下面是R語言實(shí)戰(zhàn)。 數(shù)據(jù)探索數(shù)據(jù)可以從這個網(wǎng)站免費(fèi)下載：https://www./13-appendix.html?；蛘呒尤隥Q群免費(fèi)下載。 library(dplyr) # 數(shù)據(jù)操作
library(ggplot2) # 可視化
library(lme4) # 多水平模型
library(lmerTest) # 計(jì)算P值
library(performance) # 計(jì)算模型表現(xiàn)
data <- read.csv('datasets/heck2011.csv') # 加載數(shù)據(jù)
簡單介紹下這個數(shù)據(jù)，一共有6871行，各個變量的含義如下，其中math是因變量，建模目的是用其他變量預(yù)測學(xué)生的數(shù)學(xué)成績（math），或者叫探索和學(xué)生成績有關(guān)系的變量： schcode:學(xué)校id，一共有419所學(xué)校Rid:每個學(xué)校內(nèi)的學(xué)生id，每個學(xué)校內(nèi)都是從1開始ses:衡量學(xué)校內(nèi)學(xué)生的社會經(jīng)濟(jì)地位構(gòu)成的Z-scorefemses:以性別（女性）為中心的測量學(xué)生社會經(jīng)濟(jì)地位的變量math:學(xué)生數(shù)學(xué)成績測試分?jǐn)?shù)ses_mean:以學(xué)校為單位衡量學(xué)生社會經(jīng)濟(jì)地位的均值pro4yrc:每個學(xué)校中計(jì)劃就讀四年制大學(xué)的學(xué)生比例public:學(xué)校類型，1=公立學(xué)校，0=私立學(xué)校
str(data)
## 'data.frame': 6871 obs. of 10 variables:
## $ schcode : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Rid : int 1 2 3 4 5 6 7 8 9 10 ...
## $ id : int 6701 6702 6703 6704 6705 6706 6707 6708 6709 6710 ...
## $ female : int 1 1 1 0 0 0 0 1 0 1 ...
## $ ses : num 0.586 0.304 -0.544 -0.848 0.001 -0.106 -0.33 -0.891 0.207 -0.341 ...
## $ femses : num 0.586 0.304 -0.544 0 0 0 0 -0.891 0 -0.341 ...
## $ math : num 47.1 63.6 57.7 53.9 58 ...
## $ ses_mean: num -0.267 -0.267 -0.267 -0.267 -0.267 ...
## $ pro4yrc : num 0.0833 0.0833 0.0833 0.0833 0.0833 ...
## $ public : int 0 0 0 0 0 0 0 0 0 0 ...
這個數(shù)據(jù)很明顯是有層次的，學(xué)生是分別屬于不同的學(xué)校的，所以這是一個兩個水平的數(shù)據(jù)。學(xué)生是1級水平，學(xué)校是2級水平。 不同學(xué)校之間的水平是有差異的，為了演示先選擇其中10個學(xué)校（一共有419個學(xué)校）： data_sub <- data %>%
filter(schcode <= 10)
如果不考慮學(xué)校水平的差異，探索社會經(jīng)濟(jì)地位（ses）和成績（math）之間的關(guān)系，可以畫出如下的散點(diǎn)圖和擬合線： data_sub %>%
ggplot(mapping = aes(x = ses, y = math)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, fullrange = TRUE)
## `geom_smooth()` using formula = 'y ~ x'
 如果考慮到學(xué)校之間的不同水平，按照學(xué)校分別擬合，會得到如下的圖： data_sub %>%
ggplot(mapping = aes(x = ses, y = math, colour = factor(schcode))) +
geom_point() +
geom_smooth(mapping = aes(group = schcode), method = "lm", se = FALSE, fullrange = TRUE) +
labs(colour = "schcode")
## `geom_smooth()` using formula = 'y ~ x'
 每個學(xué)校的ses的截距和斜率都是不一樣的，所以不同學(xué)校對學(xué)生成績的影響也是不同的，普通的多元線性回歸沒有考慮到這種差異，但是多水平模型可以發(fā)現(xiàn)這些差異并量化它們。 空模型先給大家演示一下最簡單的，也就是“只有隨機(jī)截距的模型”（random intercept only model），又被稱為空模型（null model）。在這種模型中，沒有任何自變量來解釋學(xué)生成績的差異，模型只考慮了學(xué)校間的差異，這些差異是由隨機(jī)截距來表示的。 null_model <- lmer(math ~ (1|schcode), data = data)
summary(null_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ (1 | schcode)
## Data: data
##
## REML criterion at convergence: 48877.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6336 -0.5732 0.1921 0.6115 5.2989
##
## Random effects:
## Groups Name Variance Std.Dev.
## schcode (Intercept) 10.64 3.262
## Residual 66.55 8.158
## Number of obs: 6871, groups: schcode, 419
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.6742 0.1883 416.0655 306.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
在lme4和lmerTest中（lmerTest可以在結(jié)果中給出P值，lme4不會給出P值，這是這兩個包的差別），擬合多水平模型的基本語法是： # DV是因變量，IV是自變量
lmer(DV ~ 1 + IV1 + IV2 + ... + IVp + (random_effect1 + random_effect2 + ... + random_effect3 | grouping_variable), data = dataset)
該數(shù)據(jù)的因變量是math，波浪號右邊的(1|schcode)表示分組變量schcode的隨機(jī)截距，空模型沒有其他自變量。 結(jié)果主要是兩部分： - Random effects：隨機(jī)效應(yīng)的估計(jì)
- Fixed effects：固定效應(yīng)的估計(jì)
先看固定效應(yīng)：截距的固定效應(yīng)為57.67，表示所有學(xué)校的平均數(shù)學(xué)成績，此時這個截距其實(shí)就是因變量的平均值，即：mean(data$math)=57.73391（微小差異，可以忽略）； 再看隨機(jī)效應(yīng)：不同學(xué)校之間成績的方差（Variance）為10.64，不同學(xué)校之間的學(xué)生成績會圍繞總體均值（也就是此時的截距57.67）波動，其標(biāo)準(zhǔn)差（Std.Dev）是3.262，同一個學(xué)校內(nèi)不同學(xué)生之間成績的方差為66.55，標(biāo)準(zhǔn)差為8.158。 組內(nèi)相關(guān)系數(shù)（intraclass correlation coefficient，ICC）是衡量多水平模型中群體間變異（即隨機(jī)效應(yīng)）與總變異（群體間變異+群體內(nèi)變異）之比的一種指標(biāo)。在多水平模型中，ICC用來評估因變量的變異有多少可以歸因于群體間的差異，而有多少是群體內(nèi)的個體差異。 若ICC為0，說明群體間的差異為0，說明數(shù)據(jù)不存在層次結(jié)構(gòu)，在本例中也就是學(xué)校之間沒有差異，可以不使用多水平模型，使用普通的多元線性回歸即可。 根據(jù)模型輸出，該例的總方差為10.64+66.55=77.19（var(data$math)），ICC=10.64/77.19=0.138；也就是說成績總變異中的13.8%是學(xué)校不同導(dǎo)致的，剩下的86.2%的變異是由同一學(xué)校內(nèi)學(xué)生之間的差異造成的。也就是說，大部分成績變異來源于學(xué)校內(nèi)部的差異，而學(xué)校之間的差異相對較小。 ICC也可以使用performance包自動計(jì)算： performance::icc(null_model)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.138
## Unadjusted ICC: 0.138
由于沒有其他自變量，所以調(diào)整的ICC和未調(diào)整的ICC是一樣的。 添加1級水平的固定效應(yīng)上面演示的空模型沒有自變量，下面我們添加一個社會經(jīng)濟(jì)狀況（ses）作為自變量，這個ses是在水平1單位（不同的學(xué)生，不同的學(xué)校屬于水平2單位）間變化的，所以這個變量產(chǎn)生的效應(yīng)屬于水平1單位的固定效應(yīng)。 像這種只有隨機(jī)截距，沒有隨機(jī)斜率（但是有自變量，不是“空模型”）的多水平模型被稱為隨機(jī)截距模型，又叫方差成分模型，相關(guān)概念的理解請參考本文開頭推薦的幾篇文章。 用lme4擬合隨機(jī)截距模型： # RMEL表示限制性最大似然估計(jì)，這也是默認(rèn)方法
ses_l1 <- lmer(math ~ ses + (1|schcode), data = data, REML = TRUE)
summary(ses_l1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + (1 | schcode)
## Data: data
##
## REML criterion at convergence: 48215.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7733 -0.5540 0.1303 0.6469 5.6908
##
## Random effects:
## Groups Name Variance Std.Dev.
## schcode (Intercept) 3.469 1.863
## Residual 62.807 7.925
## Number of obs: 6871, groups: schcode, 419
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.5960 0.1329 375.6989 433.36 <2e-16 ***
## ses 3.8739 0.1366 3914.6382 28.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ses -0.025
先看固定效應(yīng)：Intercept的估計(jì)值為57.5960，表示當(dāng)所有自變量都為0的時候，因變量的預(yù)測值；ses的估計(jì)值為3.8739，表示ses每增加一個單位，數(shù)學(xué)成績可以提高3.8739分；固定效應(yīng)的解讀和普通的多元線性回歸沒有差別。 加入自變量后，我們可以發(fā)現(xiàn)固定效應(yīng)中截距的值發(fā)生了變化，因?yàn)榇藭r的截距不再是因變量的平均值，它是一個基于當(dāng)前模型的預(yù)測值。 再看隨機(jī)效應(yīng)：schcode的方差為3.469，殘差的方差為62.807，該結(jié)果與上面的模型的結(jié)果解讀是類似的。schcode作為分組變量，它的方差越大說明這個變量對因變量的影響越大，殘差的方差表示模型無法解釋的隨機(jī)誤差部分，這個值越大說明模型無法解釋的隨機(jī)誤差越大。 未調(diào)整的ICC可以量化分層變量（這里是schcode）所能解釋的方差（變異）比例。在模型中加入ses這個自變量后，未調(diào)整的ICC=0.046： performance::icc(ses_l1)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.052
## Unadjusted ICC: 0.046
說明在考慮社會經(jīng)濟(jì)地位的影響后，數(shù)學(xué)成績差異中有4.6%是學(xué)校不同造成的?？梢钥吹郊尤?code style="color: rgb(155, 110, 35);font-size: 14px;line-height: 1.8em;letter-spacing: 0em;background-attachment: scroll;background-clip: border-box;background-color: rgb(255, 245, 227);background-image: none;background-origin: padding-box;background-position-x: 0%;background-position-y: 0%;background-repeat: no-repeat;background-size: auto;width: auto;height: auto;margin-top: 0px;margin-bottom: 0px;margin-left: 2px;margin-right: 2px;padding-top: 2px;padding-bottom: 2px;padding-left: 4px;padding-right: 4px;border-top-style: none;border-bottom-style: none;border-left-style: none;border-right-style: none;border-top-width: 3px;border-bottom-width: 3px;border-left-width: 3px;border-right-width: 3px;border-top-color: rgb(0, 0, 0);border-bottom-color: rgba(0, 0, 0, 0.4);border-left-color: rgba(0, 0, 0, 0.4);border-right-color: rgba(0, 0, 0, 0.4);border-top-left-radius: 4px;border-top-right-radius: 4px;border-bottom-right-radius: 4px;border-bottom-left-radius: 4px;overflow-wrap: break-word;font-family: 39;Operator Mono39;, Consolas, Monaco, Menlo, monospace;word-break: break-all;">ses這個變量后，不同學(xué)校能夠解釋的變異減少了（13.8%到4.6%），這個也很好理解，因?yàn)橛幸徊糠肿儺惐?code style="color: rgb(155, 110, 35);font-size: 14px;line-height: 1.8em;letter-spacing: 0em;background-attachment: scroll;background-clip: border-box;background-color: rgb(255, 245, 227);background-image: none;background-origin: padding-box;background-position-x: 0%;background-position-y: 0%;background-repeat: no-repeat;background-size: auto;width: auto;height: auto;margin-top: 0px;margin-bottom: 0px;margin-left: 2px;margin-right: 2px;padding-top: 2px;padding-bottom: 2px;padding-left: 4px;padding-right: 4px;border-top-style: none;border-bottom-style: none;border-left-style: none;border-right-style: none;border-top-width: 3px;border-bottom-width: 3px;border-left-width: 3px;border-right-width: 3px;border-top-color: rgb(0, 0, 0);border-bottom-color: rgba(0, 0, 0, 0.4);border-left-color: rgba(0, 0, 0, 0.4);border-right-color: rgba(0, 0, 0, 0.4);border-top-left-radius: 4px;border-top-right-radius: 4px;border-bottom-right-radius: 4px;border-bottom-left-radius: 4px;overflow-wrap: break-word;font-family: 39;Operator Mono39;, Consolas, Monaco, Menlo, monospace;word-break: break-all;">ses解釋了，另外殘差的變異也減小了（66.55到62.807），也是這個原因?qū)е碌摹?/p> 多水平模型肯定是要比單水平模型的擬合程度更好的，因?yàn)樗軌蚪忉尫纸M變量導(dǎo)致的變異，也就是讓不能解釋的變異更少了。下面我們擬合一個普通的多元線性回歸，并比較一下兩個模型。 f <- lm(math ~ ses, data = data)
compare_performance(f,ses_l1,metrics = "common") # 比較下兩個模型
## # Comparison of Model Performance Indices
##
## Name | Model | AIC (weights) | BIC (weights) | RMSE | R2 | R2 (adj.) | R2 (cond.) | R2 (marg.) | ICC
## --------------------------------------------------------------------------------------------------------------------------
## f | lm | 48304.0 (<.001) | 48324.5 (<.001) | 8.131 | 0.143 | 0.143 | | |
## ses_l1 | lmerModLmerTest | 48219.1 (>.999) | 48246.4 (>.999) | 7.810 | | | 0.167 | 0.121 | 0.052
結(jié)果表明，AIC、BIC、RMSE都是多水平模型更小，也就是模型表現(xiàn)更好。 多元線性回歸中常用的提取結(jié)果的方法也都適用于多水平模型，比如計(jì)算可信區(qū)間等： confint(ses_l1)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 1.575429 2.144559
## .sigma 7.789560 8.063828
## (Intercept) 57.335234 57.856673
## ses 3.596455 4.152745
添加2級水平的固定效應(yīng)前面我們添加了ses作為1級水平的自變量，以解釋學(xué)生在數(shù)學(xué)成績方面的部分差異。下面我們再添加一個在水平2單位（不同的學(xué)校）間變化的自變量，該自變量在不同的學(xué)校間是不同的，但是在同一所學(xué)校內(nèi)的所有學(xué)生中是相同的，符合條件的自變量有3個： ses_mean:以學(xué)校為單位衡量學(xué)生社會經(jīng)濟(jì)地位的均值pro4yrc:每個學(xué)校中計(jì)劃就讀四年制大學(xué)的學(xué)生比例public:學(xué)校類型，1=公立學(xué)校，0=私立學(xué)校
我們選擇public作為水平2單位的固定效應(yīng)，這個模型依然是一個隨機(jī)截距模型，或者叫方差成分模型： ses_l1_public_l2 <- lmer(math ~ ses + public + (1|schcode),
data = data, REML = TRUE)
summary(ses_l1_public_l2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + public + (1 | schcode)
## Data: data
##
## REML criterion at convergence: 48216
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7718 -0.5541 0.1309 0.6477 5.6916
##
## Random effects:
## Groups Name Variance Std.Dev.
## schcode (Intercept) 3.486 1.867
## Residual 62.807 7.925
## Number of obs: 6871, groups: schcode, 419
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.63143 0.25535 381.81733 225.693 <2e-16 ***
## ses 3.87338 0.13673 3928.37427 28.329 <2e-16 ***
## public -0.04859 0.29862 385.93649 -0.163 0.871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ses
## ses 0.013
## public -0.854 -0.031
先看固定效應(yīng)：Intercept的估計(jì)值為57.5960，表示當(dāng)所有自變量都為0的時候，因變量的預(yù)測值；ses的估計(jì)值為3.87338，表示ses每增加一個單位，數(shù)學(xué)成績可以提高3.87338分；public的估計(jì)值為-0.04859，說明相對于私立學(xué)校，公立學(xué)校在數(shù)學(xué)成績上平均降低0.04859分。固定效應(yīng)的解讀和普通的多元線性回歸沒有差別。 再看隨機(jī)效應(yīng)：結(jié)果解讀和上面一個模型的解讀類似的，就不再重復(fù)了。schcode的方差變大了一些，殘差的方差沒有變化。 理論上如果新加入的自變量能夠解釋更多的因變量變異，那么殘差的變異（方差）通常會減少，群體間（在本例中也就是學(xué)校間）的變異也會減少，因?yàn)檫@部分變異都被新加入的自變量解釋了。但是很明顯我們新加入的這個public自變量不太行，它幾乎解釋不了因變量的變異，從它的系數(shù)也可以看出來，只有-0.04859，比ses差遠(yuǎn)了。說明公立學(xué)校和私立學(xué)校對數(shù)學(xué)成績影響很?。≒值也表明這個變量沒有統(tǒng)計(jì)學(xué)意義）。 我們也可以通過計(jì)算模型的一些指標(biāo)看看這個自變量到底行不行，并且和前面的模型比較一下： compare_performance(null_model,ses_l1,ses_l1_public_l2)
## # Comparison of Model Performance Indices
##
## Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ------------------------------------------------------------------------------------------------------------------------------------------
## null_model | lmerModLmerTest | 48881.8 (<.001) | 48881.8 (<.001) | 48902.3 (<.001) | 0.138 | 0.000 | 0.138 | 7.977 | 8.158
## ses_l1 | lmerModLmerTest | 48219.1 (0.729) | 48219.1 (0.729) | 48246.4 (0.988) | 0.167 | 0.121 | 0.052 | 7.810 | 7.925
## ses_l1_public_l2 | lmerModLmerTest | 48221.1 (0.271) | 48221.1 (0.271) | 48255.2 (0.012) | 0.167 | 0.121 | 0.053 | 7.810 | 7.925
結(jié)果發(fā)現(xiàn)，加入public后，AIC、BIC還變大了，R2沒啥太大的變化，充分說明這個變量真的作用不大，可以不加。 具有隨機(jī)斜率的MLM前面我們選擇了10個學(xué)校，以展示了不同學(xué)校間數(shù)學(xué)成績（math）與社會經(jīng)濟(jì)狀況（ses）之間的關(guān)系：  從圖中可以看出，不同學(xué)校ses的截距和斜率值差異很大。例如，學(xué)校3的截距約為38，斜率較小且為正值，而學(xué)校8的截距約為55，斜率較大且為正值。 在ses_l1這個模型中，我們假定每個學(xué)校ses的斜率是相同的，只估計(jì)了隨機(jī)截距的變異，但是從圖中可以看出，其實(shí)每個學(xué)校ses的斜率都是不一樣的。也就是說，我們只估計(jì)了ses的平均效應(yīng)，卻忽略了不同學(xué)校的ses是不同的。 下面我們在模型中添加一個隨機(jī)斜率項(xiàng)來模擬學(xué)校間ses斜率的這種變異。 在lme4的語法中，只需要將想要估計(jì)的具有不同斜率的變量名字放在|前面即可： # 估計(jì)ses的隨機(jī)斜率和隨機(jī)截距
ses_l1_random <- lmer(math ~ ses + (ses|schcode),
data = data, REML = TRUE)
## boundary (singular) fit: see help('isSingular')
summary(ses_l1_random)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + (ses | schcode)
## Data: data
##
## REML criterion at convergence: 48190.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8578 -0.5553 0.1290 0.6437 5.7098
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schcode (Intercept) 3.2042 1.7900
## ses 0.7794 0.8828 -1.00
## Residual 62.5855 7.9111
## Number of obs: 6871, groups: schcode, 419
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.6959 0.1315 378.6378 438.78 <2e-16 ***
## ses 3.9602 0.1408 1450.7730 28.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ses -0.284
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
上面的公式也可以寫成math~ses+(1+ses|schcode)，1表示隨機(jī)截距項(xiàng)，這是默認(rèn)設(shè)置，也可以僅使用(ses|schcode)估計(jì)隨機(jī)截距和隨機(jī)斜率。如果你想從模型中排除隨機(jī)截距，需要寫成(0+ses|schcode)來覆蓋默認(rèn)設(shè)置。像這種既有隨機(jī)截距又有隨機(jī)斜率的模型又被稱為隨機(jī)系數(shù)模型。 先看固定效應(yīng)：Intercept的估計(jì)值為57.6959，表示當(dāng)所有自變量都為0的時候，因變量的預(yù)測值；ses的估計(jì)值為3.9602，表示ses每增加一個單位，數(shù)學(xué)成績可以提高3.9602分。固定效應(yīng)的解讀和普通的多元線性回歸沒有差別。 再看隨機(jī)效應(yīng)：schcode（Intercept）的方差為3.2042，標(biāo)準(zhǔn)差是1.7900，它衡量的是不同學(xué)校之間截距的變異。中間的ses的方差為0.7794，標(biāo)準(zhǔn)差為0.8828，它衡量的不同學(xué)校之間斜率的變異，意思是不同學(xué)校的ses的斜率圍繞總體平均斜率變化的方差為0.7794。殘差的方差為62.5855，標(biāo)準(zhǔn)差為7.9111，它衡量的是模型無法解釋的變異，可以發(fā)現(xiàn)在考慮了隨機(jī)斜率后，殘差的方差又變小了（62.807到62.5855）。Corr是-1表示隨機(jī)截距和隨機(jī)斜率的相關(guān)系數(shù)是-1，說明有些學(xué)校的截距越高，斜率就越低（結(jié)合上面的圖看，是不是有這種趨勢），也就是說：隨著平均數(shù)學(xué)成績的增加，ses與數(shù)學(xué)成績之間的關(guān)系降低。 如果想查看每個學(xué)校的截距和斜率，可以使用ranef： head(ranef(ses_l1_random))
## $schcode
## (Intercept) ses
## 1 0.9746642943 -0.4806908392
## 2 1.0450460989 -0.5154021638
## 3 -3.4842479301 1.7183824946
## 4 1.8810910566 -0.9277278791
## 省略......
## 418 -0.6233121254 0.3074088487
## 419 -1.2366360507 0.6098916565
具有交互效應(yīng)的MLM前面分別探索了ses和public對數(shù)學(xué)成績的影響，如果要探索它們的交互效應(yīng)，只需要像普通的單水平模型一樣，將交互項(xiàng)添加到公式中即可，比如將ses和public的交互項(xiàng)添加到模型中： # 添加交互效應(yīng)
crosslevel_model <- lmer(math ~ ses + public + ses:public + (ses|schcode),
data = data, REML = TRUE)
## boundary (singular) fit: see help('isSingular')
# 也可以寫成math ~ ses*public + (ses|schcode)
summary(crosslevel_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + public + ses:public + (ses | schcode)
## Data: data
##
## REML criterion at convergence: 48187.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8509 -0.5593 0.1294 0.6412 5.6998
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schcode (Intercept) 3.2144 1.7929
## ses 0.8013 0.8951 -1.00
## Residual 62.5555 7.9092
## Number of obs: 6871, groups: schcode, 419
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.72440 0.25183 382.39815 229.216 <2e-16 ***
## ses 4.42383 0.27427 1283.55623 16.130 <2e-16 ***
## public -0.02632 0.29472 387.41741 -0.089 0.9289
## ses:public -0.62520 0.31957 1363.95274 -1.956 0.0506 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ses public
## ses -0.232
## public -0.852 0.197
## ses:public 0.198 -0.858 -0.250
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
先看固定效應(yīng)（由于public是分類變量，1表示公立學(xué)校，0表示私立學(xué)校，所以在進(jìn)行回歸分析時會自動進(jìn)行啞變量編碼，以0（也就是私立學(xué)校）為參考，這里涉及一個基礎(chǔ)知識，即回歸分析中的啞變量編碼）： - 截距（Intercept）是57.72440，即：當(dāng)學(xué)校為是私立學(xué)校（public=0）且
ses也為0時的平均數(shù)學(xué)預(yù)期成績； ses的估計(jì)值為4.42383，即：對于私立學(xué)校（public=0）來說，ses每增加一個單位，數(shù)學(xué)成績會增加4.42383分；public的估計(jì)值為-0.02632，即：公立學(xué)校（public=1）相比于私立學(xué)校（public=0），平均數(shù)學(xué)成績會減少0.02632分；ses:public的估計(jì)值為-0.62520，即：ses對數(shù)學(xué)成績的影響在公立學(xué)校（public=1）比在私立學(xué)校（public=0）平均減少0.62520分。- 借助這些系數(shù)，我們可以估算
ses在公立學(xué)校（public=1）的預(yù)期斜率，即4.42383-0.62520=3.79863。因此公立學(xué)校中ses每增加一個單位，數(shù)學(xué)成績會增加3.79863分，略小于私立學(xué)校（4.42383分）。
再看隨機(jī)效應(yīng)：（和上面ses_l1_random的結(jié)果解讀基本類似，這里簡單說一下）。schcode（Intercept）的方差為3.2144，標(biāo)準(zhǔn)差是1.7929，它衡量的是不同學(xué)校之間截距的變異。中間的ses的方差為0.8013，標(biāo)準(zhǔn)差為0.8951，它衡量的不同學(xué)校之間斜率的變異，意思是不同學(xué)校的ses的斜率圍繞總體平均斜率變化的方差為0.8013。殘差的方差為62.5555，標(biāo)準(zhǔn)差為7.9092，它衡量的是模型無法解釋的變異。Corr是-1表示隨機(jī)截距和隨機(jī)斜率的相關(guān)系數(shù)是-1。 重復(fù)測量數(shù)據(jù)的MLM重復(fù)測量數(shù)據(jù)的結(jié)構(gòu)非常適合多水平模型，測量的時間點(diǎn)可以看做是1級水平，每個患者可以看做是2級水平，每個患者都包括多次測量數(shù)據(jù)，這是一個具有兩個層次的結(jié)構(gòu)。 下面用一個簡單的例子進(jìn)行演示。 使用某藥治療10個高血壓患者，分別測量每個患者治療前后的血壓，請對該數(shù)據(jù)進(jìn)行分析。 # 模擬數(shù)據(jù)
data12_1 <- data.frame(id=c(1:10,1:10),
stat=rep(c("治療前","治療后"),each=10),
bp=c(130,124,136,128,122,118,116,138,126,124,
114,110,126,116,102,100,98,122,108,106)
)
# 設(shè)置下因子水平，變成治療后vs治療前
data12_1$stat <- factor(data12_1$stat,levels = c("治療前","治療后"))
data12_1$id <- factor(data12_1$id)
str(data12_1)
## 'data.frame': 20 obs. of 3 variables:
## $ id : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ stat: Factor w/ 2 levels "治療前","治療后": 1 1 1 1 1 1 1 1 1 1 ...
## $ bp : num 130 124 136 128 122 118 116 138 126 124 ...
data12_1 # 數(shù)據(jù)長這樣
## id stat bp
## 1 1 治療前 130
## 2 2 治療前 124
## 3 3 治療前 136
## 4 4 治療前 128
## 5 5 治療前 122
## 6 6 治療前 118
## 7 7 治療前 116
## 8 8 治療前 138
## 9 9 治療前 126
## 10 10 治療前 124
## 11 1 治療后 114
## 12 2 治療后 110
## 13 3 治療后 126
## 14 4 治療后 116
## 15 5 治療后 102
## 16 6 治療后 100
## 17 7 治療后 98
## 18 8 治療后 122
## 19 9 治療后 108
## 20 10 治療后 106
先畫個圖看看不同患者治療前后的血壓值?？梢园l(fā)現(xiàn)每個患者的斜率和截距都是不一樣的： library(ggplot2)
ggplot(data12_1, aes(stat,bp))+
geom_line(aes(color=id,group = id))
 下面就是建立多水平模型即可，這里為了省事，我們默認(rèn)斜率是相等的，只考慮隨機(jī)截距： library(lmerTest)
f <- lmer(bp ~ stat + (1 | id), data = data12_1)
summary(f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: bp ~ stat + (1 | id)
## Data: data12_1
##
## REML criterion at convergence: 113.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1422 -0.5311 0.1307 0.4070 1.5714
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 63.511 7.969
## Residual 4.889 2.211
## Number of obs: 20, groups: id, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 126.2000 2.6153 9.6662 48.25 7.56e-13 ***
## stat治療后 -16.0000 0.9888 9.0000 -16.18 5.83e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
## (Intr)
## stat治療后 -0.189
結(jié)果解讀不再重復(fù)。 廣義混合效應(yīng)模型上面介紹的例子，都是假定因變量為連續(xù)分布，而在醫(yī)學(xué)和公共衛(wèi)生領(lǐng)域，許多應(yīng)變量是離散型的，例如個體的健康狀態(tài)可能與吸煙、飲酒、鍛煉等日常生活方式有關(guān)。在離散型應(yīng)變量的情形下，若數(shù)據(jù)具有層次結(jié)構(gòu)特征，則最低水平的觀察單位發(fā)生某事件的概率并不完全相互獨(dú)立，故不再服從二項(xiàng)分布或Poisson分布，而服從超二項(xiàng)（extra-binomial）分布或超Poisson（extra-Poisson）分布。因此，在擬合這種類型的模型時，結(jié)局的聚集效應(yīng)和離散型誤差的復(fù)雜分布應(yīng)考慮在模型中。假定在某試驗(yàn)中對某事件的測量為發(fā)生與不發(fā)生，即二分類的資料，若將其作為因變量，則在多水平框架內(nèi)，處理這類資料的統(tǒng)計(jì)模型一般稱為多水平廣義線性模型。 使用方法和結(jié)果解讀基本類似，以后遇到了再詳細(xì)介紹。 參考資料- Introduction to Multilevel Modelling
- A Cheatsheet for Building Multilevel Models in R
- CRAN Task View: Mixed, Multilevel, and Hierarchical Models in R
|
