# R-commands on Poisson regression for the suicide data

# We read the data

suicides=read.table("http://folk.uio.no/borgan/BGC1-2012/data/suicides.txt", header=T)

# We may compute occurrence/exposure rates for males using the glm-command:
fit.males=glm(suicides~offset(log(pyears))+factor(agegr)-1, family=poisson,data=suicides,subset=(sex==1))
summary(fit.males)

# Note that the parameter estimates from the summary(fit.males) command are the 
# logarithms of the occurrence/exposure rates from practical exercise 10

# Also note that the standard error estimates from the summary(mfit)  command equal 
# the standard errors of the logarithms of the occurrence/exposure rates given 
# as 1/sqrt(males.suicides) from practical exercise 10.


# We then fit a model with proportional rates for males and females
fit.prop=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex)-1, family=poisson,data=suicides)
summary(fit.prop)


# We may test the proportionality assumption by comparing the fit 
# with a model with separate rates for males and females
fit.both=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex) +factor(agegr):factor(sex)-1, family=poisson,data=suicides)
anova(fit.prop,fit.both,test="Chisq")

# The proportionality assumptions seems fine.