拓端tecdat|R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調查資料

原文連結：http://tecdat.cn/?p=22813

原文出處：拓端資料部落公衆号

本教程為讀者提供了使用 頻率學派的廣義線性模型（GLM）的基本介紹。具體來說，本教程重點介紹邏輯回歸在二進制結果和計數/比例結果情況下的使用，以及模型評估的方法。本教程使用教育資料例子進行模型的應用。此外，本教程還簡要示範了用R對GLM模型進行的多層次擴充。最後，還讨論了GLM架構中的更多分布和連結函數。

本教程包含以下結構。

1. 準備工作。

2. 介紹GLM。

3. 加載教育資料。

4. 資料準備。

5. 二進制（伯努利）Logistic回歸。

6. 二項式 Logistic 回歸。

7. 多層次Logistic回歸。

8. 其他族和連結函數。

本教程介紹了：

- 假設檢驗和統計推斷的基本知識。

- 回歸的基本知識。

- R語言編碼的基本知識。

- 進行繪圖和資料處理的基本知識。

廣義線性模型（GLM）簡介

對于y是連續值得情況，我們可以用這種方式處理，但當y是離散值我們用普通線性模型就不合适了，這時我們引用另外一種模型 --- Generalised Linear Models 廣義線性模型。

為了擷取GLM模型，我們列出3個條件：

1. 

，也就是y|x為指數族分布，指數族分布形式：

2. 如果我們判斷y的假設為 

，則

。

3. 自然參數和輸入x呈線性關系：

這3個條件的來由我們不讨論，我們隻知道做這樣的假設是基于“設計”的選擇，而非必然。

我們以泊松回歸為例, y服從泊松分布 

，化為指數族形式，我們可以得到

。是以

之後即為最大似然法的過程。

教育資料

本教程中使用的資料是教育資料。

該資料來源于全國性的國小教育調查。資料中的每一行都是指一個學生。結果變量留級是一個二分變量，表示一個學生在國小教育期間是否留過級。學校變量表示一個學生所在的學校。個人層面的預測因素包括。  性别（0=女性，1=男性）和學前教育（受過學前教育，0=沒有，1=有）。學校層面是學校平均SES（社會經濟地位）得分。

本教程利用教育資料試圖回答的主要研究問題是。

忽略資料的結構，性别和學前教育對學生是否留級的影響是什麼？

忽略資料的結構，學校平均SES對學生留級比例的影響是什麼？

考慮到資料的結構，性别、學前教育和學校平均SES對學生是否留級有什麼影響？

這三個問題分别用以下這些模型來回答：二進制邏輯回歸；二項邏輯回歸；多層次二進制邏輯回歸。

資料準備

加載必要的軟體包

# 如果你還沒有安裝這些包，請使用install.packages("package_name")指令。
library(lme4) # 用于多層次模型
library(tidyverse) # 用于資料處理和繪圖
           

導入資料

head(Edu)
           

資料處理

mutate(學校 = factor(學校),
         性别 = if_else(性别 == 0, "girl", "boy"),
         性别 = factor(性别, levels = c("girl", "boy")),
         受過學前教育 = if_else(受過學前教育 == 0, "no", "yes"),
         受過學前教育 = factor(受過學前教育, levels = c("no", "yes")))
           

檢查缺失的資料

summarise_each((~sum(is.na(.))
           

資料中，經濟地位變量有1066個觀測值缺失。對缺失資料的處理本身就是一個複雜的話題。為了友善起見，我們在本教程中簡單地将資料缺失的案例删除。

二進制邏輯回歸

探索資料：按性别和學前教育分類的留級數量 

group_by(性别) %>%
  summarise(是否留過級 = sum(是否留過級))
           

看來，留級的學生人數在男女之間有很大的不同，更多的男學生留級。更多沒有接受過學前教育的學生留級。這一觀察結果表明，性别和學前教育可能對留級有預測作用。

建構二進制邏輯回歸模型

R預設安裝了基礎包，其中包括運作GLM的glm函數。glm的參數與lm的參數相似：公式和資料。然而，glm需要一個額外的參數：family，它指定了結果變量的假設分布；在family中我們還需要指定連結函數。family的預設值是gaussian(link = "identity")，這導緻了一個線性模型，相當于由lm指定的模型。在二進制邏輯回歸的情況下，glm要求我們指定一個帶有logit連結的二項分布，即family = binomial(link = "logit") 。

glm(formula ,
                    family = binomial(link = "logit"))
           

解釋

從上面的總結輸出中，我們可以看到，性别對學生留級的機率有正向和顯著的預測，而學前教育則有負向和顯著的預測。具體來說，與女孩相比，男孩更有可能留級。以前上過學的學生不太可能導緻留級。

為了解釋參數估計值，我們需要對估計值進行指數化處理。

請注意，參數估計的解釋與幾率而不是機率有關。賠率的定義是。P（事件發生）/P（事件未發生）。在本分析中，假設其他一切保持不變，與女孩相比，男孩增加了54%的留級幾率；與沒有學前教育相比，假設其他一切保持不變，擁有學前教育降低了（1-0.54）%=46%的留級幾率。

參數效應的可視化

為了使參數效應的解釋更加容易，我們可以對參數效應可視化。

plot(Effects)
           

請注意，在這兩張圖中，Y刻度指的是留級的機率，而不是幾率。機率比幾率更容易解釋。每個變量的機率分數是通過假設模型中的其他變量是常數并采取其平均值來計算的。正如我們所看到的，假設一個學生有平均的學前教育，作為一個男孩比作為一個女孩有更高的留級機率（~0.16）~0.11）。同樣，假設一個學生有一個平均的性别，有學前教育的學生比沒有學前教育的學生留級的機率低（~0.11）（~0.18）。請注意，在這兩幅圖中，還包括了估計值的置信區間，以使我們對估計值的不确定性有一些了解。

請注意，平均學前教育和性别的概念可能聽起來很奇怪，因為它們是分類變量（即因素）。如果你對假設一個平均因素的想法感到奇怪，你可以指定你的預期因素水準作為參考點。

predictors = list( values=c(性别boy=0, 受過學前教育yes = 0))
           

設定性别boy = 0意味着在學前教育效應圖中，性别變量的參考水準被設定為0；學前教育yes = 0導緻0成為性别效應圖中學前教育變量的參考水準。

是以，正如上面兩幅圖所示，假設學生沒有接受過學前教育，作為男孩的留級機率（~0.20）比作為女孩的留級機率（~0.14）要高；假設學生是女性，有學前教育的留級機率（~0.09）比沒有學前教育的留級機率（~0.15）要低。

模型評估:拟合度

評價邏輯回歸模型的拟合度有不同的方法。

似然比檢驗

如果一個邏輯回歸模型與預測因子較少的模型相比，顯示出拟合度的提高，則該模型對資料有較好的拟合度。這是用似然比檢驗進行的，它将完整模型下資料的似然性與較少預測因素的模型下資料的似然性進行比較。從一個模型中删除預測變量幾乎總是會使模型的拟合度降低（即模型的對數似然率較低），但測試觀察到的模型拟合度差異是否具有統計學意義是很有用的。

#指定一個隻有`性别'變量的模型
#使用`anova()`函數來運作似然比測試
anova(ModelTest, Model, test ="Chisq")
           

我們可以看到，同時包含性别和學前教育的預測因子的模型比隻包含性别變量的模型對資料的拟合效果要好得多。請注意，這種方法也可以用來确定是否有必要包括一個或一組變量。

 AIC

Akaike資訊準則（AIC）是另一個模型選擇的衡量标準。與似然比檢驗不同，AIC的計算不僅要考慮模型的拟合度，還要考慮模型的簡單性。通過這種方式，AIC處理了模型的拟合度和複雜性之間的權衡，是以，不鼓勵過度拟合。較小的AIC是首選。

在AIC值較小的情況下，同時具有性别和學前教育預測因子的模型優于隻具有性别預測因子的模型。

正确分類率

正确分類率是另一個有用的衡量标準，可以看出模型對資料的合适程度。

#使用`predict()`函數，從拟合的模型中計算出原始資料中學生的預測機率
Pred <- if_else(Pred > 0.5, 1, 0)
ConfusionMatrix <- table(Pred, TRUE)
#正确的分類率
           

我們可以看到，該模型對所有觀測值的85.8%進行了正确分類。然而，仔細觀察可以發現，模型預測所有的觀察值都屬于 "0 "類，也就是說，所有的學生都被預測為不留級。考慮到留級變量的多數類别是0（不），該模型在分類上的表現并不比簡單地将所有觀測值配置設定到多數類别0（不）更好。

AUC（曲線下面積）

使用正确分類率的一個替代方法是曲線下面積（AUC）測量。AUC測量區分度，即測試對有目标反應和無目标反應的人進行正确分類的能力。在目前的資料中，目标變量是留級。我們從 "留級 "組和 "不留級 "組中随機抽取一名學生。預測機率較高的學生應該是 "留級 "組中的學生。AUC是随機抽出的對子的百分比。這個程式将AUC與正确分類率區分開來，因為AUC不依賴于結果變量中類的比例的變化。0.50的值意味着該模型的分類效果不比随機好。一個好的模型應該有一個遠遠高于0.50的AUC分數（最好高于0.80）。

# 計算用該模型預測類别的AUC

AUC <- performance(Pred, measure = "auc")
AUC <- [email protected][[1]]
AUC
           

AUC分數為0.60，該模型的判别能力不強。

二項式 Logistic 回歸

正如開頭提到的，邏輯回歸也可以用來為計數或比例資料模組化。二項邏輯回歸假設結果變量來自伯努利分布（這是二項分布的一個特例），其中試驗次數n為1，是以結果變量隻能是1或0。 相反，二項邏輯回歸假設目标事件的數量遵循二項分布，試驗次數n，機率q。這樣一來，二項邏輯回歸允許結果變量取任何非負整數值，是以能夠處理計數資料。

教育資料記錄了集中在學校内的個别學生的資訊。通過彙總各學校留級的學生人數，我們得到一個新的資料集，其中每一行代表一所學校，并有關于該學校留級學生的比例資訊。學校平均社會經濟地位（平均SES分數）也是在學校層面上的；是以，它可以用來預測在某個學校留級的學生的比例或數量。

轉換資料

在這個新的資料集中，留級指的是留級的學生人數；TOTAL指的是某所學校的學生總數。

探索資料

ggplot(aes(x , y)) +
  geom_smooth(method = "lm")
           

我們可以看到，留級的學生比例與學校平均社會經濟地位的反對數呈負相關。請注意，我們将變量學校平均社會經濟地位模組化為其反對數，因為在二項式回歸模型中，我們假設線性預測因子的反對數與結果（即事件比例）之間存線上性關系，而不是預測因子本身與結果之間存線上性關系。

拟合二項式Logistic回歸模型

為了拟合二項式邏輯回歸模型，我們也使用glm函數。唯一的差別是在公式中對結果變量的說明。我們需要指定目标事件的數量（留級）和非事件的數量（TOTAL-留級），并将它們包在cbind()中。

glm(cbind(是否留過級, TOTAL-是否留過級) ~ 學校平均社會經濟地位,
                  family = binomial(logit))
           

解釋

二項式回歸模型的參數解釋與二項式邏輯回歸模型相同。從上面的模型總結中我們知道，一所學校的平均SES分數與該校學生留級的幾率呈負相關。為了提高可解釋性，我們再次使用summ()函數來計算學校平均社會經濟地位的指數化系數估計。由于學校平均社會經濟地位是一個連續的變量，我們可以将指數化的學校平均社會經濟地位估計值标準化（通過将原始估計值與變量的SD相乘，然後将所得數字指數化）。

#注意，為了對二項回歸模型使用summ()函數，我們需要将結果變量作為對象。
是否留過級 <- (filter(edu, !is.na(學校平均社會經濟地位)), 是否留過級)
           

我們可以看到，随着學校平均社會經濟地位的SD增加，學生留級的幾率降低了1 - 85% = 15%。

我們可以直覺地看到學校平均社會經濟地位的效果。

plot(allEffects)
           

上面的圖表顯示了學校平均社會經濟地位對學生留級機率的預期影響。在其他因素不變的情況下，随着學校平均社會經濟地位的增加，一個學生留級的機率會降低（從0.19到0.10）。藍色陰影區域表示每個學校平均社會經濟地位值的預測值的95%置信區間。

多層次二進制邏輯回歸

前面介紹的二進制邏輯回歸模型僅限于對學生層面的預測因素的影響進行模組化；二進制邏輯回歸僅限于對學校層面的預測因素的影響進行模組化。為了同時納入學生層面和學校層面的預測因素，我們可以使用多層次模型，特别是多層次二進制邏輯回歸。

除了上述動機外，還有更多使用多層次模型的理由。例如，由于資料是在學校内分類的，來自同一學校的學生很可能比來自其他學校的學生更相似。正因為如此，在一所學校，一個學生留級的機率可能很高，而在另一所學校，則很低。此外，即使是結果（即留級）和預測變量（如性别、學前教育、學校平均社會經濟地位）之間的關系，在不同的學校也可能不同。還要注意的是，學校平均社會經濟地位變量中存在缺失值。使用多層次模型可以較好地解決這些問題。

請看下面的圖作為例子。該圖顯示了各學校留級學生的比例。我們可以看到不同學校之間的巨大差異。是以，我們可能需要多層次模型。

group_by(學校) %>%
  summarise(PROP = sum(是否留過級)/n()) %>%
  plot()
           

我們還可以通過學校來繪制性别和留級之間的關系，以了解性别和留級之間的關系是否因學校而異。

mutate(性别 = if_else(性别 == "boy", 1, 0)) %>%
  ggplot(aes(x = 性别, y = 是否留過級, color = as.factor(學校))) +
           

在上面的圖中，不同的顔色代表不同的學校。我們可以看到，不同學校的性别和留級之間的關系似乎有很大不同。

我們可以為學前教育和留級做同樣的圖。

mutate(性别 = if_else(性别 == "girl", 0, 1),
         受過學前教育 = if_else(受過學前教育 == "yes", 1, 0)) %>%
  group_by(學校) %>%
  mutate(性别 = 性别 - mean(性别),
           

學前教育和留級之間的關系在不同的學校也顯得相當不同。然而，我們也可以看到，大多數的關系都呈下降趨勢，從0（以前沒有上過學）到1（以前上過學），表明學前教育和留級之間的關系為負。

由于上述觀察結果，我們可以得出結論，在目前的資料中需要建立多層次的模型，不僅要有随機截距（學校），還可能要有性别和學前教育的随機斜率。

中心化變量

在拟合多層次模型之前，有必要采用适當的中心化方法（即均值中心化）對預測變量進行中心化，因為中心化方法對模型估計的解釋很重要。根據Enders和Tofighi（2007）的建議，我們應該對第一層次的預測因子性别和學前教育使用中心化，對第二層次的預測因子學校平均社會經濟地位使用均值中心化。

受過學前教育 = if_else(受過學前教育 == "yes", 1, 0)) %>%
  group_by(學校) %>%
  mutate(性别 = 性别 - mean(性别),
         受過學前教育 = 受過學前教育 - mean(受過學前教育)) %>%
  ungroup() %>%
           

隻有截距模型

為了指定一個多層次模型，我們使用lme4軟體包。随機斜率項和聚類項應該用|分隔。注意，我們使用了一個額外的參數指定比預設值（10000）更大的最大疊代次數。因為一個多層次模型可能需要大量的疊代來收斂。

我們首先指定一個純截距模型，以評估資料聚類結構的影響。

glmer(是否留過級 ~ 1 + (1|學校),
                             optCtrl = list(maxfun=2e5))
           

下面我們計算一下純截距模型的ICC（類内相關）。

0.33的ICC意味着結果變量的33%的變化可以被資料的聚類結構所解釋。這提供了證據表明，與非多層次模型相比，多層次模型可能會對模型的估計産生影響。是以，多層次模型的使用是必要的，也是有保證的。

完整模型

按部就班地建立一個多層次模型是很好的做法。然而，由于本文的重點不是多層次模型，我們直接從純截距模型到我們最終感興趣的全模型。在完整模型中，我們不僅包括性别、學前教育和學校平均社會經濟地位的固定效應項和一個随機截距項，還包括性别和學前教育的随機斜率項。請注意，我們指定 family = binomial(link = "logit")，因為這個模型本質上是一個二進制邏輯回歸模型。

glmer(是否留過級 ~ 性别 + 受過學前教育 + 學校平均社會經濟地位 + (1 + 性别 + 受過學前教育|學校)
           

結果（與固定效應有關）與之前二進制邏輯回歸和二項邏輯回歸模型的結果相似。在學生層面上，性别對學生留級的幾率有顯著的正向影響，而學前教育有顯著的負向影響。在學校層面上，學校地位對結果變量有顯著的負向影響。我們也來看看随機效應項的方差。

同樣，我們可以使用summ()函數來檢索指數化的系數估計值，便于解釋。

sum(Model_Full)
           

我們還可以顯示參數估計的效果。請注意，由于第一級分類變量（性别和學前教育）是中心化的，是以在模型中它們被當作連續變量，在下面的效果圖中也是如此。

plot((Model)
           

除了固定效應項之外，我們也來看看随機效應項。從之前的ICC值來看，我們知道有必要包括一個随機截距。但是，包括性别和學前教育的随機斜率的必要性就不太清楚了。為了弄清楚這一點，我們可以用似然比檢驗和AIC來判斷随機斜率的加入是否能改善模型的拟合。

glmer(是否留過級 ~ 性别 + 受過學前教育 + 學校平均社會經濟地位 + (1 + 受過學前教育|學校),
           

#拟合一個不完整的模型，剔除`受過學前教育'的随機斜率項
glmer(是否留過級 ~ 性别 + 受過學前教育 + 學校平均社會經濟地位 + (1 + 性别|學校),

           

似然比檢驗

比較完整的模型和排除了`性别'的模型 

将完整的模型與排除了 "受過學前教育 "的模型進行比較 

從所有不顯著的似然比檢驗結果（Pr（>Chisq）>0.05），我們可以得出結論，增加任何随機斜率項對模型拟合都沒有明顯的改善。

AIC

AIC #full模型
AIC#＃沒有性别的模型
AIC #＃沒有受過學前教育的模型
AIC#＃沒有随機斜率的模型
           

從AIC的結果來看，我們發現包括随機斜率項要麼沒有大幅提高AIC（用較低的AIC值表示），要麼導緻更差的AIC（即更高）。是以，我們也得出結論，沒有必要包括随機效應項。

其他族（分布）和連結函數

到目前為止，我們已經介紹了二進制和二項邏輯回歸，這兩種回歸都來自于二項家族的logit連結。然而，還有許多分布族和連結函數，我們可以在glm分析中使用。例如，為了對二進制結果進行模組化，我們還可以使用probit連結或log-log（cloglog）來代替logit連結。為了給計數資料模組化，我們也可以使用泊松回歸，它假設結果變量來自泊松分布，并使用對數作為連結函數。

參考文獻

Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01

Enders, C. K., & Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: A new look at an old issue. Psychological Methods, 12(2), 121-138. doi:10.1037/1082-989X.12.2.121

最受歡迎的見解

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