使用 merTools 中的 predictInterval() 从 glmer 对象生成边际预测置信区间

Question

我正在尝试使用 predictInterval 函数 further described here 为边际预测生成置信区间。

这里我使用了 ResourceSelection 包中的 goats 数据，其中包含已使用和可用的位置（分别编码为 1 和 0）以及感兴趣的协变量值（例如海拔、坡度等）建立一个可复制的模型。

套餐

library(lme4)
library(ResourceSelection)
library(merTools)

df 包含 10 只动物的已用和可用位置。

table(goats$ID, goats$STATUS)
       0    1
  1  1404  702
  2  1112  556
  3  1026  513
  4   634  317
  5  1272  636
  6  1456  728
  7  1394  697
  8  1468  734
  9  1608  804
  10 1302  651

下面是一个示例模型，其中为个人 (ID) 指定了随机截距。协变量使用 scale().

在模型拟合中居中和缩放

 mod <- glmer(STATUS ~ scale(ELEVATION) + scale(SLOPE) + scale(ET) + scale(HLI) + (1|ID),
             family=binomial, data = goats, verbos = 1) 
summary(mod)

我现在想预测 ELEVATION 的范围，所有其他协变量均取其平均值。因为我使用的是按比例缩放和居中的协变量，所以平均值为 0。比例的最小值和最大值 (ELEVATION) 是 -1.97056 和 2.52926，我用它们来制作下面的新预测数据。

PredDat <- data.frame(ELEVATION = seq(-1.97056, 2.52926, length.out = 1000),
                      SLOPE = 0,
                      ET = 0,
                      HLI = 0)

虽然我可以手动生成预测，但我不确定如何在大型数据集使 bootstrap 方法 (recommended here) 望而却步时估计 95% CI。是否可以在不考虑个体随机效应的情况下使用 predictInterval 函数生成边际预测和 CI？下面的代码导致错误 Error in eval(expr, envir, enclos) : object 'ID' not found，因为 PredDat 数据框中没有 ID。如果我将 ID 添加到 PredDat 数据框，代码运行正常。

Preds <- predictInterval(mod, newdata = PredDat, type = "probability")

任何有关如何从 glmer 对象生成边际预测的建议都将不胜感激。

重要的会话信息粘贴在 FYI 下面。

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

other attached packages:
[1] merTools_0.2.0          plyr_1.8.3             
[3] arm_1.8-6               MASS_7.3-45            
[5] ResourceSelection_0.2-5 lme4_1.1-10            
[7] Matrix_1.2-3            sp_1.2-1

Answer 1

这里是 merTools 的包维护者。我们实现此功能的方式并不是很简单，但是可以做到。

您需要添加一个步骤，将中位数随机效应添加到您的 data.frame。在大多数情况下，中值随机效应应该为 0，或者足够接近，它接近于您正在寻找的结果。为此，您只需稍微修改代码并使用 merTools:

中的 REquantile 函数

medEff = REquantile(mod, quantile = 0.5, 
                    groupFctr = "ID", 
                    term = "(Intercept)")

PredDat <- data.frame(ELEVATION = seq(-1.97056, 2.52926, length.out = 1000),
                      SLOPE = 0, ET = 0, HLI = 0, ID = medEff)

Preds <- predictInterval(mod, newdata = PredDat, type = "probability")

这会产生预测，但包括随机效应的不确定性，包括 0 的中值随机效应。在上面的示例中，这最终消除了 ELEVATION 变量在观察中的影响，因为中值随机效应没有被非常精确地估计。所以，这可能不是你想要的。

此外，如果您有 更复杂的随机效应规范 以及斜率和截距，那么这种方法会变得更难，因为截距的中值效应可能为 0，但它不会'不是斜坡。

如果您真的只想捕获基于固定效应及其不确定性的预测中的方差——自从构建包以来我学到的东西很常见——有很多方法可以在merTools。这不是最优雅的，但它是在 predictInterval 的幕后发生的事情，以获得固定效应预测的可变性：

PredDat <- data.frame(Intercept = 1, 
           ELEVATION = seq(-1.97056, 2.52926,length.out = 1000), 
           SLOPE = 0, ET = 0, HLI = 0)

 fe.tmp <- fixef(mod)
 vcov.tmp <- as.matrix(vcov(mod))
 n.sims <- 1000
 sigmahat <- rep(1, n.sims)

 # Make n.sims draws for each element of the fixed effects

 betaSim <- abind::abind(lapply(1:n.sims,
  function(x) mvtnorm::rmvnorm(n = 1, mean = fe.tmp, 
       sigma = sigmahat[x]*vcov.tmp, method = "chol")), along=1)
# Calculate n.sims predictions for each row in PredDat
fixed <- as.matrix(PredDat) %*% t(betaSim)
# For each row (observation) in PredDat calculate the median, upr and lwr 
Preds <- data.frame(fit = apply(fixed, 1, median), 
                upr = apply(fixed, 1, quantile, 0.9), 
                lwr = apply(fixed, 1, quantile, 0.1))
# Calculate the probability from the linear predictor
Preds <- apply(Preds, 2, invlogit)

你应该得到这样的东西：

head(Preds)

     fit       upr       lwr
1 0.1860053 0.2482220 0.1427370
2 0.1860058 0.2482226 0.1427373
3 0.1860062 0.2482231 0.1427377
4 0.1860066 0.2482237 0.1427380
5 0.1860071 0.2482242 0.1427384
6 0.1860075 0.2482248 0.1427388

这不包括与分组因素或模型本身的变化相关的观察水平的任何不确定性。

使用 merTools 中的 predictInterval() 从 glmer 对象生成边际预测置信区间

Generating marginal prediction confidence intervals from a glmer object using predictInterval() from merTools

r

predict

lme4