#Commands for computing
occurrence/exposure rates (cf. example
#Poisson regression (cf.
section 5.3.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/h06/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)
#One
should also check for other interactions, but no more are found.
#The "final model"
is sm.fit2:
summary(sm.fit2)