#R-code STAT2008 assignment 2
#u6405056 and u6380618

setwd("/Users/amyyu/Desktop/stat2008 rstudio")
Visits <- read.csv("Visits.csv", header = T)
Visits
attach(Visits)
summary(Visits)

#Question A
res.y <- (Visits$nondocco + Visits$doctorco)
res.y
summary(res.y)

hist(res.y)
boxplot(res.y)

#sex
plot(res.y ~ factor(Visits$sex), col="cyan2")

#age
hist(Visits$age)
boxplot(Visits$age, xlab="Visits$age", ylab="Age")
plot(res.y ~ Visits$age)
plot(res.y ~ factor(Visits$age), col="cyan2")

#income
hist(Visits$income)
boxplot(Visits$income, xlab="Visits$income", ylab="Income")
plot(res.y ~ Visits$income)

#insurance
#estimates for beta
plot(res.y ~ factor(Visits$insurance), col="cyan2")

#illness
hist(Visits$illness)
boxplot(Visits$illness, xlab="Visits$illness", ylab="Illness")
plot(res.y ~ Visits$illness)

#actdays
hist(Visits$actdays)
boxplot(Visits$actdays, xlab="Visits$actdays", ylab="Actdays")
plot(res.y ~ Visits$actdays)

#hscore
hist(Visits$hscore)
boxplot(Visits$hscore, xlab="Visits$hscore", ylab="Hscore")
plot(res.y ~ Visits$hscore)

#chcond
plot(res.y ~ factor(Visits$chcond), col="cyan2")

#hospadmi
hist(Visits$hospadmi)
boxplot(Visits$hospadmi, xlab="Visits$hospadmi", ylab="Hospadmi")
plot(res.y ~ Visits$hospadmi)

#hospdays
hist(Visits$hospdays)
boxplot(Visits$hospdays, xlab="Visits$hospdays", ylab="Hospdays")
plot(res.y ~ Visits$hospdays)

#prescrib
hist(Visits$prescrib)
boxplot(Visits$prescrib, xlab="Visits$prescrib", ylab="Prescrib")
plot(res.y ~ Visits$prescrib)

#nonpresc
hist(Visits$nonpresc)
boxplot(Visits$nonpresc, xlab="Visits$nonpresc", ylab="Nonpresc")
plot(res.y ~ Visits$nonpresc)


#Question B
Visits.lm <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lm, which = 1)
summary(Visits.lm)
anova(Visits.lm)


#Question C
#log(y+1)
Visits.lmt1 <- lm(log(res.y+1) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt1, which = 1, sub="Transformation: log(y+1)")
summary(Visits.lmt1)
anova(Visits.lmt1)
#sqrt(y)
Visits.lmt2 <- lm(sqrt(res.y) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt2, which = 1, sub="Transformation: sqrt(y)")
summary(Visits.lmt2)
anova(Visits.lmt2)
#y^(1/4)
Visits.lmt3 <- lm((res.y)^(1/4) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt3, which = 1, sub="Transformation: y^(1/4)")
summary(Visits.lmt3)
anova(Visits.lmt3)
#sqrt(y+1)
Visits.lmt4 <- lm(sqrt(res.y+1) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt4, which = 1, sub="Transformation: sqrt(y+1)")
summary(Visits.lmt4)
anova(Visits.lmt4)
#log(y+2)
Visits.lmt5 <- lm(log(res.y+2) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt5, which = 1, sub="Transformation: log(y+2)")
summary(Visits.lmt5)
anova(Visits.lmt5)


#Question D
require(MASS)
Visits.lmt6 <- lm((res.y+1) ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
plot(Visits.lmt6, which = 1)
boxcox(Visits.lmt6, plotit=T, lambda=seq(-6.5, 0, by=0.1))
boxcox(Visits.lmt6, plotit=T, lambda=seq(-4, -2.5, by=0.1))

y.bc <- (res.y+1)^(-3.2)
summary(y.bc)
Visits.lmboxcox <- lm(y.bc ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
summary(Visits.lmboxcox)
anova(Visits.lmboxcox)
plot(Visits.lmboxcox, which = 1)
lines(lowess(Visits.lmboxcox$fit, Visits.lmboxcox$res), lwd=2, col="red")


#Question E
d <- residuals(lm(res.y ~ factor(sex)+age+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc))
g <- residuals(lm(income ~ factor(sex)+age+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc))
plot(g, d,xlab="income residuals",
     ylab="res.y residuals",
     main="Added Variable Plot: Income",
     pch=16, cex.lab=1.5)
lines(lowess(g, d), lwd=2, col="red")
abline(lm(d~g), lwd=2, col="cyan4")

d1 <- residuals(lm(res.y ~ factor(sex)+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc))
g1 <- residuals(lm(age ~ factor(sex)+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc))
plot(g1, d1,xlab="age residuals",
     ylab="res.y residuals",
     main="Added Variable Plot: Age",
     pch=16, cex.lab=1.5)
lines(lowess(g1, d1), lwd=2, col="red")
abline(lm(d1~g1), lwd=2, col="cyan4")


#Question F
choose(4,2)
#6
#m=6, 0.05/6 = 0.008 bonferroni
insurance <- relevel(factor(insurance), ref="freepor")
mod.1 <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
confint(mod.1, level=(1-0.008))
int.1 <- confint(mod.1, level=(1-0.008))[5:7,]
int.1

insurance <- relevel(factor(insurance), ref="freerepa")
mod.2 <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
#m=6, 0.05/6 = 0.008
confint(mod.2, level=(1-0.008))
int.2 <- confint(mod.2, level=(1-0.008))[6:7,]
int.2

insurance <- relevel(factor(insurance), ref="levyplus")
mod.3 <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
#m=6, 0.05/6 = 0.008
confint(mod.3, level=(1-0.008))
int.3 <- confint(mod.3, level=(1-0.008))[7,]
int.3

pw.comps1 <- rbind(int.1, int.2, int.3)

row.names(pw.comps1) <- c("freepor to freerpa", "freepor to levyplus", "freepor to medlevy",
                         "freerpa to levyplus", "freerpa to medlevy",
                         "levyplus to medlevy")
pw.comps1


#Question G
choose(3,2)
#3
#m=3, bonferroni correction 0.05/3 = 1/60
chcond <- relevel(factor(chcond), ref="la")
mod.1 <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
confint(mod.1, level=(1-(1/60)))
int.1 <- confint(mod.1, level=(1-1/60))[11:12,]
int.1

chcond <- relevel(factor(chcond), ref="nla")
mod.2 <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
confint(mod.2, level=(1-(1/60)))
int.2 <- confint(mod.2, level=(1-1/60))[12,]
int.2

pw.comps2 <- rbind(int.1, int.2)
row.names(pw.comps2) <- c("la to nla", "la to np",
                         "nla to np")
pw.comps2


#Question H
#hypothesis test 1, beta insurance = age = sex = income = nonpresc = 0
mod.l1 <- lm(res.y ~ illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc+income+age+factor(sex)+factor(insurance))
anova(mod.l1)

#hypothesis test 2, beta insurance = age = sex = 0
mod.l2 <- lm(res.y ~ income+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc+age+factor(sex)+factor(insurance))
anova(mod.l2)

#hypothesis test 3, beta insurance = 0
mod.l3 <- lm(res.y ~ factor(sex)+age+income+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc+factor(insurance))
anova(mod.l3)

Visits.lmfinal <- lm(res.y ~ factor(sex)+age+income+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib+nonpresc)
anova(Visits.lmfinal)


#Question I
Visits.lmallinteractions <- lm(res.y ~ factor(sex)+factor(sex)*age+factor(sex)*income+factor(sex)*illness+factor(sex)*actdays+factor(sex)*hscore+factor(sex)*factor(chcond)+factor(sex)*hospadmi+factor(sex)*hospdays+factor(sex)*prescrib+factor(sex)*nonpresc)
anova(Visits.lmfinal, Visits.lmallinteractions)
anova(Visits.lmallinteractions)

#testing the order by chcond last
Visits.lmorder <- lm(res.y ~ factor(sex)+factor(sex)*age+factor(sex)*income+factor(sex)*illness+factor(sex)*actdays+factor(sex)*hscore+factor(chcond)+factor(sex)*hospadmi+factor(sex)*hospdays+factor(sex)*prescrib+factor(sex)*nonpresc+factor(sex):factor(chcond))
anova(Visits.lmorder)


#Question J
plot(fitted(Visits.lmfinal), rstandard(Visits.lmfinal), xlab="Fitted Values", ylab="Internally Studentized Residuals")
title(main="Internally Studentised Residuals vs Fitted Values", sub="Visits.lmfinal")
lines(lowess(fitted(Visits.lmfinal), rstandard(Visits.lmfinal)), lwd=2, col="red")

plot(Visits.lmfinal, which = 2)
plot(Visits.lmfinal, which = 4)

sort(cooks.distance(Visits.lmfinal), decreasing=TRUE)[1:7]
sort(abs(rstandard(Visits.lmfinal)), decreasing=TRUE)[1:7]
sort(hatvalues(Visits.lmfinal), decreasing=TRUE)[1:7]

dfbetas(Visits.lmfinal)[c(948,965,739,913),]
dffits(Visits.lmfinal)[c(948,965,739,913)]
covratio(Visits.lmfinal)[c(948,965,739,913)]


#Question K
summary(Visits.lmfinal)

Visits.lminsurance <- lm(res.y ~ factor(sex)+age+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib)

#age
plot(age,res.y,type="n", main="Response vs Age")
x <- seq(0,12, by = 0.01)
d0 <- ifelse(sex == "0", 1, 0)
d1 <- ifelse(sex == "1", 1, 0)
points((age[d1==0]), (res.y[d1==0]),
       col="cyan2", pch=16)
points((age[d1==1]), (res.y[d1==1]),
       col="orange", pch=17)
lines(x, predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="freepor", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="red")
pr <- predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="freepor", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="red")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="freerepa", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="purple")
pr <- predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="freerepa", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="purple")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="levyplus", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="green")
pr <- predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="levyplus", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="green")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="medlevy", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="yellow")
pr <- predict(Visits.lminsurance, newdata=data.frame(age=x, insurance="medlevy", sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="yellow")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
legend("topleft", legend=c("freepor", "freerepa", "levyplus", "medlevy"),
       col=c("red", "purple", "green", "yellow"), lty=1, lwd=2, cex = 0.75)

#hscore
plot(hscore,res.y,type ="n", main="Response vs Hscore")
x <- seq(0,12, by = 0.01)
d0 <- ifelse(sex == "0", 1, 0)
d1 <- ifelse(sex == "1", 1, 0)
points((hscore[d1==0]), (res.y[d1==0]),
       col="cyan2", pch=16)
points((hscore[d1==1]), (res.y[d1==1]),
       col="orange", pch=17)
lines(x, predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="freepor", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="red")
pr <- predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="freepor", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="red")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="freerepa", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="purple")
pr <- predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="freerepa", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="purple")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="levyplus", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="green")
pr <- predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="levyplus", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="green")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="medlevy", age=mean(age), sex="1",illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="yellow")
pr <- predict(Visits.lminsurance, newdata=data.frame(hscore=x, insurance="medlevy", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="yellow")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
legend("topleft", legend=c("freepor", "freerepa", "levyplus", "medlevy"),
       col=c("red", "purple", "green", "yellow"), lty=1, lwd=2, cex = 0.75)

#income
Visits.lmincome <- lm(res.y ~ factor(sex)+age+income+factor(insurance)+illness+actdays+hscore+factor(chcond)+hospadmi+hospdays+prescrib)

plot(income,res.y,type="n", main="Response vs Income")
x <- seq(0,12, by = 0.01)
d0 <- ifelse(sex == "0", 1, 0)
d1 <- ifelse(sex == "1", 1, 0)
points((income[d1==0]), (res.y[d1==0]),
       col="cyan2", pch=16)
points((income[d1==1]), (res.y[d1==1]),
       col="orange", pch=17)
lines(x, predict(Visits.lmincome, newdata=data.frame(income=x, insurance="freepor", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="red")
pr <- predict(Visits.lmincome, newdata=data.frame(income=x, insurance="freepor", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="red")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lmincome, newdata=data.frame(income=x, insurance="freerepa", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="purple")
pr <- predict(Visits.lmincome, newdata=data.frame(income=x, insurance="freerepa", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="purple")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lmincome, newdata=data.frame(income=x, insurance="levyplus", age=mean(age), sex="1",illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="green")
pr <- predict(Visits.lmincome, newdata=data.frame(income=x, insurance="levyplus", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="green")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
lines(x, predict(Visits.lmincome, newdata=data.frame(income=x, insurance="medlevy", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib))), lwd=2, col="yellow")
pr <- predict(Visits.lmincome, newdata=data.frame(income=x, insurance="medlevy", age=mean(age), sex="1", illness=mean(illness), actdays=mean(actdays), hscore=mean(hscore), chcond="np", hospadmi=mean(hospadmi), hospdays=mean(hospdays), prescrib=mean(prescrib)), interval="confidence")
lines(x, pr[,"fit"], lwd=2, col="yellow")
lines(x, pr[,"lwr"], lty=2)
lines(x, pr[,"upr"], lty=2)
legend("topleft", legend=c("freepor", "freerepa", "levyplus", "medlevy"),
       col=c("red", "purple", "green", "yellow"), lty=1, lwd=2, cex = 0.75)