---

可下载资源

---

本文分析的数据分享至会员群

---

作者

例2：一家大型HMO想知道哪些病人和医生的因素与病人的肺癌在治疗后是否得到缓解最相关，这是一项关于肺癌病人的治疗效果和生活质量的研究的一部分。

例3：一家电视台想知道时间和广告活动如何影响人们是否观看电视节目。他们对四个城市的人进行了为期六个月的抽样调查。每个月，他们都会询问人们在过去一周是否观看了某个节目。三个月后，他们在四个城市中的两个城市推出了一个新的广告活动，并继续监测人们是否观看了该节目。

within(hp, {
  已婚 <- factor(已婚, levels = 0:1, labels = c("否", "是"))
  DID <- factor(DID)
  HID <- factor(HID)
  癌症阶段 <- factor(癌症阶段)

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回归。这些可以调整非独立性，但不允许有随机效应。

第一部分告诉我们，估计值是基于自适应高斯-赫米特的似然性近似。为了避免出现不收敛的警告，我们用参数control=glmerControl(optimizer=”bobyqa”)指定不同的优化器。

下一节给我们提供了可用于比较模型的基本信息，接着是随机效应估计值。这表示对数尺度上截距的估计变化。如果有其他随机效应，比如随机斜率，它们也会出现在这里。

最后一节是固定效应估计值的表格。这些估计值代表回归系数。这些是未标准化的，而且是在对数尺度上。估计值后面是它们的标准误差（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是一种重抽样方法，就是利用有限的样本资料经由多次重复抽样，重新建立起足以代表母体样本分布的新样本。它决不是完美的，但它在概念上是直接易懂的，而且容易在代码中实现。一个缺点是，它对计算要求很高。对于大型数据集或复杂的模型，每个模型的运行需要几分钟，在成千上万的样本上进行估计，很容易需要几个小时或几天。在本页的例子中，我们使用了非常少的样本，但在实践中你会使用更多的样本。

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)
}

 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

在普通逻辑回归中，你可以保持所有预测因子不变，只改变你感兴趣的预测因子。然而，在混合效应逻辑模型中，随机效应也对结果产生影响。因此，如果你保持一切不变，那么只有当所有协变量保持不变，并且你在同一组或具有相同随机效应的一组时，结果的概率变化才是真的。

首先，让我们使用这里的符号来定义一般程序。我们通过获取 并将感兴趣的特定预测因子，比如说在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)

在这种情况下，医生之间和医院之间的截距（在对数赔率尺度上）的变化。还显示了标准差（只是方差的平方根，而不是估计方差的标准误差）。我们还得到了每个层次上的单位的数量。最后是固定效应。

看一下条件模型的分布也是很有用的，下面我们用“毛毛虫图”来做。蓝点是带有误差条的条件模型。我们对医生和医院都是这样做的。例如，对于医生来说，我们可以看到一个有点长的右尾，即极端的正值比负值多。

lattice::dotplot( (m) ))

## $DID

## $HID

## $HID