[R-sig-ME] Test of random effect in lme4

[Ben Bolker]

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