############### Exercise sheet 4

#### Exercise 1 - Compute and plot SMRs of oral cavity cancer 
## read in the dataset
oral<-read.table("oral.txt",header=TRUE)
head(oral)

## load function map.germany, cutpoints: (0.5,0.6,0.8,1.0,1.2,1.5,2.0)
source("germany.txt")
source("map.germany.txt")

## caluclate SMR
smroral<-oral$observed/oral$expected
summary(smroral)

## map SMRs
map.germany(smroral)

pdf(file="smroral.pdf",width=7,height=9)
map.germany(smroral)
dev.off()

postscript(file="smroral.ps",width=7,height=9,horizontal=FALSE)
map.germany(smroral)
dev.off()

## plot log(SMRs) vs. log(number of expected cases)
plot(log(oral$expected),log(smroral),xlab="log(expected number of oral cavity cancer cases per region)",ylab="log(SMR of oral cavity cancer cancer per region)",pch="+",col="red",main="Variability of SMRs for oral cancer")




##### Exercise 2 - Test for overdispersion when fitting a Poisson GLM
## fit a Poisson GLM including only the offset
oralglm <- glm(oral$observed ~ offset(log(oral$expected)), family="poisson")  

## extract deviance of the model
dev <- oralglm$deviance
## extract degrees of freedom of the model
df <- oralglm$df.residual

## compute measure for overdispersion
overdispersion <- dev/df
overdispersion

## chi^2 test for overdispersion
pvalue <- pchisq(dev, df,lower.tail=FALSE)
pvalue





##### Exercise 3 - Relative risk estimates by Clayton & Kaldor
## load library DCluster
library(DCluster)

## run the Clayton & Kaldor method
ckoralall<-empbaysmooth(oral$observed,oral$expected)
names(ckoralall)

## extract the Clayton & Kaldor estimates for relative risk
ckoral<-ckoralall$smthrr

## map the Clayton & Kaldor estimates for relative risk
map.germany(ckoral)

## compare SMRs and Clayton & Kaldor estimates by mapping both
par(mfrow=c(1,2))
map.germany(smroral)
map.germany(ckoral)

## compare SMRs and Clayton & Kaldor estimates by a boxplot
par(mfrow=c(1,2))
boxplot(smroral,main="SMRs",ylim=c(0,2))
boxplot(ckoral,main="CK",ylim=c(0,2))

## extract the Empirical Bayes estimates for alpha and nu
alphahat<-ckoralall$alpha
nuhat<-ckoralall$nu

## update the posterior Gamma distribution
alphapost<-oral$observed+alphahat
nupost<-oral$expected+nuhat

## compute posterior probabilities
pp<-pgamma(1,shape=alphapost,rate=nupost,lower.tail=FALSE)
summary(round(pp,digits=2))

## map posterior probabilities
map.germany(pp,cutpoints=c(0.3,0.4,0.5,0.6,0.7,0.8,0.9))

## in which regions is the posterior probability larger than 70 %?
for(i in 1:544){
if (pp[i]>0.7){
pp[i]<-1}
else {pp[i]<-0}}
map.germany(pp,cutpoints=c(0.3,0.4,0.5,0.6,0.7,0.8,0.9))