[R-sig-ME] variance structure in lme/lmer
Jarrod Hadfield
j.hadfield at ed.ac.uk
Fri Apr 9 11:04:34 CEST 2010
Hi,
I think
model=lme(above~Treatment * Range, random = list(~ at.level(Range,
"native"):Treatment-1|Pop/mom, ~ at.level(Range,
"invasive"):Treatment-1|Pop/mom), weights=varIdent(form=~1|
Treatment*Range))
should do the job.
Cheers,
Jarrod
On 8 Apr 2010, at 18:08, Alessia Guggisberg wrote:
> Hi Jarrod
>
> Thank you so much for your reply!
>
> Indeed, each population/mom belongs to one type of range (native or
> invasive) solely, but occurs in both treatment plots. I've just
> tried your suggestion and its seems to work fine:
>
> model=lmer(above~Treatment * Range * Ct0.LLL + (at.level(Range,
> "native"):Treatment-1|Pop/Mom) +
> (at.level(Range,"invasive"):Treatment-1|Pop/Mom))
>
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Mom:Pop at.level(Range, "invasive"):TreatmentControl 1.1326e
> +01 3.3654e+00
> at.level(Range, "invasive"):TreatmentHerbivory
> 9.1458e-01 9.5633e-01 0.292
> Mom:Pop at.level(Range, "native"):TreatmentControl 1.1868e
> +01 3.4451e+00
> at.level(Range, "native"):TreatmentHerbivory 1.7996e
> +00 1.3415e+00 0.396
> Pop at.level(Range, "invasive"):TreatmentControl
> 1.1083e-09 3.3291e-05
> at.level(Range, "invasive"):TreatmentHerbivory
> 2.9827e-01 5.4614e-01 0.000
> Pop at.level(Range, "native"):TreatmentControl 5.5788e
> +00 2.3619e+00
> at.level(Range, "native"):TreatmentHerbivory 1.5701e
> +00 1.2530e+00 0.839
> Residual 1.1545e+00 1.0745e+00
> Number of obs: 368, groups: Mom:Pop, 184; Pop, 42
>
> Also, the BLUPs are now 'logical' as I get zero-values for
> populations/moms that do not occur in a given range.
>
> The question remains though, whether I should additionnally account
> for heterogeneity in the innermost residuals (like 'weights' in
> lme), and if so, how I should implement that in lmer.
>
> Thanks again for your precious help; thanks to you, I've just made a
> big step today....
>
> Best,
>
> Alessia
>
>
> *************************************************************
> Alessia Guggisberg, PhD
> Postdoctoral Fellow
> Department of Botany
> University of British Columbia
> 3529-6270 University Boulevard
> Vancouver, BC, V6T 1Z4
> Canada
>
> Email: alessiag at interchange.ubc.ca
>
> *************************************************************
> ----- Original Message ----- From: "Jarrod Hadfield" <j.hadfield at ed.ac.uk
> >
> To: "Alessia Guggisberg" <alessiag at interchange.ubc.ca>
> Cc: <r-sig-mixed-models at r-project.org>
> Sent: Thursday, April 08, 2010 2:22 AM
> Subject: Re: [R-sig-ME] variance structure in lme/lmer
>
>
>> Dear Alessia,
>>
>> I think the problem is that individuals from a population all
>> belong to one type (native or invasive) (is this true?) If this is
>> the case your 4x4 covariance matrices needs to look like
>>
>> TreatmentControl:Rangenative * * 0 0
>> TreatmentHerbivory:Rangenative * * 0 0
>> TreatmentHerbivory:Rangeinvasive 0 0 * *
>> TreatmentHerbivory:Rangeinvasive 0 0 * *
>>
>> where * indicates it could be estimated and 0 not.
>>
>> To fit such a covariance matrix try
>>
>> (at.level(Range, "native"):Treatment-1|Pop)+(at.level(Range,
>> "invasive"):Treatment-1|Pop)
>>
>> in lmer. The at.level function is at the bottom of the email
>>
>> If individuals with the same mother have the same treatment then
>> the pop/mom matrix needs to be
>>
>>
>> TreatmentControl:Rangenative * 0 0 0
>> TreatmentHerbivory:Rangenative 0 * 0 0
>> TreatmentHerbivory:Rangeinvasive 0 0 * 0
>> TreatmentHerbivory:Rangeinvasive 0 0 0 *
>>
>> and you can use either
>>
>> (at.level(Range, "native"):at.level(Treatment, "Herbivory")-1|
>> Pop:mom)+(at.level(Range, "invasive"):at.level(Treatment,
>> "Herbivory")-1|Pop:mom)+(at.level(Range,
>> "native"):at.level(Treatment, "Control")-1|Pop:mom)+
>> (at.level(Range, "invasive"):at.level(Treatment, Control")-1|Pop:mom)
>>
>> or
>>
>> (at.level(subgroup, 1)|Pop:mom)+ (at.level(subgroup, 2)|Pop:mom) +
>> (at.level(subgroup, 3)|Pop:mom) + (at.level(subgroup, 4)|Pop:mom)
>>
>> I'm not sure if you can use this type of syntax in lme, but without
>> allowing different residual variances for the subgroups I would be
>> very cautious about the results.
>>
>> Cheers,
>>
>> Jarrod
>>
>>
>>
>> at.level<-function (x, level)
>> {
>> if (is.numeric(level)) {
>> M <- outer(x, levels(x)[level], "==")
>> }
>> else {
>> M <- outer(x, level, "==")
>> }
>> mode(M) <- "numeric"
>> M
>> }
>> <environment: namespace:MCMCglmm>
>>
>>
>>
>> On 8 Apr 2010, at 06:50, Alessia Guggisberg wrote:
>>
>>> Dear R-sig-mixed-models list members
>>>
>>> I apologise in advance, if my question seems trivial; I'm new to
>>> mixed-models in R.
>>>
>>> To make things easier, I'm giving you an example, followed by the
>>> questions I'm hoping to be able to address with one of you.
>>>
>>> Study: Response of different populations of plants from two
>>> different ranges (Europe vs. North America) to various stresses in
>>> a greenhouse experiment
>>> Study question: Do populations from the native range (Europe)
>>> differ in terms of response from the ones of the introduced range
>>> (North America)?
>>> Type of data: Growth and reproductive data
>>> Example of mixed model: above-ground biomass as a function of
>>> treatment (control vs. stress), range (Europe vs. North America)
>>> and a covariate (in this case, longest leaf length at t0, to
>>> account for size difference prior to stress)
>>>
>>> model=lme(above~Treatment * Range * t0.LLL, random= ~1|Pop/Mom)
>>>
>>> model=lmer(above~Treatment * Range * t0.LLL + (1|Pop/Mom))
>>>
>>> Due to the introduction history (bottleneck), I expect less
>>> genetic variation in the introduced than in the native range,
>>> hence less variation in response in the introduced than native
>>> range. Indeed, inspection of raw data and residuals indicate
>>> higher variance for native than introduced range. Similarly, I
>>> observe less variation in response in the stress than in the
>>> control plot. Obviously, these observations suggest violation of
>>> homoscedasticity. I therefore would like to test and eventually
>>> account for difference in variances between the different
>>> subgroups: (i) Europe-control, (ii) Europe-stress, (iii) North
>>> America-control, (iv) North America-stress.
>>>
>>> Following options are possible:
>>>
>>> - account for heteroscedasticity in the innermost residuals:
>>>
>>> model=lme(above~Treatment * Range * t0.LLL, random= ~1|Pop/Mom,
>>> weights=varIdent(form=~1|Treatment*Range))
>>>
>>> Random effects:
>>> Formula: ~1 | Pop
>>> (Intercept)
>>> StdDev: 0.997503
>>>
>>> Formula: ~1 | Mom %in% Pop
>>> (Intercept) Residual
>>> StdDev: 1.150242 3.684324
>>>
>>> Variance function:
>>> Structure: Different standard deviations per stratum
>>> Formula: ~1 | Treatment * Range
>>> Parameter estimates:
>>> Control*native Herbivory*native Control*invasive
>>> Herbivory*invasive
>>> 1.0000000 0.3343747 0.9324055 0.2389115
>>>
>>> Is it true, that weights cannot be incorporated in lmer?
>>>
>>> - stratify the variance for the random effects:
>>>
>>> subgroup=factor(Treatment:Range)
>>> model=lme(above~Treatment * Range * t0.LLL, random= ~subgroup-1|
>>> Pop/ Mom)
>>> ## warning message: Fewer observations than random effects in all
>>> level 2 groups
>>>
>>> model=lmer(above~Treatment * Range * t0.LLL + (0 + subgroup|Pop/
>>> Mom))
>>>
>>> Random effects:
>>> Groups Name Variance
>>> Corr
>>> Mom:Pop TreatmentControl:Rangenative 11.783704 0.000
>>> TreatmentHerbivory:Rangenative 1.714778
>>> 0.419 0.000
>>> TreatmentControl:Rangeinvasive 378
>>> -0.003 0.311 -0.001
>>> TreatmentHerbivory:Rangeinvasive 0.819599
>>> Pop TreatmentControl:Rangenative 5.578492
>>> TreatmentHerbivory:Rangenative 1.569965 1.000
>>> TreatmentControl:Rangeinvasive 0.029164
>>> 0.839 0.839
>>> TreatmentHerbivory:Rangeinvasive 0.315794
>>> 1.000 1.000 0.839
>>> Residual
>>> 1.239337
>>> Number of obs: 368, groups: Mom:Pop, 184; Pop, 42
>>>
>>> Why do I get an error message, when I run the model in lme?
>>>
>>> - another option found in R forums would be to code each subgroup
>>> with a dummy variable, as follows (I only know how to implement
>>> that method with lmer):
>>>
>>> nativeControl=as.numeric(factor(Range=="native" &
>>> Treatment=="Control"))
>>> invasiveControl=as.numeric(factor(Range=="invasive" &
>>> Treatment=="Control"))
>>> nativeHerbivory=as.numeric(factor(Range=="native" &
>>> Treatment=="Herbivory"))
>>> invasiveHerbivory=as.numeric(factor(Range=="invasive" &
>>> Treatment=="Herbivory"))
>>> model=lmer(above~Treatment*Range*t0.LLL + (0 + nativeControl|Pop/
>>> Mom) + (0 + invasiveControl|Pop/Mom) + (0 + nativeHerbivory|Pop/
>>> Mom) + (0 + invasiveHerbivory|Pop/Mom))
>>>
>>> Random effects:
>>> Groups Name Variance
>>> Mom:Pop invasiveHerbivory 0.0000e+00
>>> Mom:Pop nativeHerbivory 0.0000e+00
>>> Mom:Pop invasiveControl 9.7939e-01
>>> Mom:Pop nativeControl 1.0516e+00
>>> Pop invasiveHerbivory 0.0000e+00
>>> Pop nativeHerbivory 2.4815e-10
>>> Pop invasiveControl 1.7994e-09
>>> Pop nativeControl 9.1133e-01
>>> Residual 4.1959e+00
>>> Number of obs: 368, groups: Mom:Pop, 184; Pop, 42
>>>
>>> This last model entails less parameters than the previous one
>>> (because it doesn't calculate any correlations), but its AIC is
>>> also much higher. Other than that, what differs between the two
>>> previous models?
>>>
>>> In view of my problematic, what do you suggest me to do?
>>>
>>> I thank you in advance for any help you might be able to provide
>>> and look forward to receiving your recommendations.
>>>
>>> Best,
>>>
>>> Alessia
>>>
>>> *************************************************************
>>> Alessia Guggisberg, PhD
>>> Postdoctoral Fellow
>>> Department of Botany
>>> University of British Columbia
>>> 3529-6270 University Boulevard
>>> Vancouver, BC, V6T 1Z4
>>> Canada
>>>
>>> Email: alessiag at interchange.ubc.ca
>>>
>>> *************************************************************
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>
>>
>> --
>> The University of Edinburgh is a charitable body, registered in
>> Scotland, with registration number SC005336.
>
>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100409/b026dbbf/attachment.pl>
More information about the R-sig-mixed-models
mailing list