Test of random effect in lme4
Luciano La Sala wrote:
Dear R-list members, I am running mixed models using lme4 package. In model selection, terms were eliminated from a maximum model (with random intercept) to achieve a simpler model that retained only the significant main effects and interactions, using the Akaike information criterion. My final model includes three fixed factors plus random intercept. Then I perform a likelihood ratio test to test the significance of the random term.
In general testing significance of components in a model *after* model selection is dubious ... I agree with most of what Zuur et al say, but I am only comfortable with stepwise procedures as a relatively necessary evil for eliminating non-significant interaction terms to simplify interpretation of the remaining model.
However, because when testing on the boundary the p-value from the table is incorrect, I followed Zuur et al (2009) to get the corrected p-value by dividing the p value obtained by 2. Briefly, my best fit model consists of three main effects: Year (2006, 2007), Hatching Order (1st, 2nd, 3rd) and Sibling Competence (Present/Absent) plus NestID as random intercept. The modelled outcome is level of plasma proteins (continuous). I test the random effect (Nest ID), which has variance 2.1795e-16 and Std. Dev. 1.4763e-08 (see output). LRT yields a p-value of 0.00031 (0.00015 after dividing it by 2 as suggested). This would mean that adding a random effect Nest ID to the model is a significant improvement. However, random effect variance is near zero, which would indicate otherwise. In support of the non-significant random effect I think, coefficients and standard error are exactly the same for models with and without the RE, as seen in the outputs.
Your main problem is that the log-likelihoods returned by lm and lmer are **NOT COMPARABLE**. Sooner or later there should probably be a warning to that effect somewhere in the documentation ... You may be able to use the RLRsim package to solve your problem.
Q 1. In your opinion, should I trust this LRT with a small p-value and leave the random effect in my model, or follow the parsimony principle and eliminated it?
I would leave it in whether or not it is significant (and it's probably not). Note that as expected all the fixed effect parameters are estimated identically under lmer and lm ... the reason to drop it would be to have the convenience of not dealing with mixed effects at all.
Q 2. Is it possible, under certain conditions, to have a random effect with such low variance and still obtain a LTR p-value indicating that model fit is improved by it?
Unlikely at best, unless your response variable has a very small magnitude (e.g., you are comparing differences hummingbird weights across different diet treatments, and measuring them in units of petagrams)
Outputs for both models, with and without random effect, followed by LRT output: MIXED MODEL
full.1 <- lmer(TP10Diff~HatchOrder+Year+SibComp+(1|NestID))
Linear mixed model fit by REML
AIC BIC logLik deviance REMLdev
739.4 758.5 -362.7 738.5 725.4
Random effects:
Groups Name Variance Std.Dev.
NestID (Intercept) 2.1795e-16 1.4763e-08
Residual 4.4754e+01 6.6898e+00
Number of obs: 112, groups: NestID, 81
Fixed effects:
Estimate Std. Error t value
(Intercept) 6.959 1.078 6.453
HatchOrderSecond -1.130 2.472 -0.457
HatchOrderThird -12.483 3.514 -3.552
Year2007 7.157 1.299 5.509
SibCompPresente -2.120 2.641 -0.803
Correlation of Fixed Effects:
(Intr) HtchOS HtchOT Yr2007
HtchOrdrScn -0.219
HtchOrdrThr -0.154 0.677
Year2007 -0.676 0.016 0.020
SibCmpPrsnt 0.019 -0.816 -0.738 -0.079
MODEL WITHOUT RANDOM EFFECT
null.1 <- lm(TP10Diff~HatchOrder+Year+SibComp)
Residuals:
Min 1Q Median 3Q Max
-16.4597 -3.8812 -0.2394 4.1472 17.4203
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.959 1.078 6.453 3.28e-09 ***
HatchOrderSecond -1.130 2.472 -0.457 0.648649
HatchOrderThird -12.483 3.514 -3.552 0.000569 ***
Year2007 7.157 1.299 5.509 2.51e-07 ***
SibCompPresente -2.120 2.641 -0.803 0.424037
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.69 on 107 degrees of freedom
(155 observations deleted due to missingness)
Multiple R-squared: 0.3771, Adjusted R-squared: 0.3539
F-statistic: 16.2 on 4 and 107 DF, p-value: 2.117e-10
TEST OF SIGNIFICANCE FOR RANDOM TERM
as.numeric(2*(logLik(full.1)-logLik(null.1)))
[1] 13.01191
0.5*(1-pchisq(13.01191, 1))
[1] 0.0001547580 Thank you very much for previous assistance! Luciano
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Ben Bolker Associate professor, Biology Dep't, Univ. of Florida bolker at ufl.edu / people.biology.ufl.edu/bolker GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc