Package 'spatPomp'

Description

The spatPomp package provides facilities for inference on panel data using spatiotemporal partially-observed Markov process (SpatPOMP) models. To do so, it relies on and extends a number of facilities that the pomp package provides for inference on time series data using partially-observed Markov process (POMP) models.

The spatPomp package concerns models consisting of a collection of interacting units. The methods in spatPomp may be applicable whether or not these units correspond to spatial locations.

Data analysis using spatPomp

The first step in using spatPomp is to encode one's model(s) and data in objects of class spatPomp. This can be done via a call to the spatPomp constructor function.

Extending the pomp platform for developing inference tools

spatPomp extends to panel data the general interface to the components of POMP models provided by pomp. In doing so, it contributes to the goal of the pomp project of facilitating the development of new algorithms in an environment where they can be tested and compared on a growing body of models and datasets.

Documentation

spatPomp is described by Asfaw et al. (2020)

License

spatPomp is provided under the MIT License.

Author(s)

Kidus Asfaw, Joonha Park, Allister Ho, Edward Ionides, Aaron A. King

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

See Also

pomp

Description

An algorithm for estimating the likelihood of a spatiotemporal partially-observed Markov process model. Running abf causes the algorithm to run bootstrap replicate jobs which each yield an imperfect adapted simulation. Simulating from the "adapted filter" distribution runs into a curse of dimensionality (COD) problem, which is mitigated by keeping particles in each replicate close to each other through resampling down to one particle per replicate at each observation time point. The adapted simulations are then weighted in a way that mitigates COD by making a weak coupling assumption to get an approximate filter distribution. As a by-product, we also get an estimate of the likelihood of the data.

Usage

## S4 method for signature 'spatPomp'
abf(
  object,
  Nrep,
  Np,
  nbhd,
  params,
  tol = 1e-100,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'abfd_spatPomp'
abf(
  object,
  Nrep,
  Np,
  nbhd,
  tol,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

Value

Upon successful completion, abf() returns an object of class ‘abfd_spatPomp’ containing the algorithmic parameters used to run abf() and the estimated likelihood.

Methods

The following methods are available for such an object:

logLik

yields an estimate of the log-likelihood of the data under the model.

Author(s)

Kidus Asfaw

References

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

See Also

likelihood maximization algorithms: ienkf(), igirf, iubf, ibpf

Other likelihood evaluation algorithms: abfir(), bpfilter(), enkf(), girf()

Examples

# Complete examples are provided in the package tests
## Not run: 
# Create a simulation of a Brownian motion
b <- bm(U=2, N=5)

# Create a neighborhood function mapping a point in space-time
# to a list of neighboring points in space-time
bm_nbhd <- function(object, time, unit) {
  nbhd_list = list()
  if(time > 1 && unit > 1){
    nbhd_list = c(nbhd_list, list(c(unit-1, time-1)))
  }
  return(nbhd_list)
}

# Run ABF specified number of Monte Carlo replicates and particles per replicate
abfd_bm <- abf(b, Nrep=2, Np=10, nbhd=bm_nbhd)

# Get the likelihood estimate from ABF
logLik(abfd_bm)

## End(Not run)

Description

An algorithm for estimating the filter distribution and likelihood of a spatiotemporal partially-observed Markov process model. Running abfir causes the algorithm to run Monte Carlo replicated jobs which each carry out an adapted simulation using intermediate resampling. Adapted simulation is an easier task than filtering, since particles in each replicate remain close to each other. Intermediate resampling further assists against the curse of dimensionality (COD) problem for importance sampling. The adapted simulations are then weighted in a way that mitigates COD by making a weak coupling assumption to get an approximate filter distribution. As a by-product, we also get an approximation to the likelihood of the data.

Usage

## S4 method for signature 'spatPomp'
abfir(
  object,
  Np,
  Nrep,
  nbhd,
  Ninter,
  tol = (1e-100),
  params,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'abfird_spatPomp'
abfir(object, Np, Nrep, nbhd, Ninter, tol, params, ...)

Arguments

Value

Upon successful completion, abfir() returns an object of class ‘abfird_spatPomp’ containing the algorithmic parameters used to run abfir() and the estimated likelihood.

Methods

The following methods are available for such an object:

logLik

yields a biased estimate of the log-likelihood of the data under the model.

Author(s)

Kidus Asfaw

References

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

See Also

likelihood maximization algorithms: ienkf(), igirf(), iubf(), ibpf()

Other likelihood evaluation algorithms: abf(), bpfilter(), enkf(), girf()

Examples

# Complete examples are provided in the package tests
## Not run: 
# Create a simulation of a Brownian motion
b <- bm(U=2, N=5)

# Create a neighborhood function mapping a point in space-time
# to a list of ``neighboring points" in space-time
bm_nbhd <- function(object, time, unit) {
  nbhd_list = list()
  if(time > 1 && unit > 1){
    nbhd_list = c(nbhd_list, list(c(unit-1, time-1)))
  }
  return(nbhd_list)
}
# Run ABFIR with specified number of Monte Carlo replicates and particles
# per replicate
abfird_bm <- abfir(b,
                  Nrep = 2,
                  Np=10,
                  nbhd = bm_nbhd,
                  Ninter = length(unit_names(b)))
# Get the likelihood estimate from ABFIR
logLik(abfird_bm)

## End(Not run)

Description

Fits independent log-ARMA models for each unit, and calculates the conditional log-likelihood for each observation, as well as log-likelihood for each unit and total log-likelihood. A simple tool, but one with practical applicability, as demonstrated by King et al (2008) and Wheeler et al (2023). This function is designed for non-negative data, and adds 1 to each observation to avoid log(0).

Usage

arma_benchmark(spo, order = c(2, 0, 1))

Arguments

Author(s)

Edward L. Ionides

References

King, A. A., Ionides, E. L., Pascual, M. and Bouma, M. J. (2008). Inapparent infections and cholera dynamics. Nature 454 877-880.

Wheeler, J., Rosengart, A. L., Jiang, Z., Tan, K., Treutle, N. and Ionides, E. L. (2023). Informing policy via dynamic models: Cholera in Haiti. arxiv:2301.08979.

See Also

Other utilities: expand_params()

Examples

# Complete examples are provided in the package tests
## Not run: 
m <- he10(U = 5)
arma_benchmark(m)

## End(Not run)

Description

spatPomp objects can be recast as data frames.

Usage

## S3 method for class 'spatPomp'
as.data.frame(x, ...)

Arguments

Details

When object is a simple ‘spatPomp’ object, as(object,"data.frame") or as.data.frame(object) results in a data frame with the times, units, observables, states (if known), and interpolated covariates (if any).

Value

A ‘data.frame’ with columns for time, spatial unit and observations.

Description

Generate a class ‘spatPomp’ object representing a U-dimensional Brownian motion with spatial correlation decaying geometrically with distance around a circle. The model is defined in continuous time though in this case an Euler approximation is exact at the evaluation times.

Usage

bm(U = 5, N = 100, delta_t = 0.1)

Arguments

Value

An object of class ‘spatPomp’ representing a simulation from a U-dimensional Brownian motion

Author(s)

Edward L. Ionides

See Also

Other spatPomp model generators: bm2(), gbm(), he10(), lorenz(), measles()

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=4, N=20)
# See all the model specifications of the object
spy(b)
# Examples of methodologies applied to this model
# are provided in the tests directory

## End(Not run)

Description

Computes the exact likelihood for a model constructed using bm, using the Kalman filter. This model is useful for testing methods in a situation where an exact answer is available

Usage

bm_kalman_logLik(bm_object, params = coef(bm_object))

Arguments

Value

A numeric value for the log-likelihood.

Author(s)

Edward L. Ionides

Examples

# Further examples are provided in the tests directory
## Not run: 
b <- bm()
bm_kalman_logLik(b)

## End(Not run)

Description

An extension of bm allowing for shared or unit-specific parameters. Generate a class ‘spatPomp’ object representing a U-dimensional Brownian motion with spatial correlation decaying geometrically with distance around a circle. The model is defined in continuous time though in this case an Euler approximation is exact at the evaluation times.

Usage

bm2(
  U = 5,
  N = 100,
  delta_t = 0.1,
  unit_specific_names = "rho",
  shared_names = NULL,
  unit_params = c(rho = 0.4, sigma = 1, tau = 1, X_0 = 0)
)

Arguments

Value

An object of class ‘spatPomp’ representing a simulation from a U-dimensional Brownian motion

Author(s)

Edward L. Ionides

See Also

Other spatPomp model generators: bm(), gbm(), he10(), lorenz(), measles()

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm2(U=4, N=20,shared_names="rho",unit_specific_names=c("sigma","tau"))
# See all the model specifications of the object
spy(b)
# Examples of methodologies applied to this model
# are provided in the tests directory

## End(Not run)

Description

Computes the exact likelihood for a model constructed using bm2, using the Kalman filter. This model is useful for testing methods for models with unit-specific parameters, or method such as ibpf which require a unit-specific extension of shared parameters.

Usage

bm2_kalman_logLik(bm2_object, params = coef(bm2_object))

Arguments

Value

A numeric value for the log-likelihood.

Author(s)

Edward L. Ionides

Examples

# Further examples are provided in the tests directory
## Not run: 
b <- bm2()
bm2_kalman_logLik(b)

## End(Not run)

Description

An implementation of the block particle filter algorithm of Rebeschini and van Handel (2015), which is used to estimate the filter distribution of a spatiotemporal partially-observed Markov process. bpfilter requires a partition of the spatial units which can be provided by either the block_size or the block_list argument. The elements of the partition are called blocks. We perform resampling for each block independently based on sample weights within the block. Each resampled block only contains latent states for the spatial components within the block which allows for a “cross-pollination" of particles where the highest weighted segments of each particle are more likely to be resampled and get combined with resampled components of other particles. The method mitigates the curse of dimensionality by resampling locally.

Usage

## S4 method for signature 'missing'
bpfilter(object, ...)

## S4 method for signature 'ANY'
bpfilter(object, ...)

## S4 method for signature 'spatPomp'
bpfilter(
  object,
  Np,
  block_size,
  block_list,
  save_states,
  filter_traj,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'bpfilterd_spatPomp'
bpfilter(
  object,
  Np,
  block_size,
  block_list,
  save_states,
  filter_traj,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

Value

Upon successful completion, bpfilter() returns an object of class ‘bpfilterd_spatPomp’ containing the algorithmic parameters used to run bpfilter() and the estimated likelihood.

Details

Only one of block_size or block_list should be specified. If both or neither is provided, an error is triggered.

Methods

The following methods are available for such an object:

logLik

yields an estimate of the log-likelihood of the data under the model.

Author(s)

Kidus Asfaw

References

Rebeschini, P., & Van Handel, R. (2015). Can local particle filters beat the curse of dimensionality?. The Annals of Applied Probability, 25(5), 2809-2866.

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

See Also

likelihood maximization algorithms: ienkf(), igirf(), iubf(), ibpf()

Other likelihood evaluation algorithms: abf(), abfir(), enkf(), girf()

Examples

# Complete examples are provided in the package tests
## Not run: 
# Create a simulation of a Brownian motion
b <- bm(U=4, N=2)

# Run BPF with the specified number of units per block
bpfilterd_b1 <- bpfilter(b, Np = 10, block_size = 2)

# Run BPF with the specified partition
bpfilterd_b2 <- bpfilter(b,
                         Np = 10,
                         block_list = list(c(1,2),c(3,4)) )

# Get a likelihood estimate
logLik(bpfilterd_b2)

## End(Not run)

Description

Population and birth information about cities in England and Wales during the measles pre-vaccine era.

Details

Data includes births and population at bi-weekly observations from 40 cities and towns.

Value

a ‘data.frame’ of the 40 largest cities and towns in the UK and Wales, their latitude, longitude and mean population during the measles pre-vaccine period.

References

Dalziel, Benjamin D. et al. (2016) Persistent chaos of measles epidemics in the prevaccination United States caused by a small change in seasonal transmission patterns. PLoS Computational Biology, 12(2), e1004655. doi:10.5061/dryad.r4q34

See Also

Other datasets: measlesUK

Description

Concatenate two or more ‘pomp’ objects into a list-like ‘listie’.

Usage

## S3 method for class 'SpatPomp'
c(...)

Arguments

Details

concat applied to one or more ‘pomp’ objects or lists of ‘pomp’ objects converts the list into a ‘listie’. In particular, concat(A,B,C) is equivalent to do.call(c,unlist(list(A,B,C))).

Description

dunit_measure dunit_measure evaluates the unit measurement density of a unit's observation given the entire state

Usage

## S4 method for signature 'spatPomp'
dunit_measure(object, y, x, unit, time, params, log = TRUE, ...)

Arguments

Value

A class ‘matrix’ with the unit measurement density for spatial unit unit corresponding to the corresponding measurement in y and states in x.

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=3)
s <- states(b)[,1,drop=FALSE]
rownames(s) -> rn
dim(s) <- c(3,1,1)
dimnames(s) <- list(variable=rn, rep=NULL)
p <- coef(b); names(p) -> rnp
dim(p) <- c(length(p),1); dimnames(p) <- list(param=rnp)
o <- obs(b)[,1,drop=FALSE]
dunit_measure(b, y=o, x=s, unit=1, time=1, params=p)

## End(Not run)

Description

A function to perform filtering using the ensemble Kalman filter of Evensen, G. (1994). This function is generalized to allow for an measurement covariance matrix that varies over time. This is useful if the measurement model varies with the state.

Usage

## S4 method for signature 'spatPomp'
enkf(data, Np, ..., verbose = getOption("spatPomp_verbose", FALSE))

Arguments

Value

An object of class ‘enkfd_spatPomp’ that contains the estimate of the log likelihood (via the loglik attribute), algorithmic parameters used to run enkf(). Also included are estimated filter means, prediction means and forecasts that are generated during an enkf() run.

References

G. Evensen. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans 99, 10143–10162, 1994.

G. Evensen. Data assimilation: the ensemble Kalman filter. Springer-Verlag, 2009.

J.L. Anderson. An Ensemble Adjustment Kalman Filter for Data Assimilation. Monthly Weather Review 129, 2884–2903, 2001.

See Also

ienkf(), igirf, iubf, ibpf

Other likelihood evaluation algorithms: abf(), abfir(), bpfilter(), girf()

Examples

# Complete examples are provided in the package tests
## Not run: 
# Create a simulation of a Brownian motion
b <- bm(U=2, N=5)

# Run EnKF
enkfd_bm <- enkf(b, Np = 20)

# Get a likelihood estimate
logLik(enkfd_bm)

## End(Not run)

Description

eunit_measure evaluates the expectation of a unit's observation given the entire state

Usage

## S4 method for signature 'spatPomp'
eunit_measure(object, x, unit, time, params, Np = 1, log = FALSE)

Arguments

Value

A class ‘matrix’ with the unit expected observation for spatial unit unit corresponding to the corresponding states in x.

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=3)
s <- states(b)[,1,drop=FALSE]
rownames(s) -> rn
dim(s) <- c(3,1,1)
dimnames(s) <- list(variable=rn, rep=NULL)
p <- coef(b); names(p) -> rnp
dim(p) <- c(length(p),1); dimnames(p) <- list(param=rnp)
o <- obs(b)[,1,drop=FALSE]
eunit_measure(b, x=s, unit=2, time=1, params=p)

## End(Not run)

Description

Iterated block particle filters require shared parameters to be expanded into having a value at each unit. expand_params, contract_params and mean_by_unit provide tools for moving between representations. For a unit-specific expansion of a shared parameter, all the values for different units should be the same, and mean_by_unit ensures this by taking an average.

Usage

expand_params(params, expandedParNames, U)

contract_params(params, expandedParNames, U, average = FALSE)

mean_by_unit(params, expandedParNames, U)

Arguments

Details

These functions assume that expanded parameters have names ending in "1" through "U", where U is the number of units. Contracted parameters, meaning any parameter that is not expanded, should have a name ending in "1". This numerical suffix convention is useful for writing model-building code that allows parameters to be either expanded or contracted.

See Also

Other utilities: arma_benchmark()

Description

Generate a spatPomp object representing a U-dimensional geometric Brownian motion with spatial correlation decaying geometrically with distance around a circle. The model is defined in continuous time, but an Euler approximation is used for this numerical implementation.

Usage

gbm(U = 5, N = 100, delta_t = 0.1, IVP_values = 1, delta_obs = 1)

Arguments

Value

An object of class ‘spatPomp’ representing a simulation from a U-dimensional geometric Brownian motion

Author(s)

Kidus Asfaw

References

Asfaw, K. T. (2021). Simulation-based Inference for Partially Observed Markov Process Models with Spatial Coupling. University of Michigan Doctoral dissertation. doi:10.7302/2751

See Also

Other spatPomp model generators: bm(), bm2(), he10(), lorenz(), measles()

Examples

# Complete examples are provided in the package tests
## Not run: 
g <- gbm(U=4, N=20)
# See all the model specifications of the object
spy(g)

## End(Not run)

Description

An implementation of the algorithm of Park and Ionides (2020), following the pseudocode in Asfaw et al. (2020).

Usage

## S4 method for signature 'missing'
girf(object, ...)

## S4 method for signature 'ANY'
girf(object, ...)

## S4 method for signature 'spatPomp'
girf(
  object,
  Np,
  Ninter,
  lookahead = 1,
  Nguide,
  kind = c("bootstrap", "moment"),
  tol,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'girfd_spatPomp'
girf(
  object,
  Np,
  Ninter,
  lookahead,
  Nguide,
  kind = c("bootstrap", "moment"),
  tol,
  ...
)

Arguments

Value

Upon successful completion, girf() returns an object of class ‘girfd_spatPomp’ which contains the algorithmic parameters that were used to run girf() and the resulting log likelihood estimate.

Methods

The following methods are available for such an object:

logLik

yields an unbiased estimate of the log-likelihood of the data under the model.

Author(s)

Kidus Asfaw

References

Park, J. and Ionides, E. L. (2020) Inference on high-dimensional implicit dynamic models using a guided intermediate resampling filter. Statistics and Computing, doi:10.1007/s11222-020-09957-3

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

See Also

likelihood maximization algorithms: ienkf(), igirf(), iubf(), ibpf()

Other likelihood evaluation algorithms: abf(), abfir(), bpfilter(), enkf()

Examples

# Complete examples are provided in the package tests
## Not run: 
#
# Create a simulation of a Brownian motion
b <- bm(U=2, N=5)

# Run GIRF
girfd_bm <- girf(b,
                 Np = 10,
                 Ninter = length(unit_names(b)),
                 lookahead = 1,
                 Nguide = 10
)
# Get the likelihood estimate from GIRF
logLik(girfd_bm)

# Compare with the likelihood estimate from particle filter
pfd_bm <- pfilter(b, Np = 10)
logLik(pfd_bm)

## End(Not run)

Description

Generate a spatPomp object for measles adding spatial coupling to The model and data from He et al. (2010) with gravity transport as in Park and Ionides (2020). Other transport models may be added in future. The data in the object matches He et al. (2010). The model matches that analysis in the specific case where there is no coupling and all parameters are unit-specific.

Usage

he10(
  U = 6,
  dt = 2/365,
  Tmax = 1964,
  expandedParNames = c("alpha", "iota", "R0", "cohort", "amplitude", "gamma", "sigma",
    "sigmaSE", "rho", "psi", "g", "S_0", "E_0", "I_0"),
  basic_params = c(alpha = 1, iota = 0, R0 = 30, cohort = 0, amplitude = 0.5, gamma = 52,
    sigma = 52, mu = 0.02, sigmaSE = 0.15, rho = 0.5, psi = 0.15, g = 400, S_0 = 0.032,
    E_0 = 5e-05, I_0 = 4e-05),
  towns_selected = NULL
)

Arguments

Details

The code for this spatPomp has duplication with measles(), but in future the two models may diverge. The measles() spatPomp is a simplified situation useful for testing some methods. However, measles() does not permit unit-specific parameters, which he10() allows. Also, the structure of this spatPomp is compatible with the spatiotemporal iterated filtering algorithm ibpf(). This requires shared parameters to be represented with a value for each unit, which should be the same for each unit in a valid model instance but may vary between units while optimizing.

Value

An object of class ‘spatPomp’ representing a U-dimensional spatially coupled measles POMP model.

Relationship to published analysis

The model generator he10() differs from measles() in some details necessitated to reproduce the results of He et al (2010). The measles() model follows the decision of Park and Ionides (2020) and Ionides et al (2021) to apply the mixing exponent $\alpha_u$ to $(I_u/P_u)$ rather than just to $I_u$. he10() does this for the infections arising from individuals traveling to another town (which don't arise for the panel model of He et al (2010)). However, for infections arising within a city, in order to reproduce the results of He et al (2010), he10() uses $(I_u^{\alpha_u}/P_u)$. This is not fully documented in the text of Ionides et al (2022). Models fitted to data have $alpha_u$ close to $1$, so this issue may be negligible in practice.

Another discrepancy between the he10() code and the mathematical model written by Ionides et al (2022) arises in whether individuals traveling from $u$ to $v$ use mixing exponent $\alpha_u$ or $\alpha_v$. Ionides et al (2022) wrote $u$ but the code used implemented $v$. The implementation in he10() matches the implementation of Ionides et al (2022) and so uses $v$.

It might seem surprising that immigrant infections affect only the first term in the expression for $\mu_{SE}$ in Ionides et al (2022), and in the corresponding he10() code. This immigration term is needed in the first term to make the model of He et al (2010) a proper sub-model, when coupling is removed by setting the gravitational constant parameter equal to zero. When this constant is allowed to be positive, the role of immigrant infections transmitting to traveling individuals is anticipated to be a negligible, second-order effect which has been omitted from the model.

Note

This function goes through a typical workflow of constructing a typical spatPomp object (1-4 below). This allows the user to have a file that replicates the exercise of model building as well as function that creates a typical nonlinear model in epidemiology in case they want to test a new inference methodology. We purposely do not modularize this function because it is not an operational piece of the package and is instead useful as an example.
1. Getting a measurements data.frame with columns for times, spatial units and measurements.
2. Getting a covariates data.frame with columns for times, spatial units and covariate data.
3. Constructing model components (latent state initializer, latent state transition simulator and measurement model). Depending on the methods used, the user may have to supply a vectorfield to be integrated that represents the deterministic skeleton of the latent process.
4. Bringing all the data and model components together to form a spatPomp object via a call to spatPomp().

Author(s)

Edward L. Ionides

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

Ionides, E. L., Ning, N., and Wheeler, J. (2022). An iterated block particle filter for inference on coupled dynamic systems with shared and unit-specific parameters. Statistica Sinica, to appear. doi:10.48550/arXiv.2206.03837

Park, J. and Ionides, E. L. (2020) Inference on high-dimensional implicit dynamic models using a guided intermediate resampling filter. Statistics and Computing, doi:10.1007/s11222-020-09957-3

See Also

he10coordinates, he10measles, he10mle, he10demography

Other spatPomp model generators: bm(), bm2(), gbm(), lorenz(), measles()

Examples

# Complete examples are provided in the package tests
## Not run: 
m <- he10(U = 5)
# See all the model specifications of the object
spy(m)

## End(Not run)

Description

Longitude and lattitude for the 20 towns in England and Wales studied by He et al (2010).

Value

a ‘data.frame’ of longitude and lattitude for each town.

References

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

See Also

Other datasets he10: he10demography, he10measles, he10mle

Description

Population and birth information for some towns in England and Wales during the measles pre-vaccine era.

Details

Data are annual statistics for the 20 towns analyzed by He et al (2010).

Value

a ‘data.frame’ of with variables town, year, pop and births.

References

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

See Also

Other datasets he10: he10coordinates, he10measles, he10mle

Description

Measles case data from various cities and towns in England and Wales during the pre-vaccine era.

Details

Data are weekly case counts for the 20 towns analyzed by He et al (2010).

Value

a ‘data.frame’ of reported measles cases for 20 towns, analyzed by He et al (2010).

References

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

See Also

Other datasets he10: he10coordinates, he10demography, he10mle

Description

Maximum likelihood estimate for fitting a susceptible-exposed-infected- recovered model to the measles case report data analyzed by He et al (2010). The values are similar, but not identical, to those reported by He et al.

Value

a ‘data.frame’ containing the estimated parameters.

References

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

See Also

Other datasets he10: he10coordinates, he10demography, he10measles

Description

An iterated block particle filter, for both shared and unit-specific parameters. We require that the spatPomp has been constructed to have a unit-specific parameter "thetau" for unit u corresponding to an estimated parameter "theta", whether theta is shared or unit-specific. This permits IBPF to implement a spatiotemporal random walk to estimate theta. We require that rw.sd is positive for, and only for, all parameters of the form "thetau" if "theta" is listed in sharedParNames or unitParNames.

Usage

## S4 method for signature 'missing'
ibpf(data, ...)

## S4 method for signature 'ANY'
ibpf(data, ...)

## S4 method for signature 'spatPomp'
ibpf(
  data,
  Nbpf,
  Np,
  rw.sd,
  sharedParNames,
  unitParNames,
  cooling.type = "geometric",
  cooling.fraction.50,
  block_size,
  block_list,
  spat_regression,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'ibpfd_spatPomp'
ibpf(
  data,
  Nbpf,
  Np,
  rw.sd,
  sharedParNames,
  unitParNames,
  cooling.type = "geometric",
  cooling.fraction.50,
  block_size,
  block_list,
  spat_regression,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'bpfilterd_spatPomp'
ibpf(
  data,
  Nbpf,
  Np,
  rw.sd,
  sharedParNames,
  unitParNames,
  cooling.type = "geometric",
  cooling.fraction.50,
  block_size,
  block_list,
  spat_regression,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

ifelse(time==time[1],s,0).

ifelse(time==time[lag],s,0).

Value

Upon successful completion, ibpf returns an object of class ‘ibpfd_spatPomp’.

Methods

The following methods are available for such an object:

coef

gives the Monte Carlo estimate of the maximum likelihood.

Author(s)

Edward L. Ionides

References

Ionides, E. L., Ning, N., and Wheeler, J. (2022). An iterated block particle filter for inference on coupled dynamic systems with shared and unit-specific parameters. Statistica Sinica, to appear. doi:10.48550/arXiv.2206.03837

See Also

likelihood evaluation algorithms: girf(), enkf(), bpfilter(), abf(), abfir()

Other likelihood maximization algorithms: ienkf(), igirf(), iubf()

Examples

# Complete examples are provided in the package tests
## Not run: 
# Create a simulation of a Brownian motion, for an extended model with
# unit-specific parameters for all estimated, even if the parameter
# takes the same shared value for each unit.
U <- 4
b2 <- bm2(U=U,N=5,unit_specific_names="rho")

# Run ibpf with two blocks of two units each. estimating rho as a
# shared parameter with all other parameters being fixed.
b2_ibpf <- ibpf(b2,
  sharedParNames="rho",
  unitParNames=NULL,
  Nbpf=5,
  spat_regression=0.1,
  Np=50,
  rw.sd=do.call(rw_sd,setNames(rep(list(0.01),times=U),paste0("rho",1:U))),
  cooling.fraction.50=0.5,
  block_size=2
)

# Get a likelihood estimate
logLik(b2_ibpf)

## End(Not run)

Description

An implementation of a parameter estimation algorithm that uses the ensemble Kalman filter (Evensen, G. (1994)) to perform the filtering step in the parameter-perturbed iterated filtering scheme of Ionides et al. (2015) following the pseudocode in Asfaw, et al. (2020).

Usage

## S4 method for signature 'spatPomp'
ienkf(
  data,
  Nenkf = 1,
  rw.sd,
  cooling.type = c("geometric", "hyperbolic"),
  cooling.fraction.50,
  Np,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

ifelse(time==time[1],s,0).

ifelse(time==time[lag],s,0).

Value

Upon successful completion, ienkf returns an object of class ‘ienkfd_spatPomp’. This object contains the convergence record of the iterative algorithm with respect to the likelihood and the parameters of the model (which can be accessed using the traces attribute) as well as a final parameter estimate, which can be accessed using the coef().

Methods

The following methods are available for such an object:

coef

gives the Monte Carlo estimate of the maximum likelihood.

Author(s)

Kidus Asfaw

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

Evensen, G. (1994) Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics Journal of Geophysical Research: Oceans 99:10143–10162

Evensen, G. (2009) Data assimilation: the ensemble Kalman filter Springer-Verlag.

Anderson, J. L. (2001) An Ensemble Adjustment Kalman Filter for Data Assimilation Monthly Weather Review 129:2884–2903

See Also

likelihood evaluation algorithms: girf(), enkf(), bpfilter(), abf(), abfir()

Other likelihood maximization algorithms: ibpf(), igirf(), iubf()

Description

An implementation of a parameter estimation algorithm combining the intermediate resampling scheme of the guided intermediate resampling filter of Park and Ionides (2020) and the parameter perturbation scheme of Ionides et al. (2015) following the pseudocode in Asfaw, et al. (2020).

Usage

## S4 method for signature 'missing'
igirf(data, ...)

## S4 method for signature 'ANY'
igirf(data, ...)

## S4 method for signature 'spatPomp'
igirf(
  data,
  Ngirf,
  Np,
  rw.sd,
  cooling.type,
  cooling.fraction.50,
  Ninter,
  lookahead = 1,
  Nguide,
  kind = c("bootstrap", "moment"),
  tol = 1e-100,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

## S4 method for signature 'igirfd_spatPomp'
igirf(
  data,
  Ngirf,
  Np,
  rw.sd,
  cooling.type,
  cooling.fraction.50,
  Ninter,
  lookahead,
  Nguide,
  kind = c("bootstrap", "moment"),
  tol,
  ...,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

ifelse(time==time[1],s,0).

ifelse(time==time[lag],s,0).

Value

Upon successful completion, igirf() returns an object of class ‘igirfd_spatPomp’. This object contains the convergence record of the iterative algorithm with respect to the likelihood and the parameters of the model (which can be accessed using the traces attribute) as well as a final parameter estimate, which can be accessed using the coef(). The algorithmic parameters used to run igirf() are also included.

Methods

The following methods are available for such an object:

coef

gives the Monte Carlo maximum likelihood parameter estimate.

Author(s)

Kidus Asfaw

References

Park, J. and Ionides, E. L. (2020) Inference on high-dimensional implicit dynamic models using a guided intermediate resampling filter. Statistics and Computing, doi:10.1007/s11222-020-09957-3

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

See Also

likelihood evaluation algorithms: girf(), enkf(), bpfilter(), abf(), abfir()

Other likelihood maximization algorithms: ibpf(), ienkf(), iubf()

Examples

# Complete examples are provided in the package tests
## Not run: 
igirf(bm(U=2,N=4),Ngirf=2,
  rw.sd = rw_sd(rho=0.02,X1_0=ivp(0.02)),
  cooling.type="geometric",cooling.fraction.50=0.5,
  Np=10,Ninter=2,lookahead=1,Nguide=5)

## End(Not run)

Description

An algorithm for estimating the parameters of a spatiotemporal partially-observed Markov process. Running iubf causes the algorithm to perform a specified number of iterations of unadapted simulations with parameter perturbation and parameter resamplings. At each iteration, unadapted simulations are performed on a perturbed version of the model, in which the parameters to be estimated are subjected to random perturbations at each observation. After cycling through the data, each replicate's weight is calculated and is used to rank the bootstrap replictates. The highest ranking replicates are recycled into the next iteration. This extra variability introduced through parameter perturbation effectively smooths the likelihood surface and combats particle depletion by introducing diversity into particle population. As the iterations progress, the magnitude of the perturbations is diminished according to a user-specified cooling schedule.

Usage

## S4 method for signature 'spatPomp'
iubf(
  object,
  Nubf = 1,
  Nrep_per_param,
  Nparam,
  nbhd,
  prop,
  rw.sd,
  cooling.type = c("geometric", "hyperbolic"),
  cooling.fraction.50,
  tol = (1e-18)^17,
  verbose = getOption("spatPomp_verbose"),
  ...
)

Arguments

ifelse(time==time[1],s,0).

ifelse(time==time[lag],s,0).

Value

Upon successful completion, iubf() returns an object of class ‘iubfd_spatPomp’. This object contains the convergence record of the iterative algorithm with respect to the likelihood and the parameters of the model (which can be accessed using the traces attribute) as well as a final parameter estimate, which can be accessed using the coef(). The algorithmic parameters used to run iubf() are also included.

Methods

The following methods are available for such an object:

coef

extracts the point estimate

Author(s)

Kidus Asfaw

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

See Also

likelihood evaluation algorithms: girf(), enkf(), bpfilter(), abf(), abfir()

Other likelihood maximization algorithms: ibpf(), ienkf(), igirf()

Description

Extracts the estimated log likelihood from a fitted model

Usage

## S4 method for signature 'girfd_spatPomp'
logLik(object)

## S4 method for signature 'bpfilterd_spatPomp'
logLik(object)

## S4 method for signature 'abfd_spatPomp'
logLik(object)

## S4 method for signature 'iubfd_spatPomp'
logLik(object)

## S4 method for signature 'abfird_spatPomp'
logLik(object)

## S4 method for signature 'igirfd_spatPomp'
logLik(object)

Arguments

Value

a numeric which is the estimated log likelihood

Description

Generate a spatPomp object representing a U-dimensional stochastic Lorenz '96 process with N measurements made at times $t_n = n * delta_obs$, simulated using an Euler method with time increment delta_t.

Usage

lorenz(
  U = 5,
  N = 100,
  delta_t = 0.01,
  delta_obs = 0.5,
  regular_params = c(F = 8, sigma = 1, tau = 1)
)

Arguments

Value

An object of class ‘spatPomp’ representing a simulation from a U-dimensional Lorenz 96 model

Author(s)

Edward L. Ionides

References

Lorenz, E. N. (1996) Predictability: A problem partly solved. Proceedings of the seminar on predictability

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

See Also

Other spatPomp model generators: bm(), bm2(), gbm(), he10(), measles()

Examples

# Complete examples are provided in the package tests
## Not run: 
l <- lorenz(U=5, N=100, delta_t=0.01, delta_obs=1)
# See all the model specifications of the object
spy(l)

## End(Not run)

Description

Generate a spatPomp object for measles in the top-U most populous cities in England and Wales. The model is adapted from He et al. (2010) with gravity transport following Park and Ionides (2020). The data are from Dalziel et al (2016).

Usage

measles(
  U = 6,
  dt = 2/365,
  fixed_ivps = TRUE,
  S_0 = 0.032,
  E_0 = 5e-05,
  I_0 = 4e-05
)

Arguments

Value

An object of class ‘spatPomp’ representing a U-dimensional spatially coupled measles POMP model.

Relationship to published analysis

This model was used to generate the results of Ionides et al (2021). However, their equation (6) is not exactly correct for the Binomial Gamma infinitesimal model used in the code, as shown by Proposition 5 of Breto and Ionides, 2011. If Poisson Gamma infinitesimal increments were used (Proposition 4 of Breto and Ionides, 2011) then (6) would be correct, but the resulting unbounded increments could break the non-negativity requirement for compartment membership. The same issue arises with the description in Park and Ionides (2020), though that analysis was based on a different model implementation since the spatPomp package was not yet available.

A difference between (6) of Ionides et al (2021) and (2.1) of He et al (2010) is that in (6) the mixing exponent $\alpha$ is applied to $(I_u/P_u)$ rather than just to $I_u$. In the context of He et al (2010) this changes the parameterization but has negligible effect on the model itself since $P_u(t)$ is approximately constant and so changing its power can be compensated by a corresponding change in the transmission rate, $\beta$. In practice, models fitted to data have $alpha$ close to $1$, so this issue may be moot and this modeling mechanism may not be an effective empirical way to carry out the goal of making allowance for heterogeneous mixing.

The code here includes a cohort effect, $c$, following He et al (2010), that was not included by Ionides et al (2021). This effect leads to a non-differentiability of expected increments which is problematic for the spatPomp implementation of GIRF. For the results of Ionides et al (2021), this was set to $c=0$.

The analysis of He et al (2010), and the model generated by he10(), use weekly aggregated cases. Weekly reports were not available beyond the 20 cites studied by He et al (2010) so measles() relies on the biweekly reports used by Ionides et al (2021) and Ionides & Park (2020).

It turns out to be an important detail of the model by He et al (2010) that a delay is included between birth and entry into the susceptible compartment. He et al (2010) found a 4 year delay fits the data. This value is fixed to be the variable birth_delay in the code for measles(). The code for Ionides et al (2021) uses a 3 year delay, and the delay is not explained in the abbreviated model description. In measles() we have reverted to the 4 year delay identified by He et al (2010).

Note

This function goes through a typical workflow of constructing a typical spatPomp object (1-4 below). This allows the user to have a file that replicates the exercise of model building as well as function that creates a typical nonlinear model in epidemiology in case they want to test a new inference methodology. We purposely do not modularize this function because it is not an operational piece of the package and is instead useful as an example.
1. Getting a measurements data.frame with columns for times, spatial units and measurements.
2. Getting a covariates data.frame with columns for times, spatial units and covariate data.
3. Constructing model components (latent state initializer, latent state transition simulator and measurement model). Depending on the methods used, the user may have to supply a vectorfield to be integrated that represents the deterministic skeleton of the latent process.
4. Bringing all the data and model components together to form a spatPomp object via a call to spatPomp().

Author(s)

Edward L. Ionides

References

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

Dalziel, Benjamin D. et al. (2016) Persistent chaos of measles epidemics in the prevaccination United States caused by a small change in seasonal transmission patterns. PLoS Computational Biology, 12(2), e1004655. doi:10.5061/dryad.r4q34

Park, J. and Ionides, E. L. (2020) Inference on high-dimensional implicit dynamic models using a guided intermediate resampling filter. Statistics and Computing, doi:10.1007/s11222-020-09957-3

Breto, C. and Ionides, E.L. (2011) Compound Markov counting processes and their applications to modeling infinitesimally over-dispersed systems. Stochastic Processes and their Applications 121, 2571-2591. doi:10.1016/j.spa.2011.07.005

See Also

measles_UK, city_data_UK

Other spatPomp model generators: bm(), bm2(), gbm(), he10(), lorenz()

Examples

# Complete examples are provided in the package tests
## Not run: 
m <- measles(U = 5)
# See all the model specifications of the object
spy(m)

## End(Not run)

Description

Generate a spatPomp object for measles in the top-U most populous cities in England and Wales. The model is adapted from He et al. (2010) with gravity transport following Park and Ionides (2019). The structure of this spatPomp is designed to accommodate shared and unit-specific parameters. If carrying out spatiotemporal iterated filtering for shared parameters via ibpf, it is necessary to have a unit-specific expansion and so these parameters should be included in expandedParNames. This model and data correspond to the biweekly analysis of Park and Ionides (2020) and Ionides et al (2021). There are small differences with the weekly model and data of He et al (2010) and Ionides, Ning and Wheeler (2022).

Usage

measles2(
  U = 6,
  dt = 2/365,
  N = 391,
  expandedParNames = c("R0", "c", "A", "muIR", "muEI", "sigmaSE", "rho", "psi", "g",
    "S_0", "E_0", "I_0"),
  contractedParNames = NULL,
  simulated = FALSE,
  basic_params = c(alpha = 0.98, iota = 0.1, R0 = 30, c = 0.3, A = 0.5, muIR = 52, muEI =
    52, muD = 0.02, sigmaSE = 0.15, rho = 0.5, psi = 0.15, g = 400, S_0 = 0.032, E_0 =
    5e-05, I_0 = 4e-05)
)

Arguments

Value

An object of class ‘spatPomp’ representing a U-dimensional spatially coupled measles POMP model.

References

He, D., Ionides, E. L., and King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. doi:10.1098/rsif.2009.0151

Park, J. and Ionides, E. L. (2020) Inference on high-dimensional implicit dynamic models using a guided intermediate resampling filter. Statistics and Computing, doi:10.1007/s11222-020-09957-3

Ionides, E. L., Asfaw, K., Park, J., and King, A. A. (2021). Bagged filters for partially observed interacting systems. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1974867

Ionides, E. L., Ning, N., and Wheeler, J. (2022). An iterated block particle filter for inference on coupled dynamic systems with shared and unit-specific parameters. Statistica Sinica, to appear. doi:10.48550/arXiv.2206.03837

Examples

# Complete examples are provided in the package tests
## Not run: 
m <- measles2(U = 5)
# See all the model specifications of the object
spy(m)

## End(Not run)

Description

Measles case data from various cities and towns in England and Wales during the pre-vaccine era.

Details

Data includes bi-weekly case counts as well as births and population from 40 cities and towns.

Value

a ‘data.frame’ of the 40 largest cities and towns in the UK and Wales, their latitude, longitude and bi-weekly measles case counts, population and birthrates.

References

Dalziel, Benjamin D. et al. (2016) Persistent chaos of measles epidemics in the prevaccination United States caused by a small change in seasonal transmission patterns. PLoS Computational Biology, 12(2), e1004655. doi:10.5061/dryad.r4q34

See Also

Other datasets: city_data_UK

Description

munit_measure returns a moment-matched parameter set given an empirically calculated measurement variance and latent states. This is used in girf() and igirf() when they are run with kind='moment'.

Usage

## S4 method for signature 'spatPomp'
munit_measure(object, x, vc, unit, time, params, Np = 1)

Arguments

Value

An array with dimensions dim(array.params)[1] by dim(x)[2] by length(unit) bylength(time) representing the moment-matched parameter set(s) corresponding to the variance of the measurements, vc, and the states, x.

Author(s)

Kidus Asfaw

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=3)
s <- states(b)[,1,drop=FALSE]
rownames(s) -> rn
dim(s) <- c(3,1,1)
dimnames(s) <- list(variable=rn, rep=NULL)
p <- coef(b); names(p) -> rnp
dim(p) <- c(length(p),1); dimnames(p) <- list(param=rnp)
o <- obs(b)[,1,drop=FALSE]
array.params <- array(p,
                      dim = c(length(p),
                              length(unit_names(b)), 1, 1),
                      dimnames = list(params = rownames(p)))
vc <- c(4, 9, 16); dim(vc) <- c(length(vc), 1, 1)
munit_measure(b, x=s, vc=vc, Np=1, unit = 1, time=1, params=array.params)

## End(Not run)

Description

Visualize data in a spatPomp object or a derived class. This gives a quick view; the data can be extracted from the object to make a customized plot.

Usage

## S4 method for signature 'spatPomp'
plot(x, type = c("l", "h"), log = FALSE, plot_unit_names = TRUE, ...)

## S4 method for signature 'igirfd_spatPomp'
plot(x, params = names(coef(x)), ncol = 3)

Arguments

Value

a ggplot plot of class ‘gg’ and ‘ggplot’.

Description

Prints its argument.

Usage

## S4 method for signature 'spatPomp'
print(x)

Arguments

Value

An object of class ‘spatPomp’ is returned *invisibly*. The user is notified on the console only the class of the object.

Note

Use spy() to see model components of x instead.

Description

runit_measure simulates a unit's observation given the entire state

Usage

## S4 method for signature 'spatPomp'
runit_measure(object, x, unit, time, params, log = FALSE)

Arguments

Value

A matrix with the simulated observation corresponding to state x and unit unit with parameter set params.

Author(s)

Kidus Asfaw

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=3)
s <- states(b)[,1,drop=FALSE]
rownames(s) -> rn
dim(s) <- c(3,1,1)
dimnames(s) <- list(variable=rn, rep=NULL)
p <- coef(b); names(p) -> rnp
dim(p) <- c(length(p),1); dimnames(p) <- list(param=rnp)
o <- obs(b)[,1,drop=FALSE]
runit_measure(b, x=s, unit=2, time=1, params=p)

## End(Not run)

Description

simulate generates simulations of the latent and measurement processes.

Usage

## S4 method for signature 'spatPomp'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  format = c("spatPomps", "data.frame"),
  include.data = FALSE,
  ...
)

Arguments

Value

if format='spatPomps' and nsim=1 an object of class ‘spatPomp’ representing a simulation from the model in object is returned. If format='spatPomps' and nsim>1 a list of class ‘spatPomp’ objects is returned. If format='data.frame' then a class ‘data.frame’ object is returned.

Author(s)

Kidus Asfaw

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

Examples

# Complete examples are provided in the package tests
## Not run: 
# Get a spatPomp object
b <- bm(U=2, N=5)
# Get 2 simulations from same model as data.frame
sims <- simulate(b, nsim=2, format='data.frame')

## End(Not run)

Description

This function constructs a class ‘spatPomp’ object, encoding a spatiotemporal partially observed Markov process (SpatPOMP) model together with a uni- or multi-variate time series on a collection of units. Users will typically develop a POMP model for a single unit before embarking on a coupled SpatPOMP analysis. Consequently, we assume some familiarity with pomp and its description by King, Nguyen and Ionides (2016). The spatPomp class inherits from pomp with the additional unit structure being a defining feature of the resulting models and inference algorithms.

Usage

spatPomp(
  data,
  units,
  times,
  covar,
  t0,
  ...,
  eunit_measure,
  munit_measure,
  vunit_measure,
  dunit_measure,
  runit_measure,
  rprocess,
  rmeasure,
  dprocess,
  dmeasure,
  skeleton,
  rinit,
  rprior,
  dprior,
  unit_statenames,
  unit_accumvars,
  shared_covarnames,
  globals,
  paramnames,
  params,
  cdir,
  cfile,
  shlib.args,
  PACKAGE,
  partrans,
  compile = TRUE,
  verbose = getOption("spatPomp_verbose", FALSE)
)

Arguments

Details

One implements a SpatPOMP model by specifying some or all of its basic components, including:

rinit,

the simulator from the distribution of the latent state process at the zero-time;

rprocess,

the transition simulator of the latent state process;

dunit_measure,

the evaluator of the conditional density at a unit's measurement given the unit's latent state;

eunit_measure,

the evaluator of the expectation of a unit's measurement given the unit's latent state;

munit_measure,

the evaluator of the moment-matched parameter set given a unit's latent state and some empirical measurement variance;

vunit_measure,

the evaluator of the variance of a unit's measurement given the unit's latent state;

runit_measure,

the simulator of a unit's measurement conditional on the unit's latent state;

dprocess,

the evaluator of the density for transitions of the latent state process;

rmeasure,

the simulator of the measurements conditional on the latent state;

dmeasure,

the evaluator of the conditional density of the measurements given the latent state;

rprior,

the simulator from a prior distribution on the parameters;

dprior,

the evaluator of the prior density;

skeleton,

which computes the deterministic skeleton of the unobserved state process;

partrans,

which performs parameter transformations.

The basic structure and its rationale are described in Asfaw et al. (2020).

Each basic component is supplied via an argument of the same name to spatPomp(). The five unit-level model components must be provided via C snippets. The remaining components, whose behaviors are inherited from pomp may be furnished using C snippets, R functions, or pre-compiled native routine available in user-provided dynamically loaded libraries.

Value

An object of class ‘spatPomp’ representing observations and model components from the spatiotemporal POMP model.

Author(s)

Kidus Asfaw, Edward L. Ionides, Aaron A. King

References

Asfaw, K., Park, J., Ho, A., King, A. A., and Ionides, E. L. (2020) Partially observed Markov processes with spatial structure via the R package spatPomp. ArXiv: 2101.01157. doi:10.48550/arXiv.2101.01157

King, A. A., Nguyen, D. and Ionides, E. L. (2016) Statistical Inference for Partially Observed Markov Processes via the R Package pomp. Journal of Statistical Software, 69(12), 1–43. doi:10.18637/jss.v069.i12

Description

spatPomp_Csnippet() is used to provide snippets of C code that specify model components. It functions similarly to Csnippet() from the pomp package; in fact, the output of spatPomp_Csnippet is an object of class Csnippet. It additionally provides some arguments that allow the user to stay focused on model development in the spatiotemporal context where model size grows.

Usage

## S4 method for signature 'character'
spatPomp_Csnippet(
  code,
  method = "",
  unit_statenames,
  unit_obsnames,
  unit_covarnames,
  unit_ivpnames,
  unit_paramnames,
  unit_vfnames
)

Arguments

Value

An object of class ‘Csnippet’ which represents a model specification in C code.

Author(s)

Kidus Asfaw

Examples

# Set initial states for Brownian motion
bm_rinit <- spatPomp_Csnippet(
  method = "rinit",
  unit_statenames = c("X"),
  unit_ivpnames = c("X"),
  code = "
    for (int u = 0; u < U; u++) {
      X[u]=X_0[u];
    }
  "
)
# Skeleton for Brownian motion
bm_skel <- spatPomp_Csnippet(
  method = "skeleton",
  unit_statenames = c("X"),
  unit_vfnames = c("X"),
  code = "
      for (int u = 0 ; u < U ; u++) {
        DX[u] = 0;
      }
  "
)

Description

An S4 class to represent a spatiotemporal POMP model and data.

Slots

unit_names

A vector containing the spatial units of a spatiotemporal POMP.

unit_statenames

A vector containing the state names such that appending the unit indices to the unit statenames will result in the each unit's corresponding states.

unit_obsnames

A vector of observation types for a spatial unit.

eunit_measure

A pomp_fun representing the expected measurement for each spatial unit given its states.

dunit_measure

A pomp_fun representing the unit measurement density for each spatial unit.

runit_measure

A pomp_fun representing the unit observation simulator.

Description

unit_names outputs the contents of the unit_names slot of a spatPomp object. The order in which the units appear in the output vector determines the order in which latent states and observations for the spatial units are stored.

Usage

## S4 method for signature 'spatPomp'
unit_names(x)

Arguments

Value

A character vector with the unit names used to create the ‘spatPomp’ object.

Description

Evaluate the unit measurement model density function for each unit. This method is used primarily as part of likelihood evaluation and parameter inference algorithms.

Usage

## S4 method for signature 'spatPomp'
vec_dmeasure(object, y, x, units, times, params, log = FALSE, ...)

Arguments

Value

An array of dimension length(unit_names(object)) by dim(x)[2] by dim(x)[3] representing each unit's measurement density assessed for each replicate in x for each observation time.

Author(s)

Kidus Asfaw

Description

Simulate from the unit measurement model density function for each unit

Usage

## S4 method for signature 'spatPomp'
vec_rmeasure(object, x, times, params, ...)

Arguments

Value

An array of dimension length(unit_names(object)) by dim(x)[2] by dim(x)[3] representing each unit's simulated measurement assessed for each replicate in x for each observation time.

Author(s)

Kidus Asfaw

Description

vunit_measure evaluates the variance of a unit's observation given the entire state

Usage

## S4 method for signature 'spatPomp'
vunit_measure(object, x, unit, time, params, Np = 1)

Arguments

Value

A matrix with the unit measurement variance implied by the state, x, and the parameter set params for unit unit.

Author(s)

Kidus Asfaw

Examples

# Complete examples are provided in the package tests
## Not run: 
b <- bm(U=3)
s <- states(b)[,1,drop=FALSE]
rownames(s) -> rn
dim(s) <- c(3,1,1)
dimnames(s) <- list(variable=rn, rep=NULL)
p <- coef(b); names(p) -> rnp
dim(p) <- c(length(p),1); dimnames(p) <- list(param=rnp)
o <- obs(b)[,1,drop=FALSE]
vunit_measure(b, x=s, unit=2, time=1, params=p)

## End(Not run)