repeated measures lme question

Dear Melissa,

I would go for a different model with a different dataset. Instead of lumping the pre-treatment data I would keep them as individual observations, all having the 'A' treatment.

library(nlme)
set.seed(12345)
#create a toy dataset
dataset <- expand.grid(Year = 1995:2009, ID = 1:20, Disc = 1:5)
dataset$CO2 <- factor(ifelse(dataset$Year >= 2001 & dataset$ID >= 11, "E", "A"))
dataset$rw <- rnorm(nrow(dataset), mean = c(2, 4)[dataset$CO2])
dataset$yearCat <- factor(dataset$Year)

#create a new variable to define the interaction.
dataset$yearCatCO2 <- ifelse(dataset$CO2 == "A", "A", paste(dataset$yearCat, "E", sep = ""))

#model with interaction between year and treatment
m1 <- lme(rw ~ yearCat + yearCatCO2, random = ~1|ID/Disc, data = dataset, correlation = corAR1(form = ~ Year), method = "ML")
#model with only main effects of year and treatment
m0 <- lme(rw ~ yearCat + CO2, random = ~1|ID/Disc, data = dataset, correlation = corAR1(form = ~ Year), method = "ML")
#a likelihoodratio test for the interaction between treatment and year
anova(m1, m0)


Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Melissa Dawes
Verzonden: maandag 12 augustus 2013 10:54
Aan: r-sig-mixed-models at r-project.org
Onderwerp: [R-sig-ME] repeated measures lme question

Dear All,

I am trying to fit a repeated-measures linear mixed effects model using lme and have some doubts about the ANOVA results. My experiment consists of 20 trees, which were assigned to ambient (A) or elevated (E) CO2 for
9 years. The trees were then harvested and ring width was measured on multiple stem discs per tree, including several pre-treatment years.
Therefore, my dataset has the following variables:

CO2 (A or E treatment)
ID (code for individual tree, 20 with E and 20 with A CO2) yearCat (year of ring width, 2001-2009; set as categorical variable) Disc (individual disc with ring width measurements, a-c except for one tree with 4 discs a-d) rw (ring width measurement for individual year within disc within tree)
cov95to00 (mean ringwidth over pre-treatment years 1995 to 2000 for an individual disc)

I would like to analyze the effect of CO2 enrichment and whether it varied during the 9 years of treatment, using the pre-treatment ring width measured on individual discs as a covariable. So far, I have set up the following model, including corAR1 to account for auto-correlation of residuals between years:

rw.model <- lme(rw ~ cov95to00 + CO2 * yearCat, random = ~ 1 | ID/Disc, data = ringsyr, na.action = na.omit, method = "REML", corAR1())

However, the ANOVA output shows that year and CO2:yearCat are tested at the individual disc level, which means a very large denom DF, and I am not convinced this is correct.

numDF denDFF-value p-value
(Intercept)14128.025640.0048
cov95to001403.331490.0754
CO21182.705690.1173
yearCat8412 103.19579<.0001
CO2:yearCat84123.375340.0009

I also get huge confidence intervals (using the function intervals()) for the random effect associated with Disc or even the message "cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance", which makes me think the model structure might be incorrect.

While I am more confident in a model where measurements are averaged at the tree level (mean of all discs per tree for each year), I think a model with individual discs is better, especially because then I can use pre-treatment measurements that are specific to individual discs rather than averaged at the tree level. I get rather different ANOVA results when I average at the tree level.

I would be grateful for any suggestions about how to correctly apply lme in this situation.

Thank you,
Melissa



_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.