#########################################################
#########################################################
##########                                      #########
##########	R Scripts for MATH 204 Winter 2013  #########
##########                                      #########
#########################################################
#########################################################

#Simple Linear Regression
attach(Viscosity)
plot(pcv,viscosity,pch=19)
detach(Viscosity)

plot(Viscosity$pcv,Viscosity$viscosity,pch=19,xlab="pcv",ylab="Viscosity")

summary(Viscosity$pcv)
summary(Viscosity$viscosity)
cor.test(Viscosity$pcv,Viscosity$viscosity)

m0.visc<-lm(viscosity~pcv,data=Viscosity)
abline(coef(m0.visc))
be0.hat<-round(coef(m0.visc)[1],3)
be1.hat<-round(coef(m0.visc)[2],3)
text(40,5.9,substitute(hat(beta)[0] == b, list(b = be0.hat)),adj=0)
text(40,5.7,substitute(hat(beta)[1] == b, list(b = be1.hat)),adj=0)
summary(m0.visc)


#########################################################
#########################################################
#Coleman data
attach(Coleman)
plot(status,testscore,pch=19)
detach(Coleman)

m0.cole<-lm(testscore~status,data=Coleman)
summary(m0.cole)
abline(coef(m0.cole))

summary(aov(m0.cole))

#########################################################
#########################################################
#Hooker data
attach(Hooker)
plot(pressure,temp,pch=19)
detach(Hooker)

m0.hooker<-lm(temp~pressure,data=Hooker)
summary(m0.hooker)
abline(coef(m0.hooker))

summary(aov(m0.hooker))

#########################################################
#Polynomial
attach(Hooker)
plot(pressure,temp,pch=19)
detach(Hooker)
m1.hooker<-lm(temp~pressure+I(pressure^2),data=Hooker)
summary(m1.hooker)

#Fitted values for quadratic regression
xp<-seq(15,30,by=0.01)
m1.fit<-cbind(rep(1,length(xp)),xp,xp^2) %*% coef(m1.hooker)
lines(xp,m1.fit,col="red",lty=2)
abline(coef(m0.hooker))
legend(16,210,c("Straight Line","Quadratic"),lty=c(1,2),col=c("black","red"))
title("Hooker data")

par(mfrow=c(2,2))
attach(Hooker)
plot(pressure,temp,pch=19)
abline(coef(m0.hooker))
title("Straight-line fit")
plot(pressure,residuals(m0.hooker),pch=19,ylab="Residual",ylim=range(-1.5,1.5))
title("Residuals from straight-line fit")
abline(h=0,lty=2)
plot(pressure,temp,pch=19)
lines(xp,m1.fit,lty=2)
title("Quadratic fit")
plot(pressure,residuals(m1.hooker),pch=19,ylab="Residual",ylim=range(-1.5,1.5))
title("Residuals from quadratic fit")
abline(h=0,lty=2)
detach(Hooker)
par(mfrow=c(1,1))
anova(m0.hooker,m1.hooker)

#########################################################
#########################################################
#Multiple regression

mv.visc<-lm(viscosity~pcv+plasfib+plaspro,data=Viscosity)
summary(mv.visc)
m1.visc<-lm(viscosity~pcv,data=Viscosity)
summary(m1.visc)
m2.visc<-lm(viscosity~plasfib,data=Viscosity)
summary(m2.visc)
m3.visc<-lm(viscosity~plaspro,data=Viscosity)
summary(m3.visc)


anova(m1.visc,mv.visc)
anova(m2.visc,mv.visc)
anova(m3.visc,mv.visc)

#########################################################
#########################################################
#Creating dummy variables for factor predictors
#Diabetes data
summary(lm(glutest~group,data=Diabetes))

summary(lm(glutest~-1+group,data=Diabetes))

Diabetes.copy<-Diabetes
Diabetes.copy$group2<-factor(Diabetes.copy$group,c("Normal","Chemically Diabetic","Overt Diabetic"))
summary(lm(glutest~group2,data=Diabetes.copy))
summary(lm(glutest~-1+group2,data=Diabetes.copy))

boxplot(glutest~group2,data=Diabetes.copy)
title("Glucose level on test \n (original scale)")

Diabetes.copy$group3 <- relevel(Diabetes.copy$group, ref="Normal")
summary(lm(glutest~group3,data=Diabetes.copy))
summary(lm(glutest~-1+group3,data=Diabetes.copy))


leveneTest(glutest~group2,data=Diabetes.copy)

boxplot(log(glutest)~group2,data=Diabetes.copy)
title("Glucose level on test \n (log scale)")

leveneTest(log(glutest)~group2,data=Diabetes.copy)

Diabetes.copy$X1<-as.numeric(Diabetes.copy$group=="Normal")
Diabetes.copy$X2<-as.numeric(Diabetes.copy$group=="Chemically Diabetic")
Diabetes.copy$X3<-as.numeric(Diabetes.copy$group=="Overt Diabetic")

summary(lm(glutest~-1+X1+X2+X3,data=Diabetes.copy))

#############################################################
#Factor predictor and covariate
attach(Diabetes)
plot(log(glufast),log(glutest),pch=19)
plot(log(glufast),log(glutest),type="n")
points(log(glufast[group=="Normal"]),log(glutest[group=="Normal"]),pch=17)
points(log(glufast[group=="Chemically Diabetic"]),log(glutest[group=="Chemically Diabetic"]),pch=4)
points(log(glufast[group=="Overt Diabetic"]),log(glutest[group=="Overt Diabetic"]),pch=3)
legend(4.2,7.3,c("Normal","Chem. Diabetic","Overt Diabetic"),pch=c(17,4,3))
detach(Diabetes)

m.Diabetes<-lm(log(glutest)~log(glufast)*group2,data=Diabetes.copy)
summary(m.Diabetes)
summary(aov(m.Diabetes))

attach(Diabetes)
par(mfrow=c(2,2))
plot(log(glufast),log(glutest),type="n")
points(log(glufast[group=="Normal"]),log(glutest[group=="Normal"]),pch=17)
be0<-coef(m.Diabetes)[1]
be1<-coef(m.Diabetes)[2]
abline(be0,be1)
title("Normal")
legend(4.2,7.3,c("Normal"),pch=c(17))

plot(log(glufast),log(glutest),type="n")
points(log(glufast[group=="Chemically Diabetic"]),log(glutest[group=="Chemically Diabetic"]),pch=4)
be0<-coef(m.Diabetes)[1]+coef(m.Diabetes)[3]
be1<-coef(m.Diabetes)[2]+coef(m.Diabetes)[5]
abline(be0,be1)
title("Chemically Diabetic")
legend(4.2,7.3,c("Chem. Diabetic"),pch=c(4))

plot(log(glufast),log(glutest),type="n")
points(log(glufast[group=="Overt Diabetic"]),log(glutest[group=="Overt Diabetic"]),pch=3)
be0<-coef(m.Diabetes)[1]+coef(m.Diabetes)[4]
be1<-coef(m.Diabetes)[2]+coef(m.Diabetes)[6]
abline(be0,be1)
title("Overt Diabetic")
legend(4.2,7.3,c("Overt Diabetic"),pch=c(3))

plot(log(glufast),log(glutest),type="n")
be0<-coef(m.Diabetes)[1]
be1<-coef(m.Diabetes)[2]
abline(be0,be1)
be0<-coef(m.Diabetes)[1]+coef(m.Diabetes)[3]
be1<-coef(m.Diabetes)[2]+coef(m.Diabetes)[5]
abline(be0,be1,lty=2)
be0<-coef(m.Diabetes)[1]+coef(m.Diabetes)[4]
be1<-coef(m.Diabetes)[2]+coef(m.Diabetes)[6]
abline(be0,be1,lty=3)
legend(4.2,7.3,c("Normal","Chem. Diabetic","Overt Diabetic"),lty=c(1,2,3))
par(mfrow=c(1,1))
detach(Diabetes)

#Results separately for groups
summary(lm(log(glutest)~log(glufast),data=Diabetes,subset=(group=="Normal")))
summary(lm(log(glutest)~log(glufast),data=Diabetes,subset=(group=="Chemically Diabetic")))
summary(lm(log(glutest)~log(glufast),data=Diabetes,subset=(group=="Overt Diabetic")))

#Model selection
#Try to drop the interaction
drop1(m.Diabetes,test="F")
#Cannot drop !

m0.diabetes<-lm(log(glutest)~1,data=Diabetes.copy)
m1.diabetes<-lm(log(glutest)~group2,data=Diabetes.copy)
m2.diabetes<-lm(log(glutest)~log(glufast),data=Diabetes.copy)
m3.diabetes<-lm(log(glutest)~group2+log(glufast),data=Diabetes.copy)
m4.diabetes<-lm(log(glutest)~group2*log(glufast),data=Diabetes.copy)

#############################################################
#Unbalanced designs

m.potato<-lm(leak~potato*regime*temp,data=Potato)
summary(m.potato)

summary(aov(m.potato))

drop1(m.potato,test="F")
#Can drop 3-way interaction

m1.potato<-update(m.potato,~.-potato:regime:temp)

drop1(m1.potato,test="F")

#Can drop 2-way interaction regime:temp
m2.potato<-update(m1.potato,~.-regime:temp)
drop1(m2.potato,test="F")

#Can drop 2-way interaction potato:temp
m3.potato<-update(m2.potato,~.-potato:temp)
drop1(m3.potato,test="F")

summary(m3.potato)

summary(aov(m3.potato))

####################################################################
#Two factors, one covariate

m.Task<-lm(Errors~Group*Task*Distract,data=TaskDistract)
summary(m.Task)

summary(aov(m.Task))
drop1(m.Task,test="F")

#Can drop 3-way interaction

m1.Task<-update(m.Task,~.-Group:Task:Distract)
drop1(m1.Task,test="F")

#Maybe drop Group:Distract
m2.Task<-update(m1.Task,~.-Group:Distract)
drop1(m2.Task,test="F")

summary(m1.Task)
summary(aov(m1.Task))

m3.Task<-lm(Errors~Group+Task+Distract,data=TaskDistract)
summary(m3.Task)

###################################################################
#Prediction

#Prediction A: Subject in Group "NonSmoker", doing Task "Driving Simulation", with Distraction level 180.
#Prediction B: Subject in Group "Delayed Smoker", doing Task "Cognitive Task", with Distraction level 180.
#Prediction C: Subject in Group "Active Smoker", doing Task "Pattern Recognition", with Distraction level 100.

newTask<-data.frame(Group=c("NonSmoker","Delayed Smoker","Active Smoker"),
                    Task=c("Driving Simulation","Cognitive Task","Pattern Recognition"),
                    Distract=c(180,180,100))
predict(m1.Task,newdata=newTask,interval="confidence")
predict(m1.Task,newdata=newTask,interval="prediction")


#Prediction plus confidence interval
predict(m.Task,newdata=newTask,interval="confidence")
predict(m3.Task,newdata=newTask,interval="confidence")

#Prediction plus prediction interval
predict(m.Task,newdata=newTask,interval="prediction")
predict(m3.Task,newdata=newTask,interval="prediction")


###################################################################
#Stepwise selection using AIC

best1<-step(m.Task,k=2,direction="backward")
best2<-step(m.Task,k=2,direction="both")
formula(best1)
formula(best2)


#Stepwise selection using BIC

best3<-step(m.Task,k=log(nrow(TaskDistract)),direction="backward")
best4<-step(m.Task,k=log(nrow(TaskDistract)),direction="both")
formula(best3)
formula(best4)


###################################################################
###################################################################
#Residuals


Task.full<-lm(formula(m.Task),data=TaskDistract)
Task.best1<-lm(formula(best1),data=TaskDistract)
Task.best3<-lm(formula(best3),data=TaskDistract)

res.Task<-residuals(m.Task)
res.best1<-residuals(Task.best1)
res.best3<-residuals(Task.best3)

stres.Task<-rstandard(m.Task)
stres.best1<-rstandard(Task.best1)
stres.best3<-rstandard(Task.best3)

fit.Task<-fitted(m.Task)
fit.best1<-fitted(Task.best1)
fit.best3<-fitted(Task.best3)

par(mfrow=c(2,2))
plot(TaskDistract$Distract,res.Task,
	xlab="Distraction Level",ylab="Residuals",ylim=range(-25,25),pch=20)
plot(TaskDistract$Distract,stres.Task,
	xlab="Distraction Level",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(fit.Task,stres.Task,
	xlab="Fitted response",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(TaskDistract$Errors,fit.Task,
	xlab="Observed response",ylab="Fitted response",xlim=range(0,70),ylim=range(0,70),pch=20)
abline(0,1,lty=2)

par(mfrow=c(2,2))
plot(TaskDistract$Distract,res.best1,
	xlab="Distraction Level",ylab="Residuals",ylim=range(-25,25),pch=20)
plot(TaskDistract$Distract,stres.best1,
	xlab="Distraction Level",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(fit.best1,stres.best1,
	xlab="Fitted response",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(TaskDistract$Errors,fit.best1,
	xlab="Observed response",ylab="Fitted response",xlim=range(0,70),ylim=range(0,70),pch=20)
abline(0,1,lty=2)

par(mfrow=c(2,2))
plot(TaskDistract$Distract,res.best3,
	xlab="Distraction Level",ylab="Residuals",ylim=range(-25,25),pch=20)
plot(TaskDistract$Distract,stres.best3,
	xlab="Distraction Level",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(fit.best3,stres.best1,
	xlab="Fitted response",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(TaskDistract$Errors,fit.best3,
	xlab="Observed response",ylab="Fitted response",xlim=range(0,70),ylim=range(0,70),pch=20)
abline(0,1,lty=2)

###################################################################

summary(m.Diabetes)

par(mfrow=c(2,2))
plot(log(Diabetes$glufast),residuals(m.Diabetes),
	xlab="log(glufast)",ylab="Residuals",ylim=range(-0.5,0.5),pch=20)
plot(log(Diabetes$glufast),rstandard(m.Diabetes),
	xlab="log(glufast)",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(fitted(m.Diabetes),rstandard(m.Diabetes),
	xlab="Fitted response",ylab="Standardized Residuals",ylim=range(-3,3),pch=20)
plot(log(Diabetes$glutest),fitted(m.Diabetes),
	xlab="Observed response",ylab="Fitted response",xlim=range(5.5,7.3),ylim=range(5.5,7.3),pch=20)
abline(0,1,lty=2)
par(mfrow=c(1,1))

###################################################################
#Simulated example demonstrating the influence of outliers
set.seed(43723)
n<-25
X<-sort(runif(n,0,20))
b0<-2.0
b1<-3.0
Y<-b0+b1*X+rnorm(n)*3
cbind(X,Y)

par(mfrow=c(2,2))
plot(X,Y,pch=20)
abline(coef(lm(Y~X)))
text(0,60,substitute(hat(beta)[0] == b, list(b = round(coef(lm(Y~X))[1],3))),adj=0)
text(0,55,substitute(hat(beta)[1] == b, list(b = round(coef(lm(Y~X))[2],3))),adj=0)
title("Original data")

newY<-Y
newY[20]<-Y[20]-30
plot(X,newY,pch=20)
abline(coef(lm(newY~X)))
abline(coef(lm(Y~X)),lty=2)
points(X[20],newY[20],cex=2,col="red")

text(0,60,substitute(hat(beta)[0] == b, list(b = round(coef(lm(newY~X))[1],3))),adj=0)
text(0,55,substitute(hat(beta)[1] == b, list(b = round(coef(lm(newY~X))[2],3))),adj=0)
title("Data with outlier")

plot(X,rstandard(lm(Y~X)),ylim=range(-5,5),pch=20,ylab="Standardized residuals")
abline(h=c(-2,0,2),lty=c(2,3,2))

plot(X,rstandard(lm(newY~X)),ylim=range(-5,5),pch=20,ylab="Standardized residuals")
abline(h=c(-2,0,2),lty=c(2,3,2))

par(mfrow=c(1,1))

#######################################################################

m0.visc<-lm(viscosity~pcv,data=Viscosity)
m0.hooker<-lm(temp~pressure,data=Hooker)
m1.hooker<-lm(temp~pressure+I(pressure^2),data=Hooker)

res.visc<-residuals(m0.visc)
res.hooker0<-residuals(m0.hooker)
res.hooker1<-residuals(m1.hooker)

pdf('2013-PPplots.pdf',paper="USr",height=8.5,width=11)
par(mfrow=c(1,3))
#P-P plots
library(e1071)
probplot(res.visc);title('P-P plot for Viscosity residuals')
probplot(res.hooker0);title('P-P plot for Linear fit to Hooker data')
probplot(res.hooker1);title('P-P plot for Quadratic fit to Hooker data')
dev.off()

pdf('2013-QQplots.pdf',paper="USr",height=8.5,width=11)
par(mfrow=c(1,3))
#Q-Q plots
qqnorm(res.visc,main='Q-Q plot for Viscosity residuals')
qqline(res.visc)
qqnorm(res.hooker0,main='Q-Q plot for Linear fit to Hooker data')
qqline(res.hooker0)
qqnorm(res.hooker1,main='Q-Q plot for Quadratic fit to Hooker data')
qqline(res.hooker1)
dev.off()

par(mfrow=c(1,1))