关于数据挖掘:数据分享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）要低。

点击题目查阅往期内容

多程度模型、分层线性模型 HLM、混合效应模型钻研老师的受欢迎水平

左右滑动查看更多

01

02

03

04

模型评估: 拟合度

评估逻辑回归模型的拟合度有不同的办法。

似然比测验

如果一个逻辑回归模型与预测因子较少的模型相比，显示出拟合度的进步，则该模型对数据有较好的拟合度。这是用似然比测验进行的，它将残缺模型下数据的似然性与较少预测因素的模型下数据的似然性进行比拟。从一个模型中删除预测变量简直总是会使模型的拟合度升高（即模型的对数似然率较低），但测试察看到的模型拟合度差别是否具备统计学意义是很有用的。

# 指定一个只有 ` 性别 ' 变量的模型
#应用 `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 <- AUC@y.values[[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

点击文末 “浏览原文”

获取全文残缺材料。

本文选自《R 语言用 lme4 多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据》。

点击题目查阅往期内容

R 语言线性混合效应模型（固定效应 & 随机效应）和交互可视化 3 案例
非线性混合效应 NLME 模型反抗哮喘药物茶碱动力学钻研
生态学模仿对狭义线性混合模型 GLMM 进行功率（效用、效力、效劳）剖析 power analysis 环境监测数据
无限混合模型聚类 FMM、狭义线性回归模型 GLM 混合利用剖析威士忌市场和钻研专利申请数据
如何用潜类别混合效应模型 (Latent Class Mixed Model ,LCMM) 剖析老年痴呆年龄数据
R 语言用 lme4 多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据 R 语言 线性混合效应模型实战案例
R 语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据
R 语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言建设和可视化混合效应模型 mixed effect model
R 语言 LME4 混合效应模型钻研老师的受欢迎水平
R 语言 线性混合效应模型实战案例
R 语言用 Rshiny 摸索 lme4 狭义线性混合模型（GLMM）和线性混合模型（LMM）
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言如何解决线性混合模型中畸形拟合 (Singular fit) 的问题
基于 R 语言的 lmer 混合线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言分层线性模型案例
R 语言用 WinBUGS 软件对学术能力测验（SAT）建设分层模型
应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
SPSS 中的多层（等级）线性模型 Multilevel linear models 钻研整容手术数据
用 SPSS 预计 HLM 多层（档次）线性模型模型