############### Exercise sheet 5

#### Exercise 1 - Estimate BYM model using different hyperprioris
# load library(INLA)
library(INLA)

# source map of germany and function to plot map
source("map.germany.txt")
source("germany.txt")

# read data of oral cancer
oral<-read.table("oral.txt",header=TRUE)

# create dataframe
Y<-oral$observed
E<-oral$expected
region<-seq(1,544,1)
region.struct<-seq(1,544,1)
oraldata<-data.frame(Y,E,region,region.struct)


# estimate bym model with hyperprior 1
hypera<-Y~f(region.struct,model="besag",graph.file="germany.graph.txt",param=c(1,0.005))+f(region,model="iid",param=c(1,0.005))

mhypera<-inla(hypera,family="poisson",data=oraldata,E=E,control.predictor=list(compute=TRUE))


rrhypera<-exp(mhypera$summary.linear.predictor$mean)

# estimate bym model with hyperprior 2
hyperc<-Y~f(region.struct,model="besag",graph.file="germany.graph.txt",param=c(0.25,0.0005))+f(region,model="iid",param=c(0.25,0.00025))
mhyperc<-inla(hyperc,family="poisson",data=oraldata,E=E,control.predictor=list(compute=TRUE))
rrhyperc<-exp(mhyperc$summary.linear.predictor$mean)

# compare results when using different hyperpriors
summary(rrhypera)
summary(rrhyperc)

par(mfrow=c(1,2))
map.germany(rrhypera)
map.germany(rrhyperc)

# compare BYM results with SMRs
smroral<-oral$observed/oral$expected

par(mfrow=c(1,2))
map.germany(smroral)
map.germany(rrhyperc)


###### Exercise 2 - Compute Bayesian p-values for BYM estimates of relative risk
# estimate BYM model including option cdf=c(0)
mhypera_2<-inla(hypera,family="poisson",data=oraldata,E=E,control.predictor=list(compute=TRUE,cdf=c(0)))

# calculate Bayesian p-values from BYM results
# look at:
head(mhypera_2$summary.linear.predictor)

pps<-1-mhypera_2$summary.linear.predictor[,6]

# define cutpoints for plot
cutpoints<-c(0.3,0.4,0.5,0.6,0.7,0.8,0.9)

# plot map of Bayesian p-values
map.germany(pps,cutpoints)

par(mfrow=c(1,2))
map.germany(rrhypera)
map.germany(pps,cutpoints)


###### Exercise 3 - Ecological regression
scotland<-read.table("scotland.txt",header=TRUE)
head(scotland)

Y<-scotland$observed
E<-scotland$expected
X<-scotland$cov
region.scot<-seq(1,56,1)
region.struct.scot<-seq(1,56,1)
scotdata<-data.frame(Y,E,X,region.scot,region.struct.scot)

fscot<-Y~f(region.struct.scot,model="besag",graph.file="scotland.graph.txt",param=c(1,0.005))+f(region.scot,model="iid",param=c(1,0.005))+X

mscot<-inla(fscot,family="poisson",data=scotdata,E=E,control.predictor=list(compute=TRUE))

mscot$summary.fixed[,c(1,3,4,5)]

######## Exercise 4 - Analyse simulated data

### Simulation 1: Relative risk is 1
# generate data
Ysim1 <- rpois(544,oral$expected)

E<-oral$expected
region<-seq(1,544,1)
region.struct<-seq(1,544,1)

# estimate BYM model
datasim1<-data.frame(Ysim1,E,region,region.struct)
sim1<-Ysim1~f(region.struct,model="besag",graph.file="germany.graph.txt",param=c(0.25,0.0005))+f(region,model="iid",param=c(0.25,0.00025))
msim1<-inla(sim1,family="poisson",data=datasim1,E=E,control.predictor=list(compute=TRUE))
ssim1<-exp(msim1$summary.linear.predictor$mean)

# result
map.germany(ssim1)

### Simulation 2: Relative risk is 1 for region 1:200 and 2 for region 201:544
# generate data
Ysim2 <- rpois(544, c(rep(1,200),rep(2,344))*oral$expected)

# estimate BYM model
datasim2<-data.frame(Ysim2,E,region,region.struct)
sim2<-Ysim2~f(region.struct,model="besag",graph.file="germany.graph.txt",param=c(0.25,0.0005))+f(region,model="iid",param=c(0.25,0.00025))
msim2<-inla(sim2,family="poisson",data=datasim2,E=E,control.predictor=list(compute=TRUE))
ssim2<-exp(msim2$summary.linear.predictor$mean)

# result
map.germany(ssim2)