2.1 Observed cases

pennLC$data contains the cases in each county stratified by race, gender and age. We can obtain the number of cases in each county, Y, by aggregating the rows of pennLC$data by county and adding up the observed number of cases.

library(dplyr)
d <- group_by(pennLC$data, county) %>% summarize(Y = sum(cases))
head(d)

2.2 Expected cases

Now we calculate the expected number of cases in each county using indirect standardization. The expected counts in each county represent the total number of disease cases one would expect if the population in the county behaved the way the Pennsylvania population behaves. We can do this by using the expected() function of the SpatialEpi package. This function has three arguments, namely,

• population: a vector of population counts for each strata in each area,
• cases: a vector with the number of cases for each strata in each area,
• n.strata: number of strata considered.

Vectors population and cases have to be sorted by area first and then, within each area, the counts for all strata need to be listed in the same order. All strata need to be included in the vectors, including strata with 0 cases.

In order to obtain the expected counts we first sort the data using the order() function where we specify the order as county, race, gender and finally age.

pennLC$data <- pennLC$data[order(pennLC$data$county, pennLC$data$race, pennLC$data$gender, pennLC$data$age), ]

Then, we obtain the expected counts E in each county by calling the expected() function where we set population equal to pennLC$data$population and cases equal to pennLC$data$cases. There are 2 races, 2 genders and 4 age groups for each county, so number of strata is set to 2 x 2 x 4 = 16.

E <- expected(population = pennLC$data$population, cases = pennLC$data$cases, n.strata = 16)

Now we add the vector E to the data frame d which contains the counties ids (county) and the observed counts (Y) making sure the E elements correspond to the counties in d$county in the same order. To do that, we use match() to calculate the vector of the positions that match d$county in unique(pennLC$data$county) which are the corresponding counties of E. Then we rearrange E using that vector.

d$E <- E[match(d$county, unique(pennLC$data$county))]
head(d)

2.3 Smokers proportions

We also add to d the variable smoking which represents the proportion of smokers in each county. We add this variable using the merge() function where we specify county as the column for merging.

d <- merge(d, pennLC$smoking, by = "county")

2.4 SMRs

Finally, we compute the vector of SMRs as the ratio of the observed to the expected counts, and add it to the data frame d.

d$SMR <- d$Y/d$E

2.5 Data

The final dataset is this:

head(d)

2.6 Add data to map

We have constructed a data frame d containing the observed and expected disease counts, the smokers proportions and the SMR for each of the counties. The map of Pennsylvania counties is given by the SpatialPolygons object called map. Using this object and the data frame d we can create a SpatialPolygonsDataFrame that will allow us to make maps of the variables in d. In order to do that, we first set the row names of the data frame d equal to d$county. Then we merge the SpatialPolygons map and the data frame d matching the SpatialPolygons member Polygons ID slot values with the data frame row names (match.ID = TRUE).

library(sp)
rownames(d) <- d$county
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)
head(map@data)

4.1 Model

Let \(Y_i\) and \(E_i\) be the observed and expected number of disesase cases, respectively, and let \(\theta_i\) be the relative risk for county \(i=1,\ldots,n\). The model is specified as follows:

\[Y_i|\theta_i \sim Poisson(E_i \times \theta_i),\ i=1,\ldots,n,\]

\[\log(\theta_i) = \beta_0 + \beta_1 \times smoking_i + u_i + v_i.\]

• \(\beta_0\) is the intercept,
• \(\beta_1\) is the coefficient of the smokers proportion covariate,
• \(u_i\) is a structured spatial effect, \(u_i|\mathbf{u_{-i}} \sim N(\bar{u}_{\delta_i}, \frac{1}{\tau_u n_{\delta_i}})\) (Conditionally autoregressive model (CAR)),
• \(v_i\) is an unstructured spatial effect, \(v_i \sim N(0, 1/\tau_v)\).

4.2 Neighbourhood matrix

We create the neighbourhood matrix needed to define the spatial random effect using the poly2nb() and the nb2INLA() functions of the spdep package. First, we use poly2nb() to create a neighbours list based on areas with contiguous boundaries. Each element of the list nb represents one area and contains the indices of its neighbours. For example, nb[[2]] contains the neighbours of area 2.

Then, we use nb2INLA() to convert this list into a file with the representation of the neighbourhood matrix as required by INLA. Then we read the file using the inla.read.graph() function of INLA, and store it in the object g which we will later use for specifying the spatial disease model with INLA.

library(spdep)
library(INLA)
nb <- poly2nb(map)
head(nb)
nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")

## [[1]]
## [1] 21 28 67
## 
## [[2]]
## [1]  3  4 10 63 65
## 
## [[3]]
## [1]  2 10 16 32 33 65
## 
## [[4]]
## [1]  2 10 37 63
## 
## [[5]]
## [1]  7 11 29 31 56
## 
## [[6]]
## [1] 15 36 38 39 46 54

4.3 Inference using INLA

The model includes two random effects, namely, \(u_i\) for modelling spatial residual variation, and \(v_i\) for modelling unstructured noise. We need to include two vectors in the data that denote the indices of these random effects. We call re_u the indices vector for \(u_i\), and re_v the indices vector for \(v_i\). We set both re_u and re_v equal to \(1,\ldots,n\), where \(n\) is the number of counties. In our example, \(n\)=67 and this can be obtained with the number of rows in the data (nrow(map@data)).

map$re_u <- 1:nrow(map@data)
map$re_v <- 1:nrow(map@data)

We specify the model formula by including the response in the left-hand side, and the fixed and random effects in the right-hand side. Random effects are set using f() with parameters equal to the name of the variable and the chosen model. For \(u_i\), we use model = "besag" with neighbourhood matrix given by g. For \(v_i\) we choose model = "iid".

formula <- Y ~ smoking + f(re_u, model = "besag", graph = g) + f(re_v, model = "iid")

We fit the model by calling the inla() function and using the default priors in INLA. We specify the formula, family, data, and the expected counts, and set control.predictor equal to list(compute = TRUE) to compute the posterior means of the predictors.

res <- inla(formula, family = "poisson", data = map@data, E = E, control.predictor = list(compute = TRUE))

4.4 Results

We can inspect the results object res using summary().

summary(res)

## 
## Call:
## c("inla(formula = formula, family = \"poisson\", data = map@data, ",  "    E = E, control.predictor = list(compute = TRUE))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.8867          0.6446          0.2443          2.7756 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.3236 0.1503    -0.6212  -0.3234    -0.0279 -0.3231   0
## smoking      1.1567 0.6247    -0.0810   1.1582     2.3853  1.1619   0
## 
## Random effects:
## Name   Model
##  re_u   Besags ICAR model 
## re_v   IID model 
## 
## Model hyperparameters:
##                        mean       sd 0.025quant 0.5quant 0.975quant
## Precision for re_u    93.05    49.95      30.95    81.84     220.42
## Precision for re_v 17957.16 18118.75    1155.88 12549.22   65968.48
##                       mode
## Precision for re_u   63.63
## Precision for re_v 3103.18
## 
## Expected number of effective parameters(std dev): 18.61(4.365)
## Number of equivalent replicates : 3.60 
## 
## Marginal log-Likelihood:  -320.53 
## Posterior marginals for linear predictor and fitted values computed

We see the intercept \(\hat \beta_0=\) -0.3236 with a 95% credible interval equal to (-0.6212, -0.0279), and the coefficient of smoking is \(\hat \beta_1=\) 1.1567 with a 95% credible interval equal to (-0.0810, 2.3853). We can plot the posterior distribution of the smoking coefficient. We do this by calculating a smoothing of the marginal distribution of the coefficient with inla.smarginal() and then plot it with ggplot() of the ggplot2 package.

library(ggplot2)
marginal <- inla.smarginal(res$marginals.fixed$smoking)
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() + labs(x = expression(beta[1]), y = "Density") +
  geom_vline(xintercept = 0, col = "blue") + theme_bw()

4.5 Add results to map

The disease risk estimates and uncertainty for each of the counties are given by the mean posterior and the 95% credible intervals of \(\theta_i\), \(i=1,\ldots,n\) which are in res$summary.fitted.values. Column mean is the mean posterior and 0.025quant and 0.975quant are the 2.5 and 97.5 percentiles, respectively. We add these data to map to be able to make maps of these variables. We assign mean to the estimate of the relative risk, and 0.025quant and 0.975quant to the lower and upper limits of 95% credible intervals of the risks.

head(res$summary.fitted.values)
map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]