Is this the appropriate way to get standard errors on mixed model intercepts?

Appreciate your guidance on calculating standard errors from mixed models. There are some other questions on this subject, but mine is specifically for the random effect intercepts. I am not a statistician by training, so I am looking for an applied perspective.

It seems there are two strategies.

1. methods that trace the likelihood and/or do probability calculations like the beta method.
2. resampling methods like the bootstrap.

For calculating standard errors for fixed effects and variance components, I can find examples of both 1. and 2. and it seems both ways are valid.

For standard errors of random effect intercepts and slopes, I can't find examples. Is it the case that 1. is unavailable and I can only do 2.?

Now, here is a minimal example of bootstrapped standard errors for random effect intercepts. Is this the appropriate way to get SE on the Subject intercepts? Is there a different method you'd recommend?

And finally, is there a way to retain the names? I was careful to have my FUN argument return a named vector, (subjects "308", "309", etc.), but the result of bootMer returns a result with default names "t1*", "t2*", etc.

library(lme4)
fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)

get_intercepts <- function(model) {
# returns the Subject random effect intercepts as a named vector
int <- ranef(model)\$Subject
setNames(int[, "(Intercept)", drop = TRUE], rownames(int))
}

# test the function
get_intercepts(fm2) # test
#>         308         309         310         330         331         332
#>   1.5126648 -40.3738728 -39.1810279  24.5189244  22.9144470   9.2219759
#>         333         334         335         337         349         350
#>  17.1561243  -7.4517382   0.5787623  34.7679030 -25.7543312 -13.8650598
#>         351         352         369         370         371         372
#>   4.9159912  20.9290332   3.2586448 -26.4758468   0.9056510  12.4217547

bootMer(fm2, get_intercepts, nsim = 20, seed = 123)
#>
#> PARAMETRIC BOOTSTRAP
#>
#>
#> Call:
#> bootMer(x = fm2, FUN = get_intercepts, nsim = 20, seed = 123)
#>
#>
#> Bootstrap Statistics :
#>         original     bias    std. error
#> t1*    1.5126648  -5.042886    19.64321
#> t2*  -40.3738728  39.465677    14.76306
#> t3*  -39.1810279  39.504116    26.16501
#> t4*   24.5189244 -26.669108    23.11084
#> t5*   22.9144470 -17.033506    21.98820
#> t6*    9.2219759 -14.157263    15.48943
#> t7*   17.1561243 -16.433882    18.62850
#> t8*   -7.4517382   8.888964    24.77319
#> t9*    0.5787623   4.786413    20.75570
#> t10*  34.7679030 -42.369727    20.36696
#> t11* -25.7543312  28.378799    19.52419
#> t12* -13.8650598  16.680981    23.02667
#> t13*   4.9159912  -3.050187    21.21267
#> t14*  20.9290332 -23.140041    15.98342
#> t15*   3.2586448   1.612509    22.07355
#> t16* -26.4758468  22.991405    20.68792
#> t17*   0.9056510  -4.403906    23.93385
#> t18*  12.4217547 -10.008357    23.22543


Created on 2023-11-03 with reprex v2.0.2

• Would a random intercept satisfy the basic assumptions of bootstrapped intervals, such as that, as a statistic, they are pivotal? Commented Nov 3, 2023 at 19:27

I'm not sure that there is a statistically appropriate way to get standard errors of individual random intercepts/slopes in a mixed model. That would be treating the random effects as the equivalent of fixed effects.

Strictly, a mixed model estimates the (co)variances of the random effects, not their individual values (even though the model returns random-effect estimates for individuals). Estimates of variability in those (co)variances is typically more useful from a statistical standpoint, which I suppose you might accomplish by resampling cases in some way.

The bootMer() function you invoke is not, however, a case-resampling method. With default parameter settings, it simulates "new values of both the 'spherical' random effects u and the i.i.d. errors $$\epsilon$$, using rnorm() with parameters corresponding to the fitted model x." It's a parametric bootstrap.

So you aren't working with the original 18 individuals at all any more. For each resample, the function is simulating 18 new individual values from the model's estimate of the random-effect variance (the "'spherical' random effects u"). Notice that the sum of the "bias" in each individual intercept in bootMer() with the "original" value is very close to 0 in all cases, which is what you would expect if the resampled values came from a normal distribution with 0 mean.

• Thanks. I see your point about the parametric bootstrap. Should I do a traditional bootstrap instead? My categorical variables have unbalanced levels and some intercept estimates are likely very uncertain. It seems like a problem to present the intercepts without any estimate of standard errors. Commented Nov 3, 2023 at 20:17
• @Arthur I don't think that it's standard practice to report the individual intercept estimates (or other random-effect estimates) from a mixed model. The individual random intercepts in a mixed model aren't of direct concern. You aren't modeling these individuals; you are modeling a population from which these individuals (ideally) represent a random sample. You report the variance of the estimates, or its square root for the standard deviation among the estimates, to put this into the context of the population of interest.
– EdM
Commented Nov 3, 2023 at 21:25
• thanks again. My model is for prices of different products that exist within a hierarchy of product types. The prices are variable on the quantity sold (discount for large quantity), the region and country of sale, and the calendar date (due to inflation). It seems a perfect application of hierarchical mixed models! From this model, it also seems reasonable to provide an estimate of the price differences in (say) different countries and the corresponding standard errors. I could get these from OLS, but the large number of levels and lack of hierarchy would be a real shame. Commented Nov 3, 2023 at 21:41