How to get an "overall" p-value and effect size for a categorical factor in a mixed model (lme4)?

Question

I would like to get a p-value and an effect size of an independent categorical variable (with several levels) -- that is "overall" and not for each level separately, as is the normal output from lme4 in R. It is just like the thing people report when running an ANOVA.

How can I get this?

Answer 1

Both of the concepts you mention (p-values and effect sizes of linear mixed models) have inherent issues. With respect to effect size, quoting Doug Bates, the original author of lme4,

Assuming that one wants to define an $R^2$ measure, I think an argument could be made for treating the penalized residual sum of squares from a linear mixed model in the same way that we consider the residual sum of squares from a linear model. Or one could use just the residual sum of squares without the penalty or the minimum residual sum of squares obtainable from a given set of terms, which corresponds to an infinite precision matrix. I don't know, really. It depends on what you are trying to characterize.

For more information, you can look at this thread, this thread, and this message. Basically, the issue is that there is not an agreed upon method for the inclusion and decomposition of the variance from the random effects in the model. However, there are a few standards that are used. If you have a look at the Wiki set up for/by the r-sig-mixed-models mailing list, there are a couple of approaches listed.

One of the suggested methods looks at the correlation between the fitted and the observed values. This can be implemented in R as suggested by Jarrett Byrnes in one of those threads:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

So for example, say we estimate the following linear mixed model:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

We can calculate the effect size using the function defined above:

r2.corr.mer(fm1)
# [1] 0.0160841

A similar alternative is recommended in a paper by Ronghui Xu, referred to as $\Omega^{2}_{0}$, and can be calculated in R simply:

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

With respect to the p-values, this is a much more contentious issue (at least in the R/lme4 community). See the discussions in the questions here, here, and here among many others. Referencing the Wiki page again, there are a few approaches to test hypotheses on effects in linear mixed models. Listed from "worst to best" (according to the authors of the Wiki page which I believe includes Doug Bates as well as Ben Bolker who contributes here a lot):

• Wald Z-tests
• For balanced, nested LMMs where df can be computed: Wald t-tests
• Likelihood ratio test, either by setting up the model so that the parameter can be isolated/dropped (via anova or drop1), or via computing likelihood profiles
• MCMC or parametric bootstrap confidence intervals

They recommend the Markov chain Monte Carlo sampling approach and also list a number of possibilities to implement this from pseudo and fully Bayesian approaches, listed below.

Pseudo-Bayesian:

• Post-hoc sampling, typically (1) assuming flat priors and (2) starting from the MLE, possibly using the approximate variance-covariance estimate to choose a candidate distribution
• Via mcmcsamp (if available for your problem: i.e. LMMs with simple random effects — not GLMMs or complex random effects)
  Via pvals.fnc in the languageR package, a wrapper for mcmcsamp)
• In AD Model Builder, possibly via the glmmADMB package (use the mcmc=TRUE option) or the R2admb package (write your own model definition in AD Model Builder), or outside of R
• Via the sim function from the arm package (simulates the posterior only for the beta (fixed-effect) coefficients

Fully Bayesian approaches:

• Via the MCMCglmm package
• Using glmmBUGS (a WinBUGS wrapper/R interface)
• Using JAGS/WinBUGS/OpenBUGS etc., via the rjags/r2jags/R2WinBUGS/BRugs packages

For the sake of illustration to show what this might look like, below is an MCMCglmm estimated using the MCMCglmm package which you will see yields similar results as the above model and has some kind of Bayesian p-values:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

I hope this helps somewhat. I think the best advice for somebody starting out with linear mixed models and trying to estimate them in R is to read the Wiki faqs from where most of this information was drawn. It is an excellent resource for all sorts of mixed effects themes from basic to advanced and from modelling to plotting.

Answer 2

I use the lmerTest package. This conveniently includes an estimation of the p-value in the anova() output for my MLM analyses, but does not give an effect size for the reasons given in other posts here.

Answer 3

In regard to calculating significance (p) values, Luke (2017) Evaluating significance in linear mixed-effects models in R reports that the optimal method is either the Kenward-Roger or Satterthwaite approximation for degrees of freedom (available in R packages such as lmerTest or afex).

Abstract

Mixed-effects models are being used ever more frequently in the analysis of experimental data. However, in the lme4 package in R the standards for evaluating significance of fixed effects in these models (i.e., obtaining p-values) are somewhat vague. There are good reasons for this, but as researchers who are using these models are required in many cases to report p-values, some method for evaluating the significance of the model output is needed. This paper reports the results of simulations showing that the two most common methods for evaluating significance, using likelihood ratio tests and applying the z distribution to the Wald t values from the model output (t-as-z), are somewhat anti-conservative, especially for smaller sample sizes. Other methods for evaluating significance, including parametric bootstrapping and the Kenward-Roger and Satterthwaite approximations for degrees of freedom, were also evaluated. The results of these simulations suggest that Type 1 error rates are closest to .05 when models are fitted using REML and p-values are derived using the Kenward-Roger or Satterthwaite approximations, as these approximations both produced acceptable Type 1 error rates even for smaller samples.

(emphasis added)