R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据（上）：https://developer.aliyun.com/article/1490589

二项式 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