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)
library(glmmTMB)

Effect size measures

Distinction between unscaled effect sizes (magnitude) and scaled effect sizes [e.g. Cohen’s \(d\)] (clarity)
Recommend: get comparable effect sizes by scaling continuous predictors by their standard deviation as recommended by Schielzeth (2010) (or by 2SD as suggested by Gelman (2008)). (Or use log-transformed predictors.) (dotwhisker::dwplot() scales parameters by 2SD by default.)
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 (marginal: population-level. conditional: taking random effects into account)
overall metrics 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 r2glmm::r2beta() package with method="nsj" (Nakagawa and Schielzeth 2013; Johnson 2014; Nakagawa, Johnson, and Schielzeth 2017) (also MuMIn::r.squaredGLMM(), …)

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)))

## boundary (singular) fit: see help('isSingular')

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  fYear 0.046    0.314    0.003

piecewiseSEM::rsquared(gmod_lme4_L)

##   Response  family link   method  Marginal Conditional
## 1   shells poisson  log trigamma 0.2279262   0.2279262

MuMIn::r.squaredGLMM(gmod_lme4_L)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

## Warning: the null model is correct only if all variables used by the original
## model remain unchanged.

##                   R2m         R2c
## delta     0.021324124 0.021324124
## lognormal 0.112656506 0.112656506
## trigamma  0.001291019 0.001291019

Graphical presentation of models

dotwhisker::dwplot()
emmeans
effects
sjPlot

Owls <- transform(Owls, Nest=reorder(Nest,SiblingNegotiation))
owls_nb1 <- glmmTMB(SiblingNegotiation ~ FoodTreatment*SexParent*ArrivalTime+
                        (1|Nest)+offset(log(BroodSize)),
                    family = nbinom1, zi = ~1, data=Owls)

dotwhisker::dwplot(owls_nb1,effects="fixed")+
    geom_vline(xintercept=0,lty=2)

s1 <- sjPlot::plot_model(owls_nb1, type="pred", title="")
s1$Nest <- s1$Nest + coord_flip()
s1$BroodSize <- NULL
sjPlot::plot_grid(s1)

## Warning in sjPlot::plot_grid(s1): Not enough tags labels in list. Using letters
## instead.

library(effects)

## Loading required package: carData

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

plot(allEffects(owls_nb1))

## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect

library(emmeans)
plot(emtrends(owls_nb1, ~SexParent|FoodTreatment, var="ArrivalTime"))

plot(emmeans(owls_nb1, ~SexParent|FoodTreatment, at=list(ArrivalTime=23)))

## NOTE: Results may be misleading due to involvement in interactions

Prediction

A quick-and-dirty approach to prediction (fixed-effect variability only):

## make prediction data frame
newdat <- with(Owls,
               expand.grid(SexParent=levels(SexParent),
                           FoodTreatment=levels(FoodTreatment),
                           ArrivalTime=seq(min(ArrivalTime),
                                           max(ArrivalTime),length=31),
                           BroodSize=mean(BroodSize)))
                           
## design matrix (fixed effects)
mm <- model.matrix(delete.response(terms(owls_nb1)),newdat)
## linear predictor
newdat$pred1 <- drop(mm %*% fixef(owls_nb1)$cond)
predvar <- diag(mm %*% vcov(owls_nb1)$cond %*% t(mm))
newdat$SE <- sqrt(predvar)
linkinv <- family(owls_nb1)$linkinv
newdat$lwr <- with(newdat,linkinv(pred1-2*SE))
newdat$upr <- with(newdat,linkinv(pred1+2*SE))
newdat$SiblingNegotiation <- linkinv(newdat$pred1)

newdat$Nest <- NA
newdat$BroodSize <- 1
predfun <- function(m) predict(m,data=newdat)
bb <- bootMer(owls_nb1, FUN=predfun, nsim=100)
## ugh ...
refit(owls_nb1,simulate(owls_nb1))

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 …

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
## 
##  Random Effects:
##   Level: Subject 
##                           lower       est.      upper
## sd((Intercept))        5.486184 14.8791959 40.3541841
## sd(Days)               2.645800  4.7598609  8.5631102
## cor((Intercept),Days) -0.997924  0.8967291  0.9999938
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.2789544 0.4870368 0.6514356
## 
##  Within-group standard error:
##    lower     est.    upper 
## 25.45047 30.50069 36.55305

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. https://doi.org/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. https://doi.org/10.1002/sim.3107.

Johnson, Paul C. D. 2014. “Extension of Nakagawa & Schielzeth’s R2GLMM to Random Slopes Models.” Methods in Ecology and Evolution 5 (9): 944–46. https://doi.org/10.1111/2041-210X.12225.

Nakagawa, Shinichi, Paul C. D. Johnson, and Holger Schielzeth. 2017. “The Coefficient of Determination R2 and Intra-Class Correlation Coefficient from Generalized Linear Mixed-Effects Models Revisited and Expanded.” Journal of the Royal Society Interface 14 (134): 20170213. https://doi.org/10.1098/rsif.2017.0213.

Nakagawa, Shinichi, and Holger Schielzeth. 2013. “A General and Simple Method for Obtaining R2 from Generalized Linear Mixed-Effects Models.” Methods in Ecology and Evolution 4 (2): 133–42. https://doi.org/10.1111/j.2041-210x.2012.00261.x.

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

Schielzeth, Holger. 2010. “Simple Means to Improve the Interpretability of Regression Coefficients.” Methods in Ecology and Evolution 1: 103–13. https://doi.org/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.