96
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • 3
    $\begingroup$ 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. $\endgroup$
    – Firebug
    Mar 25, 2017 at 20:02
  • $\begingroup$ library(car) Anova(mod, test="F") $\endgroup$
    – RichardBJ
    Jul 22, 2021 at 15:04

11 Answers 11

Reset to default
81
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • 5
    $\begingroup$ 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. $\endgroup$
    – John
    Aug 5, 2012 at 2:45
  • 15
    $\begingroup$ 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 ... $\endgroup$
    – Ben Bolker
    Aug 10, 2012 at 17:07
  • 2
    $\begingroup$ try summary(m1) instead (I use this with nlme package) $\endgroup$
    – jena
    Sep 13, 2017 at 17:33
  • $\begingroup$ What about the variable "experiment"? The output is a variance, how do can know if that variable is significant? How can I know if any of the levels of that categorical variable is significant? $\endgroup$
    – skan
    Nov 27, 2020 at 15:03
  • 2
    $\begingroup$ For testing RE variances: bbolker.github.io/mixedmodels-misc/… . For testing levels of the categorical variable: under standard frequentist, you can't perform significance tests on these effects, because they are random variables rather than parameters. (Feel free to look around CrossValidated some more for answers to these questions and post a new question if you don't find anything satisfactory ...) $\endgroup$
    – Ben Bolker
    Nov 28, 2020 at 0:32
67
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ What about the variable "experiment"? The output is a variance, how do can know if that variable is significant? How can I know if any of the levels of that categorical variable is significant? $\endgroup$
    – skan
    Nov 27, 2020 at 15:03
  • $\begingroup$ see my answer to the same comment above. $\endgroup$
    – Ben Bolker
    Nov 28, 2020 at 0:32
  • $\begingroup$ worked nicely for me! $\endgroup$
    – bob
    Feb 9 at 6:15
33
$\begingroup$

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))
Improve this answer
$\endgroup$
14
  • 19
    $\begingroup$ Note that likelihood ratios are asymptotic, i.e. don't account for uncertainty in the estimate of the residual variance ... $\endgroup$
    – Ben Bolker
    Feb 17, 2012 at 0:07
  • 5
    $\begingroup$ 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? $\endgroup$
    – mguzmann
    Jul 10, 2015 at 10:24
  • 1
    $\begingroup$ @RosaMaria negative values reflect evidence in favor of the null $\endgroup$
    – Mike Lawrence
    Mar 31, 2022 at 13:02
  • 1
    $\begingroup$ @RosaMaria In my answer, the two models were constructed such that they are identical with the exception that the restricted model did not include the status predictor. In that case $H_0 : \Beta_{status} = 0$ and $H_1 : \Beta_{status} != 0$ Note also that in the intervening decade I converted to Bayes, and the “estimation” flavour therein, where one would fit/sample a single model that includes a possible effect of status and thereby yield a posterior distribution on the magnitude of such effect. The brms package in R is a great entry-point to Bayes for those already familiar with lme4 $\endgroup$
    – Mike Lawrence
    Apr 7, 2022 at 11:12
  • 1
    $\begingroup$ @RosaMaria hm, as you wrote them, the restricted and unrestricted models share the same fixed-effects structure and differ only in the random-effects structure such that the unrestricted model has by-subject variability in both intercepts and Time effects, with correlation therein, while the restricted model has only by-subject variability in intercepts. So the resulting likelihood ratio reflects evidence associated with two quantities (the correlation and the by-subject variance in time) either both being precisely zero or both being non-zero. Rarely do such multi-quantity LRs have usefulness $\endgroup$
    – Mike Lawrence
    Apr 8, 2022 at 16:47
17
$\begingroup$

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.

Improve this answer
$\endgroup$
0
11
$\begingroup$

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.

Improve this answer
$\endgroup$
9
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 5
    $\begingroup$ lme4 no longer supports this. This post could be updated to spare people having to find this out like I just did. $\endgroup$
    – timothy.s.lau
    Mar 11, 2015 at 19:39
5
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Any idea how car::Anova() avoids the sticky issues surrounding computation of p-values that Michelle links? $\endgroup$
    – Mike Lawrence
    Feb 16, 2012 at 20:29
  • $\begingroup$ 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.) $\endgroup$
    – smillig
    Feb 16, 2012 at 20:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – John
    Aug 5, 2012 at 2:17
  • $\begingroup$ 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. $\endgroup$
    – Firebug
    Mar 25, 2017 at 19:59
2
$\begingroup$

Simply loading the afex package will print the p-values in the output of the lmer function from the lme4 package (you don't need to be using afex; just load it):

library(lme4)  #for mixed model
library(afex)  #for p-values

This will automatically add a p-value column to the output of the lmer(yourmodel) for the fixed effects.

Improve this answer
$\endgroup$
1
$\begingroup$

You could use parameters::p_value() to get the p-values. I found it to be very useful.

Improve this answer
$\endgroup$
0
$\begingroup$

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")
Improve this answer
$\endgroup$
0
$\begingroup$

You can install the jtools package and use the summ function on the output model to get the p-value of the fixed effects.

library(lme4)
library(jtools)
model <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
summ(model)

FIXED EFFECTS:
---------------------------------------------------------
                      Est.   S.E.   t val.    d.f.      p
----------------- -------- ------ -------- ------- ------
(Intercept)         251.41   6.82    36.84   17.00   0.00
Days                 10.47   1.55     6.77   17.00   0.00
---------------------------------------------------------

p values calculated using Satterthwaite d.f.
```
Improve this answer
$\endgroup$

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