Repeated measures ANOVA using R

Data
 

# 16 adults followed a diet for 12 months. 
# 3 weight measurements at time 0, 6 and 12 months (weight1, weight2 and weight3 resp.)
gender = c(0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
gender = factor(gender, levels = c(0, 1), labels = c("Woman", "Man"))
contrasts(gender) = contr.helmert
weight1 = c(81.2, 76.7, 75.7, 81.2, 71.7, 71.2, 68.5, 89.8, 107.5, 105.7, 99.3, 100.7, 90.3, 105.7, 98.0, 116.6)
weight2 = c(78.0, 73.0, 73.0, 78.5, 69.9, 64.9, 63.5, 87.1, 102.1, 102.5, 97.1, 95.3, 87.5, 102.5, 93.4, 112.9)
weight3 = c(78.4, 72.1, 73.7, 78.5, 69.7, 65.4, 63.3, 85.5, 101.3, 102.2, 95.7, 93.9, 86.5, 102.5, 93.9, 113.9)
dataF = data.frame(gender, weight1, weight2, weight3)

1. Summarize data

my.summary = function(avar){
  return (paste(round(mean(avar), 1), " (", round(sd(avar), 1), ")"))
}

for(column in 2:4){
  print(my.summary(dataF[,column]))
  print(my.summary(dataF[,column][dataF$gender == "Woman"]))
  print(my.summary(dataF[,column][dataF$gender == "Man"])) 
}

2. Define time factor.

time = factor(rep("weight_measurement", 3), levels=c("weight1", "weight2", "weight3"))
time[1] = "weight1"; time[2] = "weight2"; time[3] = "weight3"
idata = data.frame(time)

3. Compute repeated measures ANOVA

mod.ok = lm(cbind(weight1, weight2, weight3) ~ gender, data=dataF)
library(car)
av.ok = Anova(mod.ok, idata=idata, idesign=~time, type = "III") 
summary(av.ok, multivariate=FALSE)

4. Residuals

allresiduals = as.data.frame(residuals(mod.ok))

5. Normality of residuals

library(moments)
lapply(allresiduals, skewness)
lapply(allresiduals, kurtosis)
lapply(allresiduals, shapiro.test)

6. Restructure of the data (3 weight columns to one)

library(afex)
dataF$id = rep(1:16)
data.long = data.frame(dataF[5], dataF[1], stack(dataF[2:4]))

6.1 Rename of columns (not absolutely necessary)

library(plyr)
data.long = rename(data.long, c("values"="weight", "ind"="time"))

7. Computation of p according to Huynh - Feldt correction
(necessary when sphericity is not assumed. The user should be aware from step 3)

rmd2 = aov_ez("id", "weight", data.long, between = c("gender"), within = c("time"), anova_table=list(correction = "HF", es = "none"))
nice(rmd2)

8. Label modification (weight1, weight2, weight3 to 1st, 2nd, 3rd)

data.long$time=revalue(data.long$time, c("weight1"="1st", "weight2"="2nd", "weight3"="3rd"))

9. Plot

par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
interaction.plot(data.long$time, data.long$gender, data.long$weight, type="b", col=c("red","blue"), legend=F,
                 lty=c(1,2), lwd=2, pch=c(18,24), xlab="Measurement",  ylab="Weight",
                 main="Weight evolution by gender", ylim = c(65, 105))

9.1 Add legend

legend("topright",c("Woman","Man"), bty="n", lty=c(1,2), lwd=2, pch=c(18,24), col=c("red","blue"), title="Gender",inset=c(-0.5,0))

10. Locate the position of labels: After the command "locator(6)" the user will select with the mouse on the plot the 6 points where labels will appear.
Woman - 1st, Woman - 2nd, Woman - 3rd, Man - 1st, Man - 2nd, Man - 3rd.

coord_for_labels = locator(6)

9.1 Position labels

counter = 1
for (agender in levels(data.long$gender)) {
  for (ameasurement in levels(data.long$time)) {
      meantoshow = mean(data.long[ which(data.long$gender==agender & data.long$time == ameasurement), ]$weight)
      text(coord_for_labels$x[counter], coord_for_labels$y[counter], round(meantoshow, 1), pos=3) 
      counter = counter + 1
    }
}

The above example is contained in the paragraph 5.4 of the book "Στατιστική ανάλυση με τη γλώσσα R" (in Greek, ISBN: 978-960-93-9445-1) published in Thessaloniki, 2017.