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

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?

• What statistics do you want exactly? You can use the anova() function to get an anova table with linear mixed models just as with linear models. Apr 24, 2014 at 15:47
• I have tried anova() but it gives me Df, Sum Sq, Mean Sq, and F value. I don't see effect size and p value. Do you have any idea about this? Apr 24, 2014 at 16:04
• By effect size, do you mean something like an equivalent to $R^{2}$? With respect to p-values, there is a long and substantial debate around their estimation and around the implementation of them in lme4. Have a look at the discussion in this question for more details. Apr 25, 2014 at 7:39
• Thanks for the link, Smilig. Does that mean that because there is a problem with p value calculation, the effect size of factor in overall is also a problem? Apr 25, 2014 at 7:54
• They are not directly related issues. However, you should keep in mind that a linear mixed model does not behave exactly like a linear model without random effects so a measure that may be appropriate for the linear model does not necessarily generalize to mixed models. Apr 25, 2014 at 8:12

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) #  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)))) #  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. • Thank very much smilig. So I might not report effect size for the overall parameters. Apr 25, 2014 at 10:40 • ... and of course the resulting$r^2\$ will be biased towards sampling units (and it goes without saying conditions) with more observations. Jul 7, 2014 at 14:57
• +6, impressively clear, comprehensive, & thoroughly annotated. Jul 7, 2014 at 16:38
• additionally, you could have a look at the afex-package and especially the mixed-function. see here Jul 7, 2014 at 20:26
• Sorry for the silly question but what does it mean high values obtained from your r2.corr.mer ? In some cases I got 0.958351 or other 0.558351 ? Thanks May 20 at 10:42

In regard to calculating significance (p) values, Luke (2016) 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 with 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.

• +1 Thanks for sharing this link. I will just briefly comment that Kenward-Roger approximation is available in the lmerTest package. Oct 13, 2016 at 20:37
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.