dac=read.table("exercise_diabetes_study.txt", header=T, sep="\t")

library("ggplot2")

p1=qplot(fbg, data=dac, geom="histogram")
p2=qplot(tg, data=dac, geom="histogram")
p3=qplot(bmi, data=dac, geom="histogram")
p4=qplot(fins, data=dac, geom="histogram",xlab="imputed fins")
p5=qplot(fins, data=d, geom="histogram")
p6=qplot(hba1c, data=dac, geom="histogram")

multiplot(p1, p2, p3, p6, cols=2)

multiplot(p5, p4, cols=2)

shapiro.test(dac$fbg)
shapiro.test(dac$fins)
shapiro.test(dac$tg)
shapiro.test(dac$bmi)
shapiro.test(dac$hba1c)


p11=ggplot(dac, aes(x=fbg, y=hba1c)) + geom_point(shape=1) + geom_smooth(method=lm, se=F)
p12=ggplot(dac, aes(x=tg, y=hba1c)) + geom_point(shape=1) + geom_smooth(method=lm, se=F)
p13=ggplot(dac, aes(x=bmi, y=hba1c)) + geom_point(shape=1) + geom_smooth(method=lm, se=F)
p14=ggplot(dac, aes(x=fins, y=hba1c)) + xlab("fins") + geom_point(shape=1) + geom_smooth(method=lm, se=F)

multiplot(p11, p12, p13, p14, cols=2)

r1=lm(hba1c~fbg, data=dac)
r2=lm(hba1c~tg, data=dac)
r3=lm(hba1c~bmi, data=dac)
r4=lm(hba1c~fins, data=dac)

l_null=lm(hba1c~1, data=dac)
l_full1=lm(hba1c~sex+Age+bmi+fbg+tg+tcho+hdl+ldl+fins, data=dac)
l_f1=step(l_null, scope=list(lower=l_null, upper=l_full1), direction="forward")
library(ggfortify)
autoplot(l_f1, which = 1:4, ncol = 2, label.size = 3)

library(car)
crPlots(l_f1)
box_lam=boxCox(l_f1,plotit=T)

lam=box_lam$x[which(box_lam$y==max(box_lam$y))]

dac_lam=dac
dac_lam$hba1c=(dac$hba1c^lam-1)/lam

l_null=lm(hba1c~1, data=dac_lam)
l_full1=lm(hba1c~sex+Age+bmi+fbg+tg+tcho+hdl+ldl+fins, data=dac_lam)
l_f1_lam=step(l_null, scope=list(lower=l_null, upper=l_full1), direction="forward")
summary(l_f1_lam)
plot(l_f1_lam)
crPlots(l_f1_lam)
autoplot(l_f1_lam, which = 1:4, ncol = 2, label.size = 3)