#    TWO WAY ANOVA in R SCRIPT FILE.
#    www_statstutor_ac_uk Community Project.
#    Sofia Maria Karadimitriou and Ellen Marshall, Sheffield University.
#     Reviewed by Jim Bull, University of Swansea.
#    Dataset: diet csv.
#    Resource: Two way ANOVA in R.

#Open the diet dataset which is saved as a csv file and call it dietR.  
#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.
dietR<-read.csv("E:\\diet.csv",header=T,sep=",")

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

#calculate the weight lost by person (difference in weight before and after the diet) and add to the dataset.
dietR$weightlost<-pre.weight-weight6weeks
#attach again your data and remove the missing values.
dietR<-na.omit(dietR)
attach(dietR)

#Carrying out the two-way ANOVA, telling R than Diet and gender are categorical using as.factor().
anova2<-aov(weightlost~as.factor(gender)*as.factor(Diet),data=dietR)

#Ask for the residuals (difference between each individual and their group mean).
#save the residuals of the anova in a separate object.
res<-residuals(anova2)

#Checking the normality of the residuals.
hist(res,main="Histogram of residuals",xlab="Residuals")

#The Levene's test for equality of variances is in the car package.  Load the additional library car.
library(car)
#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 'car' and click.  A lot of extra menus will download. Then try library(car) again.

#Carry out Levene's test.
leveneTest(weightlost~as.factor(gender)*as.factor(Diet),data=dietR)


#Note: Rstudio currently has some issues with not all commands will work.

#If the test doesn't work, compare standard deviations.

#To see the ANOVA output use summary().
summary(anova2)

#To produce an interaction plot.  First give R the labels for gender.
# The factor command uses variable<-factor(variable,c(category numbers),labels=c(category names)).
gender<-factor(gender,c(0,1),labels=c('Female','Male'))

#To generate an interaction plot for weight lost and the interaction of diet and gender 
#use the interaction.plot(independent1,independent2,dependent) command
interaction.plot(Diet, gender, weightlost,main="Mean weight lost by diet and gender",ylab = "mean of weightlost",xlab = "Type of Diet")

## Changing the type of lines ##.
#lty specifies the type of lines
#lwd specifies the thickness of the lines. If you want the lines to be different, try lwd=c(1,3)
#col specifies the colours of the lines.
interaction.plot(Diet, gender, weightlost,lty = c(1, 12),col=c(2:3),lwd = 3,ylab = "mean of weightlost",
xlab = "Type of Diet",main="Mean weight lost by diet and gender")

## Adding and editing points ##.
#Type = b adds points to the line.
#pch controls the shape of the points.
interaction.plot(Diet, gender, weightlost, main="Mean weight lost by diet and gender", type="b",pch=c(18,24),lty = c(1, 12),col=c(2:3),lwd = 3,ylab = "mean of weightlost",
xlab = "Type of Diet",main="Mean weight lost by diet and gender")

## Legend formatting ##.
#trace.label controls the name of the legend eg Gender.
#leg.bty ="o" puts a box around the legend.
#leg.bg=" " controls the background colour of the legend.
#Use ?interaction.plot to find out more.

interaction.plot(Diet, gender, weightlost, main="Mean weight lost by diet and gender",type="b",pch=c(18,24),lty = c(1, 12),col=c(2:3),lwd = 3,trace.label = "Gender",leg.bty="o", leg.bg="beige",ylab = "mean of weightlost",
xlab = "Type of Diet")

#If there are significant results in the ANOVA, post hoc test should be carried out.
#To carry out Tukey's post hoc adjustments for the pairwise comparisons.
#If the interaction is NOT significant, interpret the main effects.
#if it is significant interpret the interaction post hoc tests only.
TukeyHSD(anova2)

#-------------------------------------------------------------------.
#If you prefer to split the file and carry out separate ANOVA's by gender follow these instructions.
#Splitting the data set into two subsets.
#Save the female data into a new dataset.
females<-split(dietR,gender)[[1]]
#Save the male data into a new dataset.
males<-split(dietR,gender)[[2]]
#Do two separate ANOVAs.
fitfemale<-aov(weightlost~as.factor(Diet),data=females)
fitmale<-aov(weightlost~as.factor(Diet),data=males)
#Do two separate Tukeys.
TukeyHSD(fitfemale)
TukeyHSD(fitmale)