[R-sig-ME] Test of random effect in lme4
Anna Renwick
anna.renwick at bto.org
Thu Mar 11 16:30:31 CET 2010
There has been a lot of discussion previously whether we should remove
random effects based on LRT. The reason is that you added the random effect
based on your study design and whether it is significant or not it should
remain in there. I am not sure there is any definite rule and maybe it
depends on your study and personal view point.
Dr Anna R. Renwick
Research Ecologist
British Trust for Ornithology,
The Nunnery,
Thetford,
Norfolk,
IP24 2PU,
UK
Tel: +44 (0)1842 750050; Fax: +44 (0)1842 750030
-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Bolker
Sent: 11 March 2010 14:50
To: Luciano La Sala
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] 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
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list