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@interchange.ubc.ca

*************************************************************
  ----- Original Message ----- 
  From: Jarrod Hadfield 
  To: Alessia Guggisberg 
  Cc: r-sig-mixed-models@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@interchange.ubc.ca

    *************************************************************
    ----- Original Message ----- From: "Jarrod Hadfield" <j.hadfield@ed.ac.uk>
    To: "Alessia Guggisberg" <r-sig-mixed-models@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@interchange.ubc.ca



          *************************************************************

          [[alternative HTML version deleted]]



        < _____________________________________________

        R-sig-mixed-models@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.

	[[alternative HTML version deleted]]


