Package 'AMISforInfectiousDiseases'

Title:	Implement the AMIS Algorithm for Infectious Disease Models
Description:	Implements the Adaptive Multiple Importance Sampling (AMIS) algorithm, as described by Retkute et al. (2021, <doi:10.1214/21-AOAS1486>), to estimate key epidemiological parameters by combining outputs from a geostatistical model of infectious diseases (such as prevalence, incidence, or relative risk) with a disease transmission model. Utilising the resulting posterior distributions, the package enables forward projections at the local level.
Authors:	Evandro Konzen [aut] , Renata Retkute [aut] , Raiha Browning [aut] , Thilbault Lestang [aut], Simon Spencer [aut, cre] , University of Warwick [cph], Oxford Research Software Engineering [cph]
Maintainer:	Simon Spencer <[email protected]>
License:	MIT + file LICENSE
Version:	0.1.1
Built:	2025-01-21 08:47:14 UTC
Source:	https://github.com/drsimonspencer/amisforinfectiousdiseases-dev

Help Index

Run the AMIS algorithm to fit a transmission model to a map
Calculate summaries of weighted statistics
Produce list containing the default AMIS parameters
Wrapper function for plot.Mclust()
Plot histogram or credible interval of weighted distributions given a model fitted by amis()
Print method for object of class amis
Sample parameters from their weighted distributions given a model fitted by amis()
Summary method for object of class amis

Run the AMIS algorithm to fit a transmission model to a map

Description

This implements the AMIS algorithm as described in Retkute et al. (2021). Each iteration of the algorithm produces a set of parameters from a proposal distribution (the prior in the first iteration). For each parameter set, a simulation is run from the transmission model. Then, each preceding simulation is weighted at each location according to the distribution of prevalences (or likelihood function) at that location. A Gaussian mixture model is then fitted to the parameter samples with weights averaged over the active locations (ie locations that have yet to reach the effective sample size target). This Gaussian mixture informs the proposal for the next iteration. The algorithm continues until every location has reached the ESS target, or the maximum number of iterations is reached.

Usage

amis(
  prevalence_map,
  transmission_model,
  prior,
  amis_params = default_amis_params(),
  seed = NULL,
  output_dir = NULL,
  initial_amis_vals = NULL
)
amis(
  prevalence_map,
  transmission_model,
  prior,
  amis_params = default_amis_params(),
  seed = NULL,
  output_dir = NULL,
  initial_amis_vals = NULL
)

Arguments

`prevalence_map`	For a single timepoint, `"prevalence_map"` can be an $L \times M$ matrix or data frame containing samples from a geostatistical model, where $L$ is the number of locations and $M$ the number of samples per location. If there are multiple timepoints and/or a parametric `likelihood` function is to be used, `"prevalence_map"` must be a list with $T$ elements, one for each timepoint $t=1,\dots,T$ . Each element must itself be a list with the following objects: `data` An $L \times M$ matrix as above `likelihood` (optional) A function taking arguments: `data`: A vector of length $M_l$ , where $M_l$ is the number of samples from a geostatistical model for location $l$ or the number of likelihood parameters; `sim_prev`: A numeric value for a prevalence simulated from the transmission model; `log`: Logical indicating if calculations are to be performed on log scale (specified in `"amis_params"`, see below). The function `likelihood` must return a numeric value representing the (log-)likelihood of observing a simulated prevalence given the data from a particular location. The location names are inherited from `rownames(prevalence_map)` if `"prevalence_map"` is a matrix, and from `rownames(prevalence_map[[1]]$data)` if `"prevalence_map"` is a list. If `likelihood` is not specified, then it is assumed that the data consist of samples from a geostatistical model and a nonparametric method is used. The nonparametric method to be used is specified in `"amis_params"` using the options `breaks`, `delta`, or `sigma` (see `"amis_params"`).
`transmission_model`	A function taking arguments: `seeds`: a vector of $n$ seeds; `params`: an $n \times d$ matrix of parameter vectors; `n_tims`: number of time points. This function must return an $n \times T$ matrix of prevalences (it must be a matrix even when $T=1$ ). The vector `seeds` will be the vector of indexes of the simulated samples. If `n_samples` new samples are drawn within each iteration of the AMIS algorithm, then the vector `seeds` will be `1:n_samples` at the first iteration, `(n_samples+1):(2*n_samples)` at the second iteration, and so on.
`prior`	A list containing the functions `dprior` and `rprior` (density and random number generator, respectively). The two arguments of `dprior` must be: a $d$ -length vector of transmission model parameters; and a logical `log` to indicate whether to calculate log-density or not. The only argument of `rprior` must be a single integer $n$ that determines the number of samples to draw. `rprior` must produce an $n \times d$ matrix of parameters even when $d=1$ . Parameter names are inherited from the `colnames` of the output of `rprior`.
`amis_params`	A list containing control parameters for the AMIS algorithm (`default_amis_params()` returns the default values): `n_samples` Number of new samples drawn within each AMIS iteration. Default to `500`. `target_ess` Target effective sample size. Default to `500`. `max_iters` Maximum number of AMIS iterations. Default to `12`. `boundaries` A vector of length two with the left and right boundaries for prevalences. Default to `c(0,1)`. If, for instance, left boundary is zero and there is no right boundary, set `boundaries = c(0,Inf)`. `boundaries_param` If specified, it should be a $d \times 2$ matrix with the lower and upper boundaries for the $d$ transmission model parameters. Default to `NULL`. `log` Logical indicating if calculations are to be performed on log scale. Default to `TRUE`. `delete_induced_prior` Logical indicating whether the induced prior density is to be deleted in the update of weights. Default to `FALSE`. `mixture_samples` Number of samples used to represent the weighted parameters in the mixture fitting. `df` Degrees of freedom in the $t$ -distributions, used to yield a heavy tailed proposal. Default to `3`. `q` Parameter (between 0 and 1) controlling how the weights are calculated for active locations. Default to `0`. See Details below. `delta` Optional smoothing parameter if uniform kernel (default) is used. Default to `0.01`. `sigma` Optional smoothing parameter if Gaussian kernel is used. Default to `NULL`. `breaks` Optional vector specifying the breaks for the histogram. Default to `NULL`. For finite `boundaries`, the first and last entries of `breaks` must be equal to the left and right boundaries, respectively. For non-finite `boundaries`, ensure that the range of `breaks` includes any possible prevalence value. Uniform kernel is the default method for the density estimator of the likelihood. If `sigma` is provided, then Gaussian kernel will be used instead. If `breaks` is provided, then histogram-based method will be the nonparametric method being used. Note that if `likelihood` is provided in `prevalence_map`, then a parametric method will be implemented.
`seed`	Optional single value interpreted as an integer. It is the seed for the random number generator for the AMIS algorithm. This is not the same as the `seeds` argument passed to `"transmission_model"`.
`output_dir`	A string specifying the local directory where to save outputs after each iteration of the algorithm. At the end of the string, use the correct path separator for your machine's operating system. If the directory is specified, the outputs will be saved in a file called ‘amis_output.rds’. Default to `NULL` (i.e. outputs are not saved in a local directory).
`initial_amis_vals`	Optional list of intermittent outputs from a previous run (where at least one iteration was successful). These outputs can be saved by specifying the directory `"output_dir"`.

Details

The average weight of parameter vectors for the set of active locations at iteration $i$ $\left(A_i\right)$ has weights determined by how far the effective sample size for location $l$ $\left(\text{ESS}_l^i\right)$ is from the target $\left(\text{ESS}^R\right)$ :

$\bar{w}_j^i = \frac{\sum_{l \in A_i} \left(\text{ESS}^R-\text{ESS}_l^i\right)^q \hat{w}_{lj}^i }{ \sum_{l \in A_i} \left(\text{ESS}^R-\text{ESS}_l^i\right)^q} , \qquad q \in [0,1].$

If $q=0$ (default), the simple average of individual weights will be calculated. If $q>0$ , more weight will be assigned to locations with low ESS.

Value

A list of class amis. If the algorithm completed $I$ iterations, it simulated a total of $N = I \times$ n_samples, and therefore the list returned by amis() will contain:

seeds

An $N$ -length vector with the simulation seeds that were used.

param

An $N \times d$ matrix with the $d$ -dimensional transmission model parameters simulated by the algorithm.

simulated_prevalences

An $N \times T$ matrix with the simulated prevalences, where $T$ is the number of timepoints.

weight_matrix

An $N \times L$ , where $L$ is the number of locations.

likelihoods

A $T \times L \times N$ array with the likelihood of observing a simulated prevalence in each location at each time.

ess

An $L$ -length vector with the final effective sample size (ESS) for each location.

prevalence_map

List with the prevalence map supplied by the user.

locations_with_no_data

Vector indicating which locations have no data at any time point.

components

A list of the mixture components of all iterations, containing:

G: number of components in each iteration;
probs: the mixture weights;
Mean: the locations of the components;
Sigma: the covariance matrices of the components.

components_per_iteration

A list with the mixture components at each iteration. This object is used in plot_mixture_components().

ess_per_iteration

An $L \times I$ matrix with with the ESS for each location after each iteration.

prior_density

An $N$ -length vector with the density function evaluated at the simulated parameter values.

amis_params

List supplied by the user.

evidence

A list containing an estimate of the log model evidence and corresponding log variance of this estimate for both the full likelihood model (product over all locations), and for each location individually.

References

Retkute, R., Touloupou, P., Basanez, M. G., Hollingsworth, T. D., Spencer, S. E. (2021). Integrating geostatistical maps and infectious disease transmission models using adaptive multiple importance sampling. The Annals of Applied Statistics, 15(4), 1980-1998. doi:10.1214/21-AOAS1486.

Examples

# Define simple "transmission" model where prevalence equals first parameter
transmission_model_identity <- function(seeds, parameters, n_tims=1) {
  return(matrix(parameters[,1], ncol=1))
}

# Generate samples for prevalence map with 3 locations given by B(2,1), B(1,1)=Uniform, B(1,2). 
set.seed(123)
L <- 3    # Number of locations
M <- 500 # Number of map samples
prevalence_map <- matrix(NA, L, M)
for (l in 1:L) {
  prevalence_map[l,] <- rbeta(M, max(1,l-1), max(1,3-l))
}
rownames(prevalence_map) <- c("Here","There","Somewhere else")

# Define 2D exponential prior
rprior <- function(n) {
  params <- matrix(NA, n, 2)
  colnames(params) <- c("a","b")
  params[,1] <- rexp(n)
  params[,2] <- rexp(n)
  return(params)
}
dprior <- function(x, log=FALSE) {
  if (log) {
    return(sum(dexp(x, log=TRUE)))
  } else {
    return(prod(dexp(x)))
  }
}
prior <- list(rprior=rprior,dprior=dprior)

# Run AMIS with default control parameters
amis_params <- default_amis_params()
output <- amis(prevalence_map, transmission_model_identity, prior, amis_params, seed=1)

print(output)
summary(output)

original_par <- par(no.readonly = TRUE)
par(cex.lab=1.5, cex.main=1.5, mar=c(5,4.5,4,2)+0.1)

par(mfrow=c(1,2))
plot_mixture_components(output, what = "uncertainty", cex=3)
plot_mixture_components(output, what = "density", nlevels=200)

par(mfrow=c(3,3))
plot(output, what = "a", type="hist", locations=1:L, breaks=100)
plot(output, what = "b", type="hist", locations=1:L, breaks=100)
plot(output, what = "prev", type="hist", locations=1:L, time=1, breaks=100)

par(mar=c(5,7.5,4,2)+0.1)
par(mfrow=c(1,3))
plot(output, what=c("a","b","prev"), type="CI", locations=1:L, ylab=NA,
     cex=3, lwd=3, measure_central="median", display_location_names=TRUE)

calculate_summaries(output, what="prev", locations=1:L, alpha=0.05)

# Generate new samples from the weighted posterior distributions
new_samples <- sample_parameters(output, n_samples = 200, locations = "Here")
head(new_samples)
plot_hist <- function(column_name){
  hist(new_samples[, column_name], xlab=column_name, main=paste("Histogram of", column_name))
}
par(mfrow=c(1,3))
plot_hist("a")
plot_hist("b")
plot_hist("prevalence")

par(original_par)
# Define simple "transmission" model where prevalence equals first parameter
transmission_model_identity <- function(seeds, parameters, n_tims=1) {
  return(matrix(parameters[,1], ncol=1))
}

# Generate samples for prevalence map with 3 locations given by B(2,1), B(1,1)=Uniform, B(1,2). 
set.seed(123)
L <- 3    # Number of locations
M <- 500 # Number of map samples
prevalence_map <- matrix(NA, L, M)
for (l in 1:L) {
  prevalence_map[l,] <- rbeta(M, max(1,l-1), max(1,3-l))
}
rownames(prevalence_map) <- c("Here","There","Somewhere else")

# Define 2D exponential prior
rprior <- function(n) {
  params <- matrix(NA, n, 2)
  colnames(params) <- c("a","b")
  params[,1] <- rexp(n)
  params[,2] <- rexp(n)
  return(params)
}
dprior <- function(x, log=FALSE) {
  if (log) {
    return(sum(dexp(x, log=TRUE)))
  } else {
    return(prod(dexp(x)))
  }
}
prior <- list(rprior=rprior,dprior=dprior)

# Run AMIS with default control parameters
amis_params <- default_amis_params()
output <- amis(prevalence_map, transmission_model_identity, prior, amis_params, seed=1)

print(output)
summary(output)

original_par <- par(no.readonly = TRUE)
par(cex.lab=1.5, cex.main=1.5, mar=c(5,4.5,4,2)+0.1)

par(mfrow=c(1,2))
plot_mixture_components(output, what = "uncertainty", cex=3)
plot_mixture_components(output, what = "density", nlevels=200)

par(mfrow=c(3,3))
plot(output, what = "a", type="hist", locations=1:L, breaks=100)
plot(output, what = "b", type="hist", locations=1:L, breaks=100)
plot(output, what = "prev", type="hist", locations=1:L, time=1, breaks=100)

par(mar=c(5,7.5,4,2)+0.1)
par(mfrow=c(1,3))
plot(output, what=c("a","b","prev"), type="CI", locations=1:L, ylab=NA,
     cex=3, lwd=3, measure_central="median", display_location_names=TRUE)

calculate_summaries(output, what="prev", locations=1:L, alpha=0.05)

# Generate new samples from the weighted posterior distributions
new_samples <- sample_parameters(output, n_samples = 200, locations = "Here")
head(new_samples)
plot_hist <- function(column_name){
  hist(new_samples[, column_name], xlab=column_name, main=paste("Histogram of", column_name))
}
par(mfrow=c(1,3))
plot_hist("a")
plot_hist("b")
plot_hist("prevalence")

par(original_par)

Calculate summaries of weighted statistics

Description

Calculate summaries of weighted statistics

Usage

calculate_summaries(
  x,
  what = "prev",
  time = 1,
  locations = NULL,
  alpha = 0.05,
  exceedance_prob_threshold = 0.35
)
calculate_summaries(
  x,
  what = "prev",
  time = 1,
  locations = NULL,
  alpha = 0.05,
  exceedance_prob_threshold = 0.35
)

Arguments

`x`	The output from the function `amis()`.
`what`	What statistic should be calculated the summaries from. It must be either `"prev"` or the name of one of the model parameters. Default to `"prev"`.
`time`	Time point. Only used if `"what"` is set to `"prev"`.
`locations`	Integer vector or location names identifying locations where summaries should be calculated for. If not specified, summary statistics of all locations will be provided.
`alpha`	Numeric value between 0 and 1. Calculations are for the `(alpha/2, 1-alpha/2)`% quantiles.
`exceedance_prob_threshold`	Numeric value. Default to `0.35`, i.e. the probability that the statistic of interest (e.g. prevalence) is higher than `0.35`.

Details

For illustrative examples, see amis().

Value

A list with mean, median, and quantiles of the weighted distribution.

Produce list containing the default AMIS parameters

Description

For description of AMIS parameters, see argument amis_params in amis().

Usage

default_amis_params()
default_amis_params()

Value

List containing the default AMIS parameters.

Wrapper function for `plot.Mclust()`

Description

Wrapper function for plot.Mclust()

Usage

plot_mixture_components(
  x,
  what = "uncertainty",
  iteration = NULL,
  datapoints = "proposed",
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  ...
)
plot_mixture_components(
  x,
  what = "uncertainty",
  iteration = NULL,
  datapoints = "proposed",
  main = NULL,
  xlim = NULL,
  ylim = NULL,
  ...
)

Arguments

`x`	The output from the function `amis()`.
`what`	A string specifying the type of plot requested: `"uncertainty"` A plot of classification uncertainty (default) `"density"` A plot of estimated density `"BIC"` A plot showing BIC values used to choose the number of components
`iteration`	Integer indicating which iteration the plot should be about. If `NULL` (default), the plot will be for the final iteration. See more details in `plot.Mclust()`.
`datapoints`	A string specifying what the datapoints should represent in the plot of classification uncertainty: `"proposed"` datapoints will represent the samples simulated from the mixture. The colours indicate which mixture components the samples were simulated from. `"fitted"` datapoints will show the samples that the mixture model was fitted to, i.e. weighted samples from the previous iteration. The colour of a datapoint indicates the most likely mixture component the sample belongs to.
`main`	Title of the plot. If `NULL`, the default title will be displayed. Set to `NA` for omitting title.
`xlim`	The x limits of the plots. Default to `NULL`
`ylim`	The y limits of the plots. Default to `NULL`.
`...`	Other arguments to match the `plot.Mclust()` function.

Details

For illustrative examples, see amis().

Value

A plot for model-based clustering results.

Plot histogram or credible interval of weighted distributions given a model fitted by `amis()`

Description

Plot histogram or credible interval of weighted distributions given a model fitted by amis()

Usage

## S3 method for class 'amis'
plot(
  x,
  what = "prev",
  type = "hist",
  locations = 1,
  time = 1,
  measure_central = "mean",
  order_locations_by = NULL,
  display_location_names = FALSE,
  alpha = 0.05,
  breaks = 500,
  cex = 1,
  lwd = 1,
  xlim = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ...
)
## S3 method for class 'amis'
plot(
  x,
  what = "prev",
  type = "hist",
  locations = 1,
  time = 1,
  measure_central = "mean",
  order_locations_by = NULL,
  display_location_names = FALSE,
  alpha = 0.05,
  breaks = 500,
  cex = 1,
  lwd = 1,
  xlim = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ...
)

Arguments

`x`	The output from the function `amis()`.
`what`	What posterior distribution should be plotted. It can be `"prev"` (default) for plotting prevalences, or one of the parameter names.
`type`	Type of plot. It can be `"hist"` (default) for histogram, or `"CI"` for credible intervals
`locations`	Integer vector or location names identifying locations the plots are made for. Default to `1` (first location).
`time`	Integer index identifying the timepoint. Default to `1`.
`measure_central`	Measure of central tendency for credible interval plots. It can be `"mean"` (default) or `"median"`.
`order_locations_by`	How the credible intervals of multiple locations should be ordered. If `NULL` (default), locations are displayed according to the argument `"locations"`. Otherwise, it must be either `"prev"` or one of the parameter names, and then the locations are ranked by the corresponding measure of central tendency.
`display_location_names`	Logical indicating whether location names are to be shown or not in credible interval plots. Default to `FALSE`.
`alpha`	Numeric value between 0 and 1 indicating the endpoints of the credible intervals, which are evaluated at (alpha/2, 1-alpha/2)% quantiles. Default (`0.05`) will create 95% credible intervals.
`breaks`	Argument passed to `wtd.hist()` for histogram plots. Default to `500`.
`cex`	Argument passed to plots of credible intervals. Default to `1`.
`lwd`	Argument passed to plots of credible intervals. Default to `1`.
`xlim`	The x limits of the plots. For for credible intervals of multiple statistics (i.e. `length(what)>1)`, it must be either `NULL` or a list with the x limits for each statistic. Default to `NULL`.
`main`	Title for the plot.
`xlab`	Lable for the x axis.
`ylab`	Lable for the y axis.
`...`	Other graphical parameters passed to `wtd.hist()`.

Details

For illustrative examples, see amis().

Value

A plot.

Print method for object of class `amis`

Description

Print method for object of class amis

Usage

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

Arguments

`x`	The output from the function `amis()`.
`...`	Other arguments to match the generic `print()` function

Details

For illustrative examples, see amis().

Value

Brief description of data and model specifications used to run amis().

Sample parameters from their weighted distributions given a model fitted by `amis()`

Description

Sample parameters from their weighted distributions given a model fitted by amis()

Usage

sample_parameters(x, n_samples = 200, locations = 1)
sample_parameters(x, n_samples = 200, locations = 1)

Arguments

`x`	The output from the function `amis()`.
`n_samples`	Number of samples to draw. Default to `200`.
`locations`	Integer identifying the locations. Default to `1`.

Details

For illustrative examples, see amis().

Value

Matrix with parameter values and corresponding prevalences for each location.

Summary method for object of class `amis`

Description

Summary method for object of class amis

Usage

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

Arguments

`object`	The output from the function `amis()`.
`...`	Other arguments to match the generic `summary()` function

Details

For illustrative examples, see amis().

Value

Summary statistics of the fitted model.

Package 'AMISforInfectiousDiseases'

Help Index

Run the AMIS algorithm to fit a transmission model to a map

Description

Usage

Arguments

Details

Value

References

Examples

Calculate summaries of weighted statistics

Description

Usage

Arguments

Details

Value

Produce list containing the default AMIS parameters

Description

Usage

Value

Wrapper function for plot.Mclust()

Description

Usage

Arguments

Details

Value

Plot histogram or credible interval of weighted distributions given a model fitted by amis()

Description

Usage

Arguments

Details

Value

Print method for object of class amis

Description

Usage

Arguments

Details

Value

Sample parameters from their weighted distributions given a model fitted by amis()

Description

Usage

Arguments

Details

Value

Summary method for object of class amis

Description

Usage

Arguments

Details

Value

Wrapper function for `plot.Mclust()`

Plot histogram or credible interval of weighted distributions given a model fitted by `amis()`

Print method for object of class `amis`

Sample parameters from their weighted distributions given a model fitted by `amis()`

Summary method for object of class `amis`