What?

🙂	When fitting linear mixed models, the `glmmTMB` package is nice, because it feels like `lme4` but additionally allows for several covariance structures for the random effects.
😞	However, this only works for random terms on the G-side of the model and not for the error term (= R-side).
🤓	But there is a workaround! We can fix the error variance to be 0 and thus force all the variance into the G-side. If we do so and also add a random *“pseudo error term”* that mimics the original error term, we have created a model that essentially leads to identical results.
😎	Once this is clear, we can use all the available variance structures on our pseudo error term!

How?

Remove original error variance from R-side via dispformula = ~ 0.
Add pseudo error variance to G-side via random term.
- In the simple case, this can be done via + (1 | unit) where unit is a factor with as many levels as there are observations in the dataset.

Basically go from here…

#
#
  
StandErrMod <- glmmTMB(
  y ~ 
    fixedeffects +
    (1 | randomeffects),
  

  REML = TRUE,
  data = dat
)

…to here!

dat <- dat %>%
  mutate(unit = as.factor(1:n()))

PseudErrMod <- glmmTMB(
  y ~
    fixedeffects +
    (1 | randomeffects) +
    (1 | unit),      # Pseudo Err
  dispformula = ~ 0, # ErrVar = 0
  REML = TRUE,
  data = dat
)

Really?

Here are two examples for the pseudo-error approach based on this dataset:

dat <- agridat::mcconway.turnip %>% 
  mutate(unit = 1:n()) %>% 
  mutate_at(vars(density, unit), as.factor)

Example 1 We assume no special variance structure for the error term and thus independent & identically distributed errors. We compare two glmmTMB models - one with a standard error term and one with a pseudo-error term.
Example 2 We assume a diagonal variance structure for the factor date and thus allow for different heterogeneous error variances / heteroscedascity per group for the effect levels of date. We compare a pseudo-error-glmmTMB model with nlme::lme() model.

Example 1: iid

StandErrMod <- glmmTMB(
  yield ~ 
    gen*date*density +
    (1 | block),
  

  REML = TRUE,
  data = dat
)

PseudErrMod <- glmmTMB(
  yield ~
    gen*date*density +
    (1 | block) +
    (1 | unit),      # Pseudo Err
  dispformula = ~ 0, # ErrVar = 0
  REML = TRUE,
  data = dat
)

## Warning in fitTMB(TMBStruc): Model convergence problem; false convergence (8).
## See vignette('troubleshooting')

The variance component estimates of both models are very similar. As expected the Residual variance in the pseudo-error model is 0 and instead there is an additional variance for the random unit term.

mixedup::extract_vc(StandErrMod)

group	variance	sd	sd_2.5	sd_97.5	var_prop
block	2.812	1.677	0.635	4.431	0.227
Residual	9.591	3.097			0.773

mixedup::extract_vc(PseudErrMod)

group	variance	sd	sd_2.5	sd_97.5	var_prop
block	2.810	1.676	0.634	4.430	0.227
unit	9.592	3.097	2.519	3.808	0.773
Residual	0.000	0.000			0.000

The Likelihood/AIC values are identical for both models.

AICcmodavg::aictab(list(StandErrMod, PseudErrMod), 
                   c("StandErrMod", "PseudErrMod"), 
                   second.ord = FALSE)

Modnames	K	AIC	Delta_AIC	LL
StandErrMod	18	308.1365	0	-136.0683
PseudErrMod	18	308.1365	0	-136.0683

Example 2: diag

diag_glmmTMB <- glmmTMB(
  yield ~ 
    gen*date*density +
    (1 | block) +
    diag(date + 0 | unit),
  dispformula = ~ 0,
  REML = TRUE,
  data = dat
)

diag_lme <- lme(
  yield ~ 
    gen*date*density, 
  random  = ~ 1 | block,
  weights = varIdent(form =  ~ 1 | date),
  
  
  data    = dat
)

The variance component estimates are similar enough.

effect	grp	variance
block	(Intercept)	1.596600
unit	date21Aug1990	4.305625
unit	date28Aug1990	15.507844
Residual		0.000000

effect	grp	variance
block	(Intercept)	1.596602
Residual	21Aug1990	4.306075
Residual	28Aug1990	15.505124

Furthermore, even though we are comparing across different model classes, we also get similar model fit statistics.

logLik	AIC	BIC
-131.9636	301.9272	342.9459

logLik	AIC	BIC
-131.9636	301.9271	337.48

More Details!

Check out the GitHub issue I’ve written on this topic at the glmmTMB repository, where Ben Bolker - the author of this R-package - replied. Here are some key points:

This works only for gaussian mixed models and thus not for generalized mixed models!
As is the case in Example 1, we sometimes get a false convergence warning for the pseudo-error-model but not for the standard model.
Actually, dispformula=~0 does not fix the residual variance to be 0, but to be a small non-zero value.
- Because of this, the variance component estimates will never be exactly identical.
- At present it is set to sqrt(.Machine$double.eps), which is the squareroot of the smallest possible positive floating-point number.
- Ben Bolker commented that “One piece of low-hanging fruit would be to allow the small non-zero value to be user-settable via glmmTMBControl”.
Possible covariance structures
- Find all possible structures here
- Note that cs is heterogeneous compound symmetry and there is no homogeneous compound symmetry!
- As can be seen in example 2, we need a + 0 in the diag(date + 0 | unit), since leaving it out would by default lead to estimating not only the desired heterogeneous variances, but an additional overall variance.
- Currently, we cannot have Kronecker product / direct product / varComb() as variance structure, as confirmed in this GitHub issue.

Mooore Details!

If you are wondering how to extract the variance component estimates as I did for example 2 and you are mad that I did not show the code, click here to find the R-code of this document.
Check out the chapters on this website to see more/other uses of this approach.

Please feel free to contact us about any of this!

schmidtpaul1989@outlook.com