# How to obtain the p-value (check significance) of an effect in a lme4 mixed model?

I use lme4 in R to fit the mixed model

lmer(value~status+(1|experiment)))


where value is continuous, status and experiment are factors, and I get

Linear mixed model fit by REML
Formula: value ~ status + (1 | experiment)
AIC   BIC logLik deviance REMLdev
29.1 46.98 -9.548    5.911    19.1
Random effects:
Groups     Name        Variance Std.Dev.
experiment (Intercept) 0.065526 0.25598
Residual               0.053029 0.23028
Number of obs: 264, groups: experiment, 10

Fixed effects:
Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
(Intr) statsD
statusD -0.204
statusR -0.193  0.476


How can I know that the effect of status is significant? R reports only $t$-values and not $p$-values.

• Going by the answers provided to this question, one wonders in what actually OP is interested here: testing coefficients against a null (vanilla $t$-test one does in regular linear regression against a null $H_0:\beta=\beta_{\text{null}}$), or testing for minimization of variance (the $F$-test we get from the many types of ANOVA). Those two aim at different things. An enlightening answer, while not about mixed-effects models, is found here. – Firebug Mar 25 '17 at 20:02

There is a lot of information on this topic at the GLMM FAQ. However, in your particular case, I would suggest using

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)


because you don't need any of the stuff that lmer offers (higher speed, handling of crossed random effects, GLMMs ...). lme should give you exactly the same coefficient and variance estimates but will also compute df and p-values for you (which do make sense in a "classical" design such as you appear to have). You may also want to consider the random term ~status|experiment (allowing for variation of status effects across blocks, or equivalently including a status-by-experiment interaction). Posters above are also correct that your t statistics are so large that your p-value will definitely be <0.05, but I can imagine you would like "real" p-values.

• I don't know about this answer. lmer could just as easily report the same kinds of p-values but doesn't for valid reasons. I guess it's the comment that there are any "real" p-values here that bugs me. You could argue that you can find one possible cutoff, and that any reasonable cutoff is passed. But you can't argue there's a real p-value. – John Aug 5 '12 at 2:45
• For a classical design (balanced, nested, etc.) I think I can indeed argue that there's a real p-vaue, i.e. a probability of getting an estimate of beta of an observed magnitude or greater if the null hypothesis (beta=0) were false ... lme4 doesn't provide these denominator df, I believe, because it's harder to detect in general from an lme4 model structure when the model specified is one where some heuristic for computing a classical denominator df would work ... – Ben Bolker Aug 10 '12 at 17:07
• try summary(m1) instead (I use this with nlme package) – jena Sep 13 '17 at 17:33

You could use the package lmerTest. You just install/load it and the lmer models get extended. So e.g.

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)


would give you results with p-values. If p-values are the right indication is a little bit disputed, but if you want to have them, this is the way to get them.

If you can handle abandoning p-values (and you should), you can compute a likelihood ratio that would represent the weight of evidence for the effect of status via:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
formula = value ~ (1|experiment) + status
, REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
formula = value ~ (1|experiment)
, REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))

• Note that likelihood ratios are asymptotic, i.e. don't account for uncertainty in the estimate of the residual variance ... – Ben Bolker Feb 17 '12 at 0:07
• I am interested in your last line. What is the interpretation of the result? Are there sources I can take a look at on that? – mguzmann Jul 10 '15 at 10:24

The issue is that the calculation of p-values for these models is not trivial, see dicussion here so the authors of the lme4 package have purposely chosen not to include p-values in the output. You may find a method of calculating these, but they will not necessarily be correct.

Consider what you're asking. If you just want to know if the overall p-value for the effect of status passes some some sort of arbitrary cutoff value, like 0.05, then that's easy. First, you want to find out the overall effect. You could get that from anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)


Now you have an F value. You can take that and look it up in some F tables. Just pick the lowest possible denom. degrees of freedom. The cutoff there is going to be around 20. Your F may be larger than that but I could be wrong. Even if it's not, look at the number of degrees of freedom from a conventional ANOVA calculation here using the number of experiments you have. Sticking that value in you're down to about 5 for a cutoff. Now you easily pass it in your study. The 'true' df for your model will be something higher than that because you're modelling every data point as opposed to aggregate values that an ANOVA would model.

If you actually want an exact p-value there's no such thing unless you're willing to make a theoretical statement about it. If you read Pinheiro & Bates (2001, and perhaps some more books on the subject... see other links in these answers) and you come away with an argument for a specific df then you could use that. But you're not actually looking for an exact p-value anyway. I mention this because you therefore shouldn't report an exact p-value, only that your cutoff is passed.

You should really consider the Mike Lawrence answer because the whole idea of just sticking with a pass point for p-values as the final and most important information to extract from your data is generally misguided (but might not be in your case since we don't really have enough information to know). Mike is using a pet version of LR calculation that is interesting, but it may be hard to find a lot of documentation on it. If you look into model selection and interpretation using AIC you may like it.

Edit: This method is no longer supported in newer versions of lme4. Use the lmerTest package as suggested in this answer by pbx101.

There is a post on the R list by lme4's author for why p-values are not displayed. He suggests using MCMC samples instead, which you do using the pvals.fnc from the languageR package:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)


See http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf for an example and details.

• what package contains pvals.fnc? – John Aug 5 '12 at 2:13
• – Jeff Aug 5 '12 at 2:54
• OK, just thought it should be in the question as well as the fact that the method of calculating them used in pvals.fnc is the same one that is recommended in the post you refer to. – John Aug 5 '12 at 7:11
• I think that the code should be model=lmer(...) (or mod <- lmer(...)). – smillig Aug 5 '12 at 9:10
• lme4 no longer supports this. This post could be updated to spare people having to find this out like I just did. – timothy.s.lau Mar 11 '15 at 19:39

Are you interested in knowing if the combined effect of status has a significant effect on value? If so, you can use the Anova function in the car package (not to be confused with the anova function in base R).

dat <- data.frame(
experiment = sample(c("A","B","C","D"), 264, replace=TRUE),
status = sample(c("D","R","A"), 264, replace=TRUE),
value = runif(264)
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)


Have a look at ?Anova after loading the car package.

• Any idea how car::Anova() avoids the sticky issues surrounding computation of p-values that Michelle links? – Mike Lawrence Feb 16 '12 at 20:29
• I don't, but my guess is that it avoids the sticky issues by ignoring them! Having reread the original post, I feel that I might have misunderstood the question. If the OP wants exact p-values for the fixed effects parameters, s/he is in trouble. But if the OP just wants to know if they're significant, I think the t-values are larger than any uncertainty in how the exact p-value would be calculated. (In other words, they are significant.) – smillig Feb 16 '12 at 20:47
• I think it was definitely a good idea to redirect to an ANOVA calculation to find out the overall effect of stats but I'm not sure finessing the p-values is good. The regular anova command will give you F's. – John Aug 5 '12 at 2:17
• I think this is a bit stickier than apparent. Running ANOVA's is valid when you want to minimize variance, but from the question wording I think OP wants to establish the marginal effect of variables, i.e. test coefficients against a null. – Firebug Mar 25 '17 at 19:59

The function pvals.fnc is not longer supported by lme4. Using the package lmerTest package, it is possible to use other method to calculate the p-value, such as the Kenward-Roger's approximations

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")