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

*************************************************************
	[[alternative HTML version deleted]]


