# Commands for the suicide data

 

# Commands for computing occurrence/exposure rates (cf. example 5.4 in ABG) and

# Poisson regression (cf. section 5.2.1 in ABG) for suicides among Danish non-manual workers

 

# First we consider occurrence/exposure rates separately for each gender disregard job status

 

# Read total number of suicides in 5-year age groups:

msuicide<-c(98,150,193,175,195,206,152,119,54)

fsuicide<-c(102,107,102,112,108,124,90,69,29)

 

#Read total number of person-years in 5-year age groups

myears<-c(459933,652446,535400,448444,420206,405289,336446,254087,188587)

fyears<-c(891951,678109,424749,351765,334995,326542,264310,199731,125847)

 

#Compute occurrence/exposure rates

mrates<-msuicide/myears

frates<-fsuicide/fyears

 

#Plot occurrence/exposure rates per 10 000

plot(seq(22.5,62.5,5),mrates*10000,type="l",xlim=c(20,70),ylim=c(0,6),xlab="Age",ylab="Suicude rates per 10000")

lines(seq(22.5,62.5,5),frates*10000,lty=2)

 

 

#Then we consider Poisson regression

 

# Read full data set into R:

selvmord<- read.table("http://www.uio.no/studier/emner/matnat/math/STK4080/h08/computing/danskselvmord.txt", header=T)

 

#Define agegroup, sex and jobgroup as factors

selvmord$agegr<-factor(selvmord$agegr)

selvmord$sex<-factor(selvmord$sex)

selvmord$jobgr<-factor(selvmord$jobgr)

 

#First we show how we may reproduce occurrence/exposure rates:

sm.mfit<-glm(suicides~offset(log(pyears))+agegr-1,family=poisson,data=selvmord,subset=(sex==1))

sm.ffit<-glm(suicides~offset(log(pyears))+agegr-1,family=poisson,data=selvmord,subset=(sex==2))

cbind(mrates,exp(sm.mfit$coef))

cbind(frates,exp(sm.ffit$coef))

plot(seq(22.5,62.5,5),exp(sm.mfit$coef)*10000,type="l",xlim=c(20,70),ylim=c(0,6),xlab="Age",ylab="Suicude rates per 10000")

lines(seq(22.5,62.5,5), exp(sm.ffit$coef)*10000,lty=2)

 

#Then we fit a Poisson regression model starting with a model with only main effects:

sm.fit1<-glm(suicides~offset(log(pyears))+agegr+sex+jobgr-1, family=poisson, data=selvmord)

summary(sm.fit1)

 

#Then we fit a model with interaction between sex and job group

sm.fit2<-glm(suicides~offset(log(pyears))+agegr+sex+jobgr+sex:jobgr-1, family=poisson, data=selvmord)

anova(sm.fit1, sm.fit2, test="Chisq")

 

#One should also check for other interactions, but no more are found.

 

#The "final model" is sm.fit2:

summary(sm.fit2)