全文鏈接:http://tecdat.cn/?p=22813

最近我們被客戶要求撰寫關(guān)于混合效應廣義線性模型的研究報告，包括一些圖形和統(tǒng)計輸出。本教程為讀者提供了使用頻率學派的廣義線性模型（GLM）的基本介紹。具體來說，本教程重點介紹邏輯回歸在二元結(jié)果和計數(shù)/比例結(jié)果情況下的使用，以及模型評估的方法

本教程使用教育數(shù)據(jù)例子進行模型的應用。此外，本教程還簡要演示了用R對GLM模型進行的多層次擴展。最后，還討論了GLM框架中的更多分布和鏈接函數(shù)。

本教程包含以下結(jié)構(gòu)。

1. 準備工作。

2. 介紹GLM。

3. 加載教育數(shù)據(jù)。

4. 數(shù)據(jù)準備。

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

6. 二項式 Logistic 回歸。

7. 多層次Logistic回歸。

8. 其他族和鏈接函數(shù)。

本教程介紹了：

• 假設檢驗和統(tǒng)計推斷的基本知識。

• 回歸的基本知識。

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

• 進行繪圖和數(shù)據(jù)處理的基本知識。

廣義線性模型（GLM）簡介

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

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

1. ，也就是y|x為指數(shù)族分布，指數(shù)族分布形式：

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

1. ，則

1. 。

2. 自然參數(shù)和輸入x呈線性關(guān)系：

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

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

，化為指數(shù)族形式，我們可以得到

。所以

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

教育數(shù)據(jù)

本教程中使用的數(shù)據(jù)是教育數(shù)據(jù)。

該數(shù)據(jù)來源于全國性的小學教育調(diào)查。數(shù)據(jù)中的每一行都是指一個學生。結(jié)果變量留級是一個二分變量，表示一個學生在小學教育期間是否留過級。學校變量表示一個學生所在的學校。個人層面的預測因素包括。??性別（0=女性，1=男性）和學前教育（受過學前教育，0=沒有，1=有）。學校層面是學校平均SES（社會經(jīng)濟地位）得分。

本教程利用教育數(shù)據(jù)試圖回答的主要研究問題是。

忽略數(shù)據(jù)的結(jié)構(gòu)，性別和學前教育對學生是否留級的影響是什么？
忽略數(shù)據(jù)的結(jié)構(gòu)，學校平均SES對學生留級比例的影響是什么？
考慮到數(shù)據(jù)的結(jié)構(gòu)，性別、學前教育和學校平均SES對學生是否留級有什么影響？
這三個問題分別用以下這些模型來回答：二元邏輯回歸；二項邏輯回歸；多層次二元邏輯回歸。

數(shù)據(jù)準備

加載必要的軟件包

#?如果你還沒有安裝這些包，請使用install.packages("package_name")命令。library(lme4)?#?用于多層次模型library(tidyverse)?#?用于數(shù)據(jù)處理和繪圖

導入數(shù)據(jù)

head(Edu)

數(shù)據(jù)處理

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

檢查缺失的數(shù)據(jù)

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

數(shù)據(jù)中，經(jīng)濟地位變量有1066個觀測值缺失。對缺失數(shù)據(jù)的處理本身就是一個復雜的話題。為了方便起見，我們在本教程中簡單地將數(shù)據(jù)缺失的案例刪除。

二元邏輯回歸

探索數(shù)據(jù)：按性別和學前教育分類的留級數(shù)量?

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

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

構(gòu)建二元邏輯回歸模型

R默認安裝了基礎包，其中包括運行GLM的glm函數(shù)。glm的參數(shù)與lm的參數(shù)相似：公式和數(shù)據(jù)。然而，glm需要一個額外的參數(shù)：family，它指定了結(jié)果變量的假設分布；在family中我們還需要指定鏈接函數(shù)。family的默認值是gaussian(link = "identity")，這導致了一個線性模型，相當于由lm指定的模型。在二元邏輯回歸的情況下，glm要求我們指定一個帶有l(wèi)ogit鏈接的二項分布，即family = binomial(link = "logit") 。

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

解釋

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

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

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

參數(shù)效應的可視化

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

plot(Effects)

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

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

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

設置性別boy = 0意味著在學前教育效應圖中，性別變量的參考水平被設置為0；學前教育yes = 0導致0成為性別效應圖中學前教育變量的參考水平。

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

點擊標題查閱往期內(nèi)容

多水平模型、分層線性模型HLM、混合效應模型研究教師的受歡迎程度

左右滑動查看更多

01

02

03

04

模型評估:擬合度

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

似然比檢驗

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

#指定一個只有`性別'變量的模型#使用`anova()`函數(shù)來運行似然比測試anova(ModelTest,?Model,?test?="Chisq")

我們可以看到，同時包含性別和學前教育的預測因子的模型比只包含性別變量的模型對數(shù)據(jù)的擬合效果要好得多。請注意，這種方法也可以用來確定是否有必要包括一個或一組變量。

?AIC

Akaike信息準則（AIC）是另一個模型選擇的衡量標準。與似然比檢驗不同，AIC的計算不僅要考慮模型的擬合度，還要考慮模型的簡單性。通過這種方式，AIC處理了模型的擬合度和復雜性之間的權(quán)衡，因此，不鼓勵過度擬合。較小的AIC是首選。

在AIC值較小的情況下，同時具有性別和學前教育預測因子的模型優(yōu)于只具有性別預測因子的模型。

正確分類率

正確分類率是另一個有用的衡量標準，可以看出模型對數(shù)據(jù)的合適程度。

#使用`predict()`函數(shù)，從擬合的模型中計算出原始數(shù)據(jù)中學生的預測概率Pred?<-?if_else(Pred?>?0.5,?1,?0)ConfusionMatrix?<-?table(Pred,?TRUE)#正確的分類率

我們可以看到，該模型對所有觀測值的85.8%進行了正確分類。然而，仔細觀察可以發(fā)現(xiàn)，模型預測所有的觀察值都屬于 "0 "類，也就是說，所有的學生都被預測為不留級?？紤]到留級變量的多數(shù)類別是0（不），該模型在分類上的表現(xiàn)并不比簡單地將所有觀測值分配到多數(shù)類別0（不）更好。

AUC（曲線下面積）

使用正確分類率的一個替代方法是曲線下面積（AUC）測量。AUC測量區(qū)分度，即測試對有目標反應和無目標反應的人進行正確分類的能力。在目前的數(shù)據(jù)中，目標變量是留級。我們從 "留級 "組和 "不留級 "組中隨機抽取一名學生。預測概率較高的學生應該是 "留級 "組中的學生。AUC是隨機抽出的對子的百分比。這個程序?qū)UC與正確分類率區(qū)分開來，因為AUC不依賴于結(jié)果變量中類的比例的變化。0.50的值意味著該模型的分類效果不比隨機好。一個好的模型應該有一個遠遠高于0.50的AUC分數(shù)（最好高于0.80）。

#?計算用該模型預測類別的AUC AUC?<-?performance(Pred,?measure?=?"auc") AUC?<-?AUC@y.values[[1]]AUC

AUC分數(shù)為0.60，該模型的判別能力不強。

二項式 Logistic 回歸

正如開頭提到的，邏輯回歸也可以用來為計數(shù)或比例數(shù)據(jù)建模。二項邏輯回歸假設結(jié)果變量來自伯努利分布（這是二項分布的一個特例），其中試驗次數(shù)n為1，因此結(jié)果變量只能是1或0。相反，二項邏輯回歸假設目標事件的數(shù)量遵循二項分布，試驗次數(shù)n，概率q。這樣一來，二項邏輯回歸允許結(jié)果變量取任何非負整數(shù)值，因此能夠處理計數(shù)數(shù)據(jù)。

教育數(shù)據(jù)記錄了集中在學校內(nèi)的個別學生的信息。通過匯總各學校留級的學生人數(shù)，我們得到一個新的數(shù)據(jù)集，其中每一行代表一所學校，并有關(guān)于該學校留級學生的比例信息。學校平均社會經(jīng)濟地位（平均SES分數(shù)）也是在學校層面上的；因此，它可以用來預測在某個學校留級的學生的比例或數(shù)量。

轉(zhuǎn)換數(shù)據(jù)

在這個新的數(shù)據(jù)集中，留級指的是留級的學生人數(shù)；TOTAL指的是某所學校的學生總數(shù)。

探索數(shù)據(jù)

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

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

擬合二項式Logistic回歸模型

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

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

解釋

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

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

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

我們可以直觀地看到學校平均社會經(jīng)濟地位的效果。

plot(allEffects)

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

多層次二元邏輯回歸

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

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

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

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

我們還可以通過學校來繪制性別和留級之間的關(guān)系，以了解性別和留級之間的關(guān)系是否因?qū)W校而異。

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

在上面的圖中，不同的顏色代表不同的學校。我們可以看到，不同學校的性別和留級之間的關(guān)系似乎有很大不同。

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

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

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

由于上述觀察結(jié)果，我們可以得出結(jié)論，在目前的數(shù)據(jù)中需要建立多層次的模型，不僅要有隨機截距（學校），還可能要有性別和學前教育的隨機斜率。

中心化變量

在擬合多層次模型之前，有必要采用適當?shù)闹行幕椒ǎ淳抵行幕︻A測變量進行中心化，因為中心化方法對模型估計的解釋很重要。根據(jù)Enders和Tofighi（2007）的建議，我們應該對第一層次的預測因子性別和學前教育使用中心化，對第二層次的預測因子學校平均社會經(jīng)濟地位使用均值中心化。

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

只有截距模型

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

我們首先指定一個純截距模型，以評估數(shù)據(jù)聚類結(jié)構(gòu)的影響。

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

下面我們計算一下純截距模型的ICC（類內(nèi)相關(guān)）。

0.33的ICC意味著結(jié)果變量的33%的變化可以被數(shù)據(jù)的聚類結(jié)構(gòu)所解釋。這提供了證據(jù)表明，與非多層次模型相比，多層次模型可能會對模型的估計產(chǎn)生影響。因此，多層次模型的使用是必要的，也是有保證的。

完整模型

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

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

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

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

sum(Model_Full)

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

plot((Model)

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

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

似然比檢驗

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

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

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

AIC

AIC?#full模型AIC#＃沒有性別的模型AIC?#＃沒有受過學前教育的模型AIC#＃沒有隨機斜率的模型

從AIC的結(jié)果來看，我們發(fā)現(xiàn)包括隨機斜率項要么沒有大幅提高AIC（用較低的AIC值表示），要么導致更差的AIC（即更高）。因此，我們也得出結(jié)論，沒有必要包括隨機效應項。

其他族（分布）和鏈接函數(shù)

到目前為止，我們已經(jīng)介紹了二元和二項邏輯回歸，這兩種回歸都來自于二項家族的logit鏈接。然而，還有許多分布族和鏈接函數(shù)，我們可以在glm分析中使用。例如，為了對二元結(jié)果進行建模，我們還可以使用probit鏈接或log-log（cloglog）來代替logit鏈接。為了給計數(shù)數(shù)據(jù)建模，我們也可以使用泊松回歸，它假設結(jié)果變量來自泊松分布，并使用對數(shù)作為鏈接函數(shù)。

參考文獻

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

本文選自《R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)》。

點擊標題查閱往期內(nèi)容

R語言線性混合效應模型（固定效應&隨機效應）和交互可視化3案例
非線性混合效應 NLME模型對抗哮喘藥物茶堿動力學研究
生態(tài)學模擬對廣義線性混合模型GLMM進行功率（功效、效能、效力）分析power analysis環(huán)境監(jiān)測數(shù)據(jù)
有限混合模型聚類FMM、廣義線性回歸模型GLM混合應用分析威士忌市場和研究專利申請數(shù)據(jù)
如何用潛類別混合效應模型(Latent Class Mixed Model ,LCMM)分析老年癡呆年齡數(shù)據(jù)
R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)R語言 線性混合效應模型實戰(zhàn)案例
R語言混合效應邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言建立和可視化混合效應模型mixed effect model
R語言LME4混合效應模型研究教師的受歡迎程度
R語言 線性混合效應模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對學術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對學術(shù)能力測驗（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對學術(shù)能力測驗建立層次（分層）貝葉斯模型
SPSS中的多層（等級）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)
用SPSS估計HLM多層（層次）線性模型模型