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. – smillig Apr 24 '14 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? – user3288202 Apr 24 '14 at 16:04
  • 1
    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. – smillig Apr 25 '14 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? – user3288202 Apr 25 '14 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. – smillig Apr 25 '14 at 8:12
up vote 44 down vote accepted
+50

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.

  • Thank very much smilig. So I might not report effect size for the overall parameters. – user3288202 Apr 25 '14 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. – russellpierce Jul 7 '14 at 14:57
  • 3
    +6, impressively clear, comprehensive, & thoroughly annotated. – gung Jul 7 '14 at 16:38
  • 1
    additionally, you could have a look at the afex-package and especially the mixed-function. see here – beginneR Jul 7 '14 at 20:26

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.

  • 1
    In my case, I prefer pairwise comparison using lsmeans as it gives me all pairs of contrasts including p values. If I use lmerTest I will have to run model six times with different baselines to see all pairs of contrasts. – user3288202 Sep 5 '14 at 15:26

In regard to the question of significance, Luke (2016) Evaluating significance in linear mixed-effects models in R reports that the most optimal (certainly most conservative) test is p values based on Kenward-Roger approximation for degrees of freedom (in lmer). With large samples, p values based on the likelihood ratio (through anova()) are just as good.

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)

  • 2
    +1 Thanks for sharing this link. I will just briefly comment that Kenward-Roger approximation is available in the lmerTest package. – amoeba Oct 13 '16 at 20:37

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.