How can I best get a prediction interval for future group or nested group means with LMM (lme4)?

Question

I am building a linear mixed model of an experimental data set. The data can be grouped into individual experiments, and further each experiment consists of clusters within which a sample is observed for each level of fixed effects. Based on this data I would like to quantify the repeatability of the experiments and the stability of clusters within an experiment. Thus I am actually most interested in inference the likely range of observed means in the hypothetical of performing an additional experiment or collecting another cluster within an experiment. In my understanding, I am interested in a prediction interval for a new experiment ignoring residual and cluster variation, the uncertainty of the estimate of $\mu_t \pm 2 \sigma_{exp}$, and complementarily ignoring residual and experiment variation to estimate $\mu_t \pm 2 \sigma_{clust}$

In the searching I have done, I haven't really found something exactly like what I want. In fact in the comments to Prediction interval for lmer() mixed effects model in R I see that the predictInterval function is not able to do it.

I have tried to do this via bootstrapping with bootMer, but I'm a bit uncertain about my approach and wonder if there there may be a better way of doing this. Essentially, the method I came up with is to explicitly generate a random intercept and add it to the output of predict in my bootstrap function.

outStats <- function(model){
  var_pars <- as.data.frame(VarCorr(model))
  clust_var <- var_pars[1,]$sdcor # var_pars$clust_exp:exp_no
  pred <- predict(model, newdata=PredData, re.form=NA, allow.new.levels=TRUE)
  out = pred + rnorm(1)*clust_var
  return(out)
}
merBoot <- bootMer(lme1, outStats,  nsim = n_boot, re.form = NA, seed=seed_boot)
clustvar_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
clustvar_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
clustvar_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))

For within experiment variation between clusters, I think this also makes sense.

inStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=~(1|exp_no), allow.new.levels=TRUE)
  return(out$sim_1)
}
merBoot <- bootMer(lme1, inStats, nsim = n_boot, re.form=~(1|exp_no), seed=seed_boot) # This is sort of as expected
inexp_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
inexp_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
inexp_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))

The second approach matches the "true" prediction interval quite well, while the first generates too wide of intervals (I have not explored bootsrap convergence yet).

I'm finding it more confusing to try and ignore the cluster variation to just generate the prediction interval for a new experiment's mean values across clusters. Both approaches seem to be conservative compared to the "true" interval (maybe the sample size and problem definition cause this). The first approach (directly sampling a random intercept 4.3.c) makes more sense to me, but I'm uncertain about it. I also tried just generating a prediction interval for a new cluster in a new experment (4.3.a), and this is only slightly wider than the direct sampling approach

On a final note, it's not clear to me if bootMer and simulate generate realizations for the residual variance or if it's only the fixed effects' uncertainty and the random effects which are simulated.

Full code example

I have setup code comparing a few approaches that make sense to me with simulated data where I know what I think should come out. I find that I'm not really understanding how bootMer and simulate handle the different cases of re.form, especially with out of sample data. Thus some approaches that I think should be identical are not. I.e. re.form=NULL with an unobserved level should generate uncertainty about the random effect value for the new level, but in fact this produces the same width of interval as for the case with observed levels. (4.2.b vs 4.1.b in the code sections)

# Bootstrapping method for CI's and PI's (?) How?
# https://stats.stackexchange.com/questions/147836/prediction-interval-for-lmer-mixed-effects-model-in-r
# https://bookdown.org/dereksonderegger/571/10-mixed-effects-models.html
# See help simulate.merMod whcih suggest re.form=0 is what I want for out of sample simulation
# https://github.com/lme4/lme4/issues/388
# Confidence Intervals
# https://stats.stackexchange.com/a/344062

# 2. Model and example data
library(lme4)
library(Metrics)
library(mixtools)
library(merTools)
library(ggplot2)
library(RColorBrewer)
n_trt = 25
n_exp = 10
n_clust = 100
n_boot = 201
seed_boot = 101
exp_no <- as.factor(rep(1:n_exp, each=n_trt*n_clust))
sigma_exp <- 1
clust_exp <- as.factor(rep(rep(1:n_clust, each=n_trt),n_exp))
sigma_clust <- 0.25
treatment <- as.factor(rep(rep(1:n_trt, n_clust), n_exp))
exp_effect <- rep(rnorm(n_exp)*sigma_exp, each=n_trt*n_clust) # Experiment variance of 1/2500)
clust_effect <- rep(rnorm(n_exp*n_clust)*sigma_clust, each=n_trt) # Cluster variance of 1/10000
sigma_e <- 0.01
response <- unclass(treatment)/10 + exp_effect + clust_effect + rnorm(n_trt*n_exp*n_clust)*sigma_e # Error variance of 1/40000
dat <- data.frame(exp_no=exp_no, clust_exp=clust_exp, treatment=treatment, response=response)
lme1 <- lmer(response ~ treatment + (1 | exp_no/clust_exp), data=dat)
summary(lme1)

cint_res <- confint(lme1)
df <- as.data.frame(cint_res)
sigma_clst_upr <- df[1,2]
sigma_exp_upr <- df[2,2]

# 3. Usage of bootMer to generate CI and PI
#"
#Does bootMer ever sample the residual error? It is not clear to me. (Test with higher variance)
#Does the re.form in bootMer matter or only that in simulate?
# NULL - condition on all RE so only the uncertainty the specific effect of a level is included.
#NA - condition on nothing so all RE are resampled? (Or are they set to population mean?)

# 3.1 Confidence interval for mean
ConfData <- data.frame(treatment=unique(dat$treatment), exp_no=200, clust_exp=200)
myStats <- function(model){
  out <- predict(model, newdata=ConfData,  re.form=NULL, allow.new.levels=TRUE)
  return(out)
}
merBoot <- bootMer(lme1, myStats, nsim = n_boot, re.form = NA, seed=seed_boot) 
CI.mean = apply(merBoot$t, 2, function(x) as.numeric(mean(x, na.rm=TRUE)))
CI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
CI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
CI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
plt_dat <- data.frame(upr=CI.upper, lwr=CI.lower, median=CI.median)
plt_dat <- cbind(plt_dat, ConfData)
plt_dat$treatment <- unclass(plt_dat$treatment)

# 4.0
# "Given sigma_exp and sigma_clust the 'true' prediction intervals are"
# Normal quantiles with true values
exp_plt_true <- plt_dat
exp_plt_true$lwr <- plt_dat$median-2*sigma_exp
exp_plt_true$upr =plt_dat$median+2*sigma_exp
clust_plt_true <- plt_dat
clust_plt_true$lwr <- plt_dat$median-2*sigma_clust
clust_plt_true$upr =plt_dat$median+2*sigma_clust

# 4.1 Case of prediction interval for fully observed level
#"This should just be a standard prediction interval for fixed effects"
PredData <- data.frame(treatment=unique(dat$treatment), exp_no=1, clust_exp=1)
# 4.1.a conditional on all
inStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=NULL, allow.new.levels=TRUE)
  return(out$sim_1)
}
merBoot <- bootMer(lme1, inStats, nsim = n_boot, re.form = NULL, seed=seed_boot)

in_PI.mean = apply(merBoot$t, 2, function(x) as.numeric(mean(x, na.rm=TRUE)))
in_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
in_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
in_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
in_PI.width = in_PI.upper - in_PI.lower

# 4.1.b condition on all explicitly
inStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=~(1|exp_no/clust_exp), allow.new.levels=TRUE)
  return(out$sim_1)
}
merBoot <- bootMer(lme1, inStats, nsim = n_boot, re.form=~(1|exp_no/clust_exp)) # This is sort of as expected
inexpclst_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
inexpclst_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
inexpclst_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
inexpclst_PI.width = inexpclst_PI.upper - inexpclst_PI.lower

# 4.2 Case of unobserved cluster in observed experiment
# 4.2.a
#"As only the experiment random effect is conditioned on, the cluster variation is resampled"
inStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=~(1|exp_no), allow.new.levels=TRUE)
  return(out$sim_1)
}
merBoot <- bootMer(lme1, inStats, nsim = n_boot, re.form=~(1|exp_no), seed=seed_boot) # This is sort of as expected
inexp_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
inexp_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
inexp_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
inexp_PI.width = inexp_PI.upper - inexp_PI.lower
### 4.2.a.2 Try to condition on exp_no for a new exp
# This seems to produce roughtly what I woant although the mean is not at the population mean.
PredData <- data.frame(treatment=unique(dat$treatment), exp_no=1, clust_exp=200)
inStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=~(1|exp_no), allow.new.levels=TRUE)
  return(out$sim_1)
}
merBoot <- bootMer(lme1, inStats, nsim = n_boot, re.form=~(1|exp_no), seed=seed_boot) # This is sort of as expected
inexp0_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
inexp0_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
inexp0_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
inexp0_PI.width = inexp0_PI.upper - inexp0_PI.lower

# 4.2.b Equivalently by passing in  data with values of cluster or experiment 
# that are not in the original data set, the conditioning has no effect 
# XXX In fact though this seems to produce the same results as exp_no=1 and clust_exp=1
PredData <- data.frame(treatment=unique(dat$treatment), exp_no=1, clust_exp=200)
outStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=NULL, allow.new.levels=TRUE) # re.form shouldn't matter?
  return(out$sim_1)
}
merBoot <- bootMer(lme1, outStats, nsim = n_boot, re.form = NULL, seed=seed_boot) 
inexp2_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
inexp2_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
inexp2_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
inexp2_PI.width = inexp2_PI.upper - inexp2_PI.lower

# 4.2.c Custom resampling based on estimate of variance of cluster
# bootstrap within experiment cluster to cluster variation
# This seems to produce too wide of a PI... Why?
#"In this case I think predict predicts the population mean, and as the model is refit for each
#bootstrap sample, then the estimate of clust_var is different for each bootstrap sample, accounting
#for uncertainty in its exact value."

outStats <- function(model){
  var_pars <- as.data.frame(VarCorr(model))
  clust_var <- var_pars[1,]$sdcor # var_pars$clust_exp:exp_no
  pred <- predict(model, newdata=PredData, re.form=NA, allow.new.levels=TRUE)
  out = pred + rnorm(1)*clust_var
  return(out)
}
merBoot <- bootMer(lme1, outStats,  nsim = n_boot, re.form = NA, seed=seed_boot)
clustvar_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
clustvar_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
clustvar_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
clustvar_PI.width = clustvar_PI.upper - clustvar_PI.lower

# 4.2.d Take CI for sigma_cluster and use Normal quantiles
df <- as.data.frame(cint_res)
sigma_clst_upr <- df[1,2]
clust_sd  <- df[1,2]
clust_plt <- plt_dat
clust_plt$lwr <- plt_dat$median-2*clust_sd
clust_plt$upr =plt_dat$median+2*clust_sd
clust_plt$width = clust_plt$upr - clust_plt$lwr

# 4.2.e predictInterval estimate

new_df_in <- data.frame(treatment=unique(dat$treatment), exp_no=1, clust_exp=5)
preds_in <- predictInterval(lme1, n.sims=5000, level=0.95, stat='median', which='full', newdata=new_df_in, include.resid.var = FALSE)
preds_in <- cbind(new_df_in, preds_in)
preds_in$treatment <-unclass(preds_in$treatment)
diff <- preds_in$fit - preds_out$fit
preds_in$fit <- preds_out$fit
preds_in$lwr <- preds_in$lwr - diff 
preds_in$upr <- preds_in$upr - diff 

### 4.2. PLOTS

diff <- CI.median - inexp_PI.median
inexp_PI.lower <- inexp_PI.lower + diff
inexp_PI.upper <- inexp_PI.upper  + diff
diff <- CI.median - in_PI.median
in_PI.lower <- in_PI.lower + diff
in_PI.upper <- in_PI.upper  + diff
diff <- CI.median - inexpclst_PI.median
inexpclst_PI.lower <- inexpclst_PI.lower + diff
inexpclst_PI.upper <- inexpclst_PI.upper + diff

ggplot(NULL, aes(x=plt_dat$treatment, y=CI.median), ) +
  geom_line()+
  geom_ribbon(data=preds_in, aes(ymin=lwr, ymax=upr, color='pI-in'), fill=NA)+
  geom_ribbon(aes(ymin=in_PI.lower, ymax=in_PI.upper, color='In'), alpha=1, fill=NA)+
  geom_ribbon(aes(ymin=inexp_PI.lower, ymax=inexp_PI.upper, color='InExp'), fill=NA)+
  geom_ribbon(aes(ymin=inexp0_PI.lower, ymax=inexp0_PI.upper, color='InExp0'), fill=NA)+
  geom_ribbon(aes(ymin=inexp2_PI.lower, ymax=inexp2_PI.upper, color='InExp2'), fill=NA)+
  geom_ribbon(aes(ymin=inexpclst_PI.lower, ymax=inexpclst_PI.upper, color='InExpCls'), fill=NA)+
  geom_ribbon(aes(ymin=clustvar_PI.lower, ymax=clustvar_PI.upper, color='Clst'), fill=NA )+
  geom_ribbon(data=clust_plt, aes(ymin=lwr, ymax=upr,color='Asy.Clst'), linetype=5, fill=NA)+
  geom_ribbon(data=clust_plt_true, aes(ymin=lwr, ymax=upr,color='True.Clst'), linetype=5, fill=NA)+
  scale_color_brewer(palette="Spectral")

### END 4.2. PLOTS

# 4.3 Case of unobserved experiment but neglect the cluster variance
#"Is this possible with bootMer and simulate directly? Since cluster is nested in 
#experiment any cluster in a new experiment is by definition also new, so 
#conditioning on cluster would have no effect if resimulating the random effects."


# 4.3.a Generate an out of sample prediction interval. 
#"I expect this to be too wide as both cluster and experiment should be resampled."

PredData <- data.frame(treatment=unique(dat$treatment), exp_no=200, clust_exp=200)
outStats <- function(model){
  out <- simulate(model, newdata=PredData, re.form=NA, allow.new.levels=TRUE) # re.form shouldn't matter?
  return(out$sim_1)
}
merBoot <- bootMer(lme1, outStats, nsim = n_boot, re.form = NA, seed=seed_boot) 
out_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
out_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
out_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
out_PI.width = out_PI.upper - out_PI.lower

# 4.3.b Take CI for sigma_cluster and use Normal quantiles
df <- as.data.frame(cint_res)
sigma_exp_upr <- df[2,2]
exp_sd <- df[2,2]
exp_plt <- plt_dat
exp_plt$lwr <- plt_dat$median-2*exp_sd
exp_plt$upr =plt_dat$median+2*exp_sd
exp_plt$width = exp_plt$upr - exp_plt$lwr


# 4.3.c bootstrap experiment to experiment variation of mean values from estimated sigma_exp
#I would expect this to be correct apart from my uncertainty of how residual errors are handled.
#In this case I think predict predicts the population mean, and as the model is refit for each
#bootstrap sample, then the estimate of exp_var is different for each bootstrap sample, accounting
#for uncertainty in it's exact value.

outStats <- function(model){
  var_pars <- as.data.frame(VarCorr(model))
  exp_var <-  var_pars[2,]$sdcor # TODO (more R way to access)
  pred <- predict(model, newdata=PredData, re.form=NA, allow.new.levels=TRUE) # Fixed effect mean
  out = pred + rnorm(1)*exp_var # Sample an experiment level random intercept
  return(out)
}
merBoot <- bootMer(lme1, outStats, nsim = n_boot, re.form = NA, seed=seed_boot)
expvar_PI.median = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.5, na.rm=TRUE)))
expvar_PI.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
expvar_PI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
expvar_PI.width = expvar_PI.upper - expvar_PI.lower

# 4.3.d use predictInterval
# As I understand, this doesn't really resample the random effects, so it's not expected to be correct.
# It is strange that both seen and unseen data give more or less the same results

new_df_out <- data.frame(treatment=unique(dat$treatment), exp_no=200, clust_exp=200)
preds_out <- predictInterval(lme1, n.sims=5000, level=0.95, stat='median', which='full', newdata=new_df_out, include.resid.var = False)
preds_out <- cbind(new_df_out, preds_out)
preds_out$treatment <-unclass(preds_out$treatment) 

### 4.3 PLOTS
ggplot(NULL, aes(x=plt_dat$treatment, y=CI.median)) +
  geom_line()+
  geom_ribbon(aes(ymin=preds_out$lwr, ymax=preds_out$upr, color='pI-exp'), fill=NA) + 
  geom_ribbon(aes(ymin=out_PI.lower, ymax=out_PI.upper, color='OutExp'),  fill=NA )+
  geom_ribbon(aes(ymin=expvar_PI.lower, ymax=expvar_PI.upper, color='Exp'), fill=NA)+
  geom_ribbon(data=exp_plt, aes(ymin=lwr, ymax=upr, color='Asy.Exp'), linetype=5, fill=NA, )+
  geom_ribbon(data=exp_plt_true, aes(ymin=lwr, ymax=upr, color='True.Exp'), linetype=5, fill=NA, )+
  scale_color_brewer(palette="Dark2")
### END 4.3 PLOTS