Hello I am trying to fit a REML to some soil mineral data which has been collected over the time period 1999 - 2004. I want to know if the 19 different treatments imposed, differ in terms of their soil mineral content. A tree model of the data has shown differences between the treatments can be attributed to the Magnesium, Potassium and organic matter content of the soil, with Magnesium being the primary separating variable. I am looking at soil mineral data were collected : 99, 02, 04. In the experiment, there are 19 different treatments (treatmentcontrol, treatment6TFYM, treatment 12TFYM etc), which are replicated in 3 blocks. For the magnesium soil data, I have created the following groupedData object: magnesium<-groupedData(Mg~year|treatment, inner=~block) Where mg=magnesium Kg/ha As it is a repeated measures I was going to use an lme. I have looked at Pinherio and Bates : Mixed-Effects models in S and S-plus and I am getting slightly confused. In order to fit the lme, should I specify the data to use in the model as the grouped structure model? If so is the following command correct: Model1<-lme(mg~treatment, random=block|year, data=magnesium)? I am slightly worried that it isn't, because in model summary, instead of listing the 19 different treatments in the fixed effects section, it writes intercept (as normal), then treatment^1, treatment^2 etc. However if I don't specify the groupedData object in the model, then in the fixed effects section, it names the treatments (i.e. intercept, treatmentcontrol, treatment6TFYM. Should I be fitting the model using the whole data set rather than the groupedData object? Thank you very much for your help Emma Pilgrim ------------------------------------------------------------------------ ------ Dr Emma Pilgrim Plant Ecologist IGER North Wyke Okehampton EX20 2SB email: emma.pilgrim@bbsrc.ac.uk [[alternative HTML version deleted]]
Hello, I'm an R-newbie, but I've been learning to use lme for repeated measures experiments as well. If I understand correctly: Outcome variable: Mg (Kg/ha) Subject/grouping variable: block Condition/treatment: treatment (19 levels) Repeated factor: time (3 levels: 99, 02, 04) I think if you specify the model formula in the lme call, then the formula structure specified in the groupedData object is ignored. One suggestion for the model: Model1<-lme(mg~treatment + year + treatment:year, random=~1|block, data=magnesium) If the question of interest is the treatment:year interaction Or Model2 <- lme(mg~treatment, random=~1|block, data=magnesium) Hope this helps ... and hope the experts chime in at this point to give more guidance. Keith ------quoting original post--- Hello I am trying to fit a REML to some soil mineral data which has been collected over the time period 1999 - 2004. I want to know if the 19 different treatments imposed, differ in terms of their soil mineral content. A tree model of the data has shown differences between the treatments can be attributed to the Magnesium, Potassium and organic matter content of the soil, with Magnesium being the primary separating variable. I am looking at soil mineral data were collected : 99, 02, 04. In the experiment, there are 19 different treatments (treatmentcontrol, treatment6TFYM, treatment 12TFYM etc), which are replicated in 3 blocks. For the magnesium soil data, I have created the following groupedData object: magnesium<-groupedData(Mg~year|treatment, inner=~block) Where mg=magnesium Kg/ha As it is a repeated measures I was going to use an lme. I have looked at Pinherio and Bates : Mixed-Effects models in S and S-plus and I am getting slightly confused. In order to fit the lme, should I specify the data to use in the model as the grouped structure model? If so is the following command correct: Model1<-lme(mg~treatment, random=block|year, data=magnesium)? I am slightly worried that it isn't, because in model summary, instead of listing the 19 different treatments in the fixed effects section, it writes intercept (as normal), then treatment^1, treatment^2 etc. However if I don't specify the groupedData object in the model, then in the fixed effects section, it names the treatments (i.e. intercept, treatmentcontrol, treatment6TFYM. Should I be fitting the model using the whole data set rather than the groupedData object? Thank you very much for your help Emma Pilgrim
emma pilgrim (IGER-NW) wrote:> Hello > > I am trying to fit a REML to some soil mineral data which has been > collected over the time period 1999 - 2004. I want to know if the 19 > different treatments imposed, differ in terms of their soil mineral > content. A tree model of the data has shown differences between the > treatments can be attributed to the Magnesium, Potassium and organic > matter content of the soil, with Magnesium being the primary separating > variable. > > I am looking at soil mineral data were collected : 99, 02, 04. > > In the experiment, there are 19 different treatments (treatmentcontrol, > treatment6TFYM, treatment 12TFYM etc), which are replicated in 3 > blocks. > > For the magnesium soil data, I have created the following groupedData > object: > > magnesium<-groupedData(Mg~year|treatment, inner=~block) > Where mg=magnesium Kg/haAre you sure you want treatment to be the grouping factor?> As it is a repeated measures I was going to use an lme. I have looked > at Pinherio and Bates : Mixed-Effects models in S and S-plus and I am > getting slightly confused. In order to fit the lme, should I specify > the data to use in the model as the grouped structure model? > > If so is the following command correct: > > Model1<-lme(mg~treatment, random=block|year, data=magnesium)? > > I am slightly worried that it isn't, because in model summary, instead > of listing the 19 different treatments in the fixed effects section, it > writes intercept (as normal), then treatment^1, treatment^2 etc.This is an unfortunate side-effect of creating a groupedData object - to create plots with panels in a natural order the grouping factor is changed to an ordered factor. In your case the treatment factor will become an ordered factor and the default contrasts for an ordered factor are the polynomial contrasts. There are two ways to get around this - don't create a groupedData object or change the default contrasts using options(contrasts = c(unordered = "contr.treatment", ordered = "contr.treatment")> However if I don't specify the groupedData object in the model, then in > the fixed effects section, it names the treatments (i.e. intercept, > treatmentcontrol, treatment6TFYM.Yes.> Should I be fitting the model using the whole data set rather than the > groupedData object?Probably that is the best course.