1

I'm using the emmeans package with a negative-binomial model implemented using the glmmTMB package. I'm trying to bias adjust my backtransformed emmeans per the workflow illustrated here: https://cran.r-project.org/web/packages/emmeans/vignettes/transformations.html#bias-adj

As I've understood, I need the bias adjustment to be based on the variance of the random effects which can be obtained with VarCorr(nbinom_mod)

The problem is, I can't figure out how to extract that Std.Dev. value

library(glmmTMB) library(emmeans) library(dplyr) data(cbpp, package="lme4") nbinom_mod <- glmmTMB(incidence ~ period + (1|herd), data = cbpp, family = nbinom2) # Here's what I want to do but it fails since VarCorr doesn't return a number, it returns an object nbinom_em <- emmeans(nbinom_mod, ~ period, bias.adjust = T, sigma = VarCorr(nbinom_mod), type = "response")

So I've tried to extract data from VarCorr(nbinom_mod)and failed as such:

> class(VarCorr(nbinom_mod)) [1] "VarCorr.glmmTMB" > > typeof(VarCorr(nbinom_mod)) [1] "list" # This didn't work > unlist(VarCorr(nbinom_mod)) cond.herd 8.186695e-09 > VarCorr(nbinom_mod) Conditional model: Groups Name Std.Dev. herd (Intercept) 9.048e-05 > VarCorr(nbinom_mod)$cond $herd (Intercept) (Intercept) 8.186695e-09 attr(,"stddev") (Intercept) 9.048036e-05 attr(,"correlation") (Intercept) (Intercept) 1 attr(,"blockCode") us 1 attr(,"sc") [1] 1.626599 attr(,"useSc") [1] FALSE > VarCorr(nbinom_mod)$cond$herd (Intercept) (Intercept) 8.186695e-09 attr(,"stddev") (Intercept) 9.048036e-05 attr(,"correlation") (Intercept) (Intercept) 1 attr(,"blockCode") us 1 # Still didnt work > VarCorr(nbinom_mod)$cond$herd[1] [1] 8.186695e-09 > VarCorr(nbinom_mod)$cond$herd[2] [1] NA
Improve this question

3 Answers 3

Reset to default
3 
c(sqrt(VarCorr(nbinom_mod)$cond$herd))

should get the standard deviation associated with this random effect.

  • $cond extracts the conditional component
  • $herd extracts the component associate with this particular random effect (there could be more than one)
  • sqrt() is obvious (we need the std dev, not the variance)
  • c() drops a lot of harmless but potentially confusing extra information
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

1

Since the object could be seen as a list with different metadata stored through atrr(), you can also try with the next line of code

sapply(VarCorr(nbinom_mod)$cond, function(x) attr(x, "stddev"))

if you also just want the names of the groups you can try

names(VarCorr(nbinom_mod)$cond)

Sometimes it generates to me an error when I try to apply everything directly on VarCorr(nbinom_mod), so an alternative option is to generate an intermediate object to extract the $cond, for example

cor_random_mod <-VarCorr(nbinom_mod) names_cor <- names(cor_random_mod$cond) print(names_cor) random_sdevs <- sapply(cor_random_mod$cond, function(x) attr(x, "stddev")) print(random_sdevs)
Improve this answer

Comments

1

Just for the sake of "completeness", you can extract the different variances using insight::get_variance():

library(glmmTMB) library(emmeans) #> Welcome to emmeans. #> Caution: You lose important information if you filter this package's results. #> See '? untidy' data(cbpp, package = "lme4") nbinom_mod <- glmmTMB(incidence ~ period + (1 | herd), data = cbpp, family = nbinom2) # check for singularity performance::check_singularity(nbinom_mod) #> [1] TRUE # no random effects variances returned, because essentially singular fit insight::get_variance(nbinom_mod) #> Warning: Can't compute random effect variances. Some variance components equal #> zero. Your model may suffer from singularity (see `?lme4::isSingular` #> and `?performance::check_singularity`). #> Decrease the `tolerance` level to force the calculation of random effect #> variances, or impose priors on your random effects parameters (using #> packages like `brms` or `glmmTMB`). #> $var.fixed #> [1] 0.5567141 #> #> $var.dispersion #> [1] 0 #> #> $var.intercept #> herd #> 8.186695e-09 # change tolerance, we now get random effects variance insight::get_variance(nbinom_mod, tolerance = 1e-10) #> $var.fixed #> [1] 0.5567141 #> #> $var.random #> [1] 8.186695e-09 #> #> $var.residual #> [1] 0.7795249 #> #> $var.distribution #> [1] 0.7795249 #> #> $var.dispersion #> [1] 0 #> #> $var.intercept #> herd #> 8.186695e-09 # sepcific variance - get random effects directly insight::get_variance_random(nbinom_mod, tolerance = 1e-10) #> var.random #> 8.186695e-09 # equals: c(VarCorr(nbinom_mod)$cond$herd) #> [1] 8.186695e-09 emmeans( nbinom_mod, ~period, bias.adjust = TRUE, sigma = sqrt(insight::get_variance_random(nbinom_mod, tolerance = 1e-10)), type = "response" ) #> period response SE df asymp.LCL asymp.UCL #> 1 4.067 0.974 Inf 2.543 6.50 #> 2 1.214 0.389 Inf 0.648 2.28 #> 3 1.000 0.340 Inf 0.514 1.95 #> 4 0.538 0.235 Inf 0.229 1.27 #> #> Confidence level used: 0.95 #> Intervals are back-transformed from the log scale #> Bias adjustment applied based on sigma = 9.048e-05

Created on 2025-05-16 with reprex v2.1.1

Note that the random effects variance is close to zero, so bias adjustment has literally no effect here.

Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.