#
# AUTHORS: BRAD CARLIN AND JON WAKEFIELD, APRIL 15th, 2006.
# Function to produce maps for Ohio. If type="e" equal-spaced intervals,
# if type="q" based on quantiles; if the former then lower and upper can
# be specified.
# April 10th, 2006.
#
OhioMap <-function(data, ncol=5, figmain="", digits=5, type="e",
                   lower=NULL, upper=NULL) {
  if (is.null(lower)) lower <- min(data)
  if (is.null(upper)) upper <- max(data) 

  if (type=="q"){p <- seq(0,1,length=ncol+1)
                 br <- round(quantile(data,probs=p),2)}
  if (type=="e"){br <- round(seq(lower,upper,length=ncol+1),2)}
  shading <- gray((ncol-1):0/(ncol-1))
  data.grp <- findInterval(data,vec=br,rightmost.closed=T,all.inside=T)
  data.shad <- shading[data.grp]
  map("county", "ohio", fill=TRUE, col=data.shad)
  leg.txt<-paste("[",br[ncol],",",br[ncol+1],"]",sep="")
  for(i in (ncol-1):1){
    leg.txt<-append(leg.txt,paste("[",br[i],",",br[i+1],")",sep=""),)
  }
  leg.txt<-rev(leg.txt)
  legend(-81.4,39.4,legend=leg.txt,fill=shading,bty="n",ncol=1,cex=.8)
  title(main=figmain,cex=1.5)
  invisible()
}
#
# QUESTION 1
#
# Calculate the reference probabilities -- we have 4 age bands (we only
# look at data for 45+) and 2 gender and 2 race categories, hence 
# 2x2x4 = 16 stratum. Yagg and Nagg are the sums (over areas) of the 
# cases and populations, picked up by summing over the relevant indices.
#
library(maps)
library(maptools)
source("ohio.dat")
names(ohio)
head(ohio)
#
# The next set of code
nstrata <- 16; ncounty <- 88; q <- Yagg <- Nagg <- rep(0,nstrata);
for (j in 1:nstrata){
  indj <- j+seq(0,ncounty-1)*16
  Yagg[j] <- sum(ohio$deaths[indj]); Nagg[j] <- sum(ohio$popn[indj])
  q[j] <- Yagg[j]/Nagg[j]
}
#
# Now the expected numbers
#
Obs <- Exp <- rep(0,ncounty)
for (i in 1:88){
  indi <- (i-1)*16+seq(1,16)
  Obs[i] <- sum(ohio$deaths[indi])
  Exp[i] <- sum(ohio$popn[indi]*q)
}
SMR <- Obs/Exp
OhioMap(SMR,ncol=8,type="e",figmain="Ohio lung cancer SMRs")
ses <- sqrt(SMR/Exp)
OhioMap(ses,ncol=8,type="e",figmain="Ohio lung cancer ses of MLE")
#
#
# QUESTION 2
#
# Empirical Bayes Approach using eBayes function
#
library(SpatialEpi)
EBest <- eBayes(Obs, Exp)
EB.a = EBest$alpha #51.73209, dispersion in Gamma distribution of R.E.
EB.alpha = EBest$beta  #-0.0253608, alpha
EBest$RR  # shows RR
OhioMap(EBest$RR,ncol=8,type="e",figmain="Ohio lung cancer smoothed estimates",
        lower=0,upper=max(SMR)) # on the same scale of SMR
#
# compare SMR and eBayes RR
#
plot(SMR,EBest$RR,xlab="SMRs",ylab="EB estimates",
     xlim=c(min(SMR,EBest$RR),max(SMR,EBest$RR)),
     ylim=c(min(SMR,EBest$RR),max(SMR,EBest$RR)))
abline(0,1)
#
for (i in 1:4){
	EBpostdens(Obs[i], Exp[i], EBest$alpha, EBest$beta, lower=0, upper=1.6, main=paste("Area", i))
}
#
thresh12 <- EBpostthresh(Obs,Exp,alpha=EBest$alpha,beta=EBest$beta,rrthresh=1.2)
OhioMap(thresh12,ncol=8,type="e",figmain="Post probs of exceedence of 1.2",
        lower=0,upper=1)
#
# QUESTION 3
#
source("http://www.math.ntnu.no/inla/givemeINLA.R")
library("INLA")
#
ohio.data <-  data.frame(Counts=Obs, E=Exp, Region=1:88)
inla.mod <- inla(Counts ~ 1 + f(Region, model="iid", param = c(1, 0.026)),
data=ohio.data, family="Poisson", E=E )
summary(inla.mod)
#
inla.alpha <- inla.mod$summary.fixed[4]
inla.RE <- exp(inla.mod$summary.random$Region[,5])  # exp(V_i) = delta_i
inla.RR <- exp(inla.alpha)*inla.RE
#
# map of the RR using Poisson-Lognormal model 
#
OhioMap(inla.RR,ncol=8,type="e",lower=0,upper=max(SMR),
figmain="Ohio lung cancer smoothed estimates: poisson-Lognormal Model")
#
# compare RR using Poisson-Gamma model and Poisson-Lognormal model
#
mymin <- min(c(EBest$RR, inla.RR))
mymax <- max(c(EBest$RR, inla.RR))
plot(EBest$RR, inla.RE, xlim=c(mymin, mymax), ylim=c(mymin, mymax),
xlab="Gamma RR", ylab="Lognormal RR")
abline(0,1)
#
# QUESTION 4
#
ohioAgg = data.frame(Counts=Obs, E=Exp, Region=1:88)
ohioAgg$Region2 = ohioAgg$Region
ICAR.mod = inla(Counts ~ 1 + f(Region, model="iid", param = c(1, 0.026)) + f(Region2, model="besag",  
graph="ohio.graph", param=c(1, 0.026)), data=ohioAgg, family="Poisson", E=E )
#
summary(ICAR.mod)
#
# Extract and plot non-spatial and spatial random effects 
REnonspatial = ICAR.mod$summary.random$Region[,5]
REspatial = ICAR.mod$summary.random$Region2[,5]
#
OhioMap(REnonspatial,ncol=8,type="e",figmain="Non-spatial Random Effect")
#
OhioMap(REspatial,ncol=8,type="e",figmain="Spatial Random Effect")
#
alpha = ICAR.mod$summary.fixed[4]
fittedRR = exp(alpha + REnonspatial + REspatial)
OhioMap(fittedRR,ncol=8,type="e",figmain="Relative Risk estimated by ICAR model")