[R-sig-ME] variance structure in lme/lmer
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Apr 8 11:22:39 CEST 2010
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 11.213378
> -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