bootstrap variance component in a mixed linear model

Hi Kevin,

In the pbkrtest package there is function PBmodcomp for calculating p-values using parametric bootstrap and in the implementation of this function I have also observed the same phenomenon. I guess that what happens is simply that some of the bootstrap samples lead to "singularities" somewhere in the estimation algorithms. My solution in PBmodcomp has been to use the suppressWarnings() function...

Cheers
S?ren
 

-----Oprindelig meddelelse-----
Fra: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] P? vegne af Kevin Spring
Sendt: 5. januar 2012 00:22
Til: r-sig-mixed-models at r-project.org
Emne: [R-sig-ME] bootstrap variance component in a mixed linear model

Hi, everyone.  I don't know what is going wrong with my bootstrap.  Could anyone help me find out what is wrong?

I am trying to bootstrap the variance component in a mixed linear model.

My statistic for the bootstrap is the following:

varcomp <- function ( formula, data, indices ) {

    d <- data[indices,] #sample for boot
    fit <- lmer(formula, data=d) #linear model
    return (attr (VarCorr(fit), "sc")^2) #output random effects 
residual var }

The formula for the model is: log( y ) ~ ( 1 | a:b ) + a, where '*b*' is a
random effect nested within '*a';* which is a fixed effect.

When I run the linear model on its own it works fine, but when I try to run
the bootstrap I get warning messages of false convergence.

Warning messages:

1: In mer_finalize(ans) : false convergence (8)
2: In mer_finalize(ans) : false convergence (8)
3: In mer_finalize(ans) : false convergence (8)
4: In mer_finalize(ans) : false convergence (8)
5: In mer_finalize(ans) : false convergence (8)
6: In mer_finalize(ans) : false convergence (8)

The example data I used can be downloaded at:
http://www.mediafire.com/?77xb24rmul90k5d

The R code I actually did:

library(boot)
library(lme4)
#
#Load data
#
fp1 <- read.csv("desktop/p1f.csv")
#Set as factors
fp1$Cell.line <- as.factor(fp1$Cell.line)
fp1$DNA.extract <- as.factor(fp1$DNA.extract)
#test mixed model
lm1 <- lmer(formula=log(CNPC)~(1|Cell.line:DNA.extract)+Cell.line,data=fp1)
attr (VarCorr(lm1), "sc")
#
##Bootstrap statistic
#
varcomp <- function ( formula, data, indices ) {
    d <- data[indices,] #sample for boot
    fit <- lmer(formula, data=d) #linear model
    return (attr (VarCorr(fit), "sc")) #output random effects residual
standard deviation
}
##Bootstrap with 1000 replicates
fp1.boot <- boot ( data = fp1, statistic=varcomp, R=1000,
formula=log(CNPC)~(1|Cell.line:DNA.extract)+Cell.line)


_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models