d = read.table("../data/SAheart.data",header=TRUE,sep=",")
d$y=d$adiposity-1
#Ordinary logistic regression
fit = glm(chd~sbp+tobacco+ldl+famhist+obesity+alcohol+age,
  family=binomial,data=d)
summary(fit)
anova(fit,test="Chisq")
fit2 = glm(chd~sbp+tobacco+ldl+famhist+obesity+age,
  family=binomial,data=d)
anova(fit2,test="Chisq")
AIC(fit,fit2)
summary(fit2)
fit3 = glm(chd~sbp+tobacco+ldl+famhist+age,
  family=binomial,data=d)
anova(fit3,test="Chisq")
AIC(fit,fit2,fit3)
summary(fit3)



#Natural cubic splines
library(splines)
fit.ns = glm(chd~ns(sbp,df=4)+ns(tobacco,df=4)+ns(ldl,df=4)+
             ns(obesity,df=4)+ns(alcohol,df=4)+ns(age,df=4)+famhist,
  family=binomial,data=d)
summary(fit.ns)
anova(fit.ns,test="Chisq")
termplot(fit.ns,se=TRUE)
logLik(fit,fit.ns)
AIC(fit,fit.ns)

anova(fit.ns,test="Chisq")

fit.ns2 = glm(chd~ns(sbp,df=4)+ns(tobacco,df=4)+ns(ldl,df=4)+
             ns(obesity,df=4)+ns(age,df=4)+famhist,family=binomial,data=d)
anova(fit.ns2,test="Chisq")
AIC(fit,fit.ns,fit.ns2)
#termplot(fit.ns2,se=TRUE)

fit.ns3 = glm(chd~ns(sbp,df=4)+ns(tobacco,df=4)+ns(ldl,df=4)+
             ns(age,df=4)+famhist,family=binomial,data=d)
anova(fit.ns3,test="Chisq")
AIC(fit3,fit.ns3)
termplot(fit.ns3,se=TRUE)

fit.ns4 = glm(chd~ns(sbp,df=4)+ns(tobacco,df=4)+ldl+
             ns(age,df=4)+famhist,family=binomial,data=d)
termplot(fit.ns4,se=TRUE)
AIC(fit3,fit.ns3,fit.ns4)

fit.ns5 = glm(chd~ns(sbp,df=4)+tobacco+ldl+
                ns(age,df=4)+famhist,family=binomial,data=d)
termplot(fit.ns5,se=TRUE)
AIC(fit3,fit.ns3,fit.ns4,fit.ns5)

library(mgcv)
fit.gam = gam(chd~s(sbp)+s(tobacco)+s(ldl)+famhist+
            s(obesity)+s(alcohol)+s(age),family=binomial,data=d)
plot(fit.gam)

#cubic splines
library(splines)
fit.bs = glm(chd~bs(sbp,df=4)+bs(tobacco,df=4)+bs(ldl,df=4)+
               bs(obesity,df=4)+bs(alcohol,df=4)+bs(age,df=4)+famhist,
             family=binomial,data=d)
summary(fit.bs)
termplot(fit.bs,se=TRUE)
logLik(fit.bs)




AIC(fit,fit.bs)