#    Repeated measures ANOVA in R SCRIPT FILE.
#    www_statstutor_ac_uk Community Project.
#    Sofia Maria Karadimitriou and Ellen Marshall, Sheffield University.
#    Reviewed by Jonathan Gillard, University of Cardiff.
#    Dataset: cholesterol csv.
#    Resource: Repeated measures ANOVA in R.

#Open the cholesterol dataset which is saved as a csv file and call it cholA.  
#If your file is saved as a standard Excel file, save it as a csv file first.  
#You will need to change the command depending on where you have saved the file 
# and what you called the file e.g.the dataset'stcp-Rdataset-cholesterol'is in the D drive.
cholA<-read.csv("D:\\stcp-Rdataset-cholesterol.csv",header=T)

#Look at the data.
cholA

#Tell R we are using the cholA dataset until further notice using attach.
#This means that 'Before' can be used instead of cholA$Before.
attach(cholA)

#### Checking the assumption of normality ###.

#One of the assumptions for repeated measures is that the residuals at each time point are normally distributed.
#As R does not have this feature, check that cholesterol at each time point is normally distributed.
par(mfrow=c(2,2))
hist(Before)
hist(After4weeks)
hist(After8weeks)

#######Changing to long format for repeated measures ########.
#The file currently has the 3 time points (before, after 4 weeks and after 8 weeks).
#The data set needs to be reformatted from wide format to long format to run repeated measures.

#Tell R you want 53 rows of data and 4 columns using nrow and ncol.
#Replicate the subject ID for the 18 participants 3 times using rep(1:18,3).
#Combine the measurements for each time point into one column.
#Create a comumn to identify which time point each observation came from using rep(1,18), rep(2,18) etc.
#Replicate the margarine type for the 18 participants 3 times using rep(1:18,3).

cholB<-matrix(nrow=54,ncol=4,c(rep(1:18,3),cholA$Before,cholA$After4weeks,cholA$After8weeks,rep(1,18),rep(2,18),rep(3,18),rep(cholA$Margarine,3)))
#Give each column a label.
colnames(cholB)<-c('subject','cholesterol','time','Margarine')
#Tell R cholB is a dataset.
cholB<-data.frame(cholB)
#R assumes all numeric values are continuous so tell it that time, subject and margarine are categorical variables.
cholBsubject<-factor(cholB$subject)
cholB$time<-factor(cholB$time)
cholB$Margarine<-factor(cholB$Margarine)

#Tell R we are using the reformatted cholesterol dataset cholB until further notice using attach.
#This means that 'time' can be used instead of cholB$time.
#Look at the dataset to see if it looks ok.
cholB


#There is no procedure within standard R to carry out repeated so the package ez needs to be loaded.
#Go to 'Packages' --> 'Load packages' and choose the package ez which will be saved in your temporary files.
#Load the additional library ez.
library(ez)
#If this command does not work, you will need to go to the Packages --> Install package(s) and select the UK (London)CRAN mirror.
#Then look for the package 'ez' and click.  A lot of extra menus will download. Then try library(ez) again.
#For Rstudio, go to Tools --> Install Packages, type 'ez' and the 'Install'.
#Then select ez from the packages window.
#Note: Earlier versions of R may not allow this package to download fully.


#Run the repeated measures ANOVA, giving it a name and ask to view the output.
repeat1<-ezANOVA(data=cholB,dv=.(cholesterol),wid=.(subject),within=.(time),type=3)
repeat1

#If the ANOVA is significant,pairwise paired t-tests need to be carried out with a Bonferroni adjustment.
#See the 'Paired t-test in R' resource and script file for more details.
#The Bonferroni adjustment multiplies the t-test p-values by the number of tests being carried out.
#For example, if there are are 3 groups, there are 3 tests 1 vs 2, 1 vs 3 and 2 vs 3.
#The data needs to be formatted so the cholesterol by time point is in different columns (as in cholA).

#Tell R we are using the cholA dataset until further notice using attach.
#This means that 'Before' can be used instead of cholA$Before.
attach(cholA)

attach(cholA)

#Pairwise paired t-tests, Saving then into objects
 t1<-t.test(Before,After4weeks,paired=T)
t2<-t.test(Before,After8weeks,paired=T)
 t3<-t.test(After4weeks,After8weeks,paired=T)

#To report the full results of each test, request each using their name.
t1
t2
t3

#To summarise key results for each.
 #Report the pvalues of each test in one row.
 pvalues<-c(t1$p.value,t2$p.value,t3$p.value)
 #As multiple tests are being carried out, an adjustment to the paired t-test results needs to be carried out.
#The easiest adjustment is the Bonferroni which multiplies each p-value by the total number of pairs being tested.
#Here, time points 1 vs 2, 1 vs 3 and 2 vs 3 are being compared so there are 3 tests.
#Adjust the p-values using the bonferroni correction (round by 4 decimals).
 round(p.adjust(pvalues,'bonferroni',3),4)



