## Example taken and adapted from
## http://www.stats.gla.ac.uk/data/NERC/NERC10/Mark/Data%20and%20code%20for%20network/bayes.R

# The rates of lip cancer in 56 counties in Scotland have been analysed by Clayton and
# Kaldor (1987) and Breslow and Clayton (1993). The form of the data includes the ob-
# served and expected cases (expected numbers based on the population and its age and
# sex distribution in the county), a covariate measuring the percentage of the population
# engaged in agriculture, fishing, or forestry, and the “position” of each county expressed as
# a list of adjacent counties.


# Load the INLA package
library(INLA)
# Make the Scottish lip cancer data available for use
data(Scotland)
# Get a description of the data
#?Scotland

# Work out the path for the file of neighbourhood information
scot.neigh <- system.file("demodata/scotland.graph", package="INLA")
# Create a file for the formula, containing only a spatial terms
# and a scaled covariate
formula <- Counts~f(Region,model="besag", graph = scot.neigh,
    hyper=list(prec=list(prior="loggamma", param=c(0.5,0.0005)))) + I(X/10)
# Fit a Bayesian model using the inla() function
scot.res <- inla(formula,family="poisson",E=E,data=Scotland,control.predictor=list(compute=TRUE))

# Look at results
summary(scot.res)
# Look at the estimated marginal posterior distributions
plot(scot.res)

# Now try plotting some maps to illustrate the INLA results
# Load the SpatialEpi package; if unavailable try 2nd option below:
library(SpatialEpi)
data(scotland)

scotland.map <- scotland$spatial.polygon
# First map - the raw counts
mapvariable(Scotland$Counts,scotland.map,
    main="Counts, Lip Cancer")
# Second map - the estimated values for the linear predictor
x11()
mapvariable(scot.res$summary.linear.predictor$mean,scotland.map,
    main="Linear Predictor, Lip Cancer")
# Actual predictions, includes population sizes again
x11()
mapvariable(exp(scot.res$summary.linear.predictor$mean)*Scotland$E,scotland.map,
    main="Model Predictions, Lip Cancer")
# Plot of smoothed random effects - showing differences from rest of model
x11()
mapvariable(scot.res$summary.random$Region$mean,scotland.map,
    main="Smoothed Random Effects, Lip Cancer")