How does emmeans calculate degrees of freedom for glmmTMB?

Question

I have a model that I'd really love to use glmmTMB for (beta distribution, mixed model, heterogeneous variance that the dispersion equation can help with).

I found that the degrees of freedom for the emmeans estimates and contrasts were nowhere near what I would expect (very large, ~380). I ran a simpler model in lme using the same overall structure and allowing for heterogenous variance with a logit transformation. The degrees of freedom in the emmeans estimates and contrasts were ~8. I read that, for lme, emmeans uses a containment method to estimate the degrees of freedom which can result in underestimation of the degrees of freedom. (https://cran.r-project.org/web/packages/emmeans/vignettes/models.html)

I'm aware that I shouldn't be expecting the denominator degrees of freedom that I'd calculate using expected mean square rules and I'd like to know how emmeans and glmmtmb are figuring out the df. I really just want to make sure that, on some level, the glmmTMB model is still accounting for the fact that my experiment is split-plot or is longitudinal or has subsamples. I also want to make sure that this information is being passed on to emmeans.

Questions

1. How exactly is emmeans calculating the df from a glmmTMB model?
2. Is it reasonable? Does it properly account for split-plots or sub-sampling?
3. If not, would an acceptable workaround be to combine the df from an lme model with the estimates and standard errors of a (more or less) corresponding glmmTMB+emmeans procedure? (http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have)

Here is a worked example that shows the difference in degrees of freedom.

library(nlme)
library(glmmTMB)

# lme (straight from Pinheiro and Bates page 47)

fm1Oats <- lme( yield ~ ordered(nitro) * Variety, data = Oats,
                random =~1| Block/Variety )

# Here's the anova() output just to show that the split plot structure 
# is correctly being interpreted with regards to the denominator degrees of freedom.

anova(fm1Oats)
                       numDF denDF   F-value p-value
(Intercept)                1    45 245.14299  <.0001
ordered(nitro)             3    45  37.68562  <.0001
Variety                    2    10   1.48534  0.2724
ordered(nitro):Variety     6    45   0.30282  0.9322

emmeans(fm1Oats,~nitro|Variety)

Variety = Golden Rain:
 nitro emmean   SE df lower.CL upper.CL
   0.0   80.0 9.11  5     56.6    103.4
   0.2   98.5 9.11  5     75.1    121.9
   0.4  114.7 9.11  5     91.3    138.1
   0.6  124.8 9.11  5    101.4    148.2

Variety = Marvellous:
 nitro emmean   SE df lower.CL upper.CL
   0.0   86.7 9.11  5     63.3    110.1
   0.2  108.5 9.11  5     85.1    131.9
   0.4  117.2 9.11  5     93.8    140.6
   0.6  126.8 9.11  5    103.4    150.2

Variety = Victory:
 nitro emmean   SE df lower.CL upper.CL
   0.0   71.5 9.11  5     48.1     94.9
   0.2   89.7 9.11  5     66.3    113.1
   0.4  110.8 9.11  5     87.4    134.2
   0.6  118.5 9.11  5     95.1    141.9

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

# glmmTMB

glmm1Oats <- glmmTMB(yield ~ ordered(nitro) * Variety +(1|Block/Variety), data = Oats)
emmeans(glmm1Oats,~nitro|Variety)

Variety = Golden Rain:
 nitro emmean   SE df lower.CL upper.CL
   0.0   80.0 8.31 57     63.4     96.6
   0.2   98.5 8.31 57     81.9    115.1
   0.4  114.7 8.31 57     98.0    131.3
   0.6  124.8 8.31 57    108.2    141.5

Variety = Marvellous:
 nitro emmean   SE df lower.CL upper.CL
   0.0   86.7 8.31 57     70.0    103.3
   0.2  108.5 8.31 57     91.9    125.1
   0.4  117.2 8.31 57    100.5    133.8
   0.6  126.8 8.31 57    110.2    143.5

Variety = Victory:
 nitro emmean   SE df lower.CL upper.CL
   0.0   71.5 8.31 57     54.9     88.1
   0.2   89.7 8.31 57     73.0    106.3
   0.4  110.8 8.31 57     94.2    127.5
   0.6  118.5 8.31 57    101.9    135.1

Confidence level used: 0.95

# 5 vs. 57...is the glmmTMB version OK or does it need a modification?

Edit: aov(glmm1Oats) for d.f.

aov(glmm1Oats)
Call:
   aov(formula = glmm1Oats)

Terms:
                ordered(nitro)   Variety     Block ordered(nitro):Variety Variety:Block Residuals
Sum of Squares       20020.500  1786.361 15875.278                321.750      6013.306  7968.750
Deg. of Freedom              3         2         5                      6            10        45

Residual standard error: 13.30727
Estimated effects may be unbalanced

glmmTMB at least recognizes the right groups and d.f. at this level, but I'm not confident that the information is being used in the way I would want.

Thanks very much!

Answer 1

Looking at the code ...

> glmmTMB:::emm_basis.glmmTMB
function (object, trms, xlev, grid, component = "cond", 
    ...) 
{
    V <- as.matrix(vcov(object)[[component]])
    misc = list()
    if (family(object)$family == "gaussian") {
        dfargs = list(df = df.residual(object))
        dffun = function(k, dfargs) dfargs$df
    }
    else {
        dffun = function(k, dfargs) Inf
        dfargs = list()
    } 

    . . . (many lines omitted)

    dfargs = list(df = df.residual(object))

    . . .
}

The d.f. are in the dfargs member, which, interestingly, is defined twice, and the second time overrides the first.

Anyway, in your example:

> df.residual(glmm1Oats)
[1] 57

In turn, this is computed very simply:

> getS3method("df.residual", "glmmTMB")
function (object, ...) 
{
    nobs(object) - length(object$fit$par)
}

You can always override the default d.f., for example:

> emmeans(glmm1Oats, ~ nitro | Variety, df = 5)
Variety = Golden Rain:
 nitro emmean   SE df lower.CL upper.CL
   0.0   80.0 8.31  5     58.6    101.4
   0.2   98.5 8.31  5     77.1    119.9
   0.4  114.7 8.31  5     93.3    136.0
   0.6  124.8 8.31  5    103.5    146.2

Variety = Marvellous:
. . .

Degrees-of-freedom method: user-specified 
Confidence level used: 0.95

I am not sure how easily some more sophisticated d.f. method could be implemented. It may, for example, be possible to fix it up so the pkrbtest package would support it. One would have to ask the developer.