sigma() for mmrm - what does it exactly do?

[Generalized]

I'd like to calculate Coheh's d effect sizes over a set of emmeans (changes-from-baseline at subsequent timepoints in a non-RCT) obtained from a mmrm model.

The emmeans package has a helper function for it, e.g. eff_size(em_mmrm, sigma=...., method = "identity", edf=.....)
While I don't care of the edf (no inference will be presented), I'd like to get the pooled residual SD.

Unfortunately, I cannot share the data but I can provide some details for easier discussion.
It's a trivial model for change scores (Diff) over time (Timepoint). Each patient was examined several times.

I fit the same models using mmrm and nlme::gls() for comparisons. The two calls are equivalent.

mmrm: m_mmrm <- mmrm(Diff ~ Timepoint + us(Timepoint|PatientId), data=shifts)
nlme::gls(): m_gls <- gls(Diff ~ Timepoint, correlation = corSymm(form=~as.numeric(Timepoint)|PatientId),weights = varIdent(form =~1|Timepoint), data=shifts)

Let's check: overall coefficients:

> coef(summary(m_mmrm))
                   Estimate Std. Error       df   t value     Pr(>|t|)
(Intercept)      -0.3872661 0.11689589 67.04247 -3.312915 1.490878e-03
TimepointWeek 2  -0.2642027 0.07852287 65.05314 -3.364659 1.289003e-03
TimepointWeek 4  -0.2881937 0.11980281 67.98458 -2.405567 1.887351e-02
TimepointWeek 8  -0.5540765 0.11411082 66.92805 -4.855600 7.538807e-06
TimepointWeek 16 -0.7978211 0.11918492 65.72696 -6.693977 5.744121e-09
TimepointWeek 56 -0.9926554 0.13902561 59.76488 -7.140090 1.470701e-09

> coef(summary(m_gls))
                      Value  Std.Error   t-value      p-value
(Intercept)      -0.3872639 0.11688626 -3.313169 1.018523e-03
TimepointWeek 2  -0.2642035 0.07852278 -3.364673 8.510453e-04
TimepointWeek 4  -0.2881979 0.11980576 -2.405543 1.666520e-02
TimepointWeek 8  -0.5540805 0.11411320 -4.855534 1.810919e-06
TimepointWeek 16 -0.7978210 0.11919124 -6.693621 8.638107e-11
TimepointWeek 56 -0.9926457 0.13903559 -7.139508 5.422987e-12

... and the LS-means with Satterthwaite DoF:

# "natural" :)
> emmeans(m_mmrm, specs = ~Timepoint)
 Timepoint emmean    SE   df lower.CL upper.CL
 Week 1    -0.387 0.117 67.0   -0.621   -0.154
 Week 2    -0.651 0.132 69.5   -0.915   -0.388
 Week 4    -0.675 0.150 69.6   -0.974   -0.377
 Week 8    -0.941 0.143 69.6   -1.226   -0.657
 Week 16   -1.185 0.154 70.7   -1.492   -0.878
 Week 56   -1.380 0.154 60.7   -1.689   -1.071

# approximate (for gls) after several minutes of waiting :)
> emmeans(m_gls, specs = ~Timepoint)
Analytical Satterthwaite method not available; using appx-satterthwaite
 Timepoint emmean    SE   df lower.CL upper.CL
 Week 1    -0.387 0.117 66.6   -0.621   -0.154
 Week 2    -0.651 0.132 67.8   -0.915   -0.388
 Week 4    -0.675 0.150 69.5   -0.974   -0.377
 Week 8    -0.941 0.143 69.3   -1.226   -0.657
 Week 16   -1.185 0.154 71.6   -1.492   -0.878
 Week 56   -1.380 0.154 62.2   -1.689   -1.071

Degrees-of-freedom method: appx-satterthwaite 
Confidence level used: 0.95 

Lets also compare the residuals (they vary a little, but at this level it's negligible):

> all(sort(residuals(m_mmrm)) - sort(residuals(m_gls)) < 0.0001)
[1] TRUE

OK. Now let's find the residual SD (called eternally in R "residual standard error"...)

> sigma(m_gls)
[1] 1.18694 ~ 1.19

which is approximately close to:
> sd(residuals(m_gls)) * (357-1)/357
[1] 1.157078 ~ 1.16

which is approximately close to:
> mean(residuals(m_gls, type = "response") / residuals(m_gls, type = "pearson"))
[1] 1.148747 ~ 1.15

and we obtain the same for mmrm:
> mean(residuals(m_mmrm, type = "response") / residuals(m_mmrm, type = "pearson"))
[1] 1.148747 ~ 1.15

But sigma(m_mmrm) gives:

> sigma(m_mmrm)
[1] 1.490217 ~ 1.5

OK, not that much different, but I'm curious why both "sigmas" return different values for (asymptotically) same residuals.
I guess mmrm employs some internals from glmmTMB (I found glmmTMB:::sigma.glmmTMB).

It alters the effect sizes then:

# df.errors to make it much faster; I don't need inference this time.
> as.data.frame(eff_size(emmeans(m_gls, specs = ~Timepoint, mode="df.error"), sigma(m_gls), method = "identity", edf=60))[, 1:2]
  contrast effect.size
1   Week 1  -0.3262707
2   Week 2  -0.5488628
3   Week 4  -0.5690781
4   Week 8  -0.7930849
5  Week 16  -0.9984369
6  Week 56  -1.1625771

sigma used for effect sizes: 1.187 
Degrees-of-freedom method: inherited from df.error when re-gridding 
Confidence level used: 0.95 

> as.data.frame(eff_size(emmeans(m_mmrm, specs = ~Timepoint), sigma(m_mmrm), method = "identity", edf=60))[, 1:2]
  contrast effect.size
1   Week 1  -0.2598723
2   Week 2  -0.4371637
3   Week 4  -0.4532628
4   Week 8  -0.6316816
5  Week 16  -0.7952448
6  Week 56  -0.9259870

sigma used for effect sizes: 1.49 
Confidence level used: 0.95 

[danielinteractive]

Thanks @Generalized for the question - honestly the generic sigma() function was not on our radar so far at all!
I would recommend that for now you implement that yourself manually so you know what is going on. I can see that currently the default method does something for mmrm objects but that is probably not what you want to do for the MMRMs:

library(mmrm)
fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
sigma(fit)
# 2.537345
methods("sigma")
# sigma.default* sigma.glm*     sigma.gls*     sigma.lme*     sigma.lmList*  sigma.mlm*    
stats:::sigma.default
# function (object, use.fallback = TRUE, ...) 
# sqrt(deviance(object, ...)/(nobs(object, use.fallback = use.fallback) - 
#.   sum(!is.na(coef(object)))))

Here we have deviance() working specifically for mmrm objects, but not nobs() which will just give you the overall number of observations. I would guess the residual based formula you used in your example will work differently because the Pearson residuals will take into account time point specific standard deviation estimates.

The gls specific method uses something precomputed stored in the object:

nlme:::sigma.gls
# function (object, ...) 
# object$sigma

and one would need to find out from https://github.com/cran/nlme/blob/master/R/lme.R how exactly that is defined I guess.