repeatedANOVA.R

########################################
########################################
# 
# Repeated Measures ANOVAs in R
# 
########################################
########################################

########################################
# Construct a sample dataset
########################################

set.seed(5250)

# Say we're measuring changes in stress on
# a numerical measure (1 = low stress; 100 =
# high stress) in response to viewing
# a happy or angry image accompanied by
# either Disney or horror music
myData <- data.frame(PID = rep(seq(from = 1,
                                   to = 50, by = 1), 20),
                     stress = sample(x = 1:100,
                                     size = 1000,
                                     replace = TRUE),
                     image = sample(c("Happy", "Angry"),
                                    size = 1000,
                                    replace = TRUE),
                     music = sample(c("Disney", "Horror"),
                                    size = 1000,
                                    replace = TRUE)
)

# Convert to factors
myData <- within(myData, {
  PID   <- factor(PID)
  image <- factor(image)
  music <- factor(music)
})

# Order for tidy viewing
myData <- myData[order(myData$PID), ]


########################################
# Extract condition means
########################################

# Group outcomes by participant, music
# condition, and image condition and
# compute the mean stress score
myData.mean <- aggregate(myData$stress,
                         by = list(myData$PID, myData$music,
                                   myData$image),
                         FUN = 'mean')

# Rename the columns for ease of use
colnames(myData.mean) <- c("PID","music","image","stress")

# Order data for tidy viewing
myData.mean <- myData.mean[order(myData.mean$PID), ]


########################################
# Construct the ANOVA
########################################

# Include error term for between-participant
# error across each within-subjects independent
# variable. Asterisks represent interaction
# effects
stress.aov <- with(myData.mean,
                   aov(stress ~ music * image +
                         Error(PID / (music * image)))
)

# Summarize results of ANOVA
summary(stress.aov)


########################################
# Include between-subjects effect
########################################

set.seed(5250)

# Basically our original data set, but with
# time of day added in as a between-subjects
# variable
myData <- data.frame(PID = rep(seq(from = 1,
                                   to = 50, by = 1), 20),
                     stress = sample(x = 1:100,
                                     size = 1000,
                                     replace = TRUE),
                     image = sample(c("Happy", "Angry"),
                                    size = 1000,
                                    replace = TRUE),
                     music = sample(c("Disney", "Horror"),
                                    size = 1000,
                                    replace = TRUE),
                     time = rep(sample(c("Day", "Night"),
                                       size = 50,
                                       replace = TRUE), 2))

# Convert to factors
myData <- within(myData, {
  PID   <- factor(PID)
  image <- factor(image)
  music <- factor(music)
  time  <- factor(time)
})

# Order data
myData <- myData[order(myData$PID), ]

# Aggregate mean stress values by participant
# by condition
myData.mean <- aggregate(myData$stress,
                         by = list(myData$PID, myData$music,
                                   myData$image, myData$time),
                         FUN = 'mean')

# Name the columns for ease of use
colnames(myData.mean) <- c("PID", "music", "image",
                           "time", "stress")

# Order data
myData.mean <- myData.mean[order(myData.mean$PID), ]

# Construct aov object, including between-subjects
# time factor. We do not include this in our error
# term -- only the within-subjects factors
stress.aov <- with(myData.mean,
                   aov(stress ~ time * music *
                         image + Error(PID / (music * image))))

# Summarize findings
summary(stress.aov)


########################################
# What does significance look like?
########################################

set.seed(982)

# Basically our original data set, but with
# time of day added in as a between-subjects
# variable
myData <- data.frame(PID = rep(seq(from = 1,
                                   to = 50, by = 1), 20),
                     stress = sample(x = 1:100,
                                     size = 1000,
                                     replace = TRUE),
                     image = sample(c("Happy", "Angry"),
                                    size = 1000,
                                    replace = TRUE),
                     music = sample(c("Disney", "Horror"),
                                    size = 1000,
                                    replace = TRUE),
                     time = rep(sample(c("Day", "Night"),
                                       size = 50,
                                       replace = TRUE), 2))

# Convert to factors
myData <- within(myData, {
  PID   <- factor(PID)
  image <- factor(image)
  music <- factor(music)
  time  <- factor(time)
})

myData$music[myData$stress <= 25] <- "Disney"

myData.mean <- aggregate(myData$stress,
                         by = list(myData$PID, myData$music,
                                   myData$image),
                         FUN = 'mean')

# Name the columns for ease of use
colnames(myData.mean) <- c("PID", "music", "image", "stress")

# Remove a participant with missing data
myData.mean <- myData.mean[myData.mean$PID != 18, ]

# Construct aov object, including between-subjects
# time factor. We do not include this in our error
# term -- only the within-subjects factors
stress.aov <- with(myData.mean,
                   aov(stress ~ music * image +
                         Error(PID / (music * image))))

# Summarize findings
summary(stress.aov)