# 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)


# a)
# We 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)


# b) 

# We read the small R-function given in the exercise
expcoef=function(glmobj)
{ 
regtab=summary(glmobj)$coef 
expcoef=exp(regtab[,1]) 
lower=expcoef*exp(-1.96*regtab[,2]) 
upper=expcoef*exp(1.96*regtab[,2]) 
cbind(expcoef,lower,upper)
}

# We use the function to obtain the basline rates and the hazard rate ratio for sex (with confidence intervals)L
expcoef(fit.prop)


# c)

# We fit a model with proportional rates for sex and job status category:
fit.c=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex)+factor(jobgr)-1, family=poisson,data=suicides)

# We use a likelihood ratio test to test if job status category is significant
anova(fit.prop,fit.c,test="Chisq")

# Job status category is significant

# We inspect the results of the model with sex and job status category:
summary(fit.c)
expcoef(fit.c)


# d)

# We test for interaction between age group and sex
fit.dsex=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex)+factor(jobgr)+factor(agegr):factor(sex)-1, family=poisson,data=suicides)
anova(fit.c,fit.dsex,test="Chisq")

# We test for interaction between age group and job status category
fit.djob=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex)+factor(jobgr)+factor(agegr):factor(jobgr)-1, family=poisson,data=suicides)
anova(fit.c,fit.djob,test="Chisq")

#The interactions are not significant


# e)

# We test for interaction between sex and job status category
fit.e=glm(suicides~offset(log(pyears))+factor(agegr)+factor(sex)+factor(jobgr)+factor(sex):factor(jobgr)-1, family=poisson,data=suicides)
anova(fit.c,fit.e,test="Chisq")

# The interaction is significant


# f)
# We inspect the results of the model with interaction between sex and job status category:
summary(fit.e)
expcoef(fit.e)

# Note that the main effects for job status groups are for males

# If we define females as the reference (using the command below)
# and rerun the model fit, we the main effects apply to females.
suicides$sex=relevel(factor(suicides$sex),ref=2)