CODEWITHSUNDEEP
rlme4p-value
SofiG
SofiG

Reputation: 33

How do I find the p-value for my random effect in my linear mixed effect model?

I am running the following line of code in R:

model = lme(divedepth ~ oarea, random=~1|deployid, data=GDataTimes, method="REML")
summary(model)

and I am seeing this result:

Linear mixed-effects model fit by REML
  Data: GDataTimes 
      AIC     BIC   logLik
  2512718 2512791 -1256352

Random effects:

 Formula: ~1 | deployid
        (Intercept) Residual
StdDev:    9.426598 63.50004

Fixed effects:  divedepth ~ oarea 
                Value Std.Error     DF   t-value p-value
(Intercept) 25.549003  3.171766 225541  8.055135  0.0000
oarea2      12.619669  0.828729 225541 15.227734  0.0000
oarea3       1.095290  0.979873 225541  1.117787  0.2637
oarea4       0.852045  0.492100 225541  1.731447  0.0834
oarea5       2.441955  0.587300 225541  4.157933  0.0000

[snip]

Number of Observations: 225554
Number of Groups: 9

However, I cannot find the p-value for the random variable: deployID. How can I see this value?

Upvotes: 0

Views: 4505

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226162

As stated in the comments, there is stuff about significance tests of random effects in the GLMM FAQ. You should definitely consider:

  • why you are really interested in the p-value (it's not never of interest, but it's an unusual case)
  • the fact that the likelihood ratio test is extremely conservative for testing variance parameters (in this case it gives a p-value that's 2x too large)

Here's an example that shows that the lme() fit and the corresponding lm() model without the random effect have commensurate log-likelihoods (i.e., they're computed in a comparable way) and can be compared with anova():

Load packages and simulate data (with zero random effect variance)

library(lme4)
library(nlme)
set.seed(101)
dd <- data.frame(x = rnorm(120), f = factor(rep(1:3, 40)))
dd$y <- simulate(~ x + (1|f),
                 newdata = dd,
                 newparams = list(beta = rep(1, 2),
                                  theta = 0,
                                  sigma = 1))[[1]]

Fit models (note that you cannot compare a model fitted with REML to a model without random effects).

m1 <- lme(y ~ x , random = ~ 1 | f, data = dd, method = "ML")
m0 <- lm(y ~ x, data = dd)

Test:

anova(m1, m0)
##    Model df      AIC      BIC    logLik   Test      L.Ratio p-value
## m1     1  4 328.4261 339.5761 -160.2131                            
## m0     2  3 326.4261 334.7886 -160.2131 1 vs 2 6.622332e-08  0.9998

Here the test correctly identifies that the two models are identical and gives a p-value of 1.

If you use lme4::lmer instead of lme you have some other, more accurate (but slower) options (RLRsim and PBmodcomp packages for simulation-based tests): see the GLMM FAQ.

Upvotes: 7

Related Questions