下载全文：http://tecdat.cn/?p=22302
在混合效应逻辑回归用于建设二元后果变量的模型，其中，当数据被分组或同时存在固定和随机效应时，后果的对数几率被建模为预测变量的线性组合 （ 点击文末“浏览原文”获取残缺 代码数据 ）。
最近咱们被客户要求撰写对于混合效应逻辑回归的钻研报告，包含一些图形和统计输入。
混合效应逻辑回归的例子
例 1：一个钻研人员对 40 所不同大学的申请进行抽样调查，以钻研预测大学录取的因素。预测因素包含学生的高中 GPA、课外活动和 SAT 分数。一些学校的选择性较多或较少，所以每所学校的基准录取概率是不同的。学校层面的预测因素包含学校是公立还是私立，目前学生与老师的比例，以及学校的排名。
例 2：一家大型 HMO 想晓得哪些病人和医生的因素与病人的肺癌在医治后是否失去缓解最相干，这是一项对于肺癌病人的医治成果和生存品质的钻研的一部分。
例 3：一家电视台想晓得工夫和广告流动如何影响人们是否观看电视节目。他们对四个城市的人进行了为期六个月的抽样调查。每个月，他们都会询问人们在过来一周是否观看了某个节目。三个月后，他们在四个城市中的两个城市推出了一个新的广告流动，并持续监测人们是否观看了该节目。
数据形容
在这个例子中，咱们将应用一个模仿的数据集来探讨对于肺癌的例子。咱们收集了病人的各种后果，他们被蕴含在医生身上，而医生又被蕴含在医院里。还有一些医生层面的变量，比方咱们将在例子中应用的 “ 医生教训 ”。
within(hp, {已婚 <- factor(已婚, levels = 0:1, labels = c("否", "是"))
DID <- factor(DID)
HID <- factor(HID)
癌症阶段 <- factor(癌症阶段)当初咱们要对咱们的间断预测变量进行绘图。数据的可视化能够帮忙咱们了解散布状况，发现编码谬误（例如，咱们晓得一个变量的取值范畴是 0 到 7，但咱们在图中看到了 999），并让咱们理解变量之间的关系。例如，咱们可能看到两个预测因子高度相干，于是决定只在模型中包含一个，或者咱们可能留神到两个变量之间有曲线关系。数据可视化是一种疾速、直观的形式，能够一次性查看所有这些状况。如果你的大多数预测因子看起来都是互相独立的，数据很好。例如，如果它们是独立的，当你输出另一个预测因子时，一个预测因子的估计值不应该有太大变动（只管标准误差和显著性测验可能会有）。咱们能够通过简略地查看数据来理解所有这些信息以及判断如何建模。
ggpairs(hp[, c("IL6", "CRP", "住院工夫", "医生教训")],咱们的间断预测因子之间仿佛没有强的线性关系。让咱们看看咱们的变量在癌症阶段中的散布状况。因为住院工夫是以天为单位的，咱们能够用气泡图来钻研癌症阶段与它的关系。每个气泡的面积与具备这些数值的察看值的数量成正比。对于间断的预测因子，咱们应用小提琴图。所有的原始数据都按癌症阶段离开显示。
点击题目查阅往期内容
[](http://mp.weixin.qq.com/s?__b…)R 语言建设和可视化混合效应模型 mixed effect model
左右滑动查看更多
01
02
03
04
咱们减少了小提琴图。小提琴图只是围绕绘图轴反映的核密度图。咱们将小提琴图绘制在具备透明度的抖动点之上，这样就能够看到原始数据。
因为 IL6 和 CRP 都有偏斜散布的偏向，所以咱们在 Y 轴上应用了平方根刻度。散布看起来相当失常和对称，你依然能够看到长的右尾，即便应用了平方根刻度（留神，只有刻度被转移，数值自身没有被转换，这很重要，因为这让你看到并解释理论的分数，而不是分数的平方根）。
ggplot(stat_sum(aes(size = ..n.., group = 1)) +
scale_size_area(max_size=10)ggplot(tp, +
geom_jitter(alpha = .1) +
geom_violin(alpha = .75) +因为很难看到二元变量在连续变量的程度上如何变动，咱们能够反过来看看二元后果的每个程度上的连续变量的散布。
ggplot(tmp,
geom_boxplot() +
facet_wrap(~variable, scales="free_y")分析方法
上面是一个分析方法的列表：
- 混合效应逻辑回归，是本页面的重点。
- 混合效应 probit 回归与混合效应 logistic 回归十分类似，但它应用的是正态 CDF 而不是 logistic CDF。两者都对二元后果进行建模，能够包含固定和随机效应。
- 固定效应逻辑回归在这种状况下是无限的，因为它可能疏忽了必要的随机效应和 / 或数据中的非独立性。
- 固定效应的 probit 回归在这种状况下是无限的，因为它可能疏忽了必要的随机效应或数据中的非独立性。
- 有聚类持重标准差的 Logistic 回归。这些能够调整非独立性，但不容许有随机效应。
- 有聚类持重标准差的 Probit 回归。这些能够调整非独立性，但不容许有随机效应。
混合效应逻辑回归
上面咱们应用 glmer 命令预计混合效应逻辑回归模型，Il6、CRP 和住院工夫为患者程度的间断预测因素，癌症阶段为患者程度的分类预测因素（I、II、III 或 IV），教训为医生程度的间断预测因素，还有 DID 的随机截距，医生 ID。
# 预计模型并将后果存储在 m 中
# 输入后果，固定成果之间不相干
print(m, corr = FALSE)第一局部通知咱们，估计值是基于自适应高斯 - 赫米特的似然性近似。为了避免出现不收敛的正告，咱们用参数 control=glmerControl(optimizer=”bobyqa”)指定不同的优化器。
下一节给咱们提供了可用于比拟模型的根本信息，接着是随机效应估计值。这示意对数尺度上截距的预计变动。如果有其余随机效应，比方随机斜率，它们也会呈现在这里。最下面的局部最初是察看值的总数和第 2 级察看值的数量。在咱们的案例中，这包含病人（8,525）和医生（407）的总数。
最初一节是固定效应估计值的表格。这些估计值代表回归系数。这些是未标准化的，而且是在对数尺度上。估计值前面是它们的标准误差（SE）。系数预计的近似值可能比 SEs 的近似值稳固得更快。Wald 测验，（frac{Estimate}{SE}），依赖于渐进实践，这里指的是当最高级别的单位大小收敛到无穷大时，这些测验将呈正态分布，并由此得出 p 值（鉴于实在估计值为 0，取得察看估计值或更极其的概率）。
取得置信区间（CI）。咱们能够应用 SE 来取得粗略的区间预计。
# 带有 95%CI 的预计表
cbind(Est = fixef(m), LL = fixef(m) - 1.96 * se, UL = fixef(m) + 1.96 *
se))如果咱们须要比值比而不是对数刻度上的系数，则能够对估计值和 CI 求幂。
多层 bootstrapping（自助法）
从 GLMMs 进行推断是很简单的。除了在每个档次（尤其是最高档次）有很多观测值的状况下，假如（frac{Estimate}{SE}）是正态分布可能不精确。人们提出了各种代替办法，包含蒙特卡洛模仿、贝叶斯预计和 bootstrapping。每种办法的施行都可能很简单。咱们将重点探讨一个小的 bootstrapping 例子。
Bootstrapping 是一种重抽样办法，就是利用无限的样本材料经由多次重复抽样，从新建设起足以代表母体样本分布的新样本。它决不是完满的，但它在概念上是间接易懂的，而且容易在代码中实现。一个毛病是，它对计算要求很高。对于大型数据集或简单的模型，每个模型的运行须要几分钟，在成千上万的样本上进行预计，很容易须要几个小时或几天。在本页的例子中，咱们应用了非常少的样本，但在实践中你会应用更多的样本。
对于单层次模型，咱们能够实现简略的随机抽样，并进行替换，以进行 bootstrapping。对于多层次数据，咱们心愿以与数据生成机制雷同的形式从新取样。咱们从最高级别开始从新取样，而后逐级向下。在咱们的案例中，咱们首先将从医生那里取样，而后在每个取样的医生中，咱们将从他们的病人那里取样。要做到这一点，咱们首先须要写一个函数，在每个档次上从新取样。
sample <- function(dat, cluvar) {id <- sample(cid, size = cid * reps, replace = TRUE)
if (replace) {lapply(seq_along(cid), function(i) {cbind( ID = i, RowID = sample( (dat[, cluvar] == recid[i]),
size = length((dat[, cluvar] == recid[i]) )
})
} else {lapply(seq_along(cid), function(i) {cbind( ID = i, RowID = (dat[, cluvar] == recid[i]))
)
} )当初，咱们将从新对数据进行取样，并采取 100 次反复。同样在实践中，你可能会采取数千次。咱们设置种子，以便咱们的后果能够反复。你也很可能须要比你最终想要的更多的反复样本，因为许多样本可能不收敛，所以你不能从它们那里失去预计。
tmp <- sample(hdp, "DID", reps = 100)接下来，咱们在从新取样的数据上从新拟合模型。首先，咱们存储原始模型的估计值，咱们将用它作为自助模型的起始值。而后，咱们建设一个有 4 个节点的本地集群。接下来，咱们导出数据并在集群上加载。最初，咱们写一个函数来拟合模型并返回估计值。对 glmer()的调用被封装在 try 中，因为不是所有的模型都能在从新采样的数据上收敛。这样能够捕捉到谬误并返回，而不是进行解决。
makeCluster(4)
clusterEvalQ(cl, require(lme4))boot <- function(i) {
object <- try(glmer(缓解 ~ IL6 + CRP + 癌症阶段 + 住院工夫 subset = Replicate == i, family = binomial,
if (class(object) == "try-error")
return(object)
}当初咱们曾经有了数据、本地集群和拟合函数的设置，咱们筹备理论进行 bootstrapping 了。来自所有节点的后果被汇总回一个繁多的列表，存储在对象 res 中。一旦实现，咱们就能够敞开本地集群，终止额定的 R 实例并开释了内存。
parLapplyLB(cl, X = levels(data$Replicate), fun = boot)
# 敞开集群
stopCluster(cl)当初咱们有了 bootstrapping 法的后果，咱们能够对其进行总结。首先，咱们计算胜利收敛的模型的数量。咱们能够计算胜利的平均数，以看到收敛的比例。
# 计算胜利收敛的模型的比例
succ <- sapply(res, is.numeric)
mean(succ)## [1] 1接下来咱们把疏导后果列表转换成矩阵，而后计算每个参数的 2.5 和 97.5 百分位。最初，咱们能够将后果制成表格，包含原始估计值和标准误差、均匀疏导估计值，以及 bootstrap 的置信区间。
# 合并胜利的后果
do.call(cbind, res[success])
# 计算出 95%CI 的 2.5 和 97.5 百分位数。apply(res, 1，quantile, probs = c(0.025, 0.975)## 2.5% 97.5%
## (Intercept) -3.61982 -0.985404
## IL6 -0.08812 -0.029664
## CRP -0.04897 0.006824
## 癌症阶段 II -0.60754 -0.228019
## 癌症阶段 III -1.30217 -0.754609
## 癌症阶段 IV -2.91414 -2.002643
## 住院工夫 -0.21596 -0.046420
## 医生教训 0.06819 0.207223
## NewID.(Intercept) 2.03868 2.476366# 所有后果取整输入
round(able, 3)## Est SE BootMean 2.5% 97.5%
## (Intercept) -2.053 0.531 -2.205 -3.620 -0.985
## IL6 -0.057 0.012 -0.059 -0.088 -0.030
## CRP -0.021 0.010 -0.022 -0.049 0.007
## 癌症阶段 II -0.414 0.076 -0.417 -0.608 -0.228
## 癌症阶段 III -1.003 0.098 -1.043 -1.302 -0.755
## 癌症阶段 IV -2.337 0.158 -2.460 -2.914 -2.003
## 住院工夫 -0.121 0.034 -0.142 -0.216 -0.046
## 医生教训 0.120 0.027 0.128 0.068 0.207
## DID.(Intercept) 2.015 NA 2.263 2.039 2.476预测的概率和绘图
这些后果很适宜放在表格中或钻研文本中；然而，数字的解释可能很麻烦。图形展现有助于解释，也有助于演讲。
在一个逻辑模型中，后果通常是
- 对数几率（也叫对数），这是线性化
- 指数化的对数几率，不在线性尺度上
- 概率
对于表格来说，人们常常出现的是几率比。对于可视化来说，对数或概率比例是最常见的。每种办法都有一些长处和毛病。对数表很不便，因为它是线性化的，这意味着一个预测因素减少 1 个单位，后果就会减少一个系数单位，而且无论其余预测因素的程度如何。毛病是这个量表的可解释性不强。读者很难对对数有一个直观的了解。相同，概率是一个很好的尺度，能够直观地了解后果；然而，它们不是线性的。这意味着预测因子减少一个单位，不等于概率的恒定减少 – 概率的变动取决于为其余预测因子抉择的值。在一般逻辑回归中，你能够放弃所有预测因子不变，只扭转你感兴趣的预测因子。然而，在混合效应逻辑模型中，随机效应也对后果产生影响。因而，如果你放弃所有不变，那么只有当所有协变量放弃不变，并且你在同一组或具备雷同随机效应的一组时，后果的概率变动才是真的。
咱们将探讨一个均匀边际概率的例子。这比条件概率须要更多的工作，因为你必须为每一组计算独自的条件概率，而后将其平均化。
首先，让咱们应用这里的符号来定义个别程序。咱们通过获取 并将感兴趣的特定预测因子，比如说在 j 列，设置为常数来创立。如果咱们只关怀预测器的一个值，那就是。然而，更常见的是，咱们心愿预测因子有肯定的取值范畴，以便绘制预测概率在其范畴内的变动状况。咱们能够通过获取预测模型的察看范畴，并在该范畴内平均地抽取 k 个样本。例如，假如咱们的预测模型的范畴是 5 到 10，咱们想要 6 个样本，，所以每个样本将与前一个样本相隔 1，它们将是. 而后咱们创立不同的 k 个不同的 Xi，其中，在每种状况下，第 j 列被设置为某个常数。而后咱们计算：
这些是所有不同的线性预测因子。最初，咱们采取，这就失去，这是原始尺度上的条件期望，在咱们的例子中是概率。而后咱们能够取每个的期望值，并将其与咱们感兴趣的预测因子的值作比照。咱们还能够绘制图表，不仅显示均匀边际预测概率，而且还显示预测概率的散布。
你可能曾经留神到，这些估计值中有很多变数。咱们在应用 时，只将咱们感兴趣的预测因子放弃在一个常数，这使得所有其余预测因子都能在原始数据中取值。另外，咱们把 留在咱们的样本中，这意味着有些组的代表性比其余组要高或低。如果咱们想的话，咱们能够对所有的群体进行从新加权，使其具备等同的权重。在这个例子中，咱们抉择让所有这些货色放弃原样，是基于这样的假如：咱们的样本的确是咱们感兴趣的人群的良好代表。咱们没有试图筛选有意义的值来放弃协变量（，而是应用了咱们样本的值。这也表明，如果咱们的样本能很好地代表总体，那么均匀边际预测概率就能很好地代表咱们总体中新的随机样本的概率。
当初咱们有了一些背景和实践，咱们看看如何理论去计算这些货色。咱们失去一个住院工夫（咱们感兴趣的预测因子）的摘要，而后在其范畴内失去 100 个值，用于预测。咱们复制一份数据，这样咱们就能够固定其中一个预测因子的值，而后应用预测函数来计算预测值。默认状况下，所有的随机效应都被包含在内。
# 计算预测的概率并存储在列表中
lapply(jvalues, function(j) {predict(m, newdata = tmpdat, type = "response")当初咱们有了所有的预测概率，能够可视化它们。例如，咱们能够看一下多数不同停留时间的均匀边际预测概率。
# 几种不同工夫的均匀边际预测概率 t(sapply(pp, function(x) {c(M = mean(x), quantile(x, c(0.25, 0.75)))
# 退出工夫的数值并转换为数据框
as.data.frame(cbind(plotdat, values))
# 显示前几行
head(plotdat# 绘制均匀边际预测概率
ggplot(plotdat) + geom_line() +咱们还能够加上上限和四分位数。能够看到 50% 的预测概率所处的范畴。
ggplot(plotdat) + geom_linerange() + geom_line(size = 2)除了扭转住院工夫之外，咱们还能够对癌症阶段的每一级做同样的均匀边际预测概率。
# 计算预测的概率并存储在一个列表中
lapply(癌症阶段, function(stage) {predict(m, newdata = tmpdat, type = "response")
})
# 取得每个级别的癌症阶段的所有 j 值的平均值和四分位数
lapply(probs, function(X) {c(M=mean(x), quantile(x, c(.25, .75)))
# 放到一个数据框
do.call(rbind, plotdat2)
# 增加癌症阶段
factor(rep(levels(癌症阶段))
# 显示前几行数据
head(plotdat2)# 绘制
ggplot(plotdat2) +
geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = 癌症阶段), alpha = .15) +对于一个住院 10 天的第四期肺癌患者，其癌症失去缓解的机会看起来相当小。看起来散布也是偏斜的。咱们能够检查一下仅针对该组的预测概率分布。
ggplot(aes(Probs)) + geom_histogram() +即便应用平方根尺度，将较低的数值拉长，它依然是极其偏斜的。据估计，绝大多数人的病情缓解的概率不到 0.1。
三层混合效应逻辑回归
咱们曾经深入研究了一个带有随机截距的两级逻辑模型。这是最简略的混合效应逻辑模型。当初咱们要简要地看一下如何减少第三档次和随机斜率效应以及随机截距。
上面咱们预计一个三层逻辑模型，医生有一个随机截距，医院有一个随机截距。在这个例子中，医生被嵌套在医院内，也就是说，每个医生属于一家而且只有一家医院。另一种状况有时被称为 “ 穿插分类 ”，意思是一个医生可能属于多家医院，比方该医生的一些病人来自 A 医院，另一些来自 B 医院。在 glmer 中，你不须要指定组是嵌套还是穿插分类，R 能够依据数据计算出来。
# 输入没有固定效应之间相关性的 mod 后果
print(m3a, corr=FALSE)输入通知咱们族（二元后果的二项式）和链接函数（logit）。接着是通常的拟合指数和随机效应的方差。在这种状况下，医生之间和医院之间的截距（在对数赔率尺度上）的变动。还显示了标准差（只是方差的平方根，而不是预计方差的标准误差）。咱们还失去了每个档次上的单位的数量。最初是固定效应。
看一下条件模型的散布也是很有用的，上面咱们用“毛毛虫图”来做。蓝点是带有误差条的条件模型。咱们对医生和医院都是这样做的。例如，对于医生来说，咱们能够看到一个有点长的右尾，即极其的正值比负值多。
lattice::dotplot((m) ))## $DID## $HID咱们也能够在模型中退出随机斜率。咱们只是要为 “ 住院工夫 “ 减少一个随机斜率，这个斜率在不同的医生之间变动。就像在惯例的 R 公式中一样，咱们应用 + 运算符来 “ 增加 “ 一个效应。
# 输入没有固定效应之间相关性的 mod 后果
print(m3b, corr = FALSE)## $DID## $HID点击文末 “浏览原文”
获取全文残缺材料。
本文选自《 R 语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据》。
点击题目查阅往期内容
多程度模型、分层线性模型 HLM、混合效应模型钻研老师的受欢迎水平
R 语言 nlme、nlmer、lme4 用（非）线性混合模型 non-linear mixed model 剖析藻类数据实例
R 语言混合线性模型、多层次模型、回归模型剖析学生均匀问题 GPA 和可视化
R 语言线性混合效应模型（固定效应 & 随机效应）和交互可视化 3 案例
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 多层（档次）线性模型模型