# R commands for practical exercise 10

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


# a)
# We compute the number of suicides and the person-years for males:
suicides.aggr=aggregate(suicides,list(age=suicides$agegr,gender=suicides$sex),sum)
males.suicides=suicides.aggr[suicides.aggr$gender==1,]$suicides
males.pyears=suicides.aggr[suicides.aggr$gender==1,]$pyears

# Then we compute the occurrence/exposure rates with standard errors:
males.occexp=males.suicides/males.pyears
males.se=males.occexp/sqrt(males.suicides)

# We plot the suicide rates per 10 000 person-years with standard confidence limits:
plot(seq(22.5,62.5,5), 10000*males.occexp,type='l', ylim=c(0,6.5), 
xlab="Age",ylab="Suicide rates per 10000",main="Males")
lines(seq(22.5,62.5,5), 10000*males.occexp-1.96*males.se*10000,lty=2)
lines(seq(22.5,62.5,5), 10000*males.occexp+1.96*males.se*10000,lty=2)


# b)
# We add the log-transformed confidence limits in the plot from a:
lines(seq(22.5,62.5,5), 10000*males.occexp*exp(-1.96/sqrt(males.suicides)),lty=2,col='red')
lines(seq(22.5,62.5,5), 10000*males.occexp*exp(1.96/sqrt(males.suicides)),lty=2,col='red')


# c) 
# We repeat the commands from a and b for females
females.suicides=suicides.aggr[suicides.aggr$gender==2,]$suicides
females.pyears=suicides.aggr[suicides.aggr$gender==2,]$pyears
females.occexp=females.suicides/females.pyears
females.se=females.occexp/sqrt(females.suicides)

plot(seq(22.5,62.5,5), 10000*females.occexp,type='l', ylim=c(0,6.5), xlab="Age",ylab="Suicide rates per 10000",main="Females")
lines(seq(22.5,62.5,5), 10000*females.occexp-1.96*females.se*10000,lty=2)
lines(seq(22.5,62.5,5), 10000*females.occexp+1.96*females.se*10000,lty=2)


lines(seq(22.5,62.5,5), 10000*females.occexp*exp(-1.96/sqrt(females.suicides)),lty=2,col='red')
lines(seq(22.5,62.5,5), 10000*females.occexp*exp(1.96/sqrt(females.suicides)),lty=2,col='red')


# d)
# Rates for males and females in the same plot
plot(seq(22.5,62.5,5), 10000*males.occexp,type='l', ylim=c(0,6.5), xlab="Age",ylab="Suicide rates per 10000")
lines(seq(22.5,62.5,5), 10000*females.occexp,lty=2)
legend("topright",lty=1:2,legend=c("Males","Females"))