suncream = read.table("https://www.mas.ncl.ac.uk/~nlf8/sellafield/suncream.txt", header = TRUE)

suncream2 = stack(suncream)

head(suncream2)

boxplot(suncream2[,1]~suncream2[,2], main="", xlab="Brand of suncream", ylab="Burning index")

fit = aov(suncream2[,1]~suncream2[,2])

summary(fit)

(tukeytest = TukeyHSD(fit))

plot(tukeytest)


carrots = read.table("https://www.mas.ncl.ac.uk/~nlf8/sellafield/carrots.txt", header = TRUE)
head(carrots)
carrots$Plot = as.factor(carrots$Plot)
boxplot(carrots$Yield~carrots$Strain, xlab="Strain", ylab="Yield (kg)")

aggregate(carrots$Yield, list(carrots$Strain), FUN=mean)

aggregate(carrots$Yield, list(carrots$Strain), FUN=sd)

fit = aov(carrots$Yield~carrots$Strain+carrots$Plot)
summary(fit)

fit = aov(carrots$Yield~carrots$Strain)
(tukeytest = TukeyHSD(fit))



#read in the data from a webpage:
rats = read.table("https://www.mas.ncl.ac.uk/~nlf8/rats.txt", header = TRUE)

#Make sure R knows that the source, and level, of protein are both factors:
rats$Source = as.factor(rats$Source)
rats$Level = as.factor(rats$Level)

#Produce some summary plots - checking for Normality along the way!
par(mfrow=c(1,2))
#Boxplots for the data by source of protein:
boxplot(rats$Gains~rats$Source)
#And now be level of protein:
boxplot(rats$Gains~rats$Level)

#How about Normal Q-Q plots from last week? Easy!
par(mfrow=c(2,3))
qqnorm(rats$Gains[rats$Source==1], main="Source = Beef")
qqline(rats$Gains[rats$Source==1])
qqnorm(rats$Gains[rats$Source==2], main="Source = Cereal")
qqline(rats$Gains[rats$Source==2])
qqnorm(rats$Gains[rats$Source==3], main="Source = Pork")
qqline(rats$Gains[rats$Source==3])

qqnorm(rats$Gains[rats$Level==1], main="Level = High")
qqline(rats$Gains[rats$Level==1])
qqnorm(rats$Gains[rats$Level==2], main="Level = Low")
qqline(rats$Gains[rats$Level==2])

#Let's now produce some numerical summaries - we can use any function we like, but let's check the assumption of equal variances:
aggregate(rats$Gains, list(rats$Source), FUN = sd)
aggregate(rats$Gains, list(rats$Level), FUN = sd)

#Now let's do the ANOVA!
fit = aov(rats$Gains~rats$Source+rats$Level+rats$Source*rats$Level)
summary(fit)

#And finally, an interaction plot:
par(mfrow=c(1,1))
interaction.plot(x.factor = rats$Source, #x-axis variable
                 trace.factor = rats$Level, #variable for lines
                 response = rats$Gains, #y-axis variable
                 fun = mean, #metric to plot
                 ylab = "Weight gain (g)",
                 xlab = "Source of protein",
                 col = c("pink", "blue"),
                 lty = 1, #line type
                 lwd = 2, #line width
                 trace.label = "Level of protein")