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

[ledonret at email.unc.edu]

Quoting Douglas Bates <bates at stat.wisc.edu>:

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

Thank you for your reply, Dr. Bates.

I think I understand now - my replicates (Box1) are distinct both 
within and between treatments, so in order to have a model that 
includes Treatment as a fixed effect, family as a random effect, and 
the family by treatment interaction as a random effect, I would use,

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

which expands to

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

Cris

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