chd.data <- importData("c:\\Temp\\CHD.xls")

#Number of data points
n<-nrow(chd.data) 

#Breakdown by Category

chd.response<-table(chd.data$CHD,chd.data$AGRP)

#Column totals
chd.totals<-apply(chd.response,2,sum)

#Group Proportions
chd.grpprops<-chd.response[2,]/chd.totals
chd.ese<-sqrt(chd.response[1,]*chd.response[2,]/chd.totals^3)

#Fit the "null" GLM, where the linear predictor is constant
fit.0<-glm(CHD ~ 1,family=binomial,data=chd.data)
summary(fit.0)


#Fit the model with AGE as the only predictor
fit.1<-glm(CHD ~ AGE,family=binomial,data=chd.data)
summary(fit.1)

beta.10<-fit.1$coeff[1]
beta.11<-fit.1$coeff[2]

#Plot the inverse-linked linear predictor
x.age<-c(0:1000)/10
linear.predictor1<-beta.10+beta.11*x.age
y.predicted1<-exp(linear.predictor1)/(1+exp(linear.predictor1))

#plot(x.age,y.predicted1,type="l",xlab="AGE",ylab="Prob. of CHD",ylim=range(-0.15,1.15))
plot(x.age,y.predicted1,type="l",xlab="AGE",ylab="Prob. of CHD")


#Try alternative links
#Probit

fit.2<-glm(CHD ~ AGE,family=binomial(link = probit),data=chd.data,maxit=200)
summary(fit.2)
beta.20<-fit.2$coeff[1]
beta.21<-fit.2$coeff[2]
linear.predictor2<-beta.20+beta.21*x.age
y.predicted2<-pnorm(linear.predictor2)
lines(x.age,y.predicted2,lty=2)

#Complementary log-log
fit.3<-glm(CHD ~ AGE,family=binomial(link = cloglog),data=chd.data,maxit=200)
summary(fit.3)
beta.30<-fit.3$coeff[1]
beta.31<-fit.3$coeff[2]
linear.predictor3<-beta.30+beta.31*x.age
y.predicted3<-1-exp(-exp(linear.predictor3))
lines(x.age,y.predicted3,lty=3)

legend(0,1.0,c("logit","probit","complementary log-log"),lty=c(1,2,3))

title("Fit to CHD data for three link functions")

chd.grpmids<-c(25,32.5,37.5,42.5,47.5,52.5,57.5,65.0)
par(cex=0.75)
text(chd.grpmids,chd.grpmids*0,round(chd.grpprops,2))
par(cex=1)

points(chd.grpmids,chd.grpprops,mark=1)