从具有随机截距的线性混合模型进行模拟

我正在尝试扩展this question的答案：具体来说， how to build simulations for linear mixed effects models with a random intercept 'from scratch' （没有 simulate.merMod 或 arm ） . 我问，因为我有兴趣重新采样从拟合模型中获得的参数估计值来模拟 - 而不是预测 - 新值的响应，同时重新采样固定效应参数估计，随机效应估计和残差方差 . 我错了或如何改进流程 .

我从Ben Bolker的建议开始，并使用 mvrnorm() 从固定效应参数的采样分布中绘制新系数 .

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

v <- vcov(mod1)
b <- mod1@beta

fixed <- MASS::mvrnorm(1000, mu=b, Sigma=v)

使用 melt 重新整形重新采样的系数，并重命名列 .

library(reshape)
library(dplyr)
fixed.resampled = reshape::melt(fixed, id.vars= "refs") %>%
  dplyr::rename(species=X2)

从答案到question about plotting random effects，提取并保存 ranef 对象，在 qq 中保存截距的方差，并将截距估计，标准差和随机效果级别放在数据框中 .

randoms <- ranef(mod1, condVar = TRUE)
qq <- attr(ranef(mod1, condVar = TRUE)[[1]], "postVar")
rand.interc <- randoms$Species
df <- data.frame(Intercepts=randoms$Species[,1],
               sd.interc=2*sqrt(qq[,,1:length(qq)]),
               lev.names=rownames(rand.interc))

与固定效应系数一样，从采样分布中抽取随机效应估计值 . 每个随机效应估计从具有均值 df$Intercepts 和标准偏差 df$sd.interc 的正态分布中采样 . 创建一个空数组来保存样本 .

ranef.resampled <- array(dim=c(1000,3))

使用物种特定方法和标准偏差，从每个物种的正态分布中取样 .

for(i in 1:length(randoms$Species[,1])){
  ranef.resampled[,i] <- rnorm(1000,mean=df$Intercepts[i],sd=df$sd.interc[i])
}

重塑和重命名数据框 .

ranef.resampled <- data.frame(ranef.resampled) %>%
   dplyr::rename(setosa=X1,versicolor=X2,virginica=X3) %>%
   gather() %>%
   dplyr::rename(species=key)

出于绘图目的，创建包含系数和随机效应估计的数据帧 .

fixed.e <- data.frame(species=names(fixef(mod1)),estimate=fixef(mod1)[1:2])
ranef.e <- data.frame(species=row.names(rand.interc),estimate=rand.interc[,1])

把所有东西放在一起 .

dfLong <- fixed.resampled %>% 
          dplyr::bind_rows(ranef.resampled) %>%
          dplyr::bind_rows(fixed.e) %>%
          dplyr::bind_rows(ranef.e) %>%
          dplyr::rename(X2=species)

绘制系数和随机效应估计值 .

ggplot(dfLong, aes (value,group=X2)) +
  geom_density() + 
  geom_vline(data=dfLong,
    aes(xintercept=estimate),
    color="red",linetype="dashed") +
  facet_grid(X2~.) + 
  theme_classic()

要模拟值，首先要创建一个新数据框来保存要执行模拟的值 . 在这里，我生成了1000个样本，这些样本与 Species="setosa" 的原始数据的分布相匹配 .

set.seed(101)
newDf <- data.frame(Sepal.Width = rnorm(1000,mean(iris$Sepal.Width[1:50]), sd(iris$Sepal.Width[1:50])), Species="setosa")

然后我使用这个新数据来使用_1349850来模拟结果（这是在 mer 模型对象上调用 simulate() 时调用的 .

simulated.setosa <- data.frame(setosa=simulate(mod1, nsim=1, newdata=newDf)[1:1000,])

我想要做的是在不使用 simulate() 的情况下生成类似的预测 . 到目前为止，我采取了以下步骤 . 总和：截距，'width'的参数估计值乘以'width'的模拟值，以及'setosa.'的重采样随机效应估计值

sim.setosa <- data.frame(setosa=fixed[1:1000,1] + fixed[1:1000,2] * newDf$Sepal.Width[1:1000] + ranef.resampled[1:1000,2])

然后我使用以下内容来比较这两个模拟的输出 . 虽然我不希望它们完全匹配（重新采样系数是 simulate() 的内部，但我认为它与我尝试将模拟拼凑在一起的输出相同 .

ggplot() +
  geom_density(data=simulated.setosa,aes(setosa)) +
  geom_density(data=sim.setosa,aes(setosa),col="red") +
  theme_classic() +
  xlim(c(0,10))