[R-sig-ME] variance structure in lme/lmer
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sun Apr 11 20:43:24 CEST 2010
Dear Alessia,
I think the first model is trying to fit the model you intend, the
problem is that it is failing to converge.
It is hard for me too see from the information I have, but I think the
problem is that the Mom effects are confounded with the residuals. Am
I right in thinking there is one observation per Mom in each
treatment? If so, you cannot separate the treatment specific Mom
variances from the treatment specific residual variances.
You have four options:
a) Drop Mom from the model and interpret the weight term as the
combined effects of Mom and residual. Here you assume the variances
differ between the treatments, and the Mom effects in the two
treatments are uncorrelated
b) Have "~at.level(Range, "invasive")|Mom" (i.e without the
interaction with Treatment) and the weight term. Here you assume the
variances differ between the treatments (again measured by the weights
term), but the Mom effects have a positive covariance which is
estimated by the Mom variance component.
c) Use ASReml to fit a model where the Mom terms in the original model
defines a residual covariance matrix. This is the same as b) except
the covariances are free to go negative
You can imagine these models as trying to fit a 2x2 covariance matrix
for each range. The covariance matrix is the Mom+Residual
(co)variances for the two treatments. If we label the estimates as a)
weight in treatment 1, b) weight in treatment 2 and c) the Mom
estimate then these matrices have the form
a) = a 0
0 b
b) = a+c c
c b+c
where as ASReml would directly estimate the covariance matrix.
MCMCglmm cannot estimate c) in this case because it would not allow
the Range interaction in the residual structure as well as the
Treatment interaction. This is an over sight on my part but do not
have time in the next few months to remedy it.
Cheers,
Jarrod
Quoting Alessia Guggisberg <alessiag at interchange.ubc.ca>:
> Dear Jarrod
>
> I've tried many times to run your function in lme, but I always got
> an error message, must often of the form:
>
>> model=lme(above~Treatment*Range*t0.LLL, random =
>> list(~at.level(Range, "invasive"):Treatment-1|Pop, ~at.level(Range,
>> "native"):Treatment-1|Pop, ~at.level(Range,
>> "invasive"):Treatment-1|Mom, ~at.level(Range,
>> "native"):Treatment-1|Mom), weights=varIdent(form=~1|Treatment*Range)
>
> Erreur dans lme.formula(above ~ Treatment * Range * t0.LLL, random =
> list(~at.level(Range, :
> nlminb problem, convergence error code = 1
> message = iteration limit reached without convergence (9)
>
> or
>
>> model=lme(above~Treatment*Range*t0.LLL, random =
>> list(~at.level(Range, "invasive"):Treatment-1|Pop/Mom,
>> ~at.level(Range, "native"):Treatment-1|Pop/Mom),
>> weights=varIdent(form=~1|Treatment*Range),
>> data=ch.above.3);summary(model)
>
> Erreur dans MEEM(object, conLin, control$niterEM) :
> NA dans un appel à une fonction externe (argument 2)
> De plus : Message d'avis :
> In ncols * c(rep(1, Q), 0, 0) :
>
> I'm afraid your function cannot be combined with weights in lme.
>
> 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
> To: Alessia Guggisberg
> Cc: r-sig-mixed-models at r-project.org
> Sent: Friday, April 09, 2010 2:04 AM
> Subject: Re: [R-sig-ME] variance structure in lme/lmer
>
>
> Hi,
>
>
> I think
>
>
> model=lme(above~Treatment, random = list(~ at.level(Range,
> "native"):Treatment-1|Pop/mom, ~ at.level(Range,
> "invasive"):Treatment-1|Pop/mom), weights=varIdent(form=~1*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 br>
> 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"):Trea 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, 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" <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
>
>
>
>
>
> &nb native * 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, lev /blockquote>
> }
>
> 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 introduct 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:
>
>
>
> - ticity 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
>
>
>
> (Intercept) Residual
>
> StdDev: 1.150242 3.684324
>
>
>
> Variance function:
>
> Structure: Different standard deviations per stratum
>
> Formula: ~1 | Treatment * Range
>
> Parameter estimates:
>
> Control*native &nbs 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)
>
> ## 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 Va te>
> 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
>
> Treatm e 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 &n sp; 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 e 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 + invasiveCon veHerbivory|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
>
> Pop invasiveHerbivory 0.0000e+00
>
> Pop nativeHerbivory 2.4815e-10
>
> Pop invasiveControl 1.7994e-09
>
> Pop nativeControl 9.1133e-01
>
> Residual &nbs ; 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 yo you might be able to provide and look forward
> to receiving your recommendations.
>
>
>
> Best,
>
>
>
> Alessia
>
>
>
> *************************************************************
>
> Alessia Guggisberg, PhD
>
> Postdoctoral Fellow
>
> Department of Botany
>
> Un bia
>
> 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.
>
>
>
>
>
>
>
>
> ------------------------------------------------------------------------------
>
>
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
--
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