\( \newcommand{\bm}[1]{\boldsymbol{#1}} \)

Chapter 8 Implementing an Incidence Model for New Zealand

8.1 Exploring the NZ Shapefile

The following New Zealand shapefile was obtained from http://archive.stats.govt.nz/browse_for_stats/Maps_and_geography/Geographic-areas/digital-boundary-files.aspx#census, New Zealand 2013 (NZTM). This download actually contains several shapefiles, each with different statistical, administrative, or electoral boundaries, and consequently different spatial resolutions, depending on the number of areas. E.g.

  • Area units (replaced by SA2s under NZSAC92) - 1911 areas
  • Wards - 258 areas
  • Community board / local board - 195 areas
  • Territorial authority - 68 areas

Ideally, one would probably want to use SA2s for spatial modelling if there is sufficient data available at the level because it is a finer resolution and statistical boundaries are probably more useful than electoral boundaries. However, for the sake of computational efficiency, the wards shapefile will be used for illustration. It should be straightforward to adapt the following code to other shapefiles.

As in Chapter 7, the shapefile can be downloaded as a zip file, extracted, and then read into R manually, or we can use the following R code. In either case, the file size is 246MB.

# Download, extract, and read shapefile 
URL <- "http://www3.stats.govt.nz/digitalboundaries/census/ESRI_Census_Based_2013_NZTM.zip"
temp <- tempfile()
download.file(URL, temp)
file.rename(temp, paste0(temp, ".zip"))
temp <- paste0(temp, ".zip")
shp.NZ <- paste0("ESRI shapefile Output/2013 Digital Boundaries Generalised Clipped/WARD2013_GV_Clipped",
                 c(".dbf", ".prj", ".shp", ".shx"))
unzip(temp, exdir = tempdir(), files = shp.NZ)
map <- paste0(tempdir(), "/ESRI shapefile Output/2013 Digital Boundaries Generalised Clipped/WARD2013_GV_Clipped.shp") %>%
    readOGR()
unlink(temp)
rm(temp)

Next, we make the following modifications:

  • Reduce the complexity of the polygon detail (file size reduced from 40MB to < 1MB with no notable loss in visual quality). Note that this can take a while for R to run if the original shapefile has high detail.
  • Convert Northings and Eastings to latitude and longitude.
  • Remove distant islands (Chatham Island and Pitt Island). In practice, this could be justified since these islands are far from the mainland (>750Km) and have very few people living there. However, the practical reason for excluding them here is because they lie west of the 0/180 degree line of longitude, which creates difficulties when trying to plot the areas. Note that these two islands compose a single ward, so the modified shapefile has 257 wards.
# Simplify polygons (this may take up to several minutes)
map <- ms_simplify(map, keep = 0.03)

# Convert Northings and Eastings to latitude and longitude
map <- sp::spTransform(map, CRS("+init=epsg:4272"))

# Remove distant island
map.df <- fortify(map)
map.df$id <- as.numeric(map.df$id)
islands <- c(
    unique(map.df[which(map.df$long < 165),]$id),
    unique(map.df[which(map.df$long > 179),]$id)
)
map <- map[-islands,]
N <- length(map)    # Now 526 areas

For convenience, the modified shapefile is also provided as supplementary data (since it is also used in Section 5.2.4). The above steps can be skipped by importing the modified shapefile directly.

# Import the modified shapefile
map <- readOGR("Downloaded Data/Shapefile_NZ_Wards.shp", verbose = FALSE)
N <- length(map)

Now we define the additional maps for the ensuing visualisations, as we did for the Queensland example in Chapter 7.

#-----------------------------------------------------------------------------------------
# Define subset maps and borders
#-----------------------------------------------------------------------------------------

# Create border map from full map
map.border <- unionSpatialPolygons(map, IDs = rep(1, N))
map.border <- fortify(map.border)

# Get map projection for insets
map.proj <- proj4string(map)

# User-defined function for creating inset windows
get.inset <- function(x1, x2, y1, y2, map, map.proj){
    Inset.window <- as(raster::extent(x1, x2, y1, y2), "SpatialPolygons")
    proj4string(Inset.window) <- map.proj
    map.inset <- gIntersection(map, Inset.window, byid = TRUE, drop_lower_td = TRUE, 
        id = sapply(map@polygons, function(x) x@ID))
    return(map.inset)
}

# City inset limits (determined by trial and error)
lims <- data.frame(
    xmin = c(174, 174.5, 172),
    xmax = c(176, 175.5, 173.2),
    ymin = -c(37.5, 41.7, 44),
    ymax = -c(36, 41, 43)
)

# Inset maps for major cities
map.1 <- get.inset(174, 176, -37.5, -36, map, map.proj)         # Auckland
map.2 <- get.inset(174.5, 175.5, -41.7, -41, map, map.proj)     # Wellington
map.3 <- get.inset(172, 173.2, -44, -43, map, map.proj)         # Christchurch

# Borders for insets
map.border.1 <- unionSpatialPolygons(map.1, IDs = rep(1, length(map.1)))
map.border.2 <- unionSpatialPolygons(map.2, IDs = rep(1, length(map.2)))
map.border.3 <- unionSpatialPolygons(map.3, IDs = rep(1, length(map.3)))

# Fortify
map.df <- fortify(map)
map.1 <- fortify(map.1)
map.2 <- fortify(map.2)
map.3 <- fortify(map.3)
map.border.1 <- fortify(map.border.1)
map.border.2 <- fortify(map.border.2)
map.border.3 <- fortify(map.border.3)

# Make shapefile dataframe IDs numeric
map.df$id <- as.numeric(map.df$id) + 1  # Need to add 1 so range is 1:N
map.1$id <- as.numeric(map.1$id) + 1
map.2$id <- as.numeric(map.2$id) + 1
map.3$id <- as.numeric(map.3$id) + 1

# Quick plot
ggplot(data = map.df, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "pink", aes(fill = factor(id))) +
    coord_map() +
    scale_fill_discrete(guide = FALSE)

This is just a plot of the 257 wards coloured by their ID value to check that the map is being plotted as expected.

8.2 Exploring Simulated Incidence Data

Now let's import some simulated incidence data and explore this. We will plot the observed cases and observed SIRs. The observed SIRs are not included in the data set, but we can create them from the observed and expected cases.

# Import data
data <- read.csv("Downloaded Data/8.2 NZ Sim Incidence.csv")
head(data)
##   area_code            area_name observed  expected  covariate
## 1      3102 Onekawa-Tamatea Ward        4 10.354585 -0.4917083
## 2      5803         Cheviot Ward        4  4.262939 -0.1641180
## 3      3903        Feilding Ward        5  8.189436 -0.1674871
## 4      4403         Eastern Ward        5  5.941533 -0.5961955
## 5      1701        Pirongia Ward       10  6.245884  0.1429696
## 6      3005      Kahuranaki Ward       13 18.697521 -0.3818318
# Re-order records to match shapefile IDs
data <- data[match(as.numeric(paste(map$WARD2013)), as.numeric(data$area_code)), ]

# Observed (raw) SIR
SIR.raw <- data$observed/data$expected
SIR.raw[which(is.nan(SIR.raw))] <- 0    # No areas for this dataset

# Truncate raw SIR values
clip <- c(0.25, 4)  # Not necessary in this example since all values are within this range
SIR.raw.clipped <- SIR.raw %>% 
    replace(., which(. > clip[2]), clip[2]) %>%
    replace(., which(. < clip[1]), clip[1])

# Dataframe of the variables we need to append to the shapefile
Append <- data.frame(
    id = 1:N,
    obs = data$observed,
    SIR = SIR.raw.clipped
)

# Fill colours for observed and SIR variables
Fill.colours.obs <- c(grey(seq(1, 0, -0.2)), grey(0))
Fill.colours.SIR <- c("#2C7BB6", "#2C7BB6", "#ABD9E9", "#FFFFBF", "#FDAE61", "#D7191C",  
    "#D7191C")

# Values corresponding to the fill colours
Fill.values.obs <- c(seq(0, 150, 30), max(data$observed)) # Solid black for obs > 150
Fill.values.SIR <- c(log2(0.25), seq(log2(1/1.5), log2(1.5), length = 5), log2(4))
# Note: we use logged values for the SIR so that the colour gradient is linear

# Rescaled fill values
Fill.values.obs.r <- rescale(Fill.values.obs)
Fill.values.SIR.r <- rescale(Fill.values.SIR)

# Merge values with shapefile dataframe
dat <- inner_join(map.df, Append, by = "id")

# Set up basic components of the plot
gg.base <-  ggplot(data = NULL, aes(x = long, y = lat, group = group)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_void() + 
    theme(legend.position = "bottom") +
    guides(fill = guide_colourbar(barwidth = 15)) +
    coord_map() # Keeps x-y ratio consistent

# Add visible layers for borders and fill colours
gg.obs <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = obs), color = NA) +
    scale_fill_gradientn("", 
        colours = Fill.colours.obs, 
        values = Fill.values.obs.r,
        limits = range(Fill.values.obs),
        breaks = seq(0, 150, 50)) +
    ggtitle("Observed Cases")

gg.SIR.raw <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log2(SIR)), color = NA) +
    scale_fill_gradientn("", 
        colours = Fill.colours.SIR, 
        values = Fill.values.SIR.r,
        limits = range(Fill.values.SIR),
        breaks = seq(-2, 2), 
        labels = round(2^seq(-2, 2), 3)) +
    ggtitle("Observed SIR")

# Render the plots
grid.arrange(gg.obs, gg.SIR.raw, nrow = 1)

8.3 Fitting a Bayesian model to the Incidence Data

As with the simulated incidence data for Queensland, it is difficult to determine what spatial patterns exist from the raw data. Now it is time to fit a statistical model in order to estimate the (log) risk surface and the covariate and spatial random effects that compose it.

For this example, we will fit the BYM model (refer to Section 5.2.1). First, we compute the spatial adjacency data using the spdep package. The resulting object is a list which contains the same information as a spatial weights matrix (namely which areas are neighours of each area, and the corresponding weights), just in a different format. For reasons explained below, there is no need to compute the corresponding spatial weights matrix. So rather than use nb2mat() like we did in Sections 4.7 and 7.3 to acquire the spatiel weights matrix, we use nb2WB() instead.

# Get list of adjacency data
adj.data <- nb2WB(poly2nb(map))

There are two areas with no immediately adjacent neighbours. The nearest neighbours will be used instead. These neighbours can be determined by computing the matrix of distances between each pair of areas as demonstrated earlier (see Sections 5.2.4 and 7.1). For the sake of brevity, these two nearest neighbours have been pre-determined; the IDs are hard-coded below. The following code updates this adjacency data list.

# Which areas have no neighbours? 
no.neighbours <- which(adj.data$num == 0)

# Areas closest to these areas (based on nearest distance)
nearest.area <- c(243, 60)

# Recipricate areas with no neighbours and new neighbours for symmetry
nearest.area.r <- c(no.neighbours, nearest.area)
no.neighbours.r <- c(nearest.area, no.neighbours)

# Update adjacent neighbours 
for(i in length(no.neighbours.r):1){
    # Cumulative sum of of neighbours
    if(no.neighbours.r[i] == 1){
        cumsum <- 0
    }else{
        cumsum <- sum(adj.data$num[1:no.neighbours.r[i]])
    }
    # Add neighbour after this index
    if(cumsum == 0){
        adj.data$adj <- c(nearest.area.r[i], adj.data$adj)
    }else if(cumsum == length(adj.data$adj)){
        adj.data$adj <- c(adj.data$adj, nearest.area.r[i])
    }else{
        adj.data$adj <- c(adj.data$adj[1:cumsum], nearest.area.r[i],
            adj.data$adj[(cumsum+1):length(adj.data$adj)])
    }
    # Update number of neighbours
    adj.data$num[no.neighbours.r[i]] <- adj.data$num[no.neighbours.r[i]] + 1
}

# Update weights
adj.data$weights <- c(adj.data$weights, rep(1, length(no.neighbours.r)))

Now we are ready to fit the model. Assuming some weakly informative priors, we will fit the following model:

\[Y_i \sim \mbox{Po}\left(E_i e^{\mu_i}\right)\] where \(Y_i\) and \(E_i\) are the observed and expected counts respectively, and \[ \mu_i = \beta_0 + \beta_1 x_i + S_i + U_i \\ \beta_0 \sim \mathcal{N}(0, 100) \\ \beta_1 \sim \mathcal{N}(0, 100) \\ S_i \mid \bm{S}_{\smallsetminus i} \sim \mathcal{N}\left(\frac{\sum_j{w_{ij}s_j}}{\sum_j{w_{ij}}}, \frac{\sigma_S^2}{\sum_j{w_{ij}}}\right) \\ U_i \sim \mathcal{N}\left(0, \sigma_U^2\right) \\ \sigma_S \sim HN(0, 10) \\ \sigma_U \sim HN(0, 10) \] where the Gaussian and half-Normal distributions are parameterised in terms of mean and variance.

This model can be fit using CARBayes (with the modification that the hyperprior for \(\sigma_S\) and \(\sigma_U\) are inverse-Gamma instead of half-Normal). However, CARBayes only returns the sum \(S_i + U_i\) rather than seperate parameters, which limits the inference that can be made about the spatial random effect. Instead, we will use the software WinBUGS to fit the BYM model, which allows this additional flexibility, at the cost of some computational efficiency. Unlike CARBayes, WinBUGS uses a list of adjacency data rather than the corresponding spatial weights matrix to define the intrinsic CAR prior distribution for \(S_i\). (For a more detailed explanation of the following code, refer to the WinBUGS example provided in Appendix A.3.)

# Define the model (WinBUGS code)
model <- function(){
    for(i in 1:N){
        y[i] ~ dpois(E.mu[i])
        E.mu[i] <- E[i] * exp(mu[i])
        
        mu[i] <- beta.0 + beta.1 * x[i] + S[i] + U[i]
        U[i] ~ dnorm(0, tau.U)
    }
    beta.0 ~ dnorm(0, 0.01)
    beta.1 ~ dnorm(0, 0.01)
    
    S[1:N] ~ car.normal(adj[], weights[], num[], tau.S)
    tau.S <- pow(sigma.S,-1)
    sigma.S ~ dnorm(0, 0.1)%_%I(0,)
    
    tau.U <- pow(sigma.U, -1)
    sigma.U ~ dnorm(0, 0.1)%_%I(0,)
}

# Save model to a temporary location
fp.model <- file.path(tempdir(), "WinBUGS Model 8.3.bug")
write.model(model, fp.model)

# Fixed values (input data)
y <- data$observed      # Observed values
E <- data$expected      # Expected values
x <- data$covariate     # Covariate

# Initial values for stochastic parameters
set.seed(4321)
inits <- function() {list(
    beta.0 = 0,
    beta.1 = 0.5,
    S = rnorm(N, 0, 0.5),
    U = rep(0, N),
    sigma.S = 0.1,
    sigma.U = 1
)}

# MCMC parameters
M.burnin <- 20000       # Number of burn-in iterations (discarded)
M <- 5000               # Number of iterations retained (after thinning)
n.thin <- 5             # Thinning factor
# Fit the model using WinBUGS
MCMC <- bugs(
    data = c(list(N = N, y = y, E = E, x = x), adj.data),
    inits = inits,
    parameters.to.save = c("beta.0", "beta.1", "S", "U", "sigma.S", "sigma.U"),
    model.file = fp.model,
    n.chains = 1,
    n.burnin = M.burnin,
    n.iter = M.burnin + (M * n.thin), # Total iterations
    n.thin = n.thin,
    DIC = FALSE
)

# Store the list of parameters for convenience (ignore the other stuff in MCMC)
pars <- MCMC$sims.list

The posterior estimates of each parameter are stored in the object pars. Because this is the BYM model, we need to modify the parameters \(S\) and \(U\) to resolve the likelihood identifiability problem (refer back to Section 5.2.1). We also recreate the log-relative risk parameter.

# User-defined function to compute proportion of excess variation
get.psi <- function(S, U){
    sd.S <- sd(apply(S, 2, median))
    sd.U <- sd(apply(U, 2, median))
    psi <- sd.S / (sd.S + sd.U)
    return(psi)
}

# Compute proportion of excess variation
psi <- get.psi(pars$S, pars$U)

# Modify random effects
pars$S <- pars$S - psi * pars$U
pars$U <- pars$U + psi * pars$U

# Recreate mu (log-relative risk)
pars$mu <- pars$beta.0 + outer(pars$beta.1, x) + pars$S + pars$U

Graphical summaries of the marginal posterior distributions:

# Total number of sub-plots
gg.post <- vector("list", 7)

# Plot for beta_0
dat <- data.frame(val = pars$beta.0)
gg.post[[1]] <- ggplot(dat, aes(x = val)) + 
    geom_density(fill = "black", alpha = 0.5, color = NA) +
    scale_x_continuous(expression(beta[0])) +
    scale_y_continuous("")

# Plot for beta_1
dat <- data.frame(val = pars$beta.1)
gg.post[[2]] <- ggplot(dat, aes(x = val)) + 
    geom_density(fill = "black", alpha = 0.5, color = NA) +
    scale_x_continuous(expression(beta[1])) +
    scale_y_continuous("")

# Plot for S
dat <- pars$S
mean <- apply(dat, 2, mean)
CI.lower <- apply(dat, 2, quantile, probs = 0.025)
CI.upper <- apply(dat, 2, quantile, probs = 0.975)
rank <- order(order(mean))
dat <- data.frame(id = rank, mean, CI.lower, CI.upper)
gg.post[[3]] <- ggplot(dat, aes(x = id, y = mean)) + 
    geom_pointrange(aes(ymin = CI.lower, ymax = CI.upper), 
        alpha = 0.3, fatten = 0.1, colour = "grey50") +
    geom_point(size = 1) +
    scale_x_continuous("Ranked Ward", breaks = seq(0, N, 50)) +
    scale_y_continuous(expression(S)) +
    coord_flip()

# Plot for U
dat <- pars$U
mean <- apply(dat, 2, mean)
CI.lower <- apply(dat, 2, quantile, probs = 0.025)
CI.upper <- apply(dat, 2, quantile, probs = 0.975)
rank <- order(order(mean))
dat <- data.frame(id = rank, mean, CI.lower, CI.upper)
gg.post[[4]] <- ggplot(dat, aes(x = id, y = mean)) + 
    geom_pointrange(aes(ymin = CI.lower, ymax = CI.upper), 
        alpha = 0.3, fatten = 0.1, colour = "grey50") +
    geom_point(size = 1) +
    scale_x_continuous("Ranked Ward", breaks = seq(0, N, 50)) +
    scale_y_continuous(expression(U)) +
    coord_flip()

# Plot for relative risk (SIR)
SIR.labels <- c(0.125, 0.25, 0.5, 1, 2, 4, 8) # Ratio scale labels
dat <- pars$mu
mean <- apply(dat, 2, mean)
CI.lower <- apply(dat, 2, quantile, probs = 0.025)
CI.upper <- apply(dat, 2, quantile, probs = 0.975)
rank <- order(order(mean))
dat <- data.frame(id = rank, mean, CI.lower, CI.upper)
gg.post[[5]] <- ggplot(dat, aes(x = id, y = mean)) + 
    geom_pointrange(aes(ymin = CI.lower, ymax = CI.upper), 
        alpha = 0.3, fatten = 0.1, colour = "grey50") +
    geom_point(size = 1) +
    scale_x_continuous("Ranked Ward", breaks = seq(0, N, 50)) +
    scale_y_continuous("SIR", breaks = log(SIR.labels), labels = SIR.labels) +
    coord_flip(ylim = log(c(0.1, 10)))

# Plot for sigma.S
dat <- data.frame(val = pars$sigma.S)
gg.post[[6]] <- ggplot(dat, aes(x = val)) + 
    geom_density(fill = "black", alpha = 0.5, color = NA) +
    scale_x_continuous(expression(sigma[S])) +
    scale_y_continuous("")

# Plot for sigma.U
dat <- data.frame(val = pars$sigma.U)
gg.post[[7]] <- ggplot(dat, aes(x = val)) + 
    geom_density(fill = "black", alpha = 0.5, color = NA) +
    scale_x_continuous(expression(sigma[U])) +
    scale_y_continuous("")

# View combined plots
grid.arrange(gg.post[[1]], gg.post[[2]], gg.post[[3]], gg.post[[4]], 
    gg.post[[5]], gg.post[[6]], gg.post[[7]], ncol = 3
)

Now we decompose the log-risk surface (\(\mu_i\)) into its components: the structured and unstructured spatial random effects, \(S_i\) and \(U_i\) respectively, and the covariate effect, \(\beta_1 x_i\). For clarity, we combine the intercept with the structured spatial random effect.

log.SIR <- as.numeric(apply(pars$mu, 2, median))

# For clarity, truncate extreme SIR values
clip <- c(0.25, 4)  # Not necessary in this example since all values are within this range
SIR.clipped <- exp(log.SIR) %>% 
    replace(., which(. > clip[2]), clip[2]) %>%
    replace(., which(. < clip[1]), clip[1])

# Dataframe of the parameters we need to append to the shapefile
Append <- data.frame(
    id = 1:N,
    # int = median(pars$beta.0) + rnorm(N, 0, 0.001),
    S = apply(pars$S, 2, median) + median(pars$beta.0),
    U = apply(pars$U, 2, median),
    cov = median(pars$beta.1) * x,
    log.SIR = log.SIR,
    SIR = SIR.clipped
)

# Merge values with shapefile dataframe
dat <- inner_join(map.df, Append, by = "id")

gg.base <-  ggplot(data = NULL, aes(x = long, y = lat, group = group)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme_void() + 
    theme(legend.position = "bottom") +
    guides(fill = guide_colourbar(barwidth = 8)) +
    coord_map() # Keeps x-y ratio consistent

# Fill colours and values for SIR
Fill.colours.SIR <- c("#2C7BB6", "#ABD9E9", "#FFFFBF", "#FDAE61", "#D7191C")
Fill.values.SIR <- seq(log2(1/1.5), log2(1.5), length = 5)

# Fill colours and values for other parameters
Fill.colours.linear <- c("#00bbcc", "#66ffff", "white", "#d580ff", "#8600b3")
Lims <- Append[,2:5] %>% range %>% abs
Lims <- Lims[Lims %>% which.max]
Fill.values.linear <- seq(-Lims, Lims, length = 5)

# Range of (truncated) ratio scale values and breaks
Lims <- log2(range(Append$SIR))
if(Lims[2] > Fill.values.SIR[5]){
    Fill.values.SIR <- c(Fill.values.SIR, Lims[2])
    Fill.colours.SIR <- c(Fill.colours.SIR, Fill.colours.SIR[5])
}
if(Lims[1] < Fill.values.SIR[1]){
    Fill.values.SIR <- c(Lims[1], Fill.values.SIR)
    Fill.colours.SIR <- c(Fill.colours.SIR[1], Fill.colours.SIR)
}

# Rescaled fill values
Values.1 <- rescale(Fill.values.linear, from = range(Append$S))
Values.2 <- rescale(Fill.values.linear, from = range(Append$U))
Values.3 <- rescale(Fill.values.linear, from = range(Append$cov))
Values.4 <- rescale(Fill.values.linear, from = range(Append$log.SIR))
Values.5 <- rescale(Fill.values.SIR, from = range(log2(Append$SIR)))

Breaks.SIR <- c(0.5, 1/1.5, 1, 1.5, 2)

gg.1 <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = S), color = NA) +
    scale_fill_gradientn("", colours = Fill.colours.linear, values = Values.1,
        breaks = seq(-10, 10, 0.25)) +
    ggtitle("Structured SRE")

gg.2 <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = U), color = NA) +
    scale_fill_gradientn("", colours = Fill.colours.linear, values = Values.2,
        breaks = seq(-10, 10, 0.25)) +
    ggtitle("Unstructured SRE")

gg.3 <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = cov), color = NA) +
    scale_fill_gradientn("", colours = Fill.colours.linear, values = Values.3,
        breaks = seq(-10, 10, 0.25)) +
    ggtitle("Covariate Effect")

gg.4 <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log.SIR), color = NA) +
    scale_fill_gradientn("", colours = Fill.colours.linear, values = Values.4,
        breaks = seq(-10, 10, 0.5)) +
    ggtitle("Log-SIR")

gg.5 <- gg.base + 
    geom_polygon(data = map.border, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log2(SIR)), color = NA) +
    scale_fill_gradientn("", colours = Fill.colours.SIR, values = Values.5,
        breaks = log2(Breaks.SIR), labels = round(Breaks.SIR, 3)) +
    ggtitle("SIR") +
    annotate("rect",
        xmin = lims$xmin - 0.01, xmax = lims$xmax + 0.01,
        ymin = lims$ymin - 0.01, ymax = lims$ymax + 0.01,
        colour = "black", fill = NA)

# View combined plots
grid.arrange(
    arrangeGrob(gg.1, gg.2, gg.3, gg.4, gg.5, nrow = 1)
)

Here are the city insets for the SIR estimates.

# Create inset maps
gg.inset <- vector("list", 3)
gg.base.inset <- gg.base + guides(fill = FALSE, alpha = FALSE) +
    scale_fill_gradientn("", 
        colours = Fill.colours.SIR, 
        values = Values.5,
        limits = range(log2(Append$SIR))
    )

# Auckland
dat <- inner_join(map.1, Append, by = "id")
gg.inset[[1]] <- gg.base.inset + 
    geom_polygon(data = map.border.1, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log2(SIR)), color = NA) +
    ggtitle("Auckland")

# Wellington
dat <- inner_join(map.2, Append, by = "id")
gg.inset[[2]] <- gg.base.inset + 
    geom_polygon(data = map.border.2, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log2(SIR)), color = NA) +
    ggtitle("Wellington")

# Christchurch
dat <- inner_join(map.3, Append, by = "id")
gg.inset[[3]] <- gg.base.inset + 
    geom_polygon(data = map.border.3, fill = "grey30", color = "black", size = 0.8) +
    geom_polygon(data = dat, aes(fill = log2(SIR)), color = NA) +
    ggtitle("Christchurch")

# Render the plots
grid.arrange(
    arrangeGrob(gg.inset[[2]], gg.inset[[3]], ncol = 1, heights = c(7, 11)),
    arrangeGrob(gg.inset[[1]]), 
    ncol = 2, widths = c(1, 2)
)

Like the Queensland example in Section 7.3, we plot the observed and smoothed SIRs against the ranked SA2s to ascertain the effect of smoothing.

# Truncate observed and smoothed SIR values
clip <- c(0.2, 10)     # Wider limits to show major deviations from raw SIR
raw.SIR <- SIR.raw %>%
    replace(., which(. > clip[2]), clip[2]) %>%
    replace(., which(. < clip[1]), clip[1])
SIR <- exp(apply(pars$mu, 2, median)) %>% 
    replace(., which(. > clip[2]), clip[2]) %>%
    replace(., which(. < clip[1]), clip[1])

# Fill colours and values
Fill.colours <- c("#2C7BB6", "#2C7BB6", "#ABD9E9", "#FFFFBF", "#FDAE61", "#D7191C", 
    "#D7191C")
Fill.values <- c(log2(1/4.5), seq(log2(1/1.5), log2(1.5), length = 5), log2(4.5))
Fill.values.r <- rescale(Fill.values, na.rm = TRUE)

Breaks.SIR <- c(0.25, 0.5, 1/1.5, 1, 1.5, 2, 4)

dat <- data.frame(
    id = 1:N,
    SIR = SIR,
    raw.SIR = raw.SIR
)
dat$rank <- order(order(dat$raw.SIR))

ggplot(dat) +
    geom_hline(yintercept = 1, size = 1, colour = grey(0.3)) +
    geom_point(aes(x = rank, y = raw.SIR, colour = log2(raw.SIR), shape = "Observed"), 
        size = 3) +
    geom_point(aes(x = rank, y = SIR, shape = "Smoothed"), size = 1.5) +
    scale_x_continuous("Ranked SA2", expand = c(0.01, 0.01)) +
    scale_y_log10("", breaks = Breaks.SIR, 
        labels = round(Breaks.SIR, 3),
        expand = c(0.01, 0.01)) +
    theme(panel.grid.minor.y = element_blank()) + 
    scale_colour_gradientn("",
        colours = Fill.colours,
        values = Fill.values.r,
        limits = range(Fill.values),
        breaks = log2(Breaks.SIR),
        labels = round(Breaks.SIR, 3)
    ) +
    scale_shape_manual("",
        breaks = c("Observed", "Smoothed"),
        values = c("Observed" = 17, "Smoothed" = 16)
    ) +
    guides(
        colour = guide_colourbar(barheight = unit(8, "cm"))
    )

Note that there appears to be less smoothing occurring here than in the Queensland example.