Hi,
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@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.
