#Let's try regression on the heptathlon data!
final_scores = c(6955, 6649, 6599, 6576, 6480, 6452, 6414, 6392,6383,6324, 6319,6300, 6267,6232,6198,6109,6107,6084,6030,6013,6013,5947,5935,5873,5856,5846,5688,5338, 5319)

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

#Re-scale the running events and store new dataset in "heptathlon2"
heptathlon2 = heptathlon
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))
par(mfrow=c(1,1))
plot(pca$x[,1], pca$x[,2], xlab="First PC", ylab="Second PC")
text(pca$x[,1], pca$x[,2], labels=rownames(heptathlon2))

attach(heptathlon2)

#Multiple linear regression:
fit.ml = lm(final_scores~H100M+R200M+R800M+SP+JT+LJ+HJ)

summary(fit.ml)

#Regression on just first two PCs:
fit.pc = lm(final_scores~pca$x[,1]+pca$x[,2])
summary(fit.pc)

predicted.lm = -215.8637+147.0022*H100M+101.6977*R200M+12.9245*R800M+66.3866*SP+19.6023*JT+269.3657*LJ+1345.9647*HJ
predicted.pc = 6162-201.325*pca$x[,1]-133.614*pca$x[,2]

par(mfrow=c(2,2))

plot(final_scores, predicted.lm, xlab="Actual score", ylab="Estimated from regression on all variables",pch=19, col=rainbow(29))
abline(0,1)

plot(final_scores, predicted.pc, xlab="Actual score", ylab="Estimated from regression on PCs",pch=19, col=rainbow(29))
abline(0,1)

plot(predicted.pc~predicted.lm, xlab="Regression on original variables", ylab="Regression on first two PCs", pch=19, col=rainbow(29))
abline(0,1)


###Further regression
albumin = seq(0, 20, 2)
peak.height = c(0, 80.33, 135.7, 207, 243.3, 303, 329, 370.3, 406, 431.7, 444.3)
plot(albumin, peak.height, main="Average peak height versus Albumin", xlab="Albumin", ylab="Average peak height", pch=16)

albumin.sq = albumin^2

fit1 = lm(peak.height~albumin)
fit2 = lm(peak.height~albumin+albumin.sq)

lines(fitted(fit2)~albumin, col="blue", lwd=2, lty=2)

leukaemia = read.table("https://www.mas.ncl.ac.uk/~nlf8/sellafield/leukaemia.txt", header = TRUE)
attach(leukaemia)
fit3 = glm(REMISS~CELL+SMEAR+INFIL+LI+BLAST+TEMP, family = binomial)
summary(fit3)

fit4 = glm(REMISS~LI, family = binomial)
summary(fit4)

plot(LI, REMISS, pch=19,lwd=3)
curve(exp(fit4$coef[1]+fit4$coef[2]*x)/(1+exp(fit4$coef[1]+fit4$coef[2]*x)), from=min(LI), to = max(LI), col="blue", lty=2, lwd=2,add=TRUE)

###Problem session in further regression###
hp = read.table("https://www.mas.ncl.ac.uk/~nlf8/sellafield/potter.txt", header=TRUE)
attach(hp)
fit1 = glm(Potter~Age+Sex, family=binomial)