R help - Oct 2009 - Reposting various problems with two-way anova, lme, etc.

Hi,

I posted the message below last week, but no answers, so I'm giving it
another attempt in case somebody who would be able to help might have missed
it and it has now dropped off the end of the list of mails.

I am fairly new to R and still trying to figure out how it all works, and I
have run into a few issues. I apologize in advance if my questions are a bit
basic, I'm also no statistics wizard, so part of my problem my be a more
fundamental lack of knowledge in the field.

I have a dataset that looks something like this:

Week Subj   Group Readout
0    1      A     352.2
1    1      A     366.6
2    1      A     377.3
3    1      A     389.8
4    1      A     392.5
0    2      B     344.7
1    2      B     360.9
.
.
.
.

So basically I have a number of subjects that are divided into different
groups receiving different interventions/treatments. Observations on these
subjects are made on 5 occasions (Week 0-4). I would like to see if there is
difference between the "treatment" groups and if the observations that
we
are making change in the individual groups over time. According to my very
limited statistics training that means that I would have to do a two-way
anova with repeated measures because the same subjects are observed on the
different weeks.

Now in R I can do something like this:

MyData$fWeek <- factor(MyData$Week) #Convert week to factor
m1 <- aov(Readout ~ Group*fWeek + Error(Subj/fWeek), data=MyData)

My first naive question is: Is the Error(...) term correct for the design
that I describe? I am not quite sure if I understand the syntax correctly.

In any case, it actually seems to work fine, but I can't figure out how to
do post hoc testing on this. TukeyHSD does not work for the aovlist which is
returned, so I am kind of stuck. Is there a simple way to do the post hoc
test based on the aovlist?

I have been reading various questions on the list and I think that I have
understood that I should be using lme from the nlme package and this is
where I run into some problems. As I understand it, the lme command that I
would need would look something like this:

m1 <- lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek, data=MyData)

Now this actually fails with an error message like:

Error in MEEM(object, conLin, control$niterEM) :
  Singularity in backsolve at level 0, block 1

The reason (I believe) is that I have some weeks where there are no
observations for a specific group. In this case, one of the groups had a lot
of drop-out and at Week 4, there are no subjects left in one of the groups
and that seems to be causing the problem. I can get it to run by excluding
the last week with something like:

m1 <- lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek,
data=MyData[which(MyData$Week < 4), ])

My next question is: Is there another way around that? I would still like to
run the analysis for the remaining groups on the last time point and I
believe that it should somehow be included into the entire analysis. I have
also managed to find a few postings with similar problems, but no real
solutions, so anything helpful comments would be much appreciated.

My final issue is how do I do the post hoc testing on the model. I
understand (I think) that I should use the glht function from the multcomp
package. For Instance, I could do something like:

summary(glht(m1,linfct=c("GroupB:fWeek1 - GroupC:fWeek1 =
0","GroupB:fWeek2
- GroupC:fWeek2 = 0")))

And that actually works fine. My problem is that Group A and fWeek 0 seems
to have turned into the (intercept), so I can't write something like:

summary(glht(m1,linfct=c("GroupB:fWeek0 - GroupB:fWeek1 = 0")))

To check for changes between baseline and week 1 in Group B because I get
the error message:

Error in chrlinfct2matrix(linfct, names(beta)) :
  variable(s) ‘GroupB:fWeek0’ not found

When I look at summary(m1) the error message makes sense, because (as I
said) I think they are listed as the (intercept). However, that is a bit
inconvenient and I was wondering if there is a different way to to so that I
don't have that problem. How do you actually control what the intercept is?
In my specific case, I don't think of Group A as baseline or vehicle or
placebo treatment.

I apologize for my long email. To sum up my questions:
1. Is my formulation of the Error(...) term or the value for the random
variable in lme correct given the design that I have described?
2. Is there a way to do post hoc testing on the aovlist object that comes
out of the anova with repested measures
3. Is there a way to make lme not fail when there are group/time
combinations with no entries
4. How do you control what is the (intercept) in the model returned by the
lme function and is there a way to still be able to refer to all groups and
timepoints in there without referring to intercept?

I would appreciate any help and comments and do feel free to educate me if I
have got it all wrong....it is quite likely.

Thanks,
Michael

	[[alternative HTML version deleted]]

Hi Michael,
>> How do you control what is the (intercept) in the model returned by the
>> lme function and is there a way to still be able to refer to all groups
>> and 
>> timepoints in there without referring to intercept?
Here is some general help. The intercept is controlled by the contrasts that
you used when fitting your model. By default these are treatment (i.e.
Dunnett) contrasts, but you can change this.

##
?options
options("contrasts")
?contrasts
options(contrasts=c("contr.sum", "contr.poly"))
options("contrasts")

You can design your own contrasts or you can use those provided by
"base" R:

##
?contr.treatment

(There is also contr.sdif in MASS, which is a very useful type for
environmental studies) See
http://users.monash.edu.au/~murray/stats/BIO4200/LinearModels.pdf, which
covers the main ready-made types in R. For ANOVA, sum-to-zero contrasts are
commonly used.

If you do use treatment contrasts you can (usually) change your reference
level "on-the-fly" or you can set it when you make your factor.

##
?factor
?levels
?reorder
?relevel

If you are using the multcomp package then look at:

?contrMat

An example of how to use it is given under the glht function.

Regards, Mark.


Michael Schacht Hansen wrote:
> 
> Hi,
> 
> I posted the message below last week, but no answers, so I'm giving it
> another attempt in case somebody who would be able to help might have
> missed
> it and it has now dropped off the end of the list of mails.
> 
> I am fairly new to R and still trying to figure out how it all works, and
> I
> have run into a few issues. I apologize in advance if my questions are a
> bit
> basic, I'm also no statistics wizard, so part of my problem my be a
more
> fundamental lack of knowledge in the field.
> 
> I have a dataset that looks something like this:
> 
> Week Subj   Group Readout
> 0    1      A     352.2
> 1    1      A     366.6
> 2    1      A     377.3
> 3    1      A     389.8
> 4    1      A     392.5
> 0    2      B     344.7
> 1    2      B     360.9
> .
> .
> .
> .
> 
> So basically I have a number of subjects that are divided into different
> groups receiving different interventions/treatments. Observations on these
> subjects are made on 5 occasions (Week 0-4). I would like to see if there
> is
> difference between the "treatment" groups and if the observations
that we
> are making change in the individual groups over time. According to my very
> limited statistics training that means that I would have to do a two-way
> anova with repeated measures because the same subjects are observed on the
> different weeks.
> 
> Now in R I can do something like this:
> 
> MyData$fWeek <- factor(MyData$Week) #Convert week to factor
> m1 <- aov(Readout ~ Group*fWeek + Error(Subj/fWeek), data=MyData)
> 
> My first naive question is: Is the Error(...) term correct for the design
> that I describe? I am not quite sure if I understand the syntax correctly.
> 
> In any case, it actually seems to work fine, but I can't figure out how
to
> do post hoc testing on this. TukeyHSD does not work for the aovlist which
> is
> returned, so I am kind of stuck. Is there a simple way to do the post hoc
> test based on the aovlist?
> 
> I have been reading various questions on the list and I think that I have
> understood that I should be using lme from the nlme package and this is
> where I run into some problems. As I understand it, the lme command that I
> would need would look something like this:
> 
> m1 <- lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek, data=MyData)
> 
> Now this actually fails with an error message like:
> 
> Error in MEEM(object, conLin, control$niterEM) :
>   Singularity in backsolve at level 0, block 1
> 
> The reason (I believe) is that I have some weeks where there are no
> observations for a specific group. In this case, one of the groups had a
> lot
> of drop-out and at Week 4, there are no subjects left in one of the groups
> and that seems to be causing the problem. I can get it to run by excluding
> the last week with something like:
> 
> m1 <- lme(Readout ~ Group*fWeek,random=~1|Subj/fWeek,
> data=MyData[which(MyData$Week < 4), ])
> 
> My next question is: Is there another way around that? I would still like
> to
> run the analysis for the remaining groups on the last time point and I
> believe that it should somehow be included into the entire analysis. I
> have
> also managed to find a few postings with similar problems, but no real
> solutions, so anything helpful comments would be much appreciated.
> 
> My final issue is how do I do the post hoc testing on the model. I
> understand (I think) that I should use the glht function from the multcomp
> package. For Instance, I could do something like:
> 
> summary(glht(m1,linfct=c("GroupB:fWeek1 - GroupC:fWeek1 >
0","GroupB:fWeek2
> - GroupC:fWeek2 = 0")))
> 
> And that actually works fine. My problem is that Group A and fWeek 0 seems
> to have turned into the (intercept), so I can't write something like:
> 
> summary(glht(m1,linfct=c("GroupB:fWeek0 - GroupB:fWeek1 = 0")))
> 
> To check for changes between baseline and week 1 in Group B because I get
> the error message:
> 
> Error in chrlinfct2matrix(linfct, names(beta)) :
>   variable(s) ?GroupB:fWeek0? not found
> 
> When I look at summary(m1) the error message makes sense, because (as I
> said) I think they are listed as the (intercept). However, that is a bit
> inconvenient and I was wondering if there is a different way to to so that
> I
> don't have that problem. How do you actually control what the intercept
> is?
> In my specific case, I don't think of Group A as baseline or vehicle or
> placebo treatment.
> 
> I apologize for my long email. To sum up my questions:
> 1. Is my formulation of the Error(...) term or the value for the random
> variable in lme correct given the design that I have described?
> 2. Is there a way to do post hoc testing on the aovlist object that comes
> out of the anova with repested measures
> 3. Is there a way to make lme not fail when there are group/time
> combinations with no entries
> 4. How do you control what is the (intercept) in the model returned by the
> lme function and is there a way to still be able to refer to all groups
> and
> timepoints in there without referring to intercept?
> 
> I would appreciate any help and comments and do feel free to educate me if
> I
> have got it all wrong....it is quite likely.
> 
> Thanks,
> Michael
> 
> 	[[alternative HTML version deleted]]
> 
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 
-- 
View this message in context:
http://www.nabble.com/Reposting-various-problems-with-two-way-anova%2C-lme%2C-etc.-tp25963667p25972403.html
Sent from the R help mailing list archive at Nabble.com.