####Analysis of Coronary Heart Disease Data
##############################################


#Age              Smokers Non-smokers
#group        Deaths Person-years Deaths Person-years
#35 ? 44            32 52407 2 18790
#45 ? 54          104 43248 12 10673
#55 ? 64          206 28612 28 5710
#65 ? 74          186 12663 28 2585
#75 ? 84            102 5317 31 1462




deaths <- c(32,2,104,12,206,28,186,28,102,31)

population <- c(52407, 18790, 43248, 10673, 28612, 5710, 12663, 2585, 5317, 1462)

smoke <-gl(2,1,10, labels=c("Yes", "No"))

age <- gl(5,2,10,labels=c("35--44","45--54","55--64", "65--74", "75--84"))

chddata=data.frame(deaths,population,smoke,age)

chddata


rate= deaths*100000/population

plot(age[smoke=="No"], rate[smoke=="No"], xlab="Age",
ylab="Deaths per 100.000 of population",lty=0,ylab=c(0,2500))
points(age[smoke=="Yes"], rate[smoke=="Yes"], pch="A")
points(age[smoke=="No"], rate[smoke=="No"], pch="B")
legend("topleft",c("smokers","non-smokers"),pch=c("A","B"))


age <- as.numeric(age)

age

smoke <- ifelse(smoke=="Yes",1,0)


smoke

agesq <- age^{2}

agesq

agesm <-ifelse(smoke==0, age, 0)

agesm

populationl <- log(population)




fit1 <- glm(deaths~offset(populationl)+smoke+age+agesq+agesm, family=poisson)

rate.ratio <- exp(fit1$coef[-1])

rate.ratio

summary(fit1)


res.pear <- residuals(fit1, type="pearson")

res.dev  <- residuals(fit1, type="deviance")

predict.fit <- predict(fit1, type="response")

cbind(age, smoke, deaths, predict.fit, res.pear, res.dev)

devian.fit <- sum(res.dev^{2})

1-pchisq(devian.fit, df=10-5)

pear.fit <- sum(res.pear^{2})

1-pchisq(pear.fit, df=10-5)

par(mfrow=c(2,2))
plot(res.pear,main="Pearson Residuals")
abline(h=0)
hist(res.pear,main="Histogram of Pearson Residuals")
plot(res.dev,main="Deviance Residuals")
abline(h=0)
hist(res.dev,main="Histogram of Deviance Residuals")