[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