GreenBlueBar.gif GreenBlueBar.gif

R Code for Chapter 18

Repeated-Measures Analysis of Variance


For this chapter I am going to give two examples. As I said in the text, repeated measures designs for R can be very messy. I will follow those with code that you might want to have in the future if you need to run repeated measures designs, but it is code that I don't think that you want to study too hard for this chapter. Think of it as an added resource.

SOMETHING WORTH POINTING OUT!! When I first wrote this code I would read in some data, use attach(), run the analysys, read more data, attach() again, etc. That was a disaster. A variable (e.g. dv) that any idiot knew had 80 values, only had 48 values. And this happened repeatedly. So I gave up on attach() and tried with(). That didn't work either. So I went to appending the name of the data frame to each variable name. It is a nuisance, but at least it works.

The first set of code comes from Exercise 18.10 and relates to a study by St. Lawrence et al. (1995) on behavioral skills training. There was only one group, so Group will not be a variable in our analysis. Each subject was measured four times, represented in four columns. For this analysis I have created the necessary variables directly instead of using what is called the reshape() function. It is a perfectly reasonable approach, and may in fact be preferable to the reshape approach. What I have done is to read the data in 4 columns corresponding to the different testing times. I then create a long variable called "dv" by concatenating the data from the four trials. But then I needed to add a variable named "time" denoting to which trial each score belongs. If you tell R to print out dv and time, you will see what I have done.These steps are shown below.


### Data from St. Lawrence et al. (1995)  Exercise 18.10
### NOTE!! Using attach() creates very many problems, so I am not attaching and, instead,
### use names like dataLong$Group for much of what follows. Messy, but ....

data.BST <- read.table("http://www.uvm.edu/~dhowell/fundamentals9/DataFiles/Ex18-10.dat", header = TRUE)
dv <- c(data.BST$Pretest, data.BST$Posttest, data.BST$FU6, data.BST$FU12)
time <- rep(1:4, each = 10)
subject <- rep(1:10, 4)
time <- factor(time)
subject <- factor(subject)
cat("\nTrial Means \n")
tapply(dv, time, mean)
cat("\nSubject Means \n")
tapply(dv, subject, mean)
BSTmodel <- aov(dv ~ time + Error(subject/time))
print(summary(BSTmodel))
rm(dv, time, subject, BSTmodel, data.BST)

________________________________________
Trial Means 
> tapply(dv, time, mean)
   1    2    3    4 
19.7 20.7 24.0 18.1 
> cat("\nSubject Means \n")

Subject Means 
> tapply(dv, subject, mean)
    1     2     3     4     5     6     7     8     9    10 
14.00 19.00 39.50 28.00 26.75 16.50 24.25 18.50  4.00 15.75 
> BSTmodel <- aov(dv ~ time + Error(subject/time))
> print(summary(BSTmodel))

Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  9   3318   368.7               

Error: subject:time
          Df Sum Sq Mean Sq F value Pr(>F)
time       3  186.3   62.09   1.042   0.39
Residuals 27 1609.0   59.59               


In the next example I have taken the data from Ex18.12 and added it to the data for Ex18.10. This gives me a group effect as well as a time effect. In addition, I have used the reshape() function to show a second way of creating a long file. I would suggest that until you become quite familiar with R you build the file directly as I did in the previous example rather than using the reshape function. I think that you will find that easier.


### Ex18.14
dataWide <- read.table("http://www.uvm.edu/~dhowell/fundamentals9/DataFiles/Ex18-14.dat", header = TRUE)
head(dataWide)
dataWide$Group <- factor(dataWide$Group)                       
dataLong <- reshape(data = dataWide, varying = 2:5, v.names = "dv", timevar
                   = "Trial", idvar = "Subject", ids = 1:20, direction = "long")
### This is a way of going from a wide file to a long one.
dataLong$Group <- factor(dataLong$Group, levels = 1:2, labels = c("BST", "Control"))
dataLong$Trial <- factor(dataLong$Trial)
cat("\nTrial Means\n")
tapply(dataLong$dv, dataLong$Trial, mean)
cat("\nSubject Means","\n")
tapply( dataLong$dv, dataLong$Subject, mean)

# Formula and calculation for long format
BST.aov <- aov(dataLong$dv ~ dataLong$Group * dataLong$Trial + Error(dataLong$Subject/dataLong$Group))

# Present the summary table (ANOVA source table)
print(summary(BST.aov))
rm(dataWide, dataLong, BST.aov)



Notice that the "aov()" function now contains three effects, one for Group, one for Trial, and one for the interaction. Here the error term is variability of subjects within groups. That is because the groups contain different subjects. The output is shown below.


Error: Subject
      Df Sum Sq Mean Sq
Group  1  98.17   98.17

Error: Subject:Group
      Df Sum Sq Mean Sq
Group  1  95.24   95.24

Error: Within
            Df Sum Sq Mean Sq F value Pr(>F)  
Group        1   1298  1298.1   5.291 0.0244 *
Trial        3    962   320.8   1.308 0.2789  
Group:Trial  3   1736   578.8   2.359 0.0789 .
Residuals   70  17175   245.4                 

The following code is presented without explanation. I simply shows the way to set up various repeated measures designs. Set it aside for future reference if you need it.


### Various Repeated Measures Designs

###One Within Subject Variable   No Between Subject Variable

data1 <- read.table("http://www.uvm.edu/~dhowell/fundamentals9/DataFiles/Tab18-1.dat", header = TRUE)
data1$person <- factor(1:25)
data1Long <- reshape(data = data1, varying = 1:5, v.names = "depression", timevar = "time", idvar = "person", direction = "long")
data1Long$time <- factor(data1Long$time)
data1Long$person <- factor(data1Long$person)
model1 <- aov(data1Long$depression ~ data1Long$time + Error(data1Long$person))
summary(model1)


### Now one between and one within
### Using King data
data2 <- read.table("http://www.uvm.edu/~dhowell/methods8/DataFiles/Tab14-4.dat", header = TRUE)
data2$person <- factor(1:24)
data2$Group <- factor(data2$Group)
data2Long <- reshape(data = data2, varying = 2:7, v.names = "outcome", timevar
                     = "time", idvar = "person", ids = "person", direction = "long")
data2Long$time <- factor(data2Long$time)
my.aov <- aov(data2Long$outcome ~ (data2Long$Group*data2Long$time) + Error(data2Long$person/(data2Long$time)))

###  THIS WORKS TOO  --  SO DOES THE NEXT

############## This was modified from one written by Joshua Wiley, in the Psychology Department at UCLA.
### Handles 2 Between and 1 Within
### Howell Stat Methods in Psych -- Table 14.7 ###
### Repeated Measures ANOVA with 2 Between and 1 Within variables
### Read in data, convert to 'long' format, and factor()

dat <- read.table("http://www.uvm.edu/~dhowell/methods7/DataFiles/Tab14-7.dat",
                  header = TRUE)
head(dat)                   
dat$subject <- factor(1:40)
dataLong <- reshape(data = dat, varying = 4:7, v.names = "outcome", timevar
                   = "time", idvar = "subject", ids = "subj", direction = "long")

dataLong$Condition <- factor(dataLong$Condition, levels = 1:2, labels = c("BST", "Control"))
dataLong$Sex <- factor(dataLong$Sex, levels = 1:2, labels = c("Male","Female"))
dataLong$time <- factor(dataLong$time)
str(dataLong)
# Actual formula and calculation
my.aov <- aov(outcome ~ (Condition*Sex*time) + Error(subject/(time)),data = dataLong)

# Present the summary table (ANOVA source table)
summary(my.aov)


# Two within subject variables, one between
# Data from Bouton & Schwartzentruber (1985)
# Methods8, p. 486
rm(list = ls())
dataBouton <- read.table("http://www.uvm.edu/~dhowell/methods8/DataFiles/Tab14-11.dat", header = TRUE)
subject <- factor(c(1:24))
## I should use "reshape," but I can't get it to work. This is pretty easy.

Phase <- factor(rep(1:2, each = 24, times = 4))
Cycle <- factor(rep(1:4, each = 48))
Group = factor(rep(1:3, each = 8,times = 8))
dv <- c(dataBouton$C1P1, dataBouton$C1P2, dataBouton$C2P1, dataBouton$C2P2, dataBouton$C3P1, dataBouton$C3P2, dataBouton$C4P1, dataBouton$C4P2)  
   # That sure was a pain in ght neck, but "with() didn't work.
Subj <- factor(rep(1:24, times = 8))

cat("Means and sd by Group \n")
tapply(dv, Group, mean); tapply(dv, Group, sd)
cat("\n Means and sd by Cycle\n")
tapply(dv, Cycle, mean); tapply(dv, Cycle, sd)
cat("\n Means and sd by Phase\n")
tapply(dv, Phase, mean); tapply(dv, Phase, sd)

options(contrasts = c("contr.sum","contr.poly"))
model1 <- aov(dv ~ (Group*Cycle*Phase) + Error(Subj/(Cycle*Phase)), contrasts = contr.sum)
print(summary(model1))
print(coefficients(model1))     # This will give mean deviations)


GreenBlueBar.gif GreenBlueBar.gif dch: