# 9 Spatial modeling of geostatistical data. Malaria in The Gambia

In this chapter we show how to fit a geostatistical model to predict malaria prevalence in The Gambia using the stochastic partial differential equation (SPDE) approach and the R-INLA package . We use data of malaria prevalence in children obtained at 65 villages in The Gambia which are contained in the geoR package , and high-resolution environmental covariates downloaded with the raster package . We show how to create a triangulated mesh that covers The Gambia, the projection matrix and the data stacks to fit the model. Then we show how to manipulate the results to obtain the malaria prevalence predictions, and 95% credible intervals denoting uncertainty. We also show how to compute exceedance probabilities of prevalence being greater than a given threshold value of interest for policymaking. Results are shown by means of interactive maps created with the leaflet package .

## 9.1 Data

First, we load the geoR package and attach the data gambia which contains information about malaria prevalence in children obtained at 65 villages in The Gambia.

library(geoR)
data(gambia)

Next we inspect the data and see it is a data frame with 2035 observations and the following 8 variables:

• x: x coordinate of the village (UTM),
• y: y coordinate of the village (UTM),
• pos: presence (1) or absence (0) of malaria in a blood sample taken from the child,
• age: age of the child in days,
• netuse: indicator variable denoting whether the child regularly sleeps under a bed net,
• treated: indicator variable denoting whether the bed net is treated,
• green: satellite-derived measure of the greenness of vegetation in the vicinity of the village,
• phc: indicator variable denoting the presence or absence of a health center in the village.
head(gambia)
##           x       y pos  age netuse treated green phc
## 1850 349631 1458055   1 1783      0       0 40.85   1
## 1851 349631 1458055   0  404      1       0 40.85   1
## 1852 349631 1458055   0  452      1       0 40.85   1
## 1853 349631 1458055   1  566      1       0 40.85   1
## 1854 349631 1458055   0  598      1       0 40.85   1
## 1855 349631 1458055   1  590      1       0 40.85   1

## 9.2 Data preparation

Data in gambia are given at an individual level. Here, we do the analysis at the village level by aggregating the malaria tests by village. We create a data frame called d with columns containing, for each village, the longitude and latitude, the number of malaria tests, the number of positive tests, the prevalence, and the altitude.

### 9.2.1 Prevalence

We can see that gambia has 2035 rows and the matrix of the unique coordinates has 65 rows. This indicates that 2035 malaria tests were conducted at 65 locations.

dim(gambia)
##  2035    8
dim(unique(gambia[, c("x", "y")]))
##  65  2

We create a data frame called d containing, for each village, the coordinates (x, y), the total number of tests performed (total), the number of positive tests (positive), and the malaria prevalence (prev). In data gambia, column pos indicates the tests results. Positive tests have pos equal to 1 and negative tests have pos equal to 0. Therefore, we can calculate the number of positive tests in each village by adding up the elements in gambia$pos. Then we calculate the prevalence in each village by calculating the proportion of positive tests (number of positive results divided by the total number of tests in each village). We can create the data frame d using the dplyr package as follows: library(dplyr) d <- group_by(gambia, x, y) %>% summarize( total = n(), positive = sum(pos), prev = positive / total ) head(d) ## # A tibble: 6 x 5 ## # Groups: x  ## x y total positive prev ## <dbl> <dbl> <int> <dbl> <dbl> ## 1 349631. 1458055 33 17 0.515 ## 2 358543. 1460112 63 19 0.302 ## 3 360308. 1460026 17 7 0.412 ## 4 363796. 1496919 24 8 0.333 ## 5 366400. 1460248 26 10 0.385 ## 6 366688. 1463002 18 7 0.389 An alternative to dplyr to calculate the number of positive tests in each village is to use the aggregate() function. We can use aggregate() to obtain the total number of tests (total) and positive results (positive) in each village. Then we can calculate the prevalence prev by dividing positive by total, and create the data frame d with these vectors and the x and y coordinates. total <- aggregate( gambia$pos,
by = list(gambia$x, gambia$y),
FUN = length
)
positive <- aggregate(
gambia$pos, by = list(gambia$x, gambia$y), FUN = sum ) prev <- positive$x / total$x d <- data.frame( x = total$Group.1,
y = total$Group.2, total = total$x,
positive = positive$x, prev = prev ) ### 9.2.2 Transforming coordinates Now we can plot the malaria prevalence in a map created with the leaflet() function of the leaflet package. leaflet() expects data to be specified in geographic coordinates (longitude/latitude). However, the coordinates in the data are in UTM format (Easting/Northing). We transform the UTM coordinates in the data to geographic coordinates using the spTransform() function of the sp package . First, we create a SpatialPoints object called sps where we specify the projection of The Gambia, that is, UTM zone 28. Then, we transform the UTM coordinates in sps to geographic coordinates using spTransform() where we set CRS to CRS("+proj=longlat +datum=WGS84"). library(sp) library(rgdal) sps <- SpatialPoints(d[, c("x", "y")], proj4string = CRS("+proj=utm +zone=28") ) spst <- spTransform(sps, CRS("+proj=longlat +datum=WGS84")) Finally, we add the longitude and latitude variables to the data frame d. d[, c("long", "lat")] <- coordinates(spst) head(d) ## # A tibble: 6 x 7 ## # Groups: x  ## x y total positive prev long lat ## <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> ## 1 349631. 1458055 33 17 0.515 -16.4 13.2 ## 2 358543. 1460112 63 19 0.302 -16.3 13.2 ## 3 360308. 1460026 17 7 0.412 -16.3 13.2 ## 4 363796. 1496919 24 8 0.333 -16.3 13.5 ## 5 366400. 1460248 26 10 0.385 -16.2 13.2 ## 6 366688. 1463002 18 7 0.389 -16.2 13.2 ### 9.2.3 Mapping prevalence Now we use leaflet() to create a map with the locations of the villages and the malaria prevalence (Figure 9.1). We use addCircles() to put circles on the map, and color circles according to the value of the prevalences. We choose a palette function with colors from viridis and four equal intervals from values 0 to 1. We use the base map given by providers$CartoDB.Positron so that point colors can be distinguised from the background map. We also add a scale bar using addScaleBar().

library(leaflet)
library(viridis)

pal <- colorBin("viridis", bins = c(0, 0.25, 0.5, 0.75, 1))
leaflet(d) %>%
addProviderTiles(providers$CartoDB.Positron) %>% addCircles(lng = ~long, lat = ~lat, color = ~ pal(prev)) %>% addLegend("bottomright", pal = pal, values = ~prev, title = "Prev." ) %>% addScaleBar(position = c("bottomleft")) FIGURE 9.1: Malaria prevalence in The Gambia. ### 9.2.4 Environmental covariates We model malaria prevalence using a covariate that indicates the altitude in The Gambia. This covariate can be obtained with the getData() function of the raster package. This package can be used to obtain geographic data from anywhere in the world. In order to get the altitude values in The Gambia, we need to call getData() with the three following arguments: • name of the data equal to "alt", • country equal to the 3 letters of the International Organization for Standardization (ISO) code of The Gambia (GMB), • mask equal to TRUE so the neighboring countries are set to NA. library(raster) r <- getData(name = "alt", country = "GMB", mask = TRUE) We make a map with the altitude raster using the addRasterImage() function of leaflet. To create this map, we use a palette function pal with the values of the raster (values(r)) and specifying that the NA values are transparent (Figure 9.2). pal <- colorNumeric("viridis", values(r), na.color = "transparent" ) leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>%
addRasterImage(r, colors = pal, opacity = 0.5) %>%
pal = pal, values = values(r),
title = "Altitude"
) %>%
addScaleBar(position = c("bottomleft"))

FIGURE 9.2: Map of altitude in The Gambia.

Now we add the altitude values to the data frame d to be able to use it as a covariate in the model. To do that, we get the altitude values at the village locations using the extract() function of raster. The first argument of this function is the altitude raster (r). The second argument is a two-column matrix with the coordinates where we want to know the values, that is, the coordinates of the villages given by d[, c("long", "lat")]. We assign the altitude vector to the column alt of the data frame d.

mesh$n ##  673 plot(mesh) points(coo, col = "red") FIGURE 9.3: Triangulated mesh used to build the SPDE model. ### 9.3.3 Building the SPDE model on the mesh Then, we use the inla.spde2.matern() function to build the SPDE model. spde <- inla.spde2.matern(mesh = mesh, alpha = 2, constr = TRUE) Here, we set constr = TRUE to impose an integrate-to-zero constraint. Here alpha is a parameter related to the smoothness parameter of the process, namely, $$\alpha=\nu+d/2$$. In this example, we set the smoothness parameter $$\nu$$ equal to 1 and in the spatial case $$d=2$$ so alpha=1+2/2=2. ### 9.3.4 Index set Now we generate the index set for the SPDE model. We do this with the function inla.spde.make.index() where we specify the name of the effect (s) and the number of vertices in the SPDE model (spde$n.spde). This creates a list with vector s equal to 1:spde$n.spde, and vectors s.group and s.repl that have all elements equal to 1s and length given by the number of mesh vertices. indexs <- inla.spde.make.index("s", spde$n.spde)
lengths(indexs)
##       s s.group  s.repl
##     673     673     673

### 9.3.5 Projection matrix

We need to build a projection matrix $$A$$ that projects the spatially continuous Gaussian random field from the observations to the mesh nodes. The projection matrix can be built with the inla.spde.make.A() function passing the mesh and the coordinates.

A <- inla.spde.make.A(mesh = mesh, loc = coo)

### 9.3.6 Prediction data

Here we specify the locations where we wish to predict the prevalence. We set the prediction locations to the locations of the raster of the covariate altitude. We can get the coordinates of the raster r with the function rasterToPoints() of the raster package. This function returns a matrix with the coordinates and values of the raster that do not have NA values. We see there are 12964 points.

dp <- rasterToPoints(r)
dim(dp)
##  12964     3

In this example, we use fewer prediction points so the computation is faster. We can lower the resolution of the raster by using the aggregate() function of raster. The arguments of the function are

• x: raster object,
• fact: aggregation factor expressed as number of cells in each direction (horizontally and vertically),
• fun: function used to aggregate values.

We specify fact = 5 to aggregate 5 cells in each direction, and fun = mean to compute the mean of the cell values. We call coop to the matrix of coordinates with the prediction locations.

ra <- aggregate(r, fact = 5, fun = mean)

dp <- rasterToPoints(ra)
dim(dp)
##  627   3
coop <- dp[, c("x", "y")]

We also construct the matrix that projects the spatially continuous Gaussian random field from the prediction locations to the mesh nodes.

Ap <- inla.spde.make.A(mesh = mesh, loc = coop)

### 9.3.7 Stack with data for estimation and prediction

Now we use the inla.stack() function to organize data, effects, and projection matrices. We use the following arguments:

• tag: string to identify the data,
• data: list of data vectors,
• A: list of projection matrices,
• effects: list with fixed and random effects.

We construct a stack called stk.e with data for estimation and we tag it with the string "est". The fixed effects are the intercept (b0) and a covariate (altitude). The random effect is the spatial Gaussian random field (s). Therefore, in effects we pass a list with a data.frame with the fixed effects, and a list s containing the indices of the SPDE object (indexs). A is set to a list where the second element is A, the projection matrix for the random effects, and the first element is 1 to indicate the fixed effects are mapped one-to-one to the response. In data we specify the response vector and the number of trials. We also construct a stack for prediction that called stk.p. This stack has tag equal to "pred", the response vector is set to NA, and the data is specified at the prediction locations. Finally, we put stk.e and stk.p together in a full stack stk.full.

# stack for estimation stk.e
stk.e <- inla.stack(
tag = "est",
data = list(y = d$positive, numtrials = d$total),
A = list(1, A),
effects = list(data.frame(b0 = 1, altitude = d$alt), s = indexs) ) # stack for prediction stk.p stk.p <- inla.stack( tag = "pred", data = list(y = NA, numtrials = NA), A = list(1, Ap), effects = list(data.frame(b0 = 1, altitude = dp[, 3]), s = indexs ) ) # stk.full has stk.e and stk.p stk.full <- inla.stack(stk.e, stk.p) ### 9.3.8 Model formula 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. In the formula, we remove the intercept (adding 0) and add it as a covariate term (adding b0), so all the covariate terms can be captured in the projection matrix. formula <- y ~ 0 + b0 + altitude + f(s, model = spde) ### 9.3.9inla() call We fit the model by calling inla() and using the default priors in R-INLA. We specify the formula, family, data, and options. In control.predictor we set compute = TRUE to compute the posteriors of the predictions. We set link=1 to compute the fitted values (res$summary.fitted.values and res$marginals.fitted.values) with the same link function as the family specified in the model. res <- inla(formula, family = "binomial", Ntrials = numtrials, control.family = list(link = "logit"), data = inla.stack.data(stk.full), control.predictor = list( compute = TRUE, link = 1, A = inla.stack.A(stk.full) ) ) ## 9.4 Mapping malaria prevalence Now we map the malaria prevalence predictions using leaflet. The mean prevalence and lower and upper limits of 95% credible intervals are in the data frame res$summary.fitted.values. The rows of res$summary.fitted.values that correspond to the prediction locations can be obtained by selecting the indices of the stack stk.full that are tagged with tag = "pred". We can obtain these indices by using inla.stack.index() passing stk.full and tag = "pred". index <- inla.stack.index(stack = stk.full, tag = "pred")$data

We create vectors with the mean prevalence and lower and upper limits of 95% credible intervals with the values of the columns "mean", "0.025quant" and "0.975quant" and the rows given by index.

prev_mean <- res$summary.fitted.values[index, "mean"] prev_ll <- res$summary.fitted.values[index, "0.025quant"]
prev_ul <- res$summary.fitted.values[index, "0.975quant"] Now we create a map with the predicted prevalence (prev_mean) at the prediction locations coop using the addCircles() function of leaflet. We use a palette function created with colorNumeric(). We could set the domain to the range of values we plot (prev_mean), but we decide to use the interval [0, 1] because we will use the same palette to plot the mean prevalence and the lower and upper limits of the 95% credible intervals (Figure 9.4). pal <- colorNumeric("viridis", c(0, 1), na.color = "transparent") leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>%
lng = coop[, 1], lat = coop[, 2],
color = pal(prev_mean)
) %>%
pal = pal, values = prev_mean,
title = "Prev."
) %>%
addScaleBar(position = c("bottomleft"))

FIGURE 9.4: Map of malaria prevalence predictions in The Gambia created with leaflet() and addCircles().

Instead of showing the prevalence predictions at points, we can also plot them using a raster. Here, the prediction locations coop are not on a regular grid; therefore, we need to create a raster with the predicted values using the rasterize() function of raster. This function is useful for transferring the values of a vector at some locations to raster cells. We transfer the predicted values prev_mean from the locations coop to the raster ra that we used to get the prediction locations. We use the following arguments:

• x = coop: coordinates where we made the predictions,
• y = ra: raster where we transfer the values,
• field = prev_mean: values to be transferred (prevalence predictions in locations coop),
• fun = mean: to assign the mean of the values to cells that have more than one point.
r_prev_mean <- rasterize(
x = coop, y = ra, field = prev_mean,
fun = mean
)

Now we plot the raster with the prevalence values in a map using the addRasterImage() function of leaflet (Figure 9.5).

pal <- colorNumeric("viridis", c(0, 1), na.color = "transparent")

leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>% addRasterImage(r_prev_mean, colors = pal, opacity = 0.5) %>% addLegend("bottomright", pal = pal, values = values(r_prev_mean), title = "Prev." ) %>% addScaleBar(position = c("bottomleft")) FIGURE 9.5: Map of malaria prevalence predictions in The Gambia created with leaflet() and addRasterImage(). We can follow the same approach to create maps with the lower and upper limits of the prevalence predictions. First, we create rasters with the lower and upper limits of the predictions. Then we make the maps with the same palette function we used to create the map of mean prevalence predictions (Figures 9.6 and 9.7). r_prev_ll <- rasterize( x = coop, y = ra, field = prev_ll, fun = mean ) leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>%
addRasterImage(r_prev_ll, colors = pal, opacity = 0.5) %>%
pal = pal,
values = values(r_prev_ll), title = "LL"
) %>%
addScaleBar(position = c("bottomleft"))

FIGURE 9.6: Map of lower limit of 95% CI of malaria prevalence in The Gambia.

r_prev_ul <- rasterize(
x = coop, y = ra, field = prev_ul,
fun = mean
)

leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>% addRasterImage(r_prev_ul, colors = pal, opacity = 0.5) %>% addLegend("bottomright", pal = pal, values = values(r_prev_ul), title = "UL" ) %>% addScaleBar(position = c("bottomleft")) FIGURE 9.7: Map of upper limit of 95% CI of malaria prevalence in The Gambia. ## 9.5 Mapping exceedance probabilities We can also calculate the exceedance probabilities of malaria prevalence being greater than a given threshold value that is of interest for policymaking. For example, we can be interested in knowing what are the probabilities that malaria prevalence is greater than 20%. Let $$p_i$$ be the malaria prevalence at location $$\boldsymbol{x}_i$$. The probability that the malaria prevalence $$p_i$$ is greater than a value $$c$$ can be written as $$P(p_i > c)$$. This probability can be calculated by substracting $$P(p_i \leq c)$$ to 1 as $P(p_i > c) = 1 - P(p_i \leq c).$ In R-INLA, $$P(p_i \leq c)$$ can be calculated using the inla.pmarginal() function with arguments on the posterior distribution of the predictions and the threshold value. Then, the exceedance probability $$P(p_i > c)$$ can be calculated as 1 - inla.pmarginal(q = c, marginal = marg) where marg is the marginal distribution of the predictions, and c is the threshold value. In our example, we can calculate the probabilites that malaria prevalence exceeds 20% as follows. First, we obtain the posterior marginals of the predictions for each location. These marginals are in the list object res$marginals.fitted.values[index] where index is the vector of the indices of the stack stk.full corresponding to the predictions. In the previous section, we obtained these indices by using the inla.stack.index() function and specifying tag = "pred".

index <- inla.stack.index(stack = stk.full, tag = "pred")$data The first element of the list, res$marginals.fitted.values[index][], contains the marginal distribution of the prevalence prediction corresponding to the first location. The probability that malaria prevalence exceeds 20% at this location is given by

marg <- res$marginals.fitted.values[index][] 1 - inla.pmarginal(q = 0.20, marginal = marg) ##  0.6931 To compute the exceedance probabilities for all the prediction locations, we can use the sapply() function with two arguments. The first argument denotes the marginal distributions of the predictions (res$marginals.fitted.values[index]) and the second argument denotes the function to compute the exceedance probabilities (1- inla.pmarginal()). Then the sapply() function returns a vector of the same length as the list res$marginals.fitted.values[index], where each element is the result of applying the function 1- inla.pmarginal() to the corresponding element of the list of marginals. excprob <- sapply(res$marginals.fitted.values[index],
FUN = function(marg){1-inla.pmarginal(q = 0.20, marginal = marg)})

head(excprob)
## fitted.APredictor.066 fitted.APredictor.067
##                0.6931                0.7213
## fitted.APredictor.068 fitted.APredictor.069
##                0.7742                0.7467
## fitted.APredictor.070 fitted.APredictor.071
##                0.6634                0.6994

Finally, we can make maps of the exceedance probabilities by creating a raster with the exceedance probabilities with rasterize() and then using leaflet() to create the map.

r_excprob <- rasterize(
x = coop, y = ra, field = excprob,
fun = mean
)
pal <- colorNumeric("viridis", c(0, 1), na.color = "transparent")

leaflet() %>%
addScaleBar(position = c("bottomleft"))