[R-sig-ME] variance structure in lme/lmer
Alessia Guggisberg
alessiag at interchange.ubc.ca
Thu Apr 8 19:08:02 CEST 2010
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.
>
More information about the R-sig-mixed-models
mailing list