additional topics

Licensed under the Creative Commons attribution-noncommercial license. Please share & remix noncommercially, mentioning its origin.

library(lme4)
library(nlme)
library(dotwhisker)
library(broom)
library(broom.mixed)
library(cowplot)
library(r2glmm)

List of possible topics

effect size measures/\(R^2\)
type III sums of squares (mostly philosophical)
zero-inflation
temporal autocorrelation

type III sums of squares

There is a long-standing argument about the principle of marginality: when does it make sense to interpret/test the significance of a main effect in a model where that effect is also involved in an interaction?

main effect parameter value/significance is always evaluated for a particular value of all of the other parameters
when the interaction is present, that implies that the value/significance changes depending on where you test it.
for an interaction with a continuous variable, the parameter will always be significant eventually
the “average” condition might not make biological sense (what is halfway between male and female?)
if marginality is violated, parameter effects/SSQ depend on how the model is parameterized/where the baseline (zero value) of each parameter lies

the bottom line

need to choose a sensible “model center” (Schielzeth 2010)
center continuous variables by their mean, or by the minimum, or some sensible reference point
center factor variables by using sum-to-zero contrasts (weights equally across levels), or see “Weights” in ?emmeans::emmeans
type-3 ANOVAs are available from lmerTest, car::Anova, afex …

Effect size measures

I would generally recommend finding standardized effect sizes by scaling continuous predictors by their standard deviation as recommended by Schielzeth (2010) (or by 2SD as suggested by Gelman (2008))
R-squared measures:
- generalized linear models: have to decide whether to include sampling variance or not, and what measure of variability to use
- mixed models: have to decide how to account for/include variability at different levels
overall metric like \(\eta^2\):

Unfortunately this is currently not possible. I believe that most of these problems are also discussed in a recent Psych Methods paper (Rights and Sterba 2018) … The fact that calculating a global measure of model fit (such as R2) is already riddled with complications and that no simple single number can be found, should be a hint that doing so for a subset of the model parameters (i.e., main-effects or interactions) is even more difficult. Given this, I would not recommend to try finding a measure of standardized effect sizes for mixed models. (Singmann 2018)

see the relevant section of the glmm FAQ for packages etc. …

Example: using the r2beta() package with method="nsj" [nakagawa_general_2013;@johnson_extension_2014]

library(lme4)
load("../data/gopherdat2.RData")
Gdat <- transform(Gdat,fYear=factor(year))
gmod_lme4_L <- glmer(shells~prev+fYear+(1|Site),
                     offset=log(Area),
                     family=poisson,data=Gdat,
                     control=glmerControl(optimizer="bobyqa",
                                          check.conv.grad=.makeCC("warning",0.05)))
r2beta(gmod_lme4_L,method="nsj")

##      Effect   Rsq upper.CL lower.CL
## 1     Model 0.224    0.522    0.064
## 2      prev 0.201    0.471    0.019
## 3 fYear2005 0.046    0.278    0.000
## 4 fYear2006 0.016    0.213    0.000

Temporal correlation

Easier than spatial correlation!
Diagnose (without trying to incorporate it in the model) via nlme::ACF
use correlation= argument in lme or MASS::glmmPQL
temporal correlation options in: glmmTMB, brms, MCMCglmm, INLA …

m1 <- lme(Reaction~Days,
          random=~Days|Subject,
          sleepstudy)
plot(ACF(m1),alpha=0.05)

## this appears to fit, but makes the intercept/slope model singular
m1A <- update(m1,
              correlation=corAR1(),
              control=lmeControl(opt="optim"))
intervals(m1A)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                  lower      est.     upper
## (Intercept) 238.713471 252.24315 265.77283
## Days          7.441343  10.46701  13.49268
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Subject 
##                            lower      est.      upper
## sd((Intercept))        5.4839676 14.879196 40.3704921
## sd(Days)               2.6437661  4.759861  8.5696975
## cor((Intercept),Days) -0.9980981  0.896729  0.9999944
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.2790912 0.4870368 0.6513502
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##    lower     est.    upper 
## 25.45216 30.50069 36.55062

plot(ACF(m1A,resType="normalized"),alpha=0.05)

Spatial correlation

use smooth terms (see E. Gonzalez’s course …)
glmmPQL, spaMM, INLA, … see Dormann et al. (2007) (now a bit out of date)

Multivariate models

via implicit melting: MCMCglmm
via melting: lme4 etc.

References

Dormann, Carsten F., Jana M. McPherson, Miguel B. Araújo, Roger Bivand, Janine Bolliger, Gudrun Carl, Richard G. Davies, et al. 2007. “Methods to Account for Spatial Autocorrelation in the Analysis of Species Distributional Data: A Review.” Ecography 30 (5): 609–28. doi:10.1111/j.2007.0906-7590.05171.x.

Gelman, Andrew. 2008. “Scaling Regression Inputs by Dividing by Two Standard Deviations.” Statistics in Medicine 27 (15): 2865–73. doi:10.1002/sim.3107.

Rights, Jason D., and Sonya K. Sterba. 2018. “Quantifying Explained Variance in Multilevel Models: An Integrative Framework for Defining R-Squared Measures.” Psychological Methods. doi:10.1037%2Fmet0000184.

Schielzeth, Holger. 2010. “Simple Means to Improve the Interpretability of Regression Coefficients.” Methods in Ecology and Evolution 1: 103–13. doi:10.1111/j.2041-210X.2010.00012.x.

Singmann, Henrik. 2018. “Compute Effect Sizes for Mixed() Objects.” Afex: Analysis of Factorial EXperiments. https://afex.singmann.science/forums/topic/compute-effect-sizes-for-mixed-objects.