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)
dotwhisker::dwplot()
scales parameters by 2SD by default.)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
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
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))
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
?emmeans::emmeans
lmerTest
, car::Anova
, afex
…nlme::ACF
correlation=
argument in lme
or MASS::glmmPQL
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)
glmmPQL
, spaMM
, INLA
, … see Dormann et al. (2007) (now a bit out of date)MCMCglmm
lme4
etc.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.