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

[Luciano La Sala]

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. 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. 

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? 

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?    

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