Type:	Package
Title:	Disaggregation Modelling
Version:	0.4.0
Description:	Fits disaggregation regression models using 'TMB' ('Template Model Builder'). When the response data are aggregated to polygon level but the predictor variables are at a higher resolution, these models can be useful. Regression models with spatial random fields. The package is described in detail in Nandi et al. (2023) <doi:10.18637/jss.v106.i11>.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.2
Imports:	splancs, Matrix, stats, TMB, dplyr, ggplot2, cowplot, rSPDE, sparseMVN, fmesher, tidyterra, terra, sf, utils
Suggests:	testthat, knitr, rmarkdown, SpatialEpi
LinkingTo:	TMB, RcppEigen
SystemRequirements:	GNU make
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2024-10-02 16:23:34 UTC; simon
Author:	Anita Nandi [aut], Tim Lucas [aut, cre], Rohan Arambepola [aut], Andre Python [aut], Simon Smart [ctb]
Maintainer:	Tim Lucas <timcdlucas@gmail.com>
Repository:	CRAN
Date/Publication:	2024-10-02 20:40:03 UTC

Function to fit the disaggregation model

Description

Function to fit the disaggregation model

Usage

as.disag_data(
  polygon_shapefile,
  shapefile_names,
  covariate_rasters,
  polygon_data,
  covariate_data,
  aggregation_pixels,
  coords_for_fit,
  coords_for_prediction,
  start_end_index,
  mesh = NULL,
  coordsForFit = NULL,
  coordsForPrediction = NULL,
  startendindex = NULL
)

Arguments

Value

A list is returned of class disag_data. The functions summary, print and plot can be used on disag_data. The list of class disag_data contains:

Build mesh for disaggregaton model

Description

build_mesh function takes a sf object and mesh arguments to build an appropriate mesh for the spatial field.

Usage

build_mesh(shapes, mesh_args = NULL, mesh.args = NULL)

Arguments

Details

The mesh is created by finding a tight boundary around the polygon data, and creating a fine mesh within the boundary and a coarser mesh outside. This speeds up computation time by only having a very fine mesh within the area of interest and having a small region outside with a coarser mesh to avoid edge effects.

Six mesh parameters can be specified as arguments: convex, concave and resolution, to control the boundary of the inner mesh, and max.edge, cutoff and offset, to control the mesh itself, with the names meaning the same as used by the fmesher functions fm_nonconvex_hull_inla and fm_mesh_2d.

Defaults are: pars <- list(convex = -0.01, concave = -0.5, resolution = 300, max.edge = c(3.0, 8), cutoff = 0.4, offset = c(1, 15)).

Value

An inla.mesh object

Examples

## Not run: 
polygons <- list()
for(i in 1:14) {
  row <- ceiling(i/10)
  col <- ifelse(i %% 10 != 0, i %% 10, 10)
  xmin = 2*(col - 1); xmax = 2*col; ymin = 2*(row - 1); ymax = 2*row
polygons[[i]] <- list(cbind(c(xmin, xmax, xmax, xmin, xmin),
                            c(ymax, ymax, ymin, ymin, ymax)))
}

polys <- lapply(polygons, sf::st_polygon)
response_df <- data.frame(area_id = 1:100,
                          response = runif(100, min = 0, max = 10))
spdf <- sf::st_sf(polys, response_df)

my_mesh <- build_mesh(spdf)

## End(Not run)

Fit the disaggregation model

Description

disag_model function takes a disag_data object created by prepare_data and performs a Bayesian disaggregation fit.

Usage

disag_model(
  data,
  priors = NULL,
  family = "gaussian",
  link = "identity",
  iterations = 100,
  field = TRUE,
  iid = TRUE,
  hess_control_parscale = NULL,
  hess_control_ndeps = 1e-04,
  silent = TRUE
)

Arguments

Details

The model definition

The disaggregation model makes predictions at the pixel level:

link(pred_i) = \beta_0 + \beta X + GP(s_i) + u_i

And then aggregates these predictions to the polygon level using the weighted sum (via the aggregation raster, agg_i):

cases_j = \sum_{i \epsilon j} pred_i \times agg_i

rate_j = \frac{\sum_{i \epsilon j} pred_i \times agg_i}{\sum_{i \epsilon j} agg_i}

The different likelihood correspond to slightly different models (y_j is the response count data):

• Gaussian: If \sigma is the dispersion of the pixel data, \sigma_j is the dispersion of the polygon data, where \sigma_j = \sigma \sqrt{\sum agg_i^2} / \sum agg_i

  dnorm(y_j/\sum agg_i, rate_j, \sigma_j)

  - predicts incidence rate.

• Binomial: For a survey in polygon j, y_j is the number positive and N_j is the number tested.

  dbinom(y_j, N_j, rate_j)

  - predicts prevalence rate.

• Poisson:

  dpois(y_j, cases_j)

  - predicts incidence count.

Specify priors for the regression parameters, field and iid effect as a single list. Hyperpriors for the field are given as penalised complexity priors you specify \rho_{min} and \rho_{prob} for the range of the field where P(\rho < \rho_{min}) = \rho_{prob}, and \sigma_{min} and \sigma_{prob} for the variation of the field where P(\sigma > \sigma_{min}) = \sigma_{prob}. Also, specify pc priors for the iid effect

The family and link arguments are used to specify the likelihood and link function respectively. The likelihood function can be one of gaussian, poisson or binomial. The link function can be one of logit, log or identity. These are specified as strings.

The field and iid effect can be turned on or off via the field and iid logical flags. Both are default TRUE.

The iterations argument specifies the maximum number of iterations the model can run for to find an optimal point.

The silent argument can be used to publish/suppress verbose output. Default TRUE.

Value

A list is returned of class disag_model. The functions summary, print and plot can be used on disag_model. The list of class disag_model contains:

References

Nanda et al. (2023) disaggregation: An R Package for Bayesian Spatial Disaggregation Modeling. <doi:10.18637/jss.v106.i11>

Examples

## Not run: 
polygons <- list()
n_polygon_per_side <- 10
n_polygons <- n_polygon_per_side * n_polygon_per_side
n_pixels_per_side <- n_polygon_per_side * 2

for(i in 1:n_polygons) {
  row <- ceiling(i/n_polygon_per_side)
  col <- ifelse(i %% n_polygon_per_side != 0, i %% n_polygon_per_side, n_polygon_per_side)
  xmin = 2*(col - 1); xmax = 2*col; ymin = 2*(row - 1); ymax = 2*row
  polygons[[i]] <- list(cbind(c(xmin, xmax, xmax, xmin, xmin),
                              c(ymax, ymax, ymin, ymin, ymax)))
}

polys <- lapply(polygons,sf::st_polygon)
N <- floor(runif(n_polygons, min = 1, max = 100))
response_df <- data.frame(area_id = 1:n_polygons, response = runif(n_polygons, min = 0, max = 1000))

spdf <- sf::st_sf(response_df, geometry = polys)

# Create raster stack
r <- terra::rast(ncol=n_pixels_per_side, nrow=n_pixels_per_side)
terra::ext(r) <- terra::ext(spdf)
r[] <- sapply(1:terra::ncell(r), function(x){
rnorm(1, ifelse(x %% n_pixels_per_side != 0, x %% n_pixels_per_side, n_pixels_per_side), 3))}
r2 <- terra::rast(ncol=n_pixels_per_side, nrow=n_pixels_per_side)
terra::ext(r2) <- terra::ext(spdf)
r2[] <- sapply(1:terra::ncell(r), function(x) rnorm(1, ceiling(x/n_pixels_per_side), 3))
cov_stack <- c(r, r2)
names(cov_stack) <- c('layer1', 'layer2')

test_data <- prepare_data(polygon_shapefile = spdf,
                          covariate_rasters = cov_stack)

 result <- disag_model(test_data, iterations = 2)
 
## End(Not run)

Roxygen commands

Description

Roxygen commands

Usage

dummy()

Get a SpatRaster of covariates from a folder containing .tif files

Description

Looks in a specified folder for raster files. Returns a multi-layered SpatRaster of the rasters cropped to the extent specified by the shape parameter.

Usage

getCovariateRasters(directory, file_pattern = ".tif$", shape)

Arguments

Value

A multi-layered SpatRaster of the raster files in the directory

Examples

## Not run: 
  getCovariateRasters('/home/rasters', '.tif$', shape)
 
## End(Not run)

Extract polygon id and response data into a data.frame from a sf object

Description

Returns a data.frame with a row for each polygon in the sf object and columns: area_id, response and N, containing the id of the polygon, the values of the response for that polygon, and the sample size respectively. If the data is not survey data (the sample size does not exist), this column will contain NAs.

Usage

getPolygonData(
  shape,
  id_var = "area_id",
  response_var = "response",
  sample_size_var = NULL
)

Arguments

Value

A data.frame with a row for each polygon in the sf object and columns: area_id, response and N, containing the id of the polygon, the values of the response for that polygon, and the sample size respectively. If the data is not survey data (the sample size does not exist), this column will contain NAs.

Examples

{
polygons <- list()
for(i in 1:100) {
  row <- ceiling(i/10)
  col <- ifelse(i %% 10 != 0, i %% 10, 10)
  xmin = 2*(col - 1); xmax = 2*col; ymin = 2*(row - 1); ymax = 2*row
  polygons[[i]] <- list(cbind(c(xmin, xmax, xmax, xmin, xmin),
                              c(ymax, ymax, ymin, ymin, ymax)))
}

polys <- lapply(polygons,sf::st_polygon)
response_df <- data.frame(area_id = 1:100, response = runif(100, min = 0, max = 10))
spdf <- sf::st_sf(response_df, geometry = polys)

 getPolygonData(spdf, id_var = 'area_id', response_var = 'response')
}

Function to match pixels to their corresponding polygon

Description

From the covariate data and polygon data, the function matches the polygon id between the two to find which pixels from the covariate data are contained in each of the polygons.

Usage

getStartendindex(covariates, polygon_data, id_var = "area_id")

Arguments

Details

Takes a data.frame containing the covariate data with a polygon id column and one column for each covariate, and another data.frame containing polygon data with a polygon id, response and sample size column (as returned by getPolygonData function).

Returns a matrix with two columns and one row for each polygon. The first column is the index of the first row in covariate data that corresponds to that polygon, the second column is the index of the last row in covariate data that corresponds to that polygon.

Value

A matrix with two columns and one row for each polygon. The first column is the index of the first row in covariate data that corresponds to that polygon, the second column is the index of the last row in covariate data that corresponds to that polygon.

Examples

{
 covs <- data.frame(area_id = c(1, 1, 1, 2, 2, 3, 3, 3, 3), response = c(3, 9, 5, 2, 3, 6, 7, 3, 5))
 response <- data.frame(area_id = c(1, 2, 3), response = c(4, 7, 2), N = c(NA, NA, NA))
 getStartendindex(covs, response, 'area_id')
}

Create the TMB model object for the disaggregation model

Description

make_model_object function takes a disag_data object created by prepare_data and creates a TMB model object to be used in fitting.

Usage

make_model_object(
  data,
  priors = NULL,
  family = "gaussian",
  link = "identity",
  field = TRUE,
  iid = TRUE,
  silent = TRUE
)

Arguments

Details

The model definition

The disaggregation model make predictions at the pixel level:

link(pred_i) = \beta_0 + \beta X + GP(s_i) + u_i

And then aggregates these predictions to the polygon level using the weighted sum (via the aggregation raster, agg_i):

cases_j = \sum_{i \epsilon j} pred_i \times agg_i

rate_j = \frac{\sum_{i \epsilon j} pred_i \times agg_i}{\sum_{i \epsilon j} agg_i}

The different likelihood correspond to slightly different models (y_j is the response count data):

• Gaussian: If \sigma is the dispersion of the pixel data, \sigma_j is the dispersion of the polygon data, where \sigma_j = \sigma \sqrt{\sum agg_i^2} / \sum agg_i

  dnorm(y_j/\sum agg_i, rate_j, \sigma_j)

  - predicts incidence rate.

• Binomial: For a survey in polygon j, y_j is the number positive and N_j is the number tested.

  dbinom(y_j, N_j, rate_j)

  - predicts prevalence rate.

• Poisson:

  dpois(y_j, cases_j)

  - predicts incidence count.

Specify priors for the regression parameters, field and iid effect as a single named list. Hyperpriors for the field are given as penalised complexity priors you specify \rho_{min} and \rho_{prob} for the range of the field where P(\rho < \rho_{min}) = \rho_{prob}, and \sigma_{min} and \sigma_{prob} for the variation of the field where P(\sigma > \sigma_{min}) = \sigma_{prob}. Also, specify pc priors for the iid effect.

The precise names and default values for these priors are:

• priormean_intercept: 0

• priorsd_intercept: 10.0

• priormean_slope: 0.0

• priorsd_slope: 0.5

• prior_rho_min: A third the length of the diagonal of the bounding box.

• prior_rho_prob: 0.1

• prior_sigma_max: sd(response/mean(response))

• prior_sigma_prob: 0.1

• prior_iideffect_sd_max: 0.1

• prior_iideffect_sd_prob: 0.01

The family and link arguments are used to specify the likelihood and link function respectively. The likelihood function can be one of gaussian, poisson or binomial. The link function can be one of logit, log or identity. These are specified as strings.

The field and iid effect can be turned on or off via the field and iid logical flags. Both are default TRUE.

The iterations argument specifies the maximum number of iterations the model can run for to find an optimal point.

The silent argument can be used to publish/supress verbose output. Default TRUE.

Value

The TMB model object returned by MakeADFun.

Examples

## Not run: 
polygons <- list()
n_polygon_per_side <- 10
n_polygons <- n_polygon_per_side * n_polygon_per_side
n_pixels_per_side <- n_polygon_per_side * 2

for(i in 1:n_polygons) {
  row <- ceiling(i/n_polygon_per_side)
  col <- ifelse(i %% n_polygon_per_side != 0, i %% n_polygon_per_side, n_polygon_per_side)
  xmin = 2*(col - 1); xmax = 2*col; ymin = 2*(row - 1); ymax = 2*row
  polygons[[i]] <- list(cbind(c(xmin, xmax, xmax, xmin, xmin),
                              c(ymax, ymax, ymin, ymin, ymax)))
}

polys <- lapply(polygons,sf::st_polygon)
N <- floor(runif(n_polygons, min = 1, max = 100))
response_df <- data.frame(area_id = 1:n_polygons, response = runif(n_polygons, min = 0, max = 1000))

spdf <- sf::st_sf(response_df, geometry = polys)

# Create raster stack
r <- terra::rast(ncol=n_pixels_per_side, nrow=n_pixels_per_side)
terra::ext(r) <- terra::ext(spdf)
r[] <- sapply(1:terra::ncell(r), function(x){
rnorm(1, ifelse(x %% n_pixels_per_side != 0, x %% n_pixels_per_side, n_pixels_per_side), 3))}
r2 <- terra::rast(ncol=n_pixels_per_side, nrow=n_pixels_per_side)
terra::ext(r2) <- terra::ext(spdf)
r2[] <- sapply(1:terra::ncell(r), function(x) rnorm(1, ceiling(x/n_pixels_per_side), 3))
cov_stack <- c(r, r2)
names(cov_stack) <- c('layer1', 'layer2')

test_data <- prepare_data(polygon_shapefile = spdf,
                          covariate_rasters = cov_stack)

 result <- make_model_object(test_data)
 
## End(Not run)

Plot input data for disaggregation

Description

Plotting function for class disag_data (the input data for disaggregation).

Usage

## S3 method for class 'disag_data'
plot(x, which = c(1, 2, 3), ...)

Arguments

Details

Produces three plots: polygon response data, covariate rasters and INLA mesh.

Value

A list of three plots: the polygon plot (ggplot), covariate plot (spplot) and INLA mesh plot (ggplot)

Plot results of fitted model

Description

Plotting function for class disag_model (the result of the disaggregation fitting).

Usage

## S3 method for class 'disag_model'
plot(x, include_iid = FALSE, ...)

Arguments

Details

Produces two plots: results of the fixed effects and in-sample observed vs predicted plot.

Value

A list of two ggplot plots: results of the fixed effects and an in-sample observed vs predicted plot

Plot mean and uncertainty predictions from the disaggregation model results

Description

Plotting function for class disag_prediction (the mean and uncertainty predictions of the disaggregation fitting).

Usage

## S3 method for class 'disag_prediction'
plot(x, ...)

Arguments

Details

Produces raster plots of the mean prediction, and the lower and upper confidence intervals.

Value

A list of plots of rasters from the prediction: mean prediction, lower CI and upper CI.

Convert results of the model ready for plotting

Description

Convert results of the model ready for plotting

Usage

plot_disag_model_data(x)

Arguments

Value

A list that contains:

Predict mean and uncertainty from the disaggregation model result

Description

predict.disag_model function takes a disag_model object created by disaggregation::disag_model and predicts mean and uncertainty maps.

Usage

## S3 method for class 'disag_model'
predict(
  object,
  new_data = NULL,
  predict_iid = FALSE,
  N = 100,
  CI = 0.95,
  newdata = NULL,
  ...
)

Arguments

Details

To predict over a different spatial extent to that used in the model, a SpatRaster covering the region to make predictions over is passed to the argument new_data. If this is not given predictions are made over the data used in the fit.

The predict_iid logical flag should be set to TRUE if the results of the iid effect from the model are to be used in the prediction.

For the uncertainty calculations, the number of the realisations and the size of the credible interval to be calculated are given by the arguments N and CI respectively.

Value

An object of class disag_prediction which consists of a list of two objects:

Examples

## Not run: 
predict(fit_result)

## End(Not run)

Function to predict mean from the model result

Description

predict_model function takes a disag_model object created by disaggregation::disag_model and predicts mean maps.

Usage

predict_model(
  model_output,
  new_data = NULL,
  predict_iid = FALSE,
  newdata = NULL
)

Arguments

Details

Function returns rasters of the mean predictions as well as the covariate and field contributions to the linear predictor.

To predict over a different spatial extent to that used in the model, a SpatRaster covering the region to make predictions over is passed to the argument new_data. If this is not given predictions are made over the data used in the fit.

The predict_iid logical flag should be set to TRUE if the results of the iid effect from the model are to be used in the prediction.

Value

The mean prediction, which is a list of:

• prediction Raster of mean predictions based.

• field Raster of the field component of the linear predictor.

• iid Raster of the iid component of the linear predictor.

• covariates Raster of the covariate component of the linear predictor.

Examples

## Not run: 
predict_model(result)

## End(Not run)

Function to predict uncertainty from the model result

Description

predict_uncertainty function takes a disag_model object created by disaggregation::disag_model and predicts upper and lower credible interval maps.

Usage

predict_uncertainty(
  model_output,
  new_data = NULL,
  predict_iid = FALSE,
  N = 100,
  CI = 0.95,
  newdata = NULL
)

Arguments

Details

Function returns a SpatRaster of the realisations as well as the upper and lower credible interval rasters.

To predict over a different spatial extent to that used in the model, a SpatRaster covering the region to make predictions over is passed to the argument new_data. If this is not given predictions are made over the data used in the fit.

The predict_iid logical flag should be set to TRUE if the results of the iid effect from the model are to be used in the prediction.

The number of the realisations and the size of the credible interval to be calculated. are given by the arguments N and CI respectively.

Value

The uncertainty prediction, which is a list of:

• realisations SpatRaster of realisations of predictions. Number of realisations defined by argument N.

• predictions_ci SpatRaster of the upper and lower credible intervals. Defined by argument CI.

Examples

## Not run: 
predict_uncertainty(result)

## End(Not run)

Prepare data for disaggregation modelling

Description

prepare_data function is used to extract all the data required for fitting a disaggregation model. Designed to be used in the disaggregation::disag_model function.

Usage

prepare_data(
  polygon_shapefile,
  covariate_rasters,
  aggregation_raster = NULL,
  id_var = "area_id",
  response_var = "response",
  sample_size_var = NULL,
  mesh_args = NULL,
  na_action = FALSE,
  make_mesh = TRUE,
  mesh.args = NULL,
  na.action = NULL,
  makeMesh = NULL,
  ncores = NULL
)

Arguments

Details

Takes a sf object with the response data and a SpatRaster of covariates.

Extract the values of the covariates (as well as the aggregation raster, if given) at each pixel within the polygons (parallelExtract function). This is done in parallel and n.cores argument is used to set the number of cores to use for covariate extraction. This can be the number of covariates used in the model.

The aggregation raster defines how the pixels within each polygon are aggregated. The disaggregation model performs a weighted sum of the pixel prediction, weighted by the pixel values in the aggregation raster. For disease incidence rate you use the population raster to aggregate pixel incidence rate by summing the number of cases (rate weighted by population). If no aggregation raster is provided a uniform distribution is assumed, i.e. the pixel predictions are aggregated to polygon level by summing the pixel values.

Makes a matrix that contains the start and end pixel index for each polygon. Builds an INLA mesh to use for the spatial field (getStartendindex function).

The mesh.args argument allows you to supply a list of INLA mesh parameters to control the mesh used for the spatial field (build_mesh function).

The na.action flag is automatically off. If there are any NAs in the response or covariate data within the polygons the prepare_data method will error. Ideally the NAs in the data would be dealt with beforehand, however, setting na.action = TRUE will automatically deal with NAs. It removes any polygons that have NAs as a response, sets any aggregation pixels with NA to zero and sets covariate NAs pixels to the median value for the that covariate.

Value

A list is returned of class disag_data. The functions summary, print and plot can be used on disag_data. The list of class disag_data contains:

Examples

polygons <- list()
for(i in 1:100) {
  row <- ceiling(i/10)
  col <- ifelse(i %% 10 != 0, i %% 10, 10)
  xmin = 2*(col - 1); xmax = 2*col; ymin = 2*(row - 1); ymax = 2*row
  polygons[[i]] <- list(cbind(c(xmin, xmax, xmax, xmin, xmin),
                              c(ymax, ymax, ymin, ymin, ymax)))
}

polys <- lapply(polygons,sf::st_polygon)
response_df <- data.frame(area_id = 1:100, response = runif(100, min = 0, max = 10))
spdf <- sf::st_sf(response_df,geometry=polys)

r <- terra::rast(nrow=20,ncol=20)
terra::ext(r) <- terra::ext(spdf)
r[] <- sapply(1:terra::ncell(r), function(x) rnorm(1, ifelse(x %% 20 != 0, x %% 20, 20), 3))

r2 <- terra::rast(nrow=20,ncol=20)
terra::ext(r2) <- terra::ext(spdf)
r2[] <- sapply(1:terra::ncell(r), function(x) rnorm(1, ceiling(x/10), 3))
cov_rasters <- c(r, r2)

test_data <- prepare_data(polygon_shapefile = spdf,
                          covariate_rasters = cov_rasters)

Print function for disaggregation input data

Description

Function that prints the input data from the disaggregation model.

Usage

## S3 method for class 'disag_data'
print(x, ...)

Arguments

Details

Prints the number of polyons and pixels, the number of pixels in the largest and smallest polygons and summaries of the covariates.

Print function for disaggregation fit result.

Description

Function that prints the result of the fit from the disaggregation model.

Usage

## S3 method for class 'disag_model'
print(x, ...)

Arguments

Details

Prints the negative log likelihood, model parameters and calculates metrics from in-sample performance.

Print function for disaggregation prediction

Description

Function that prints the prediction from the disaggregation model.

Usage

## S3 method for class 'disag_prediction'
print(x, ...)

Arguments

Details

Prints the number of polyons and pixels, the number of pixels in the largest and smallest polygons and summaries of the covariates.

Summary function for disaggregation input data

Description

Function that summarizes the input data from the disaggregation model.

Usage

## S3 method for class 'disag_data'
summary(object, ...)

Arguments

Details

Prints the number of polyons and pixels, the number of pixels in the largest and smallest polygons and summaries of the covariates.

Value

A list of the number of polyons, the number of covariates and summaries of the covariates.

Summary function for disaggregation fit result

Description

Function that summarises the result of the fit from the disaggregation model.

Usage

## S3 method for class 'disag_model'
summary(object, ...)

Arguments

Details

Prints the negative log likelihood, model parameters and calculates metrics from in-sample performance.

Value

A list of the model parameters, negative log likelihood and metrics from in-sample performance.

Summary function for disaggregation prediction

Description

Function that summarizes the prediction from the disaggregation model.

Usage

## S3 method for class 'disag_prediction'
summary(object, ...)

Arguments

Details

Prints the number of polyons and pixels, the number of pixels in the largest and smallest polygons and summaries of the covariates.

Value

A list of the number of polyons, the number of covariates and summaries of the covariates.