How to retrieve standard errors from random effects in nlme?

Question

I have conducted a linear mixed effect model with the nlme package in R.

lmm.reg.slope <- lme(V1~ V2+V3+V4+V5+V6+V7+V8+V9+V10, data = data, 
    random = V2+V3+V4+V5+V6+V7+V8+V9+V10|regions, method = 'ML', 
    control = lmeControl(opt = "optim", msMaxIter=1000, 
    maxIter = 1000, msMaxEval = 1000))

When I look at the random effect, I get the difference between the individual intercept and the global intercept. I remember it as the standard deviation of the random effect could also be retrieve if I somehow changed it to a data frame or similar. But it does not work.

ranef(lmm.reg.slope)

      (Intercept)
AK    9.815204e-09
NY   -6.132803e-09
MIN   2.393367e-08
WIS  -1.884604e-08
CA    1.469633e-08
WAS  -2.454771e-09
MAS  -1.397460e-09
CT    7.225472e-09
FL   -1.694695e-08
IL   -1.233468e-09
OH   -2.637688e-08
IO    7.647110e-09
TX   -2.296820e-09
AZ    1.448242e-08
NC   -2.484795e-09
SC    3.697730e-10

How can I retrieve the standard error related to the random effects for each state?

Edit: I simply want to do this "We tested for differences in effect sizes among ecoregions using Tukey–Kramer post hoc analysis for multiple comparisons in the package 'emmeans'", as seen in this article. I thought by getting the random effect and the standard error that would be possible?

E.g

emmeans(lmm.reg.slope, pairwise ~ states$Intercept

contrast                      estimate     SE     df     t-ratio     p-value
AK$Intercept - NY$Intercept   0.97831      2.22   1      0.288       0.0137
AK$Intercept - MIN$Intercept  0.01038      0.96   1      0.01        0.5101
...

Answer 1

To test whether the intercept differs among 50 states, consider using likelihood ratio test between one model without state as either random or fixed effects and another model with random intercepts by state, illustrated https://rpubs.com/DKCH2020/578881.

If you do pairwise comparisons between each state pair, I agree with Roland's comment at How to retrieve standard errors from random effects in nlme? that it requires modeling state as fixed effects, meaning that there will be 50 - 1 = 49 coefficients estimated representing the difference in V1 between each state and a reference state. As a result, there could be 50*49/2 = 1225 pairs of differences, which makes the results intangible. Instead, you can designate one single state (e.g. your major study area) as the reference to reduce the number of pairwise comparisons to 49. In that case, you can use adjusted CIs to visualize the state comparisons. See Wright, T., Klein, M., & Wieczorek, J. (2019). A primer on visualizations for comparing populations, including the issue of overlapping confidence intervals. The American Statistician, 73(2), 165–178. https://doi.org/10.1080/00031305.2017.1392359.

To retrieve standard errors of standard deviations or variances of random effects, see https://stats.oarc.ucla.edu/r/faq/how-can-i-calculate-standard-errors-for-variance-components-from-mixed-models/. Because standard deviations and variances of random effects have skewed sample distribution, however, standard errors for these statistics can be misleading. Instead, reporting profile confidence intervals is encouraged. See https://github.com/lme4/lme4/issues/497.