You probably can do this with lme function, but I don't know that for
sure. "aov" (included in the "stats" package), with a call
to the
"Error" function how I generally analyze data obtained from a
repeated measure design. For a very good description of how Error
and aov work together, you should read Baron's “Notes on the use of R
for psychology experiments and questionnaires” (follow the link
below) The section beginning on page 36 addresses your question.
Generally speaking, the aov function expects you to describe your
model formula similarly to how you would with a call to the "lm"
function. For example, Y~X1*X2 describes a two way ANOVA--your
dependent variable Y modeled as explanatory variables X1 and X2,
including the interaction term:
model.aov<-aov(Y~X1*X2, data=DataFile)
If X1 and X2 are within-subject variables, then the above formula
would be written as follows:
model.aov<-aov(Y~X1*X2+Error(Subject/X1*X2),data=DataFile)
The call "summary(model.aov)" will display a summary of the model.
Depending on your experience working with model formulae in R, there
are several steps in the repeated-measure ANOVA procedure that can
come back to haunt you, if you're not careful (e.g., checking the
assumption of sphericity, normally distributed errors, homogeneity of
variance, etc., etc.), so make sure you're covered on these before
you believe your results.
Here's a link to the document I referred to: http://cran.r-
project.org/doc/contrib/Baron-rpsych.pdf
Kyle H. Ambert
Graduate Student, Dept. Behavioral Neuroscience
Oregon Health & Science University
ambertk@ohsu.edu
On Apr 8, 2007, at 7:55 PM, Scott Norton wrote:
> Hi,
> I have what I believe is a repeated-measures dataset that I'm
> trying to analyze using lme(). This is *not* homework, but an
> exercise in my trying to self-teach myself repeated-measure ANOVA
> for other *real* datasets that I have and that are extremely
> similar to the following design.
>
> I'm fairly sure the dataset described below would work with lme()
> -- but it'd be great if anybody can confirm that after I describe
> the dataset below)
>
> The study involves measuring the effect of a drug on blood
> pressure. There were 16 patients in all and 6 replicate measures
> per patient of their blood pressure on one week (one measure per
> day). Two weeks later, a drug was introduced to 8 randomly selected
> patients in such a way that I had equal representation of the 4 age
> groups among the two treatment groups. Then, another two weeks
> later, 6 replicate measures per patient (per day) of blood pressure
> was retaken. So each patient had 12 total measures whether they
> were in the treatment group or in the control group (6 reads (R1-
> R6) in the baseline-week and another 6 reads (R1-R6) in the post-
> treatment week).
>
> So,
> Background: 16 patients
> Response measure: Blood pressure
> Fixed Factor: 4 Age groups
> Fixed Factor: Drug vs. NoDrug
> Random factor: Day of the read (i.e. 6 replicate reads (R1-R6) at
> the baseline time, and 6 replicate reads (R1-R6) after the drug has
> had time to take effect)
> Random Factor: Subjects 1-16
>
> Patient AgeGroup BP(Blood Pressure) Read (replicate
> reads) Pre/PostTreatmentWeek Group
> 1 20-29 83
> R1
> pre Treat
> 2 20-29 81
> R1
> pre Control
> 3 20-29 74
> R1
> pre Treat
> 4 20-29 85
> R1
> pre Control
> 5 30-39 82
> R1
> pre Treat
> … … … … … …
> 3 20-29 74
> R2
> pre Treat
> … … … … … …
> 1 20-29 83
> R1
> post Treat
> 2 20-29 82
> R1
> post Control
> 3 20-29 86
> R1
> post Treat
> 4 20-29 84
> R1
> post Control
> … … … … … …
>
> I'm trying to do an analysis of variance to decide whether there is
> a measurable change in blood pressure between the Treat and Control
> groups.
>
> Another issue is that some of the 16 patients didn't get all 6
> replicate reads in their pre/post treatment weeks, so I need to
> include the na.omit function.
>
> What I think I'm having the most trouble with is the repeated reads
> (R1 through R6) in the pre/post treatment weeks. I'm fairly sure
> this is a random variable -- their order or identify (R1 in pre-
> treatment week has no relation to R1 in the post-treatment week,
> etc). By placing Read as a random variable, am I covering myself
> there?
> If I execute:
>
>> summary(lme(BP ~ Group, random = ~ 1 | Patient, data = bloodpress,
>> na.action=na.omit))
>
> I get a result, but I'm not sure it's correct -- do I need to tell
> the model about the Read factor (R1-R6 in pre/post weeks)?
>
> I'm really trying to set the right form of the lme() function call
> to decide
> 1) if there is a statistical difference between the Treat/Control
> groups and,
> 2) if one takes into account AgeGroup, is there a statistical
> difference between Treat/Control Groups, and finally
> 3) if I don't see a statistical difference, can someone recommend
> an R function that might solve the supplemental question, "given
> the noise in day-to-day blood pressure reads, and given that I
> wanted to have enough statistical power to observe a say, a 5%
> benefit in blood pressure, how many additional reads or patients I
> would need."
>
> Basically, is lme() the proper function, and can someone offer any
> pointers on what my call to this function should look like to make
> the above to determinations?
>
> Thanks!
> -Scott
>
> ______________________________________________
> R-help@stat.math.ethz.ch 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.
>
[[alternative HTML version deleted]]