[R-sig-ME] Testing signif. of random effects by bootstrapping

[Luciano La Sala]

Dear list members, 

I am using lmer function (lme4 version 32) to fit some quite simple mixed
models. Dependent variable is chick's weight, fixed effects are HatchOrder
(3 levels), ClutchSize (3 levels), Year (2 levels) and significant two-way
interaction terms. I model Nest as a random intercept (chicks nested in
nest!). After fitting the full model I need to test the significance of the
random effect. In doing so, I've found to approaches based on Faraway's book
(Extending the linear model with R, 2006), namely likelihood ratio test and
parametric bootstrap approach. The later, I understand, would be used to
obtain a more accurate p value.  

Before presenting the outputs, here go my main questions. First off: 

1. I have fitted models where the RE should be removed according to a LRT (p
= 0.4284092), while after bootstrapping (1000 repetitions) the same model I
obtain p = 0 indicating a significant improvement of model fit by the RE.
Two opposite results. Is this normal?  

2. The bootstrap method yields p-values that are either 0 or 1 (no decimal
places), which seems a little odd to me. How can I ask R to provide, say, 4
significant digits for the bootstrapping p-value?

Following is a summary of the outputs from the model in question, with and
without the random effect, and bootstrapping right after them: 

FULL MODEL

> FULL <-
lmer(Weight~HatchOrder+ClutchSize+Year+HatchOrder*Year+ClutchSize*Year+SibCo
mp+(1|NestID))

Linear mixed model fit by REML 

   AIC   BIC logLik deviance REMLdev
 562.2 598.3 -268.1    551.1   536.2

Random effects:
 Groups   Name        Variance Std.Dev.
 NestID   (Intercept) 4.1513   2.0375  
 Residual             3.4385   1.8543  
Number of obs: 118, groups: NestID, 87

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               40.44583    0.79525   50.86
HatchOrderSecond          -0.05734    1.08618   -0.05
HatchOrderThird           -4.57739    2.01458   -2.27
ClutchSizeTwo             -0.68492    1.12169   -0.61
ClutchSizeThree           -0.79372    1.08271   -0.73
Year2007                  -0.35917    1.12465   -0.32
SibCompPresente           -1.67799    0.94786   -1.77
HatchOrderSecond:Year2007 -0.71230    1.00113   -0.71
HatchOrderThird:Year2007  -1.59001    1.87684   -0.85
ClutchSizeTwo:Year2007    -0.25690    1.66292   -0.15
ClutchSizeThree:Year2007  -2.14999    1.46019   -1.47


MODEL WITHOUT RANDOM INTERCEPT TERM: 

> REDUCED
<-lm(Weight~HatchOrder+ClutchSize+Year+HatchOrder*Year+ClutchSize*Year+SibCo
mp) 

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5723 -1.6628  0.3156  1.5917  6.4361 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               40.44583    0.80129  50.476   <2e-16 ***
HatchOrderSecond           0.04877    1.22858   0.040   0.9684    
HatchOrderThird           -4.17110    2.05523  -2.030   0.0449 *  
ClutchSizeTwo             -0.82366    1.11288  -0.740   0.4608    
ClutchSizeThree           -0.91305    1.07716  -0.848   0.3985    
Year2007                  -0.35917    1.13319  -0.317   0.7519    
SibCompPresente           -1.94834    1.06655  -1.827   0.0705 .  
HatchOrderSecond:Year2007 -0.52951    1.28902  -0.411   0.6821    
HatchOrderThird:Year2007  -1.49689    2.16455  -0.692   0.4907    
ClutchSizeTwo:Year2007    -0.24930    1.64309  -0.152   0.8797    
ClutchSizeThree:Year2007  -2.12061    1.45500  -1.457   0.1479    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 2.776 on 107 degrees of freedom
  (149 observations deleted due to missingness)
Multiple R-squared: 0.512,      Adjusted R-squared: 0.4664 
F-statistic: 11.23 on 10 and 107 DF,  p-value: 6.46e-13

LIKELIHOOD RATIO TEST 
as.numeric(2*(logLik(FULL)-logLik(REDUCED)))
[1] 28.0247

pchisq(28.0247,1,lower=FALSE) 
[1] 1.197768e-07

NOW THE BOOTSTRAPPING APPROACH

y <- simulate(REDUCED)
lrstat <- numeric(1000)
for(i in 1:1000){
y <- unlist(simulate(REDUCED))
REDUCED <-
lm(TMT10~HatchOrder+ClutchSize+Year+HatchOrder*Year+ClutchSize*Year+SibComp)
FULL <-
lmer(TMT10~HatchOrder+ClutchSize+Year+HatchOrder*Year+ClutchSize*Year+SibCom
p+(1|NestID))
lrstat[i] <- as.numeric(2*(logLik(FULL)-logLik(REDUCED)))}

mean(lrstat>28.0247)
[1] 0

3. Final question. Although LRT and bootstrap p-values seem to indicate the
same thing here, I would like to be able to get a p-value with some more
digits!   

Thank you very much in advance List! 

Luciano