从具有随机截距的多级模型生成预测模拟

Question

我正在努力理解如何在 R 中使用具有一组随机截距的多级线性回归模型为新数据生成预测模拟。按照 this text 第 146-147 页的示例，我可以针对没有随机效应的简单线性模型执行此任务。我无法解决的问题是如何扩展设置以适应添加到该模型的因子的随机截距。

我将使用 iris 和一些假数据来显示我遇到困难的地方。我将从一个简单的线性模型开始：

mod0 <- lm(Sepal.Length ~ Sepal.Width, data = iris)

现在让我们使用该模型为 250 个新案例生成 1,000 个预测模拟。我将从弥补这些情况开始：

set.seed(20912)
fakeiris <- data.frame(Sepal.Length = rnorm(250, mean(iris$Sepal.Length), sd(iris$Sepal.Length)),
                       Sepal.Width = rnorm(250, mean(iris$Sepal.Length), sd(iris$Sepal.Length)),
                       Species = sample(as.character(unique(iris$Species)), 250, replace = TRUE),
                       stringsAsFactors=FALSE)

按照上述文本中的示例，这是我为这 250 个新案例中的每一个案例进行 1,000 次预测模拟的方法：

library(arm)
n.sims = 1000  # set number of simulations
n.tilde = nrow(fakeiris)  # set number of cases to simulate
X.tilde <- cbind(rep(1, n.tilde), fakeiris[,"Sepal.Width"])  # create matrix of predictors describing those cases; need column of 1s to multiply by intercept
sim.fakeiris <- sim(mod0, n.sims)  # draw the simulated coefficients
y.tilde <- array(NA, c(n.sims, n.tilde))  # build an array to hold results
for (s in 1:n.sims) { y.tilde[s,] <- rnorm(n.tilde, X.tilde %*% sim.fakeiris@coef[s,], sim.fakeiris@sigma[s]) }  # use matrix multiplication to fill that array

效果很好，现在我们可以做 colMeans(y.tilde) 之类的事情来检查这些模拟的中心趋势，并 cor(colMeans(y.tilde), fakeiris$Sepal.Length) 将它们与 [=66 的（假）观察值进行比较=].

现在让我们尝试扩展那个简单的模型，在这个模型中我们假设截距因观察组而异——这里是物种。我将使用 lme4 包中的 lmer() 来估计与该描述相匹配的简单 multilevel/hierarchical 模型：

library(lme4)
mod1 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris)

好的，行得通，但是现在呢？我运行:

sim.fakeiris.lmer <- sim(mod1, n.sims)

当我使用 str() 检查结果时，我看到它是 class sim.merMod 的对象，具有三个组件：

@fixedef，一个 1,000 x 2 矩阵，带有固定效应的模拟系数（截距和 Sepal.Width）
@ranef，一个 1,000 x 3 矩阵，带有随机效应（三个物种）的模拟系数
@sigma，一个长度为 1,000 的向量，包含与每个模拟相关的西格玛

我不知道如何将用于简单线性模型的矩阵构造和乘法扩展到这种情况，这增加了另一个维度。我查看了文本，但我只能找到一个示例（第 272-275 页），用于单个组（此处为物种）中的单个案例。我打算执行的真实世界任务涉及运行模拟这样的 256 个新案例（职业足球比赛）均匀分布在 32 个组（主队）中。如果您能提供任何帮助，我将不胜感激。

附录。愚蠢的是，我在发布之前没有查看 lme4 中 simulate.merMod() 的详细信息。我现在有了。看起来应该可以解决问题，但是当我运行 simulate(mod0, nsim = 1000, newdata = fakeiris) 时，结果只有 150 行。这些值看起来合理，但 fakeiris 中有 250 行（个案）。那 150 是从哪里来的？

Answer 1

这可能会有所帮助：它不使用 sim()，而是使用 mvrnorm() 从固定效应参数的抽样分布中提取新系数，使用一些内部机制(setBeta0) 重新分配固定效应系数的内部值。随机效应系数的内部值由 simulate.merMod 使用默认参数 re.form=NA 自动重采样。但是，残差方差未重新采样——它在整个模拟过程中保持固定，这不是 100% 真实的。

在您的用例中，您将指定 newdata=fakeiris。

library(lme4)
mod1 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris)
simfun <- function(object,n=1,newdata=NULL,...) {
    v <- vcov(object)
    b <- fixef(object)
    betapars <- MASS::mvrnorm(n,mu=b,Sigma=v)
    npred <- if (is.null(newdata)) {
                 length(predict(object))
             } else nrow(newdata)
    res <- matrix(NA,npred,n)
    for (i in 1:n) {
        mod1@pp$setBeta0(betapars[i,])
        res[,i] <- simulate(mod1,newdata=newdata,...)[[1]]
    }
    return(res)
}
ss <- simfun(mod1,100)

Answer 2

一种可能是使用 merTools 包中的 predictInterval 函数。该软件包即将提交给 CRAN，但可以从 GitHub、

下载当前的开发版本

    install.packages("devtools")
    devtools::install_github("jknowles/merTools")

要获得 100 次模拟的中位数和 95% 可信区间：

    mod1 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris)

    out <- predictInterval(mod1, newdata=fakeiris, level=0.95,
                           n.sims=100, stat="median")

默认情况下，predictInterval 包括残差，但您可以使用以下命令关闭该功能：

    out2 <- predictInterval(mod1, newdata=fakeiris, level=0.95,
                           n.sims=100, stat="median", 
                           include.resid.var=FALSE)

希望对您有所帮助！

从具有随机截距的多级模型生成预测模拟

Generating predictive simulations from a multilevel model with random intercepts

simulation

r

bayesian

prediction

lme4