Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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

Thanks a lot

share|improve this question

6 Answers 6

up vote 27 down vote accepted

There is a lot of information on this topic at http://glmm.wikidot.com/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.

share|improve this answer
1  
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
2  
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

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
)
#compute a model where the effect of status is 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))
share|improve this answer
4  
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

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.

share|improve this answer

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.

share|improve this answer
    
what package contains pvals.fnc? –  John Aug 5 '12 at 2:13
    
1  
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
2  
I think that the code should be model=lmer(...) (or mod <- lmer(...)). –  smillig Aug 5 '12 at 9:10
    
good call on both; post edited –  Jeff Aug 5 '12 at 13:50

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.

share|improve this answer
    
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
1  
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

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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