cc 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

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)

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

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

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?

the bottom line

Temporal correlation

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

Multivariate models

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.