[R-sig-ME] Specifying Crossed vs. Nested Factors with Lmer and bootstrapping downstream statistics

[Douglas Bates]

On Wed, Sep 9, 2009 at 9:58 AM, Andrew Dolman <andydolman at gmail.com> wrote:
> andydolman at gmail.com
>
>
> 2009/9/9 <ledonret at email.unc.edu>
>
>> Dear List,
>>
>> I have been using R for a couple of years now, but I am very new to the
>> lme4 package. I am currently trying to use lmer to analyze a mixed model,
>> and also trying to bootstrap some downstream statistics that utilize the
>> variance components I am obtaining from this model. However, I'm not sure
>> that I am specifying the model correctly, and while I have written a
>> function that works with the boot package for iterations of up to about 30,
>> it will not work for greater iterations. I was hoping you or others could
>> shed some light on these issues!
>>
>> *I am using lme4 version 0.99375-28 with Mac OS X version 10.5.7
>>
>> My design is this: I created 9 Families, and I treated these families with
>> two diets. Each family by diet combination was replicated 9 times
>> (randomized and interspersed). In my model, I am treating Family as a random
>> effect, and Diet as a fixed effect. I am mostly concerned with determining
>> the variance of the Family by Diet interaction. From previous queries, I
>> gathered that the correct specification would be this (disregarding
>> replicate)
>>
>> lmer(SVL ~ Treatment + (1|Family) + (1|Family:Treatment), dataset)
>>
>> or equivalently
>>
>> lmer(SVL ~ Treatment + (1|Family/Treatment), dataset)
>>
>> where "SVL" is my trait. It's a little confusing to me, because in the
>> second formula it appears that Treatment is nested in Family, instead of
>> being completely crossed with family, but I tried both ways and they gave me
>> the same answer, so I was satisfied.
>
>
>
> Both these specifications are for nested models, for crossed you want
>
> lmer(SVL ~ Treatment + (1|Family) + (1|Treatment), dataset)
>
>
> But you've only got 2 levels of Treatment so it's probably not suitable to
> treat it as random, particularly as you want to use the variance components.
>
>
>
>
>> However, I'd also like to add replicate (Box1) as a nested factor. The way
>> I'd intuitively write that,
>>
>> lmer(SVL ~ Treatment + (1|Family/Treatment/Box1), dataset)
>>
>> should mean that replicate is completely cross with Family, instead of
>> being nested within. The way I tried to deal with this was to make sure each
>> replicate had a unique identity... is that sufficient?
>>
>>
>
> To nest Box in e.g. Treatment you'd do
>
> lmer(SVL ~ Treatment + (1|Family) + (1|Treatment/Box1), dataset)
>
> if that's what you want to do.

Actually it would be

lmer(SVL ~ Treatment + (1|Family) + (1|Treatment:Box1), dataset)

so that you don't try to model Treatment as both a fixed effect and a
random effect.  A term like (1|Treatment/Box1) expands to
(1|Treatment) + (1|Treatment:Box1)

In fact, if the levels of Box1 are defined so that they are distinct
both within and between treatments you can write the model as

lmer(SVL ~ Treatment + (1|Family) + (1|Box1), dataset)

The nesting of Box1 within Treatment can be discovered from the data
as long as the levels of Box1 are defined as described above.

>
>
>
>
>
>> My second issue concerned bootstrapping and some errors that I receive when
>> trying to use this with an R>30. The function that I am attempting to use is
>> a function that I wrote to calculate heritability (I am trying to bootstrap
>> my estimate of the heritability of the Diet by Treatment interaction).
>>
>> #The function
>>
>>> fun<-function(data,i){
>>>
>> + d<-data[i,]
>> + M<-lmer(SVL~Treatment+(1|Family/Box1/Family),d)
>>
>
> family is nested in family here - probably not at all what you want
>
>
>
>
>> + Mv<-VarCorr(M)
>> + varcomps<-c(unlist(lapply(Mv,diag))) #Variance of random factors
>> + resid<-(attr(Mv,"sc")^2) #Residual Variance
>> + VarTot<-sum(varcomps,resid) #Total variance
>> + VarFamTrt<-Mv$`Treatment:Family`[1,1] #Variance for Family by Treatment
>> + return((2*VarFamTrt)/VarTot)} #Heritability
>> #bootstrapping
>>
>>> Fam.boot<-boot(data,statistic=fun,R=10)
>>>
>> This will work up to about R=30, but beyond that, I get this error message,
>>
>> Warning message:
>> In mer_finalize(ans) : singular convergence (7)
>>
>> I realize this is probably a problem that isn't specific to lmer, but it is
>> the first time that I have run into it, so I'm at a loss for what to do.
>>
>> Thank you so much in advance, for any insight you can provide!
>>
>> Best, Cristina Ledon-Rettig
>> UNC-Chapel Hill
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>