2
$\begingroup$

I want to get confidence intervals around modelled data from a lmer model. I found that Bootmer is the way to go. There seem to be 3 ways to do this:

1.parametrically resampling both the “spherical” random effects u and the i.i.d. errors ϵ (use.u = FALSE, default, seems te lead to relatively large CI)

2.treating the random effects as fixed and parametrically resampling the i.i.d. errors (use.u = TRUE, relatively small CI)

3.treating the random effects as fixed and semi-parametrically resampling the i.i.d. errors from the distribution of residuals.

I did not find anywere when to use which specification. I am not interested in prediction, only in the modeled data (e.g. I have a model with Condition, time and time x Condition as predictors and I want the modelled data with CI of the effect per condition at the different timepoints.)

Thanks in advance for your help!

$\endgroup$
3
  • $\begingroup$ Welcome to CV! From the question it is clear your problem is statistical in nature (+1), but the title suggests a software issue (which would be off-topic). I therefore recommend changing the title to avoid close votes. $\endgroup$ Jul 15, 2019 at 14:11
  • $\begingroup$ Try using the confint function the stats package (pre-loaded), which accepts lmer objects. $\endgroup$ Jul 15, 2019 at 14:39
  • $\begingroup$ Thanks for the suggestion. However, confint function gives confidence interval for the model parameters, while I am looking for the confidence intervals of the modeled data. $\endgroup$
    – Chris H.
    Jul 16, 2019 at 7:33

1 Answer 1

0
$\begingroup$

The nice thing about bootstrap is it's (mostly) quite intuitive. So the question of whether to set resample the random effects really depends on whether you would view that future cases also would resample the random effects.

  1. Perhaps the most traditional definition of a random factor is one where the categories of our random factor are just a randomly selected subset from a broader population of possibilities and that our interest is in generalising to that broader population. For example, when the random factor models the participants in your study: you are probably not that interested in these participants, so much as in people in general. A rerun of the experiment would resample the participants, and your bootstrap should too. To resample, we set use.u=FALSE.

  2. On the other hand, say we have a categorical predictor with a fixed set of categories that will always stay the same. We might model that as a fixed factor. But in some cases, we might enter it as a random factor, modelling all the categories as one big bell-curve blob of possibilities and drawing on the advantages of random effects, such as shrinkage. This can be particularly advantageous if there are a lot of categories, or some categories have only a few datapoints. Since the categories are fixed and will never change, they should not change in the bootstrap either. Do not resample, set use.u=TRUE.

These descriptions sound like the traditional definitions of a random and fixed effects, respectively. So, it's interesting to see that the default standard errors output by fixed and mixed models (at least, on a simple example) output roughly the same numbers as the two use.u possibilities. The code at the bottom produces this output (each is the standard error of the intercept from a different approach):

 use.u.TRUE use.u.FALSE fixed.model mixed.model 
 0.06300992  0.34607292  0.06618366  0.34576152

One interesting side-note from this is that the default way we generate inferences on random effects may give extremely conservative inferences in case 2 above.

For cases where we have a mixture of types 1 and 2 above, it might be nice to be able to switch use.u differently per random effect.

Warning

It's really important not to set use.u=TRUE unless you're sure that's what you want. As you noted, the SEs (and CIs) can be wildly different and that's because the question you are asking is completely different.

For example, in an analysis I just ran, a number of replicates were taken from each beaker of "stuff". I model the beakers as a random factor. The replicates from each beaker all come out roughly the same, so the residual standard deviation is actually tiny: most of the sample-to-sample variation is encoded in the beaker random effect. If we don't model that variation beaker to beaker, you would think there was almost no noise! And consequently, your inferences would be WAY too optimistic, based on the assumption that things always come out the same.

Code

library(lme4)

n.per.group <- c(100,100)
nsim <- 500

set.seed(1)

num.groups <- length(n.per.group)
u <- rnorm(num.groups)
data <- data.frame(
  y         = unlist(mapply(rnorm,     n=n.per.group, mean=u,         SIMPLIFY=FALSE)),
  group.num = unlist(mapply(rep,   times=n.per.group, x=1:num.groups, SIMPLIFY=FALSE))
)

data$group.num <- factor(data$group.num)

# This sets up the fixed model so that the intercept models the grand mean
#  (i.e. same as in the mixed model)
contrasts(data$group.num) <- contr.sum(levels(data$group.num))

fixed.model <- lm(  y~   group.num,  data=data)
mixed.model <- lmer(y~(1|group.num), data=data)

standard.errors <- c(
  use.u.TRUE  = sd(bootMer(mixed.model, fixef, use.u=TRUE,  nsim=nsim)$t),
  use.u.FALSE = sd(bootMer(mixed.model, fixef, use.u=FALSE, nsim=nsim)$t),
  
  fixed.model = summary(fixed.model)$coefficients["(Intercept)", "Std. Error"],
  mixed.model = summary(mixed.model)$coefficients["(Intercept)", "Std. Error"]
)
print(standard.errors)

```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.