[R-sig-ME] bootstrap variance component in a mixed linear model

[Douglas Bates]

On Wed, Jan 4, 2012 at 5:22 PM, Kevin Spring <kevinjspring at gmail.com> wrote:
> 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)
>

Those warnings are coming from the optimizer used in lme4 (the same
one as in the R function nlminb).  Finding reliable code for
optimizing a nonlinear function subject to box constraints (in this
case, standard deviations being greater than zero) is not easy.
Eventually I gave up and created an implementation of the Nelder-Mead
simplex method (the version used in the optim function in R doesn't
handle constraints) which can be a bit slower but is also more
reliable.  This is used in the development version of lme4 called
lme4Eigen. (I know - "not ANOTHER development version".  I have a
rather severe "best is the enemy of the good" problem.)

For some GLMMs the version of glmer in lme4Eigen is both faster and
more reliable.

Ben has set a deadline of January for releasing what is now lme4Eigen
as lme4.  I hope he means the end of January and not the beginning of
January.

We would appreciate feedback on lme4Eigen.  Unfortunately you need to
compile it to be able to use it on Windows or on Mac OS X (sigh).  On
Windows there is a problem with compiling the 64-bit version which
traces back to something in Rcpp.  It may be possible to make a small
change in Rcpp to bypass that.  On Mac OS X the compilation of
RcppEigen, on which lme4Eigen depends, croaks because the decade-old
compiler that Apple still uses has bugs.  Apparently Apple is going to
replace gcc-4.2.1 with clang in a new release of Xcode so that logjam
may be freed up too.


> 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)
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models