#################################################################
### CODE FOR EXPLORATORY MULTIVARIATE DATA ANALYSIS PRACTICAL ###
#################################################################

# PROBLEM SESSION 1

# QUESTION: A HIDDEN OUTLIER?
#Reads the data into R from a webpage:
experiment = read.table("https://www.mas.ncl.ac.uk/~nlf8/cpact/experiment.txt")
plot(ts(experiment[,1]),main="Plot of acidity values",xlab="Time",ylab="pH")

# PROBLEM SESSION 2

#QUESTION: SCATTERPLOT MATRICES
#Reads the data into R from a webpage:
heptathlon = read.table("https://www.mas.ncl.ac.uk/~nlf8/cpact/heptathlon.txt", header=TRUE)
#Look at the head of the dataset:
head(heptathlon)

#QUESTION: CORRELATION MATRICES
cor(heptathlon)

# PROBLEM SESSION 3

#QUESTION: PCA ON HEPTATHLON DATASET
#First, copy the original dataset into a new dataframe:
heptathlon2 = heptathlon

#Now re-scale, as we saw in the lecture:
heptathlon2$H100M = max(heptathlon2$H100M) - heptathlon2$H100M
heptathlon2$R200M = max(heptathlon2$R200M) - heptathlon2$R200M
heptathlon2$R800M = max(heptathlon2$R800M) - heptathlon2$R800M

#Perform the PCA:
(pca = prcomp(x=heptathlon2, scale=TRUE))

#Produce a 'score plot':
plot(pca$x[,1], pca$x[,2], xlab="First PC", ylab="Second PC")
text(pca$x[,1], pca$x[,2], labels=rownames(heptathlon2))

# PROBLEM SESSION 4

# QUESTION: PCA ON US ARRESTS DATASET
USA = read.table("https://www.mas.ncl.ac.uk/~nlf8/cpact/USArrests.txt", header=TRUE)



###############################################
### CODE FOR ANALYSIS OF VARIANCE PRACTICAL ###               
###############################################

# PROBLEM SESSION 1

# QUESTION: ONE-WAY ANOVA
#Typing the data into R manually:
yield = c(4.85, 4.89, 5.60, 5.08, 3.42, 9.43, 6.77, 5.39, 7.39, 4.67, 5.49, 5.96, 3.71, 5.51, 5.50, 4.90, 6.00, 5.81, 4.68)
catalyst = c(rep("A", 5), rep("B",4), rep("C",5), rep("D",5))

fit = aov(yield~catalyst)
summary(fit)
(tukeytest = TukeyHSD(fit))
plot(tukeytest)

# QUESTION: RANDOMISED BLOCKS ANOVA
#Read the suncream data into R from a webpage:
suncream = read.table("https://www.mas.ncl.ac.uk/~nlf8/cpact/suncream.txt", header=TRUE)

# PROBLEM SESSION 2

# QUESTION: FACTORIAL EXPERIMENTS

# Read the dataset in from a webpage:
rats = read.table("https://www.mas.ncl.ac.uk/~nlf8/cpact/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")

# PROBLEM SESSION 3

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