Search code examples
rregressionemmeansstatistical-testglmmtmb

How to extract Std.Dev from VarCorr glmmTMB


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

Solution

  • 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