Partially observed Markov processes with spatial structure via the R package spatPomp
Kidus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward Ionides
SStatistical Inference for Spatiotemporal PartiallyObserved Markov Processes via the R PackagespatPomp
Kidus Asfaw
Univ. of Michigan
Joonha Park
Univ. of Kansas
Allister Ho
Univ. of Michigan
Aaron A. King
Univ. of Michigan
Edward Ionides
Univ. of Michigan
Abstract
We consider inference for a class of nonlinear stochastic processes with latent dynamicvariables and spatial structure. The spatial structure takes the form of a finite collectionof spatial units that are dynamically coupled. We assume that the latent processes havea Markovian structure and that unit-specific noisy measurements are made. A modelof this form is called a spatiotemporal partially observed Markov process (SpatPOMP).The R package spatPomp provides an environment for implementing SpatPOMP models,analyzing data, and developing new inference approaches. We describe the spatPomp im-plementations of some methods with scaling properties suited to SpatPOMP models. Wedemonstrate the package on a simple Gaussian system and on a nontrivial epidemiologi-cal model for measles transmission within and between cities. We show how to constructuser-specified SpatPOMP models within spatPomp .This document is provided under the Creative Commons Attribution License. Keywords : Markov processes, hidden Markov model, state space model, stochastic dynam-ical system, maximum likelihood, plug-and-play, spatiotemporal data, mechanistic model,sequential Monte Carlo, R .
1. Introduction
A partially observed Markov process (POMP) model consists of incomplete and noisy mea-surements of a latent Markov process. A POMP model in which the latent process has spatialas well as temporal structure is called a spatiotemporal POMP or SpatPOMP. Many biolog-ical, social and physical systems have the spatiotemporal structure, dynamic stochasticityand imperfect observability that characterize SpatPOMP models. The spatial structure ofSpatPOMPs adds complexity to the problems of likelihood estimation, parameter inferenceand model selection for nonlinear and non-Gaussian systems. The objective of the spatPomp package is to facilitate model development and data analysis in the context of the general classof SpatPOMP models, enabling scientists to separate the scientific task of model developmentfrom the statistical task of providing inference tools. Thus, spatPomp brings together generalpurpose methods for carrying out Monte Carlo statistical inference for such systems. Moregenerally, spatPomp provides an abstract representation for specifying SpatPOMP models.This ensures that SpatPOMP models formulated with the package can be investigated usinga range of methods, and that new methods can be readily tested on a range of models. In its a r X i v : . [ s t a t . M E ] J a n spatPomp : Spatiotemporal Partially Observed Markov Processes in R current manifestation, spatPomp is appropriate for data analysis with a moderate number ofspatial units (up to around 100 spatial units) having nonlinear and non-Gaussian dynamics.In particular, spatPomp is not targeted at very large spatiotemporal systems such as thosethat arise in geophysical data assimilation (Anderson et al. spatPomp , but avariety of alternative methods and software are available in this case (Wikle et al. et al. et al. spatPomp package builds on the pomp package described by King et al. (2016). Math-ematically, a SpatPOMP model is also a POMP model, and this property is reflected in theobject-oriented design of spatPomp : The package is implemented using S4 classes (Chambers1998; Genolini 2008; Wickham 2019) and the basic class ‘ spatPomp ’ extends the class ‘ pomp ’provided by pomp . Among other things, this allows us to test new methods against exten-sively tested methods in low-dimensional settings, use existing convenience functions, andapply methods for POMP models on SpatPOMP models. However, standard Monte Carlostatistical inference methods for nonlinear POMP models break down with increasing spatialdimension (Bengtsson et al. curse of dimensionality . Therefore,effective inference approaches must, in practice, take advantage of the special structure ofSpatPOMP models.A SpatPOMP model is characterized by the transition density for the latent Markov processand unit-specific measurement densities . Once these elements are specified, calculating andsimulating from all joint and conditional densities are well defined operations. However,different statistical methods vary in the operations they require. Some methods require onlysimulation from the transition density whereas others require evaluation of this density. Somemethods avoid working with the model directly, replacing it by an approximation, such as alinearization. For a given model, some operations may be considerably easier to implementand so it is useful to classify inference methods according to the operations on which theydepend. In particular, an algorithm is said to be plug-and-play if it utilizes simulation of thelatent process but not evaluation of transition densities (Bretó et al. et al. et al. et al. pomp have provedcapable for state-of-the-art inference problems for POMP models (e.g., King et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. spatPomp package provides a general environment for methods with and without the We use the term “density” in this article to encompass both the continuous and discrete cases. Thus,when latent variables or measured quantities are discrete, one can replace “probability density function” with“probability mass function”. idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides spatPomp ’ inthe spatPomp package. Section 3 introduces simulation and several spatiotemporal filteringmethods currently implemented in spatPomp . Section 4 introduces some parameter estima-tion algorithms currently implemented in spatPomp which build upon these simulation andfiltering techniques. Section 5 constructs a simple linear Gaussian SpatPOMP model and usesthis example to illustrate the statistical methodology. Section 6 discusses the construction ofspatially structured compartment models for population dynamics, in the context of coupledmeasles dynamics in UK cities; this demonstrates the kind of nonlinear stochastic systemprimarily motivating the development of spatPomp . Finally, Section 7 discusses extensionsand applications of spatPomp .
2. SpatPOMP models and their representation in spatPomp
To set up notation for SpatPOMP models, suppose there are U units labeled 1 : U = { , , . . . , U } .Let T be the set of times at which the latent dynamic process is defined. The Spat-POMP framework and the spatPomp package permit the timescale T to be a discrete orcontinuous subset of R . In either case, T must contain a set of N observation times, t N = { t n : n = 1 , . . . , N } , where t < t < · · · < t N . In the examples below, we take T = [ t , t N ], with t ≤ t being the time at which we begin modeling the dynamics of thelatent Markov process. The unobserved, latent Markov process is written as { X ( t ; θ ) : t ∈ T } where X ( t ; θ ) = (cid:0) X ( t ; θ ) , . . . , X U ( t ; θ ) (cid:1) , with X u ( t ; θ ) taking values in some space X and θ a D -dimensional real-valued parameter which we write as θ = θ D ∈ R D . The spatPomp package supposes that X is a space of continuous or discrete vectors in R dim( X ) . For somepurposes it is adequate to consider the latent process only at the finite collection of timesat which it is observed. We write X n = X ( t n ; θ ), noting that we sometimes choose tosuppress the dependence on θ . We can also write X n = ( X ,n , . . . , X U,n ) = X U,n . Theinitial value X = X ( t ; θ ) may be stochastic or deterministic. The observable process { Y n = Y U,n , n ∈ N } takes values in Y U . spatPomp requires that Y = R dim( Y ) , althoughthe observations can be discrete or continuous, and it is permitted that dim( Y ) > Y u,n is condition-ally independent of { Y ˜ u, ˜ n , (˜ u, ˜ n ) = ( u, n ) } given X u,n . We suppose the existence of a jointdensity f X N , Y N of X U, N and Y U, N with respect to some reference measure. We usethe same subscript notation to denote other marginal and conditional densities. The data, y ∗ U, N = ( y ∗ , . . . , y ∗ N ), are modeled as a realization of this observation process. Spatially ortemporally dependent measurement errors can be modeled by adding suitable latent variables.The SpatPOMP structure defined above is equivalent to the following factorization of thejoint density in terms of the initial density, f X ( x ; θ ), together with the conditional transitionprobability density, f X n | X n − ( x n | x n − ; θ ), and the unit-specific measurement densities, f Y u,n | X u,n ( y u,n | x u,n ; θ ) spatPomp : Spatiotemporal Partially Observed Markov Processes in R Method Argument to Mathematical terminology spatPomp()dunit_measure dunit_measure
Evaluate f Y u,n | X u,n ( y u,n | x u,n ; θ ) runit_measure runit_measure Simulate from f Y u,n | X u,n ( y u,n | x u,n ; θ ) eunit_measure eunit_measure Evaluate e u,n ( x, θ ) = E [ Y u,n | X u,n = x ; θ ] vunit_measure vunit_measure Evaluate v u,n ( x, θ ) = Var[ Y u,n | X u,n = x ; θ ] munit_measure munit_measure m u,n ( x, V, θ ) = ψ solves v u,n ( x, ψ ) = V , e u,n ( x, ψ ) = e u,n ( x, θ ) rprocess rprocess Simulate from f X n | X n − ( x n | x n − ; θ ) dprocess dprocess Evaluate f X n | X n − ( x n | x n − ; θ ) rmeasure rmeasure Simulate from f Y n | X n ( y n | x n ; θ ) dmeasure dmeasure Evaluate f Y n | X n ( y n | x n ; θ ) rprior rprior Simulate from the prior distribution π ( θ ) dprior dprior Evaluate the prior density π ( θ ) rinit rinit Simulate from f X ( x ; θ ) timezero t0 t time times t N obs data y ∗ N states — x N coef params θ Table 1 – Model component methods for class ‘ spatPomp ’ objects and their translation intomathematical notation for SpatPOMP models. For example, the rprocess method is setusing the rprocess argument to the spatPomp constructor function.for 1 ≤ n ≤ N , 1 ≤ u ≤ U : f X N , Y N ( x N , y N ; θ ) = f X ( x ; θ ) N Y n =1 f X n | X n − ( x n | x n − ; θ ) U Y u =1 f Y u,n | X u,n ( y u,n | x u,n ; θ ) . This formalism allows the transition density, f X n | X n − , and unit-specific measurement density, f Y u,n | X u,n , to depend on n and u , so we allow for the possibility of temporally and spatiallyinhomogeneous models. A SpatPOMP model is represented in spatPomp by an S spatPomp ’. Slotsin this object encode the components of the SpatPOMP model, and can be filled or changedusing the constructor function spatPomp() and various other convenience functions. Methodsfor the class ‘ spatPomp ’ use these components to carry out computations on the model.Table 1 gives the mathematical notation corresponding to the elementary methods that canbe executed on a class ‘ spatPomp ’ object.Class ‘ spatPomp ’ inherits from the class ‘ pomp ’ defined by the pomp package. One of themain ways in which spatPomp extends pomp is the addition of unit-level specifications of themeasurement model. This reflects the modeling assumption that measurements are carried outindependently in both space and time, conditional on the current value of the latent processwhich is known as the state of the dynamic system. There are five unit-level functionalities of idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides spatPomp ’ objects: dunit_measure , runit_measure , eunit_measure , vunit_measure and munit_measure . Each functionality corresponds to an S spatPomp() constructor function of the same name. (See Table 1 for details).Only functionalities that are required to run an algorithm of interest need to be supplied inadvance. dunit_measure evaluates the probability density of the measurement of a spatialunit given its latent state vector, whereas runit_measure simulates from this conditionaldistribution. Given the latent state, eunit_measure and vunit_measure give the expectationand variance of the measurement, respectively. They are used by the ensemble Kalman filter(EnKF, Section 3.3) and iterated EnKF (Section 4.2). munit_measure returns a parametervector corresponding to given moments (mean and variance), used by one of the options fora guided particle filter (Section 3.1). Specification of the initial condition X = X ( t ; θ ) of a SpatPOMP model is similar tothat of a POMP model, and is carried out using the rinit argument of the spatPomp() constructor function. The initial distribution f X ( x ; θ ) may sometimes be a known propertyof the system but in general it must be inferred. If the transition density for the dynamicmodel, f X n | X n − ( x n | x n − ; θ ), does not depend on time and possesses a unique stationarydistribution, it may be natural to set f X ( x ; θ ) to be this stationary distribution. When noclear scientifically motivated choice of f X ( x ; θ ) exists, one can treat X as a componentof the parameter set to be estimated. In this case, f X ( x ; θ ) concentrates at a point whichdepends on θ . Scientifically, one may be interested in the impact of a vector-valued covariate process {Z ( t ) } on the latent dynamic system. For instance, the transition density, f X n | X n − , and the mea-surement density, f Y n | X n , can depend on this observed process. The covariate process mightbe shared between units, or might take unit-specific values. spatPomp allows modeling andinference for spatiotemporal processes that depend on an observed covariate process. If such acovariate exists the user can supply a class ‘ data.frame ’ object to the covar argument of the spatPomp() constructor function. This data.frame will have a column for time, spatial unit,and each of the covariates. If any of the variables in the covariates data.frame is commonamong all units the user must supply the variable names as class ‘ character ’ vectors to the shared_covarnames argument of the spatPomp() constructor function. Otherwise, all co-variates are assumed to be unit-specific, meaning that they generally take on different valuesfor each unit. spatPomp manages the task of presenting interpolated values of the covariatesto the elementary model functions at the time they are called. An example implementing aSpatPOMP model with covariates is presented in Section 6. One of the requirements for creating a new class ‘ spatPomp ’ object is a class ‘ data.frame ’ ob-ject containing observations for each spatial unit at each time. This long-format data.frame is supplied to the data argument of the spatPomp() constructor function. The package does spatPomp : Spatiotemporal Partially Observed Markov Processes in R not presuppose data at all units at all observation times. Missing data in some or all of theobservations for any given observation times are preserved, which allows the user to describea measurement model that handles missingness. The name of the data column containingobservation times is supplied to the times argument; the name of the column containing theunit names is supplied to the units argument. The t0 argument, for which the user suppliesthe initial time of the dynamics in the time unit of the measurements, is the last of the fourrequired arguments for the spatPomp() constructor. The resulting class ‘ spatPomp ’ objectstores the unit names as a vector in a slot called unit_names . Internally, only the positionsof the names in the unit_names vector are used to keep track of units. In other words, theunits are labeled 1 through U just as in our mathematical notation, and these labels can beused to construct measurement or process models that differ between units. The text la-bels in unit_names are re-attached for user-readable outputs like simulation plots and class‘ data.frame ’ outputs.Another important argument to the spatPomp() constructor is unit_statenames . This argu-ment expects a class ‘ character ’ vector of length dim( X ) containing names of the componentsof any unit’s latent state at a given time. For example, to implement a SpatPOMP modelstudying the population dynamics between frogs and dragonflies in neighboring spatial units,the user can provide unit_statenames = c(‘F’,‘D’) . The package then expects that unit-level model components will use ‘F’ and ‘D’ in their instructions.Other name vectors, unit_obsnames and unit_covarnames , are used internally to keep trackof data and covariates that have corresponding values for each unit. These need not besupplied in the spatPomp() constructor, however, since the data and covariate data.frame objects provided by the user implicitly supply them. The spatPomp package extends the C snippet facility in pomp which allows users to specifythe model components in Table 1 using inline C code. The package is therefore able to performcomputationally expensive calculations in C while outputting results in higher-level R objects.The code below illustrates the creation of a C snippet unit measurement density using the Csnippet() function.
R> example_dunit_measure <- Csnippet("+ // define measurement standard deviation for first unit+ double sd0 = sd*1.5;+ // Quantities u, Y, X, sd and lik are available in the context+ // in which this snippet is executed.+ if (u==0) lik = dnorm(Y, X, sd0, give_log);+ else lik = dnorm(Y, X, sd, give_log);+ "+ )
The example C snippet above illustrates a key difference in the practical use of the five unit-level model components in Table 1 compared to the components inherited from pomp thataccess the entire state and measurement vectors. All the filtering and inference algorithms in spatPomp assume the conditional independence of the spatial unit measurements given the idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides f Y n | X n ( y n | x n ; θ ) = U Y u =1 f Y u,n | X u,n ( y u,n | x u,n ; θ )When we specify the unit-level model components we can thus assume that the segments ofthe measurement and state vectors for the current unit are passed in during the executionof the unit-level model component. This allows the user to declare the unit-level modelcomponents by using the unit_statenames and unit_obsnames without having to specifythe names of the relevant components of the latent state and observation. For instance, theexample C snippet uses X instead of X1 , X2 , etc. Similarly, it uses Y instead of Y1 , Y2 , etc.The variable u , which takes a value between and U-1 , is passed in to each unit-level modelcomponent. This allows the user to specify heterogeneity in the unit-level model componentsacross units. Since C uses zero-based numbering, a user interested in introducing modelheterogeneity for a unit must find the position of the unit in the unit_names slot and subtractone to get the corresponding value of u . The user can then use standard conditioning logic tospecify the heterogeneity. For instance, when the example_dunit_measure C snippet aboveis executed for spatial unit 1, the u is passed in as , Y will have the value of the measurementfrom unit 1 (or Y1 ) and X will have the value of the latent state in unit 1 (or X1 ). Thisexample C snippet is coded so that the first unit has a measurement error inflated relativeto that of the other units. As in pomp , the give_log variable is a logical flag indicatingwhether the desired output is on log scale: in this example, it is passed to the log argumentof dnorm . Setting give_log=TRUE allows the calling algorithm to maintain accuracy whenmeasurement densities become very small; the same convention is followed by all base R probability distribution functions.Not all of the model components need to be supplied for any specific computation. In par-ticular, plug-and-play methodology by definition never uses dprocess . An empty dprocess slot in a class ‘ spatPomp ’ object is therefore acceptable unless a non-plug-and-play algorithmis attempted. In the package, the data and corresponding measurement times and units areconsidered necessary parts of a class ‘ spatPomp ’ object whilst specific values of the parametersand latent states are not. The construction of a new class ‘ spatPomp ’ object is illustrated in Section 6. To providesome examples of class ‘ spatPomp ’ objects, the spatPomp package includes functions bm() , lorenz() and measles() . These functions create class ‘ spatPomp ’ objects with user-specifieddimensions for a correlated Brownian motion model, the Lorenz-96 atmospheric model (Lorenz1996), and a spatiotemporal measles SEIR model, respectively. These examples can be usedto better understand the components of class ‘ spatPomp ’ objects as well as to test filteringand inference algorithms for future development.For instance, we can create four correlated Brownian motions each with ten time steps asfollows. The correlation structure and other model details are discussed in Section 5. R> U <- 4; N <- 10R> bm4 <- bm(U = U, N = N) spatPomp : Spatiotemporal Partially Observed Markov Processes in R The above code results in the creation of the bm4 object of class ‘ spatPomp ’ with simulateddata. This is done by bringing together pre-specified C snippets and adapting them to afour-dimensional process. One can inspect many aspects of bm4 , some of which are listed inTable 1, by using the corresponding accessor functions: R> obs(bm4)R> unit_names(bm4)R> states(bm4)R> as.data.frame(bm4)R> plot(bm4)R> timezero(bm4)R> time(bm4)R> coef(bm4)R> rinit(bm4)
The measles() example is described in Section 6 to demonstrate user-specified model con-struction.
3. Simulation and filtering: Elementary SpatPOMP methods
Once the user has encoded one or more SpatPOMP models as objects of class ‘ spatPomp ’,the package provides a variety of algorithms to assist with data analysis. Methods can bedivided into two categories: elementary operations that investigate a model with a fixed set ofparameter values, and inference operations that estimate parameters. This section considersthe first of these categories.A basic operation on a SpatPOMP model is to simulate a stochastic realization of the latentprocess and the resulting data. This requires specifications of rprocess and rmeasure . Ap-plying the simulate function to an object of class ‘ spatPomp ’ by default returns another objectof class ‘ spatPomp ’, within which the data y ∗ N have been replaced by a stochastic realizationof Y N . The corresponding realization of X N is accessible via the states slot, and the params slot is filled with the value of θ used in the simulation. Optionally, simulate can bemade to return a class ‘ data.frame ’ object by supplying the argument format=‘data.frame’ in the call to simulate . Section 6 provides an example of constructing a class ‘ spatPomp ’object and simulating from it.Evaluating the conditional distribution of latent process variables given currently availabledata is an elementary operation called filtering . Filtering also provides an evaluation of thelikelihood function for a fixed parameter vector. The curse of dimensionality associated withspatiotemporal models can make filtering for SpatPOMP models computationally challenging.A widely used time-series filtering technique is the particle filter, available as pfilter in the pomp package. However, most particle filter algorithms scale poorly with dimension (Bengts-son et al. et al. spatPomp package contains im-plementations of five such algorithms, four of which are described below. Spatiotemporal dataanalysis using mechanistic models is a nascent topic, and future methodological developmentsare anticipated. Since the mission of spatPomp is to be a home for such analyses, the packagedevelopers welcome contributions and collaborations to further expand the functionality of idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides spatPomp . In the remainder of this section, we describe and discuss some of the filteringmethods currently implemented in the package. The guided intermediate resampling filter (GIRF, Park and Ionides 2020) is an extension ofthe auxiliary particle filter (APF, Pitt and Shepard 1999). GIRF is appropriate for moder-ately high-dimensional SpatPOMP models with a continuous-time latent process. All particlefilters compute importance weights for proposed particles and carry out resampling to focuscomputational effort on particles consistent with the data (see reviews by Arulampalam et al. et al. pomp , the pfilter function is discussed by King et al. (2016). GIRF combines two techniques for improvedscaling of particle filters: the use of a guide function and intermediate resampling.The guide function steers particles using importance weights that anticipate upcoming obser-vations. Future measurements are considered up to a lookahead horizon, L . APF correspondsto a lookahead horizon L = 2, and a basic particle filter has L = 1. Values L ≤ S sub-intervals, and carriesout reweighting and resampling on each sub-interval. Perhaps surprisingly, intermediate re-sampling can facilitate some otherwise intractable importance sampling problems (Del Moraland Murray 2015). APF and the basic particle filter correspond to S = 1, whereas choosing S = U gives favorable scaling properties (Park and Ionides 2020).In Algorithm 1 the F , G and P superscripts indicate filtered, guide and proposal particles,respectively. The goal for the pseudocode in Algorithm 1, and subsequent algorithms inthis paper, is a succinct description of the logic of the procedure rather than a completerecipe for efficient coding. The code in spatPomp takes advantage of memory overwriting andvectorization opportunities that are not represented in this pseudocode.We call the guide in Algorithm 1 a bootstrap guide function since it is based on resamplingthe Monte Carlo residuals calculated in step 5. Another option of a guide function in girf is the simulated moment guide function developed by Park and Ionides (2020) which usesthe eunit_measure , vunit_measure and munit_measure model components together withsimulations to calculate the guide. The expectation of Monte Carlo likelihood estimates doesnot depend on the guide function, so an inexact guide approximation may lead to loss ofnumerical efficiency but does not affect the consistency of the procedure.The intermediate resampling is represented in Algorithm 1 by the loop of s = 1 , . . . , S instep 6. The intermediate times are defined by t n,s = t n + ( t n +1 − t n ) · s (cid:14) S and we write X n,s = X ( t n,s ). The resampling weights (step 12) are defined in terms of guide functionevaluations g P,jn,s . The only requirement for the guide function to achieve unbiased estimatesis that it satisfies g F,j , = 1 and g P,jN − ,S = f Y N | X N (cid:0) y ∗ N | X F,jN − ,S ; θ (cid:1) , which is the case inAlgorithm 1. The particular guide function calculated in step 11 evaluates particles using aprediction centered on a function µ ( x , s, t ; θ ) ≈ E [ X ( t ) | X ( s ) = x ; θ ] . We call µ ( x , s, t ; θ ) a deterministic trajectory associated with X ( t ). For a continuous timeSpatPOMP model, this trajectory is typically the solution to a system of differential equations0 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Algorithm 1: girf(P, Np = J , Ninter = S , Nguide = K , Lookahead = L ) , using nota-tion from Table 1 where P is a class ‘ spatPomp ’ object with definitions for rprocess , dunit_measure , rinit , skeleton , obs and coef . input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); evaluator for f Y u,n | X u,n ( y u,n | x u,n ; θ ), and µ ( x , s, t ; θ ); data, y ∗ N ; parameter, θ ; number ofparticles, J ; number of guide simulations, K ; number of intermediate timesteps, S ; number of lookahead lags, L . initialize: simulate X F,j , ∼ f X ( · ; θ ) and set g F,j , = 1 for j in 1 : J for n in 0 : N − do sequence of guide forecast times, L = ( n + 1) : min( n + L, N ) guide simulations, X G,j,k L ∼ f X L | X n (cid:0) · | X F,jn, ; θ (cid:1) for j in 1 : J , k in 1 : K guide residuals, (cid:15) j,k ,‘ = X G,j,k‘ − µ (cid:0) X F,jn , t n , t ‘ ; θ (cid:1) for j in 1 : J , k in 1 : K , ‘ in L for s in 1 : S do prediction simulations, X P,jn,s ∼ f X n,s | X n,s − (cid:0) · | X F,jn,s − ; θ (cid:1) for j in 1 : J deterministic trajectory, µ P,jn,s,‘ = µ (cid:0) X P,jn,s , t n,s , t ‘ ; θ (cid:1) for j in 1 : J , ‘ in L pseudo guide simulations, ˆ X j,kn,s,‘ = µ P,jn,s,‘ + (cid:15) j,ks − ,‘ − (cid:15) j,ks − ,n +1 + q t n +1 − t n,s t n +1 − t n, (cid:15) j,ks − ,n +1 for j in 1 : J , k in 1 : K , ‘ in L discount factor, η n,s,‘ = 1 − ( t n + ‘ − t n,s ) / { ( t n + ‘ − t max( n + ‘ − L, ) · (1 + L =1 ) } g P,jn,s = Y ‘ in L U Y u =1 " K K X k =1 f Y u,‘ | X u,‘ (cid:16) y ∗ u,‘ | ˆ X j,ku,n,s,‘ ; θ (cid:17) η n,s,‘ for j in 1 : J for j in 1 : J , w jn,s = f Y n | X n (cid:0) y n | X F,jn,s − ; θ (cid:1) g P,jn,s . g F,jn,s − if s = 1 and n = 0 g P,jn,s . g F,jn,s − else log likelihood component, c n,s = log (cid:16) J − P Jq =1 w qn,s (cid:17) normalized weights, ˜ w jn,s = w jn,s . P Jq =1 w qn,s for j in 1 : J select resample indices, r J with P [ r j = q ] = ˜ w qn,s for j in 1 : J X F,jn,s = X P,r j n,s , g F,jn,s = g P,r j n,s , (cid:15) j,ks,‘ = (cid:15) r j ,ks − ,‘ for j in 1 : J , k in 1 : K , ‘ in L end set X F,jn +1 , = X F,jn,S and g Fn +1 , ,j = g Fn,S,j for j in 1 : J endoutput: log likelihood, λ GIRF = P N − n =0 P Ss =1 c n,s , and filter particles, X F, JN, complexity: O (cid:0) J LU N ( K + S ) (cid:1) idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides skeleton (Tong 1990). In spatPomp , the argument to skeleton is a map or vector field which is numerically solved to obtain µ ( x , s, t ; θ ). It canbe specified using complied C code via a C snippet argument to spatPomp , as demonstratedin Section 6. The forecast spread around this point prediction is given by the simulatedbootstrap residuals constructed in step 5. The adapted bagged filter (Ionides et al. bagging , ( b ootstrap agg regat ing ), since a basic particle filter is also called abootstrap filter. The adapted distribution is the conditional distribution of the latent processgiven its current value and the subsequent observation (Johansen and Doucet 2008). In theadapted bagged filter, each bootstrap replicate makes a Monte Carlo approximation to a drawfrom the adapted distribution. Thus, in the pseudocode of Algorithm 2, X A,i N is a MonteCarlo sample targeting the adapted sampling distribution, f X ( x ; θ ) N Y n =1 f X n | Y n , X n − ( x n | y ∗ n , x n − ; θ ) . (1)Each adapted simulation replicate is constructed by importance sampling using proposalparticles { X P,i,jn } . The ensemble of adapted simulation replicates are then weighted usingdata in a spatiotemporal neighborhood of each observation to obtain a locally combined MonteCarlo sample targeting the filter distribution, with some approximation error due to the finitespatiotemporal neighborhood used. This local aggregation of the bootstrap replicates alsoprovides an evaluation of the likelihood function.On a given bootstrap replicate i at a given time n , all the adapted proposal particles X P,i, Jn instep 3 are necessarily close to each other in state space because they share the parent particle X A,in − . This reduces imbalance in the adapted weights in step 5, which helps to battle the curseof dimensionality that afflicts importance sampling. The combination of the replicates for thefilter estimate in step 11 is carried out using only weights in a spatiotemporal neighborhood,thus avoiding the curse of dimensionality. For any point ( u, n ), the neighborhood B u,n shouldbe specified as a subset of A u,n = { (˜ u, ˜ n ) : ˜ n < n or (˜ u < u and ˜ n = n ) } . If the model has amixing property, meaning that conditioning on the observations in the neighborhood B u,n isnegligibly different from conditioning on the full set A u,n , then the approximation involvedin this localization is adequate.Steps 1 through 7 do not involve interaction between replicates and therefore iteration over i can be carried out in parallel. If a parallel backend has been set up by the user, the abf method will parallelize computations over the replicates using multiple cores. The user canregister a parallel backend using the doParallel package (Wallig and Weston 2020, 2019)prior to calling abf . R> library("doParallel")R> registerDoParallel(3)
The neighborhood is supplied via the nbhd argument to abf as a function which takes a pointin space-time, ( u, n ), and returns a list of points in space-time which correspond to B u,n . Anexample with B u,n = { ( u − , n ) , ( u, n − } follows.2 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Algorithm 2: abf(P, replicates = I , Np = J , nbhd= B u,n ) , using notation from Ta-ble 1 where P is a class ‘ spatPomp ’ object with definitions for rprocess , dunit_measure , rinit , obs and coef . input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); evaluator for f Y u,n | X u,n ( y u,n | x u,n ; θ ); data, y ∗ N ; parameter, θ ; number of particles perreplicate, J ; number of replicates, I ; neighborhood structure, B u,n initialize adapted simulation, X A,i ∼ f X ( · ; θ ) for i in 1 : I for n in 1 : N do proposals, X P,i,jn ∼ f X n | X n − (cid:0) · | X A,in − ; θ (cid:1) for i in 1 : I , j in 1 : J w i,ju,n = f Y u,n | X u,n (cid:0) y ∗ u,n | X P,i,ju,n ; θ (cid:1) for u in 1 : U , i in 1 : I , j in 1 : J adapted resampling weights, w A,i,jn = Q Uu =1 w i,ju,n for u in 1 : U , i in 1 : I , j in 1 : J set X A,in = X P,i,jn with probability w A,i,jn (cid:16)P Jq =1 w A,i,qn (cid:17) − for i in 1 : I w P,i,ju,n = n − Y ˜ n =1 J J X q =1 Y (˜ u, ˜ n ) ∈ B u,n w i,q ˜ u, ˜ n Y (˜ u,n ) ∈ B u,n w i,j ˜ u,n for u in 1 : U , i in 1 : I , j in 1 : J end filter weights, w F,i,ju,n = w i,ju,n w P,i,ju,n P I p =1 P Jq =1 w P,p,qu,n for u in 1 : U , n in 1 : N , i in 1 : I , j in 1 : J conditional log likelihood, λ u,n = log (cid:16)P I i =1 P Jj =1 w F,i,ju,n (cid:17) for u in 1 : U , n in 1 : N set X F,ju,n = X P,i,ku,n with probability w F,i,ku,n e − λ u,n for u in 1 : U , n in 1 : N , j in 1 : J output: filter particles, X F, Jn , for n in 1 : N ; log likelihood, λ ABF = P Nn =1 P Uu =1 λ u,n complexity: O ( I J U N ) R> example_nbhd <- function(object, unit, time){+ nbhd_list = list()+ if(time>1) nbhd_list <- c(nbhd_list, list(c(unit, time-1)))+ if(unit>1) nbhd_list <- c(nbhd_list, list(c(unit-1, time)))+ return(nbhd_list)+ }
ABF can be combined with the guided intermediate resampling technique used by GIRF togive an algorithm called ABF-IR (Ionides et al. abfir . Algorithm 3 is an implementation of the ensemble Kalman filter (EnKF)(Evensen 1994;Evensen and van Leeuwen 1996). The EnKF makes a Gaussian approximation to assimilateMonte Carlo simulations from a state prediction model with data observed in the correspond-ing time step. The EnKF has two steps: the prediction and the filtering steps; the latteris known the “analysis” step in the geophysical data-assimilation literature. The predictionstep advances Monte Carlo particles to the next observation time step by using simulationsfrom a possibly non-linear transition model. In the filtering step, the sample estimate ofthe state covariance matrix and the linear measurement model are used to make a Gaussianapproximation to the conditional distribution of the state given the data. idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Algorithm 3: enkf(P, Np = J ) , using notation from Table 1 where P is a class ‘ spatPomp ’object with definitions for rprocess , eunit_measure , vunit_measure , rinit , coef and obs . input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); evaluator for e u ( X u,n , θ )and v u ( X u,n , θ ); parameter, θ ; data, y ∗ N ; number of particles, J . initialize filter particles, X F,j ∼ f X ( · ; θ ) for j in 1 : J for n in 1 : N do prediction ensemble, X P,jn ∼ f X n | X n − (cid:0) · | X F,jn − ; θ (cid:1) for j in 1 : J centered prediction ensemble, ˜ X P,jn = X P,jn − J P Jq =1 X P,qn for j in 1 : J forecast ensemble, ˆ Y jn = e u ( X P,ju,n , θ ) for j in 1 : J forecast mean, Y n = J P Jj =1 ˆ Y jn centered forecast ensemble, ˜ Y jn = ˆ Y jn − Y n for j in 1 : J forecast measurement variance, R u, ˜ u = u, ˜ u J P Jj =1 v u (cid:0) X P,ju,n , θ (cid:1) for u, ˜ u in 1 : U forecast estimated covariance, Σ Y = J − P Jj =1 ( ˜ Y jn )( ˜ Y jn ) T + R prediction and forecast sample covariance, Σ XY = J − P Jj =1 ( ˜ X P,jn )( ˜ Y jn ) T Kalman gain, K = Σ XY Σ − Y artificial measurement noise, (cid:15) jn ∼ Normal( , R ) for j in 1 : J errors, r jn = ˆ Y jn − y ∗ n for j in 1 : J filter update, X F,jn = X P,jn + K (cid:0) r jn + (cid:15) jn (cid:1) for j in 1 : J λ n = log (cid:2) φ (cid:0) y ∗ n ; Y n , Σ Y (cid:1)(cid:3) where φ ( · ; µ , Σ) is the Normal( µ , Σ) density. endoutput: filter sample, X F, Jn , for n in 1 : N ; log likelihood estimate, λ EnKF = P Nn =1 λ n complexity: O ( J U N )4 spatPomp : Spatiotemporal Partially Observed Markov Processes in R In step 8 of Algorithm 3, the conditional variance of the measurement at the current timestep is approximated by constructing a diagonal covariance matrix whose diagonal elementsare the sample average of the theoretical unit measurement variances at each unit. This iswritten using an indicator function u, ˜ u which takes value 1 if u = ˜ u and 0 otherwise. The vunit_measure model component aids in this step whereas eunit_measure specifies how wecan construct forecast data (step 5) that can be used to later update our prediction particlesin step 14. In step 12 we add artificial measurement error to arrive at a consistent samplecovariance for the filtering step (Evensen 1994; Evensen and van Leeuwen 1996), writingNormal( µ , Σ) for independent draws from a multivariate normal random variable with mean µ and variance matrix Σ.EnKF achieves good dimensional scaling relative to a particle filter (by avoiding the resam-pling step) at the expense of a Gaussian approximation in the filtering update rule. Addinghierarchical layers to the model representation can help to make the EnKF approximationapplicable in non-Gaussian contexts (Katzfuss et al. spatPomp primarily for situations with relatively low dimension, this implementation does not engagein regularization issues required when the dimension of the observation space exceeds thenumber of particles in the ensemble.Our EnKF implementation supposes we have access to the functionse u,n ( x, θ ) = E (cid:2) Y u,n | X u,n = x ; θ (cid:3) v u ( x, θ ) = Var (cid:0) Y u,n | X u,n = x ; θ (cid:1) If the measurements are unbiased, e u,n ( x u , θ ) will simply extract the measured components of x u . The measurements in a SpatPOMP do not necessarily correspond to specific componentsof the state vector, and enkf permits arbitrary relationships subject to the constraint that theuser can provide the necessary e u,n ( x, θ ) and v u ( x, θ ) functions. These functions can be definedduring the construction of a class ‘ spatPomp ’ object by supplying C snippets to the arguments eunit_measure and vunit_measure respectively. For common choices of measurement model,such as Gaussian or negative binomial, e u,n and v u,n are readily available. In general, thefunctional forms of e u,n and v u,n may depend on u and n . Users can specify more generalfunctional forms in spatPomp since the variables u and t are defined for the C snippets.Similarly, both the mathematical notation in Section 2 and the spatPomp implementationpermit arbitrary dependence on covariate time series. Algorithm 4 is an implementation of the block particle filter (BPF Rebeschini and van Handel2015), also called the factored particle filter (Ng et al. B , . . . , B K . Each unit is placed in exactly one block. BPF generatesproposal particles by simulating from the joint latent process across all blocks, exactly as theparticle filter does. However, the resampling in the filtering step is carried out independentlyfor each block, using weights corresponding only to the measurements in the block. Differentproposal particles may be successful for different blocks, and the block resampling allows thefilter particles to paste together these successful proposed blocks. BPF supposes that spatialcoupling conditional on the data occurs primarily within blocks, and is negligible betweenblocks.The user has a choice of specifying the blocks using either the block_list argument or idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Algorithm 4: bpfilter(P, Np = J , block_list = B ) using notation from Table 1 where P is a class ‘ spatPomp ’ object with definitions for rprocess , dunit_measure , rinit , obs , coef . input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); number of particles, J ;evaluator for f Y u,n | X u,n ( y u,n | x u,n ; θ ); data, y ∗ N ; parameter, θ ; blocks, B K ; initialization, X F,j ∼ f X ( · ; θ ) for j in 1 : J for n in 1 : N do prediction, X P,jn ∼ f X n | X n − (cid:0) · | X F,jn − ; θ (cid:1) for j in 1 : J block weights, w jk,n = Y u ∈B k f Y u,n | X u,n (cid:0) y ∗ u,n | X P,ju,n ; θ (cid:1) for j in 1 : J , k in 1 : K resampling indices, r jk,n with P h r jk,n = i i = ˜ w ik,n . P Jq =1 w qk,n for j in 1 : J , k in 1 : K resample, X F,j B k ,n = X P,r j,k B k ,n for j in 1 : J , k in 1 : K endoutput: log likelihood, λ BPF ( θ ) = P Nn =1 P Kk =1 log (cid:16) J P Jj =1 w jk,n (cid:17) , filter particles X F, J N complexity: O ( J U N ) block_size , but not both. block_list takes a class ‘ list ’ object where each entry is avector representing the units in a block. block_size takes an integer and evenly partitions1 : U into blocks of size approximately block_size . For example, if there are 4 units, executing bpfilter with block_size=2 is equivalent to setting block_list=list(c(1,2),c(3,4)) .
4. Inference for SpatPOMP models
We focus on iterated filtering methods (Ionides et al. spatPomp . The main idea of iterated filtering isto extend a POMP model to include dynamic parameter perturbations. Repeated filtering,with parameter perturbations of decreasing magnitude, approaches the maximum likelihoodestimate. Here, we present iterated versions of GIRF and EnKF, known as IGIRF (Park andIonides 2020) and IEnKF (Li et al. et al. spatPomp provides a platformfor developing and testing new methods, in addition to new models and data analysis.
Algorithm 5 describes igirf() , the spatPomp implementation of IGIRF. This algorithmcarries out the IF2 algorithm of Ionides et al. (2015) with filtering carried out by GIRF,therefore its implementation combines the mif2 function in pomp with girf (Algorithm 1).For Algorithm 5, we unclutter the pseudocode by using a subscript and superscript notationfor free indices, meaning subscripts and superscripts for which a value is not explicitly specifiedin the code. We use the convention that a free subscript or superscript is evaluated for allvalues in its range, leading to an implicit ‘for’ loop. This does not hold for capitalizedsubscripts and superscripts, which describe the purpose of a Monte Carlo particle, matchingusage in Algorithm 1.6 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Algorithm 5: igirf(P, params = θ , Ngirf = M , Np = J , Ninter = S , Nguide = K ,Lookahead = L ), rw.sd = σ N, D , cooling.fraction.50 = a using notation from Ta-ble 1 where P is a class ‘ spatPomp ’ object with definitions for rprocess , dunit_measure , skeleton , rinit and obs input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); evaluator for f Y u,n | X u,n ( y u,n | x u,n ; θ ), and µ ( x , s, t ; θ ); data, y ∗ N ; starting parameter, θ ;iterations, M ; particles, J ; lookahead lags, L ; intermediate timesteps, S ; randomwalk intensities, σ N, D ; cooling fraction in 50 iterations, a . note: free indices are implicit ‘for’ loops, calculated for j in 1 : J , k in 1 : K , ‘ in ( n + 1) : min( n + L, N ), u in 1 : U , d, d in 1 : D . initialize parameters, Θ F, ,jN − ,S = θ for m in 1 : M do initialize parameters, Θ F,m,j , ∼ Normal (cid:0) Θ F,m − ,jN − ,S , a m/ Σ ivp (cid:1) for (cid:2) Σ ivp (cid:3) d,d = σ ,d d,d initialize filter particles, simulate X F,j , ∼ f X (cid:16) · ; Θ F,m,j , (cid:17) and set g F,jn, = 1 for n in 0 : N − do guide simulations, X G,j,k‘ ∼ f X ‘ | X n (cid:0) · | X F,jn, ; Θ F,m,jn, (cid:1) guide residuals, (cid:15) j,k ,‘ = X G,j,k‘ − µ (cid:0) X F,jn, , t n , t ‘ ; Θ F,m,jn, (cid:1) for s in 1 : S do perturb parameters, Θ P,m,jn,s ∼ Normal (cid:0) Θ F,m,jn,s − , a m/ Σ n (cid:1) for (cid:2) Σ n (cid:3) d,d = σ n,d d,d /S prediction simulations, X P,jn,s ∼ f X n,s | X n,s − (cid:0) · | X F,jn,s − ; Θ P,m,jn,s (cid:1) deterministic trajectory, µ P,jn,s,‘ = µ (cid:0) X P,jn,s , t n,s , t ‘ ; Θ P,m,jn,s (cid:1) pseudo guide simulations,ˆ X j,kn,s,‘ = µ P,jn,s,‘ + (cid:15) j,ks − ,‘ − (cid:15) j,ks − ,n +1 + q t n +1 − t n,s t n +1 − t n, (cid:15) j,ks − ,n +1 discount factor, η n,s,‘ = 1 − ( t n + ‘ − t n,s ) / { ( t n + ‘ − t max( n + ‘ − L, ) · (1 + L =1 ) } g P,jn,s = min( n + L,N ) Y ‘ = n +1 U Y u =1 " K K X k =1 f Y u,‘ | X u,‘ (cid:16) y ∗ u,‘ | ˆ X j,ku,n,s,‘ ; Θ P,m,jn,s (cid:17) η n,s,‘ w jn,s = f Y n | X n (cid:0) y n | X F,jn,s − ; Θ F,m,jn,s − (cid:1) g P,jn,s . g F,jn,s − if s = 1, n = 0 g P,jn,s . g F,jn,s − else normalized weights, ˜ w jn,s = w jn,s . P Jq =1 w qn,s resampling indices, r J with P [ r j = q ] = ˜ w qn,s set X F,jn,s = X P,r j n,s , g F,jn,s = g P,r j n,s , (cid:15) j,ks,‘ = (cid:15) r j ,ks − ,‘ , Θ F,m,jn,s = Θ
P,m,r j n,s end end endoutput: Iterated GIRF parameter swarm, Θ
F,M, JN − ,S complexity: O (cid:0) M J LU N ( K + S ) (cid:1) idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides P,m,jn,s gives a perturbed parameter vector for θ corresponding to particle j on iteration m at the s th intermediate time between n and n + 1. The perturbations inAlgorithm 5 are taken to follow a multivariate normal distribution, with a diagonal covariancematrix scaled by σ n,d . Normal perturbations are not theoretically required, but this is acommon choice in practice. The igirf function permits perturbations to be carried out on atransformed scale, specified using the partrans argument, to accommodate situations wherenormally distributed perturbations are more natural on the log or logistic scale, or any otheruser-specified scale. For regular parameters, i.e. parameters that are not related to the initialconditions of the dynamics, it may be appropriate to set the perturbation scale independentof n . If parameters are transformed so that a unit scale is relevant, for example using alogarithmic transform for non-negative parameters, a simple choice such as σ n,d = 0 .
02 maybe effective. Initial value parameters (IVPs) are those that determine only the latent stateat time t , and these should be perturbed only at the beginning of each iteration m . Thematrix σ N, D can be constructed using the rw.sd function, which simplifies the constructionfor regular parameters and IVPs. The cooling.fraction.50 argument takes the fractionof rw.sd by which to perturb the parameters after 50 iterations of igirf . If set to 0.5, forinstance, the default behavior is to lower the perturbation standard deviation geometricallyso that it is halved by the 50 th iteration of igirf . Algorithm 6 is an implementation of the iterated ensemble Kalman filter (IEnKF) whichextends the IF2 approach for parameter estimation by replacing a particle filter with anensemble Kalman filter. As described in Section 4.1, we employ a free index notation wherebysuperscripts and subscripts that are not otherwise specified have an implicit ‘for’ loop.We note a caveat in using IEnKF. If the forecast mean ˆ Y is not dependent on a parametercomponent, that component of the parameter is not updated by the Kalman gain on average.For example, in Brownian motion, the forecast ˆ Y is independent of the measurement varianceparameter τ , and so IEnKF is ineffective in estimating τ . By contrast, for geometric Brownianmotion, which is obtained by exponentiating Brownian motion, IEnKF can estimate τ becausehigh values of τ lead to higher values of ˆ Y on average. In this case, if the average forecastis different from the observed data, the τ parameter gets updated accordingly to reduce theerror. Therefore, IEnKF may need to be coupled with other parameter estimation methods(such as IGIRF) to estimate parameters that do not affect the forecast mean. Objects of class ‘ spatPomp ’ inherit methods for inference from class ‘ pomp ’ objects imple-mented in the pomp package. As discussed earlier, IF2 (Ionides et al. spatPomp using the approximateBayesian computing (ABC) method inherited from class ‘ pomp ’, abc() . ABC has previouslybeen used for spatiotemporal inference (Brown et al. spatPomp : Spatiotemporal Partially Observed Markov Processes in R Algorithm 6: ienkf(P, params = θ , Nenkf = M , Np = J , cooling.fraction.50 = a ,rw.sd = σ N, D ) , using notation from Table 1 where P is a class ‘ spatPomp ’ object withdefinitions for rprocess , eunit_measure , vunit_measure , rinit , and obs . input: simulator for f X n | X n − ( x n | x n − ; θ ) and f X ( x ; θ ); evaluator for e u,n ( x, θ ) andv u,n ( x, θ ); data, y ∗ N ; number of particles, J ; number of iterations, M ; startingparameter, θ ; random walk intensities, σ N, D ; cooling fraction in 50 iterations, a . note: free indices are implicit ‘for’ loops, calculated for j in 1 : J , u and ˜ u in 1 : U , d and d in 1 : D . initialize parameters, Θ F, ,jN = θ for m in 1 : M do initialize parameters, Θ F,m,j ∼ Normal (cid:0) Θ F,m − ,jN , a m/ Σ (cid:1) for (cid:2) Σ n (cid:3) d,d = σ n,d d,d initialize filter particles, simulate X F,j ∼ f X (cid:16) · ; Θ F,m,j (cid:17) . for n in 1 : N do perturb parameters, Θ P,m,jn ∼ Normal (cid:0) Θ F,m,jn − , a m/ Σ n (cid:1) prediction ensemble, X P,jn ∼ f X n | X n − (cid:0) · | X F,jn − ; Θ P,m,jn (cid:1) process and parameter ensemble, Z P,jn = X P,jn Θ P,m,jn ! centered process and parameter ensemble, ˜ Z P,jn = Z P,jn − J P Jq =1 Z P,qn forecast ensemble, ˆ Y ju,n = e u ( X P,ju,n , Θ P,m,jn ) centered forecast ensemble, ˜ Y jn = ˆ Y jn − J P Jq =1 ˆ Y qn forecast measurement variance, R u, ˜ u = u, ˜ u J P Jj =1 v u ( X P,ju,n , Θ P,m,jn ) forecast sample covariance, Σ Y = J − P Jj =1 ( ˜ Y jn )( ˜ Y jn ) T + R prediction and forecast sample covariance, Σ ZY = J − P Jj =1 ( ˜ Z P,jn )( ˜ Y jn ) T Kalman gain: K = Σ ZY Σ − Y artificial measurement noise, (cid:15) jn ∼ Normal( , R ) errors, r jn = ˆ Y jn − y ∗ n filter update: Z F,jn = X F,jn Θ F,m,jn ! = Z P,jn + K (cid:0) r jn + (cid:15) jn (cid:1) end end set θ M = J P Jj =1 Θ F,M,jN output:
Monte Carlo maximum likelihood estimate, θ M . complexity: O ( M J U N ) idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides
5. Demonstrating data analysis tools on a toy model
We illustrate key capabilities of spatPomp using a model for correlated Brownian motions.This allows us to demonstrate a data analysis in a simple context where we can compare resultswith a standard particle filter as well as validate all methods against the exact solutions whichare analytically available. Here we defer the details of model construction by using a modelpre-specified within the package. Section 6 proceeds to develop a model exemplifying thekinds of nonlinear, non-Gaussian spatiotemporal dynamics of moderately high dimension,which are the target of spatPomp . Consider spatial units 1 , . . . , U located evenly around acircle, with dist( u, ˜ u ) being the circle distance,dist( u, ˜ u ) = min (cid:0) | u − ˜ u | , | u − ˜ u + U | , | u − ˜ u − U | (cid:1) . We investigate a SpatPOMP where the latent process is a U -dimensional Brownian motion X ( t ) having correlation that decays with distance. Specifically, dX u ( t ) = U X ˜ u =1 ρ dist( u, ˜ u ) dW ˜ u ( t ) , where W ( t ) , . . . , W U ( t ) are independent Brownian motions with infinitesimal variance σ , and | ρ | < 1. Using the notation in Section 2, we suppose our measurement model for discrete-timeobservations of the latent process is Y u,n = X u,n + η u,n where η u,n iid ∼ Normal(0 , τ ). The model specification is completed by specifying the initialconditions, { X u (0) , u ∈ U } . A class ‘ spatPomp ’ object which simulates from this modelfor U = 10 with N = 50 evenly-spaced observations that are one unit time apart can besimulated using bm() and plotted using plot() , yielding the plot in Figure 1 as follows: R> bm10 <- bm(U=10, N=50)R> plot(bm10)
Such plots can help the user qualitatively assess dynamics within and between the units. plot() visualizes the results of coercing bm10 into a class ‘ data.frame ’ object by using the as.data.frame(bm10) method. More customized plots can thus be created by using themany plotting options in R for class ‘ data.frame ’ objects. A detailed description of thecomponents of the bm10 object can be obtained by invoking the spy() method from pomp asfollows (the output is suppressed to conserve space): R> spy(bm10) bm10 contains all the necessary model components for likelihood estimation using the algo-rithms discussed in Section 3. The standard particle filter, GIRF, ABF, EnKF and BPF can0 spatPomp : Spatiotemporal Partially Observed Markov Processes in R U1U2U3U4U5U6U7U8U9U10 10 20 30 40 50 time un i t −505101520 Y Figure 1 – Output of the plot() method on a class ‘ spatPomp ’ object representing a simulationfrom a 10-dimensional correlated Brownian motions model with 50 observations that are oneunit time apart.be run to estimate the likelihood of the data at a given parameter set. Here, we use theparameter set that was used to simulate bm10 and show one likelihood evaluation for eachmethod.
R> theta <- coef(bm10)R> thetarho sigma tau X1_0 X2_0 X3_0 X4_0 X5_0 X6_0 X7_0 X8_00.4 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0X9_0 X10_00.0 0.0R> logLik(pfilter(bm10, params=theta, Np=1000))[1] -989.9714R> logLik(girf(bm10, params=theta, Np=100, Nguide=10, Ninter=10, lookahead=1))[1] -1334.892R> logLik(abf(bm10, params=theta, Nrep=100, Np=10))[1] -986.674R> logLik(enkf(bm10, params=theta, Np=1000))[1] -937.3424R> logLik(bpfilter(bm10, params=theta, Np=1000, block_size=2))[1] -955.7845
We see considerable differences in these initial log likelihood estimates. These might beexplained by Monte Carlo variability or bias, and additional investigation is needed to makean assessment. Both Monte Carlo uncertainty and bias are typical of spatiotemporal filteringapplications because of two main factors. First, approximations that give a filtering methodscalability have non-negligible impact, primarily bias, on the resulting likelihood estimates.Second, in high-dimensional filtering problems, a Monte Carlo filter can be expected to havenon-negligible uncertainty in the likelihood estimate even for methods designed to be scalable. idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Number of spatial units R M SE A −130−120−110−100−90 100908070605040302010 Number of spatial units
Log − li k e li hood pe r un i t B Method
ABF BPF EnKF GIRF Particle Filter
Figure 2 – A: RMSE of log likelihood estimates from ABF, BPF, EnKF, GIRF and particlefilter on correlated Brownian motions of various dimensions. For a given dimension, we runeach method 10 times and calculate the RMSE against the exact log likelihood (obtained viaKalman filter). Error bars represent the variability of the RMSE. B: Log likelihood estimatesfrom ABF, BPF, EnKF, GIRF and particle filter compared to the exact likelihood obtainedvia Kalman filter (in black).2 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Table 2 – Comparison of computational resources of the filtering algorithmsMethod Resources(core-minutes) Particles(per replicate) Replicates Guideparticles LookaheadParticle Filter 311.03 30000 - - -ABF 1588.56 50 1500 - -GIRF 3740.68 750 - 40 1EnKF 126.00 10000 - - -BPF 177.81 10000 - - -This variability translates to a negative bias in log likelihood estimates, even for methods thatare unbiased for the likelihood in the natural scale, due to Jensen’s inequality. Overall, all thefilters we study have negative bias because they make probabilistic forecasts which involvesome approximation to the true forecast distribution under the postulated model. The loglikelihood is a proper scoring rule for forecasts, meaning that the exact probabilistic forecasthas a higher expected log likelihood than an approximation, if the model is correctly specified(Gneiting et al. U and evaluate the likelihood many times for each algorithmto quantify the Monte Carlo variability. We then develop our model for increasing U . As U grows the relative performances of the algorithms can vary, so we evaluate the likelihoodusing all possible methods with several repetitions for a fixed U . On this toy problem withanalytically evaluable likelihood, we can compare all likelihood evaluations with the exactlikelihood computed using the Kalman filter. As can be seen from Figure 2, as the difficultyof the problem increases through the increase in the value of U , we quickly enter the regimein which the particle filter does worse than the methods designed specifically for SpatPOMPmodels. ABF trades off a slowly growing bias for a reduced variance by conditioning on a localneighborhood; GIRF reduces variance by using guide functions at intermediate time points toguide prediction particles towards regions of high probability, but this can be computationallycostly; BPF approximates the joint filter distribution by resampling independently betweenblocks of units; EnKF uses a Gaussian-inspired update rule to improve computational effi-ciency, and on this Gaussian problem the Gaussian approximation made by EnKF is valid,leading to strong performance. In general, since each filtering method has its strengths andlimitations, it is worthwhile on a new problem to try them all.Users will also need to keep in mind considerations about computational resources used up bythe available algorithms. Computing resources used by each algorithm for Figure 2 are givenin Table 2. Each algorithm was allowed to use 36 CPU cores to evaluate all the likelihoodsand the algorithmic settings were fixed as shown in the table. The time-complexity of GIRF isquadratic in U , due to the intermediate time step loop shown in the pseudocode in Section 3.1,whereas the other algorithms scale linearly with U for a fixed algorithmic setting. In addition,GIRF is less scalable than the other filter methods designed for SpatPOMP models. However,a positive feature of GIRF is that it shares with PF the property that it targets the exact idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides U × N matrix, the user should supply more memory if U and/or N are very large. BPF and EnKF generally run the quickest and require the leastmemory. However, the Gaussian and independent blocks assumptions, respectively, of thetwo algorithms must be reasonable to obtain likelihood estimates with low bias. The correlated Brownian motions example also serves to illustrate parameter inference usingIGIRF. Suppose we have data from the correlated 10-dimensional Brownian motions modeldiscussed above. We are interested in estimating the model parameters σ , τ , ρ . The initialconditions, { X u (0) , u ∈ U } , can be considered to be known such that these parameters willnot undergo perturbations in IGIRF.We must construct a starting parameter set for our search. R> start_params <- c(rho = 0.8, sigma = 0.4, tau = 0.2,+ X1_0 = 0, X2_0 = 0, X3_0 = 0, X4_0 = 0, X5_0 = 0,+ X6_0 = 0, X7_0 = 0, X8_0 = 0, X9_0 = 0, X10_0 = 0)
We can now run igirf() . Note that we set the parameter perturbation standard deviationto zero for the initial conditions, which allows us to only estimate our parameters of interest.
R> igirf_out <- igirf(+ bm10,+ params=start_params,+ Ngirf=30,+ Np=1000,+ Ninter=10,+ lookahead=1,+ Nguide=50,+ rw.sd=rw.sd(rho=0.02, sigma=0.02, tau=0.02,+ X1_0=0, X2_0=0, X3_0=0, X4_0=0,+ X5_0=0, X6_0=0, X7_0=0, X8_0=0, X9_0=0, X10_0=0),+ cooling.type = "geometric",+ cooling.fraction.50=0.5+ )
The output of igirf() is an object of class ‘ igirfd_spatpomp ’. We can view the finalparameter estimate and obtain a likelihood evaluation at this estimate.
R> coef(igirf_out)[c('rho','sigma','tau')]rho sigma tau0.4643952 0.9335466 1.1129127 spatPomp : Spatiotemporal Partially Observed Markov Processes in R sigma tauloglik rho iteration v a l ue Figure 3 – The output of the plot() method on the object of class ‘ igirfd_spatPomp ’ thatencodes our model for correlated Brownian motions produces convergence traces for ρ , σ and τ , and the corresponding log likelihoods. Over 30 iterations igirf() has allowed us to getwithin a neighborhood of the maximum likelihood. R> logLik(igirf_out)[1] -945.166
To get a more accurate likelihood evaluation at the final estimate, the user can run thefiltering algorithms with greater computational effort. Since our model is linear and Gaussian,the maximum likelihood estimate of our model and the likelihood at this estimate can befound analytically. The maximum log likelihood is -933.4. An enkf run at our igirf() parameter estimate yields a log likelihood estimate of -939.04. This shortfall is a reminderthat Monte Carlo optimization algorithms should usually be replicated, and may be best usedwith inference methodology that accommodates Monte Carlo error, as discussed in Section 5.3.A useful diagnostic of the parameter search is the record of improvement of our parameterestimates during the course of an igirf() run. Each iteration within an igirf run providesa parameter estimate and a likelihood evaluation at that estimate. The plot method fora class ‘ igirfd_spatPomp ’ object shows the convergence record of parameter estimates andtheir likelihood evaluations.
R> plot(igirf_out, params = c("rho", "sigma", "tau"), ncol = 2)
As shown in Figure 3, igirf() has allowed us to explore the parameter space and climb signif-icantly up the likelihood surface to within a small neighborhood of the maximum likelihood.The run took 18.85 minutes on one CPU core for this example with 10 spatial units. Forlarger models, one may require starting multiple searches of the parameter space at variousstarting points by using parallel runs of igirf() on a larger machine with multiple cores.
Proper interpretation of a parameter estimate requires uncertainty estimates. For instance,we may be interested in estimating confidence intervals for the coupling parameters of spa- idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides −930−925−920 0.2 0.3 0.4 0.5 r p r o f il e l og li k e li hood Figure 4 – An approximate 95% profile likelihood confidence interval for ρ in the 10-dimensional Brownian motion example. For each fixed value of ρ we maximize the likelihoodover all values of σ and τ . This allows us to construct the corresponding Monte Carlo con-fidence interval. The dashed blue line corresponds to the value of ρ that generated the datawhereas the vertical red lines indicate the bounds of the interval.tiotemporal models. These are parameters that influence the strength of the dependencebetween the latent dynamics in different spatial units. In our correlated Brownian motionsexample, ρ plays this role. The dependence between any two units is moderated by thedistance between the units and the value of ρ .We can often estimate confidence intervals for parameters like τ and σ which drive the dy-namics of each spatial unit. However, coupling parameters can be hard to detect becauseany signal can be overwhelmed by the inevitably high variance estimates of high-dimensionalmodels. Full-information inference methods like IGIRF which are able to mitigate high vari-ance issues in the filtering step can allow us to extract what limited information is availableon coupling parameters like ρ . In Figure 4, we have used igirf() to create a profile likelihoodfor ρ with a 95% confidence interval constructed to adjust for Monte Carlo error using themethod of Ionides et al. (2017).
6. A spatiotemporal model of measles transmission
Compartment models for population dynamics divide up the population into categories (called compartments ) which are modeled as homogeneous. The rate of flow of individuals betweena pair of compartments may depend on the count of individuals in other compartments.Compartment models have widespread scientific applications, especially in the biological andhealth sciences (Bretó et al. et al. (2010). We use this example to demonstrate how to construct spatiotemporalcompartment models in spatPomp . The measles() function in spatPomp constructs such6 spatPomp : Spatiotemporal Partially Observed Markov Processes in R an object, and here we consider the key steps in this construction. Beyond the examples inthe pomp and spatPomp packages, previous analyses using pomp with published open-sourcecode provide additional examples of compartment models (Marino et al. T = [0 , T ]. The measles model divides the population of each city into suscepti-ble ( S ), exposed ( E ), infectious ( I ), and recovered/removed ( R ) compartments. The la-tent state process is { X ( t ; θ ) , t ∈ T } = { (cid:0) X ( t ; θ ) , . . . , X U ( t ; θ ) (cid:1) , t ∈ T } with X u ( t ; θ ) = (cid:0) S u ( t ; θ ) , E u ( t ; θ ) , I u ( t ; θ ) , R u ( t ; θ ) (cid:1) . The number of individuals in each compartment for city u at time t are denoted by S u ( t ), E u ( t ), I u ( t ), and R u ( t ). The population dynamics aredescribed by the following set of stochastic differential equations: dS u ( t ) = dN BS,u ( t ) − dN SE,u ( t ) − dN SD,u ( t ) dE u ( t ) = dN SE,u ( t ) − dN EI,u ( t ) − dN ED,u ( t ) dI u ( t ) = dN EI,u ( t ) − dN IR,u ( t ) − dN ID,u ( t ) for u = 1 , . . . , U .Here, N SE,u ( t ), N EI,u ( t ), and N IR,u ( t ) denote the cumulative number of transitions, betweenthe compartments identified by the subscripts, up to time t in city u . When these aremodeled as integer-valued, the system of differential equations has step function solutions.The recruitment of susceptible individuals into city u is denoted by the counting process N BS,u ( t ). Here, B denotes a source of individuals that can enter the susceptible population,primarily modeling births. Each compartment also has an outflow, written as a transition to D , primarily representing death, which occurs at a constant per-capita rate µ . The number ofrecovered individuals R u ( t ) in city u is defined implicitly given the known census population P u ( t ) = S u ( t ) + E u ( t ) + I u ( t ) + R u ( t ). R u ( t ) plays no direct role in the dynamics, beyondaccounting for individuals not in any of the other classes. A continuous time latent processmodel is defined as the limit of the Euler scheme in Box 1 as the Euler time incrementapproaches zero. We use tildes to distinguish the numerical solution from the continuoustime model. The scheme involves initializing numerical analogues ˜ S u (0), ˜ E u (0), ˜ I u (0), ˜ R u (0) = P u (0) − ˜ S u (0) − ˜ E u (0) − ˜ I u (0) and executing the one-step transitions in Box 1 at time incrementsof δ until time T . In the limit as δ approaches zero, this results in a model with infinitesimalmean and variance given by E [ N SE,u ( t + dt ) − N SE,u ( t )] ≈ µ SE,u ( t ) S u ( t ) dt + o ( dt ) V [ N SE,u ( t + dt ) − N SE,u ( t )] ≈ (cid:2) µ SE,u ( t ) S u ( t ) + µ SE,u ( t ) S u ( t ) σ SE (cid:3) dt + o ( dt ) , where µ SE,u ( t ) = β ( t ) I u P u + X ˜ u = u v u ˜ u P u (cid:26) I ˜ u P ˜ u − I u P u (cid:27) . We use an integrated noise process with independent Gamma distributed increments thatwe use to model extrademographic stochasticity on the rate of transition from susceptibleclasses to exposed classes, following Bretó et al. (2009). Extrademographic stochasticitypermits overdispersion (McCullagh and Nelder 1989) which is often appropriate for stochasticmodels of discrete populations. The ‘ ≈ ’ in the above two approximations is a consequence ofextrademographic noise, and as σ SE becomes small it approaches equality (Bretó and Ionides2011). Here, β ( t ) denotes the seasonal transmission coefficient (He et al. v u ˜ u denotesthe number of travelers from city u to ˜ u . Here, v u ˜ u is constructed using the gravity model ofXia et al. (2004), v u ˜ u = G · dist¯ P · P u · P ˜ u dist( u, ˜ u ) , idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Pseudocode for one rprocess Euler timestep for the measles spatPomp
For u in 1 : U :i. Draw unbiased, independent, multiplicative noise for each µ SE,u ( s m +1 ),∆Γ Q Q ,u ∼ Gamma( δσ Q Q , σ Q Q ).Define ∆Γ Q i Q j ,u = δ , for ( i, j ) = (1 , S u ( s m ), ˜ E u ( s m ) and ˜ I u ( s m ): (cid:16) ∆ ˜ N Q Q ,u , ∆ ˜ N Q Q ,u , ˜ S u ( s m ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u (cid:17) ∼ Multinomial (cid:16) ˜ S u ( s m ) , p Q Q ,u , p Q Q ,u , − p Q Q ,u − p Q Q ,u (cid:17) ; (cid:16) ∆ ˜ N Q Q ,u , ∆ ˜ N Q Q ,u , ˜ E u ( s m ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u (cid:17) ∼ Multinomial (cid:16) ˜ E u ( s m ) , p Q Q ,u , p Q Q ,u , − p Q Q ,u − p Q Q ,u (cid:17) ; (cid:16) ∆ ˜ N Q Q ,u , ∆ ˜ N Q Q ,u , ˜ I u ( s m ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u (cid:17) ∼ Multinomial (cid:16) ˜ I u ( s m ) , p Q Q ,u , p Q Q ,u , − p Q Q ,u − p Q Q ,u (cid:17) , where p Q i Q j ,u = (cid:0) − exp( P k µ Q i Q k ,u ( s m +1 )∆Γ Q i Q k ,u ) (cid:1) µ Q i Q j ,u ( s m +1 )∆Γ Q i Q j ,u P k µ Q i Q k ,u ( s m +1 )∆Γ Q i Q k ,u iii. New entries into susceptible class via birth:∆ ˜ N BQ ,u ∼ Poisson( µ BQ ,u ( s m +1 ) · δ )iv. Update compartments by the one-step transitions:˜ S u ( s m +1 ) = ˜ S u ( t n ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u + ∆ ˜ N BQ ,u ˜ E u ( s m +1 ) = ˜ E u ( t n ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u + ∆ ˜ N Q Q ,u ˜ I u ( s m +1 ) = ˜ I u ( t n ) − ∆ ˜ N Q Q ,u − ∆ ˜ N Q Q ,u + ∆ ˜ N Q Q ,u ˜ R u ( s m +1 ) = P ( s m +1 ) − ˜ S u ( s m +1 ) − ˜ E u ( s m +1 ) − ˜ I u ( s m +1 )Box 1 – An Euler increment between times s m = mδ to s m +1 = s m + δ of an Euler schemewhose limit as δ approaches zero is our continuous-time measles latent process model. Fornotational convenience, Q , Q , Q , Q and Q represent susceptible (S), exposed (E), infec-tious (I), recovered (R) and natural death (D) statuses, respectively. We keep track of changesto ˜ S u , ˜ E u , ˜ I u and ˜ R u , the numerical analogues of S u , E u , I u and R u in our mathematicalmodel, by updating each compartment using dynamic rates and our population covariate.The dynamics are coupled via µ SE,u terms that incorporate travel of infectives from otherunits. Here, Gamma( α, β ) is the gamma distribution with mean αβ and variance αβ . Moreinformation about the gamma, multinomial and Poisson distributions can be found in Casellaand Berger (1990). The instructions in this box are encoded in the measles_rprocess C snippet in the following subsection.8 spatPomp : Spatiotemporal Partially Observed Markov Processes in R where dist( u, ˜ u ) denotes the distance between city u and city ˜ u , P u is the average populationfor city u across time, ¯ P is the average population across cities, and dist is the averagedistance between a randomly chosen pair of cities. In this version of the model, we model v u ˜ u as fixed through time and symmetric between any two arbitrary cities, though a naturalextension would allow for temporal variation and asymmetric movement between two cities.The gravitation constant G is scaled with respect to ¯ P and dist.The discrete set of observation times is t N = { t n , n = 1 , . . . , N } . The observations forcity u are bi-weekly new case counts. The observation process { Y n = Y U,n , n ∈ N } canbe written { Y n = cases U,n , n ∈ N } . We denote the number of true transitions fromcompartment I to compartment R accumulated between an observation time and some time t before the next observation time to be C u ( t ) = N IR,u ( t ) − N IR,u ( b t c ), where b t c is the greatestelement of t N that is less than t .Our measurement model assumes that a certain fraction, ρ , called the reporting probability,of the transitions from the infectious compartment to the recovered compartment were, onaverage, counted as reported cases. Our measurement model is:cases u,n | C u ( t n ) = c ∼ Normal( ρ c, ρ (1 − ρ ) c + ( ψ ρ c ) ) , where ψ is an overdispersion parameter that allows us to have measurement variance thatis greater than the variance of the binomial distribution with number of trials c and successprobability ρ . The construction of class ‘ spatPomp ’ objects is similar to the construction of class ‘ pomp ’ ob-jects discussed by King et al. (2016). Here, we focus on the distinctive features of SpatPOMPmodels.Suppose for our example below that we have bi-weekly measles case counts from U = 5cities in England as reported by Dalziel et al. (2016) in the object measles_cases of class‘ data.frame ’. The first few rows of this data are shown here. We see the column correspondingto time is called year and is measured in years (two weeks is equivalent to 0.038 years). year city cases1950.000 LONDON 961950.000 BIRMINGHAM 1791950.000 LIVERPOOL 5331950.000 MANCHESTER 221950.000 LEEDS 171950.038 LONDON 60 We can construct a spatPomp object by supplying three minimal requirements in additionto our data above: the column names corresponding to the spatial units and times of eachobservation ( ‘city’ and ‘year’ in this case) and the time at which the latent dynamics aresupposed to begin. Here we set this to two weeks before the first recorded observations. idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides R> measles5 <- spatPomp(+ data=measles_cases,+ units='city',+ times='year',+ t0=min(measles_cases$year)-1/26)
We can successively add each model component to measles5 with a call to spatPomp() on measles5 with the argument for each component. To avoid repetition, we will construct allof our model components and supply them all at once in a later call to spatPomp() .First, we consider covariates. Suppose that we have covariate information for each city ateach observation time in a class ‘ data.frame ’ object called measles_covar . In this case, wehave census population and lagged birthrate data. We consider lagged birthrate because weassume children enter the susceptible pool when they are old enough to go to school. year city pop lag_birthrate1950.000 LONDON 3389306.0 70571.231950.000 BIRMINGHAM 1117892.5 24117.231950.000 LIVERPOOL 802064.9 19662.961950.000 MANCHESTER 704468.0 15705.461950.000 LEEDS 509658.5 10808.731950.038 LONDON 3388407.4 70265.20 If covariate information is not reported at the same frequency as the measurement data, spatPomp will linearly interpolate the covariate data, as is the default behavior in pomp .For ease of access, the spatial unit names are mapped to the entries 1 , . . . , U . The mappingfor each unit can be found by extracting the unit’s position in:
R> unit_names(measles)
We now move from preparing our covariates to writing our model components. We shall use C snippets to specify our model components due to the computational advantages discussedin Section 2.5. spatPomp compiles the C snippets when building the class ‘ spatPomp ’ object.Before coding up our model components let us specify some global variables in C that will beaccessible to all model components. The globals argument to a spatPomp() call can be usedto supply these. A global argument that is automatically created based on the units columnof our observed data is the U variable, which encodes the number of spatial units. Since themovement matrix (cid:0) v u, ˜ u (cid:1) u, ˜ u ∈ U is calculable up to the parameter G . We can then define atwo-dimensional C array, called v_by_g that supplies each (cid:0) v u, ˜ u G (cid:1) u, ˜ u ∈ U in a C snippet called measles_globals . R> measles_globals <- Csnippet("+ const double v_by_g[5][5] = {+ {0,2.205,0.865,0.836,0.599},+ {2.205,0,0.665,0.657,0.375},+ {0.865,0.665,0,1.118,0.378}, spatPomp : Spatiotemporal Partially Observed Markov Processes in R + {0.836,0.657,1.118,0,0.580},+ {0.599,0.375,0.378,0.580,0}+ };+ ") We now construct a C snippet for initializing the latent process at time t , which correspondsto t_0 above. This involves drawing from f X ( x ; θ ). The parameter vector θ includes initialproportions of the population in each of the four compartments for each city. The namesof these initial value parameters (IVPs) will be passed in alongside other parameters to the paramnames argument of the spatPomp() constructor with names S1_0 , . . . , S5_0 , E1_0 , . . . , E5_0 , I1_0 , . . . , I5_0 and
R1_0 , . . . , R5_0 . We can use spatPomp_Csnippet() to assign thelatent states S u (0), E u (0), I u (0), and R u (0) to the numbers implied by the correspondingIVPs. We do this via the unit_statenames argument of spatPomp_Csnippet() , which, inour example below, receives the vector c("S", "E", "I", "R", "C", "W") . This functionrecognizes that all S u (0) , u ∈ U are stored contiguously in a C array (with names S1 , . . . , S5 ) and gives us access to S1 via S[0] . S2 , . . . , S5 can then be accessed as S[1] , . . . , S[U-1] .Similarly, it provides us access to E1 , I1 and R1 via E[0] , I[0] and
R[0] .The unit_ivpnames argument of spatPomp_Csnippet() serves a similar purpose. If theuser provides paramnames to the spatPomp() constructor that includes IVP names storedcontiguously, their corresponding values are also stored in a C array that can be traversedeasily. Setting unit_ivpnames = c("S") then gives us access to the initial value parameterscorresponding to the susceptible class for all units (i.e. S1_0 , . . . ,
S5_0 ) via
S_0[0] , . . . ,
S_0[U-1]
Finally, the unit_covarnames argument of spatPomp_Csnippet() similarly allows us to haveaccess to pop1 , which is the population covariate for our first city, via pop[0] . The populationsof other cities can then be found via pop[1] , . . . , pop[U-1]
These arguments to spatPomp_Csnippet() allow us to have a code argument that focuseson specifying the model component. Here, we are able to write a few lines relating the latentstates for each city at the initial time to the population in each city and the IVPs.
R> measles_rinit <- spatPomp_Csnippet(+ unit_statenames = c("S", "E", "I", "R", "C", "W"),+ unit_ivpnames = c("S", "E", "I", "R"),+ unit_covarnames = c("pop"),+ code = "+ for (int u=0; u measles_rprocess <- spatPomp_Csnippet(+ unit_statenames = c("S", "E", "I", "R", "C", "W"),+ unit_covarnames = c("pop", "lag_birthrate"),+ code = "+ double beta, br, seas, foi, dw, births;+ double rate[6], trans[6];+ int u,v;++ // school term-time seasonality+ t = (t-floor(t))*365.25;+ if ((t>=7&&t<=100) || (t>=115&&t<=199) ||+ (t>=252&&t<=300) || (t>=308&&t<=356))+ seas = 1.0+amplitude*0.2411/0.7589;+ else+ seas = 1.0-amplitude;++ // transmission rate+ beta = R0*(gamma+mu)*seas;++ for (u= 0 ; u < U ; u++) {+ br = lag_birthrate[u];++ // expected force of infection+ foi = (I[u])/pop[u];+ for (v=0; v < U ; v++) {+ if(v != u)+ foi += g * v_by_g[u][v] * (I[v]/pop[v] -+ I[u]/pop[u]) / pop[u];+ }++ // white noise (extrademographic stochasticity)+ dw = rgammawn(sigmaSE,dt);+ rate[0] = beta*foi*dw/dt; // stochastic force of infection+ rate[1] = mu; // natural S death+ rate[2] = sigma; // rate of ending of latent stage spatPomp : Spatiotemporal Partially Observed Markov Processes in R + rate[3] = mu; // natural E death+ rate[4] = gamma; // recovery+ rate[5] = mu; // natural I death++ // Poisson births+ births = rpois(br*dt);++ // transitions between classes+ reulermultinom(2,S[u],&rate[0],dt,&trans[0]);+ reulermultinom(2,E[u],&rate[2],dt,&trans[2]);+ reulermultinom(2,I[u],&rate[4],dt,&trans[4]);++ S[u] += births - trans[0] - trans[1];+ E[u] += trans[0] - trans[2] - trans[3];+ I[u] += trans[2] - trans[4] - trans[5];+ R[u] = pop[u] - S[u] - E[u] - I[u];+ W[u] += (dw - dt)/sigmaSE; // standardized i.i.d. white noise+ C[u] += trans[4]; // true incidence+ }+ "+ ) The measurement model is chosen to allow for overdispersion relative to the binomial distri-bution with success probability ρ . Here, we show the C snippet defining the unit measurementmodel. The lik variable is pre-defined and is set to the evaluation of the unit measurementdensity in either the log or natural scale depending on the value of give_log . R> measles_dunit_measure <- spatPomp_Csnippet(+ code = "+ double m= rho*C;+ double v = m*(1.0-rho+psi*psi*m);+ lik = dnorm(cases,m,sqrt(v),give_log);+ "+ ) spatPomp will then multiply the unit measurement densities over u ∈ U to compute themeasurement density at each time. The user may rather directly supply dmeasure that returnsthe product of unit-specific measurement densities. We do so and store the resulting C snippetin measles_dmeasure , but do not show the code here. This may be used, for instance, to run pfilter in pomp . We use Csnippet() since the argument to theThe runit_measure argument of the spatPomp() constructor can be supplied a C snippetfor generating data for a point in time and space given the latent state at that point. This isuseful for simulating data from a model. We construct such a C snippet here. idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides R> measles_runit_measure <- Csnippet("+ double cases;+ double m= rho*C;+ double v = m*(1.0-rho+psi*psi*m);+ cases = rnorm(m,sqrt(v));+ if (cases > 0.0) cases = nearbyint(cases);+ else cases = 0.0;+ ")
We construct the corresponding rmeasure and store it in the measles_rmeasure variable.To run the methods EnKF, IEnKF, GIRF and IGIRF, we must supply more specifica-tions about the measurement model. The first two require eunit_measure whereas thelast two require skeleton and additionally eunit_measure , vunit_measure , munit_measure when kind=‘moment’ is the desired kind of guide function for GIRF. As was the case with dunit_measure and runit_measure , the C snippets for eunit_measure , vunit_measure and munit_measure can be written assuming that the unit statenames and the u and t variableshave been pre-defined. Within the C snippet for eunit_measure , a variable named ey isdefined which should be coded to compute the quantity E [ Y u,n | X u,n ] that eunit_measure is tasked to obtain. Similarly, since vunit_measure computes a unit measurement variancegiven the parameter set and the unit states, a variable named vc is pre-defined and shouldtake the value of the computed variance. Finally, munit_measure returns a moment-matchedparameter set given the existing parameter set, the unit states, and an empirically com-puted variance. Variables with the names of the parameters prefixed by M_ (e.g. M_tau ) arepre-defined and assigned to the existing parameter values. The user need only change theparameters that would take on a new value after moment-matching.For our measles example, eunit_measure multiplies the latent modeled cases by the reportingrate.
R> measles_eunit_measure <- spatPomp_Csnippet(+ code = "+ ey = rho*C;+ "+ )vunit_measure computes the variance of the unit observation given the unit states andparameter set.
R> measles_vunit_measure <- spatPomp_Csnippet(+ code = "+ double m = rho*C;+ vc = m*(1.0-rho+psi*psi*m);+ "+ )munit_measure computes a moment-matched size parameter given an empirically calculatedvariance.4 spatPomp : Spatiotemporal Partially Observed Markov Processes in R R> measles_munit_measure <- spatPomp_Csnippet(+ code = "+ double binomial_var;+ double m;+ m = rho*C;+ binomial_var = rho*(1-rho)*C;+ if(vc > binomial_var) M_psi = sqrt(vc - binomial_var)/m;+ "+ )
The skeleton model component allows the user to specify a system of differential equations,also called a vector field, which can be numerically solved to evaluate a deterministic tra-jectory of the latent process at requested times (King et al. spatPomp_Csnippet() provides an argument called unit_vfnames which provides pointers to vector field values forthe corresponding states. The time derivatives for the susceptible classes for our five spatialunits,
DS1 , . . . , DS5 can then be assigned using
DS[0] , . . . , DS[U-1] . R> measles_skel <- spatPomp_Csnippet(+ unit_statenames = c("S", "E", "I", "R", "C", "W"),+ unit_vfnames = c("S", "E", "I", "R", "C", "W"),+ unit_covarnames = c("pop", "lag_birthrate"),+ code = "+ double beta, br, seas, foi;+ int u,v;++ // term-time seasonality+ t = (t-floor(t))*365.25;+ if ((t>=7&&t<=100) || (t>=115&&t<=199) ||+ (t>=252&&t<=300) || (t>=308&&t<=356))+ seas = 1.0+amplitude*0.2411/0.7589;+ else+ seas = 1.0-amplitude;++ // transmission rate+ beta = R0*(gamma+mu)*seas;++ // deterministic skeleton for each unit+ for (u = 0 ; u < U ; u++) {+ br = lag_birthrate[u];+ foi = I[u]/pop[u];+ for (v=0; v < U ; v++) {+ if(v != u)+ foi+=g*v_by_g[u][v]*(I[v]/pop[v]-I[u]/pop[u])/pop[u];+ }++ DS[u] = br - (beta*foi + mu)*S[u]; idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides + DE[u] = beta*foi*S[u] - (sigma+mu)*E[u];+ DI[u] = sigma*E[u] - (gamma+mu)*I[u];+ DR[u] = gamma*I[u] - mu*R[u];+ DW[u] = 0;+ DC[u] = gamma*I[u];+ }+ "+ ) Finally we declare the names of states and parameters. This will allow the compilation of themodel components which refer to these names.
R> measles_unit_statenames <- c('S','E','I','R','C','W')R> measles_covarnames <- paste0(rep(c("pop","lag_birthrate"),each=U),1:U)R> measles_statenames <- paste0(rep(measles_unit_statenames,each=U),1:U)
As discussed above, some unit_statenames may be used to keep track of accumulations ofother unit_statenames over an observation time period. The spatPomp() constructor pro-vides an argument called unit_accumvars to handle this behavior. Among other things, thisextends pomp ’s feature of resetting such variables to zero at the beginning of a measurementperiod.A parameter can sometimes be classified as an initial value parameter (IVP) that determinesonly the initial condition, or a regular parameters (RP) that contributes to the process ormeasurement model throughout the observed time interval. This classification, when it exists,can be helpful since there are inferential consequences. Precision on estimates of IVPs maynot grow with increasing number, N , of observations, whereas for RPs we expect increasinginformation with increasing N . R> measles_IVPnames <- paste0(measles_statenames[1:(4*U)],"_0")R> measles_RPnames <- c("R0","amplitude","gamma","sigma","mu",+ "sigmaSE","rho","psi","g")R> measles_paramnames <- c(measles_RPnames,measles_IVPnames)
The pieces of the SpatPOMP are now combined by a call to spatPomp . R> measles5_full <- spatPomp(+ data = measles5,+ covar = measles_covar,+ unit_statenames = measles_unit_statenames,+ unit_accumvars = c("C","W"),+ paramnames = measles_paramnames,+ rinit = measles_rinit,+ rprocess = euler(measles_rprocess, delta.t=2/365),+ skeleton=vectorfield(measles_skel),+ dunit_measure = measles_dunit_measure, spatPomp : Spatiotemporal Partially Observed Markov Processes in R + eunit_measure = measles_eunit_measure,+ vunit_measure = measles_vunit_measure,+ munit_measure = measles_munit_measure,+ runit_measure = measles_runit_measure,+ dmeasure = measles_dmeasure,+ rmeasure = measles_rmeasure,+ globals = measles_globals+ ) Suppose we wanted to simulate data from our model for measles dynamics in the U =5 citiesand that we have a parameter set m5_params at which we are simulating. We can compare oursimulations to the data using the code below and the plot() method on the class ‘ spatPomp ’objects resulting from the simulation and the measles5_full object (which includes the trueobservations) respectively. For epidemiological settings, it helps to set the argument log=TRUE of the plot() method to focus more on seasonal trends and less on spikes in case counts. Theresulting plots are shown in Figure 5. This figure may indicate room for improvement inthe current parameter vector or model structure. As discussed before, such plots can becustomized by working directly with the class ‘ data.frame ’ output of as.data.frame() . R> m5_paramsR0 amplitude gamma sigma mu sigmaSE rho5.68e+01 5.54e-01 3.04e+01 2.89e+01 2.00e-02 2.00e-02 4.88e-01psi g S1_0 S2_0 S3_0 S4_0 S5_01.16e-01 1.00e+02 2.97e-02 2.97e-02 2.97e-02 2.97e-02 2.97e-02E1_0 E2_0 E3_0 E4_0 E5_0 I1_0 I2_05.17e-05 5.17e-05 5.17e-05 5.17e-05 5.17e-05 5.14e-05 5.14e-05I3_0 I4_0 I5_0 R1_0 R2_0 R3_0 R4_05.14e-05 5.14e-05 5.14e-05 9.70e-01 9.70e-01 9.70e-01 9.70e-01R5_09.70e-01R> m5_sim <- simulate(measles5_full, params=m5_params)
7. Conclusion
The spatPomp package is designed to be both a tool for data analysis based on SpatPOMPmodels and a principled computational framework for the ongoing development of inferencealgorithms. The model specification language provided by spatPomp is very general, andimplementing a SpatPOMP model in spatPomp makes a wide range of inference algorithmsavailable. These two features facilitate objective comparison of alternative models and meth-ods.As a development platform, spatPomp is particularly convenient for implementing algorithmswith the plug-and-play property, since models will typically be defined by their rprocess idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides LONDONBIRMINGHAMLIVERPOOLMANCHESTERLEEDS time A time B log(cases+1) Figure 5 – A: Bi-weekly observed measles case counts in the five largest cities in England.B: Simulations from the measles SEIR model encoded in the class spatPomp object called measles5_full . The figure indicates that the parameter vector and/or the model structureof our SEIR model need to be altered to get patterns similar to those observed in the data.simulator, together with rmeasure and often dunit_measure , but can accommodate infer-ence methods based on other model components such as dprocess if they are available. Asan open-source project, the package readily supports expansion, and the authors invite com-munity participation in the spatPomp project in the form of additional inference algorithms,improvements and extensions of existing algorithms, additional model/data examples, docu-mentation contributions and improvements, bug reports, and feature requests.Complex models and large datasets can challenge computational resources. With this in mind,key components of the spatPomp package are written in C , and spatPomp provides facilitiesfor users to write models either in R or, for the acceleration that typically proves necessary inapplications, in C . Multi-processor computing also becomes necessary for ambitious projects.The two most common computationally intensive tasks are the assessment of Monte Carlovariability and the investigation of the roles of starting values and other algorithmic settingson optimization routines. These analyses require only embarrassingly parallel computationsand need no special discussion here.Practical modeling and inference for spatiotemporal partially observed systems, capable ofhandling scientifically motivated nonlinear, non-stationary stochastic models, is the last openproblem of the challenges raised by Bjørnstad and Grenfell (2001). Recent studies haveunderscored the need for deeper analyses of spatially coupled dynamics (Dalziel et al. et al. et al. et al. spatPomp package addresses these challenges by combining access to modern8 spatPomp : Spatiotemporal Partially Observed Markov Processes in R algorithmic developments with a suitable framework for model specification. The capabilityto carry out statistically efficient inference for general spatiotemporal systems will promote thedevelopment, criticism, refinement and validation of new spatiotemporal models. Nonlinearinteracting systems are hard to understand intuitively even when there are relatively fewunits. Even the single-unit case, corresponding to a low-dimensional nonlinear stochasticdynamic system with a low-dimensional observation process, has rich mathematical theory.Statistically efficient inference for this low-dimensional case was not generally available beforethe recent development of iterated filtering and particle Markov Chain Monte Carlo methods,and application of these methods has been assisted by their implementations in pomp . Weanticipate there is much to be gained scientifically by carrying out modeling and inferencefor spatiotemporal processes with relatively few spatial units but nevertheless surpassing thecapabilities of previous software. Facilitating this task is the primary goal of spatPomp . Acknowledgments
This work was supported by National Science Foundation grants DMS-1761603 and DMS-1646108, and National Institutes of Health grants 1-U54-GM111274 and 1-U01-GM110712.
References
Anderson J, Hoar T, Raeder K, Liu H, Collins N, Torn R, Avellano A (2009). “The data as-similation research testbed: A community facility.”
Bulletin of the American MeteorologicalSociety , (9), 1283–1296. doi:10.1175/2009BAMS2618.1 .Arulampalam MS, Maskell S, Gordon N, Clapp T (2002). “A tutorial on particle filters for on-line nonlinear, non-Gaussian Bayesian tracking.” IEEE Transactions on Signal Processing , , 174–188. doi:10.1109/78.978374 .Bakker KM, Martinez-Bakker ME, Helm B, Stevenson TJ (2016). “Digital epidemiologyreveals global childhood disease seasonality and the effects of immunization.” Proceedings ofthe National Academy of Sciences , (24), 6689–6694. doi:10.1073/pnas.1523941113 .Becker AD, Birger RB, Teillant A, Gastanaduy PA, Wallace GS, Grenfell BT (2016). “Esti-mating enhanced prevaccination measles transmission hotspots in the context of cross-scaledynamics.” Proceedings of the National Academy of Sciences , (51), 14595–14600. doi:10.1073/pnas.1604976113 .Becker AD, Wesolowski A, Bjørnstad ON, Grenfell BT (2019). “Long-term dynamics ofmeasles in London: Titrating the impact of wars, the 1918 pandemic, and vaccination.” PLoS Computational Biology , (9), e1007305. doi:10.1371/journal.pcbi.1007305 .Bengtsson T, Bickel P, Li B (2008). “Curse-of-dimensionality revisited: Collapse of the particlefilter in very large scale systems.” In T Speed, D Nolan (eds.), Probability and Statistics:Essays in Honor of David A. Freedman , pp. 316–334. Institute of Mathematical Statistics,Beachwood, OH. doi:10.1214/193940307000000518 . idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Journal of the American Statistical Association , , 440–451. doi:10.1198/jasa.2011.ap10323 .Bjørnstad ON, Grenfell BT (2001). “Noisy clockwork: Time series analysis of populationfluctuations in animals.” Science , , 638–643. doi:10.1126/science.1062226 .Blackwood JC, Cummings DAT, Broutin H, Iamsirithaworn S, Rohani P (2013a). “De-ciphering the impacts of vaccination and immunity on pertussis epidemiology in Thai-land.” Proceedings of the National Academy of Sciences of the USA , , 9595–9600. doi:10.1073/pnas.1220908110 .Blackwood JC, Streicker DG, Altizer S, Rohani P (2013b). “Resolving the roles of immunity,pathogenesis, and immigration for rabies persistence in vampire bats.” Proceedings of theNational Academy of Sciences of the USA . doi:10.1073/pnas.1308817110 .Blake IM, Martin R, Goel A, Khetsuriani N, Everts J, Wolff C, Wassilak S, Aylward RB,Grassly NC (2014). “The role of older children and adults in wild poliovirus transmission.” Proceedings of the National Academy of Sciences of the USA , (29), 10604–10609. .Bretó C (2014). “On idiosyncratic stochasticity of financial leverage effects.” Statistics &Probability Letters , , 20–26. doi:http://dx.xdoi.org/10.1016/j.spl.2014.04.003 .Bretó C, He D, Ionides EL, King AA (2009). “Time series analysis via mechanistic models.” Annals of Applied Statistics , , 319–348. doi:10.1214/08-AOAS201 .Bretó C, Ionides EL (2011). “Compound Markov counting processes and their applicationsto modeling infinitesimally over-dispersed systems.” Stochastic Processes and their Appli-cations , , 2571–2591. doi:10.1016/j.spa.2011.07.005 .Brown GD, Porter AT, Oleson JJ, Hinman JA (2018). “Approximate Bayesian computationfor spatial SEIR(S) epidemic models.” Spatial and Spatio-temporal Epidemiology , , 27–37. doi:10.1016/j.sste.2017.11.001 .Buhnerkempe MG, Prager KC, Strelioff CC, Greig DJ, Laake JL, Melin SR, DeLong RL,Gulland F, Lloyd-Smith JO (2017). “Detecting signals of chronic shedding to explainpathogen persistence: Leptospira interrogans in California sea lions.” Journal of AnimalEcology , (3), 460–472. doi:10.1111/1365-2656.12656 .Cappello C, De Iaco S, Posa D (2020). “covatest: An R Package for Selecting a Class ofSpace-Time Covariance Functions.” Journal of Statistical Software , (1), 1–42. doi:10.18637/jss.v094.i01 .Casella G, Berger RL (1990). Statistical Inference . Wadsworth, Pacific Grove.Chambers JM (1998).
Programming with Data: A Guide to the S Language . Springer Science& Business Media.0 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Dalziel BD, Bjørnstad ON, van Panhuis WG, Burke DS, Metcalf CJE, Grenfell BT (2016).“Persistent chaos of measles epidemics in the prevaccination United States caused by a smallchange in seasonal transmission patterns.”
PLoS Computational Biology , (2), e1004655. doi:10.1371/journal.pcbi.1004655 .Del Moral P, Murray LM (2015). “Sequential Monte Carlo with highly informative observa-tions.” Journal on Uncertainty Quantification , , 969–997. doi:10.1137/15M1011214 .Doucet A, Johansen A (2011). “A tutorial on particle filtering and smoothing: Fifteenyears later.” In D Crisan, B Rozovsky (eds.), Oxford Handbook of Nonlinear Filtering .Oxford University Press. URL https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/johansen/publications/dj11.pdf .Earn DJ, He D, Loeb MB, Fonseca K, Lee BE, Dushoff J (2012). “Effects of school closureon incidence of pandemic influenza in Alberta, Canada.”
Annals of Internal Medicine , ,173–181. doi:10.7326/0003-4819-156-3-201202070-00005 . .Evensen G (1994). “Sequential data assimilation with a nonlinear quasi-geostrophic modelusing Monte Carlo methods to forecast error statistics.” Journal of Geophysical Research:Oceans , (C5), 10143–10162. doi:10.1029/94JC00572 .Evensen G, van Leeuwen PJ (1996). “Assimilation of geostat altimeter data for the AgulhasCurrent using the ensemble Kalman filter with a quasigeostrophic model.” Monthly WeatherReview , , 58–96. doi:10.1175/1520-0493(1996)124<0085:AOGADF>2.0.CO;2 .Genolini C (2008). “A (not so) short introduction to S4.” Technical report , The R-Projectfor Statistical Computing. URL http://christophe.genolini.free.fr/webTutorial/S4tutorialV0-5en.pdf .Gneiting T, Balabdaoui F, Raftery AE (2007). “Probabilistic forecasts, calibration and sharp-ness.”
Journal of the Royal Statistical Society, Series B (Statistical Methodology) , , 243–268. doi:10.1111/j.1467-9868.2007.00587.x .He D, Dushoff J, Day T, Ma J, Earn DJD (2013). “Inferring the causes of the three waves ofthe 1918 influenza pandemic in England and Wales.” Proceedings of the Royal Society ofLondon, Series B , , 20131345. doi:10.1098/rspb.2013.1345 .He D, Ionides EL, King AA (2010). “Plug-and-play inference for disease dynamics: Measles inlarge and small towns as a case study.” Journal of the Royal Society Interface , , 271–283. doi:10.1098/rsif.2009.0151 .Ionides EL, Asfaw K, Park J, King AA (2020). “Bagged filters for partially observed spa-tiotemporal systems.” arXiv:2002.05211v2 .Ionides EL, Breto C, Park J, Smith RA, King AA (2017). “Monte Carlo profile confidenceintervals for dynamic systems.” Journal of the Royal Society Interface , , 1–10. doi:10.1098/rsif.2017.0126 .Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA (2015). “Inference for dynamic andlatent variable models via iterated, perturbed Bayes maps.” Proceedings of the NationalAcademy of Sciences of the USA , (3), 719—-724. doi:10.1073/pnas.1410597112 . idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Statistics & Proba-bility Letters , , 1498–1504. doi:10.1016/j.spl.2008.01.032 .Kain MP, Childs ML, Becker AD, Mordecai EA (2020). “Chopping the tail: How preventingsuperspreading can help to maintain COVID-19 control.” MedRxiv . doi:10.1101/2020.06.30.20143115 .Kantas N, Doucet A, Singh SS, Maciejowski J, Chopin N, et al. (2015). “On particle methodsfor parameter estimation in state-space models.” Statistical Science , (3), 328–351. doi:10.1214/14-STS511 .Katzfuss M, Stroud JR, Wikle CK (2020). “Ensemble Kalman methods for high-dimensionalhierarchical dynamic space-time models.” Journal of the American Statistical Association , (530), 866–885. doi:10.1080/01621459.2019.1592753 .King AA, Ionides EL, Pascual M, Bouma MJ (2008). “Inapparent infections and choleradynamics.” Nature , , 877–880. doi:10.1038/nature07084 .King AA, Nguyen D, Ionides EL (2016). “Statistical inference for partially observed Markovprocesses via the R Package pomp.” Journal of Statistical Software , , 1–43. doi:10.18637/jss.v069.i12 .Lau MS, Becker AD, Korevaar HM, Caudron Q, Shaw DJ, Metcalf CJE, Bjørnstad ON,Grenfell BT (2020). “A competing-risks model explains hierarchical spatial coupling ofmeasles epidemics en route to national elimination.” Nature Ecology & Evolution , pp. 1–6. doi:10.1038/s41559-020-1186-6 .Lee EC, Chao DL, Lemaitre JC, Matrajt L, Pasetto D, Perez-Saez J, Finger F, Rinaldo A,Sugimoto JD, Halloran ME, et al. (2020). “Achieving coordinated national immunity andcholera elimination in Haiti through vaccination: a modelling study.”
The Lancet GlobalHealth , (8), e1081–e1089. doi:10.1016/S2214-109X(20)30310-7 .Li R, Pei S, Chen B, Song Y, Zhang T, Yang W, Shaman J (2020). “Substantial undocumentedinfection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2).” Science , (6490), 489–493. doi:10.1126/science.abb3221 .Lorenz EN (1996). “Predictability: A problem partly solved.” Proceedings of the Seminar onPredictability , , 1–18.Marino JA, Peacor SD, Bunnell D, Vanderploeg HA, Pothoven SA, Elgin AK, Bence JR, JiaoJ, Ionides EL (2019). “Evaluating consumptive and nonconsumptive predator effects onprey density using field time-series data.” Ecology , (3), e02583. doi:10.1002/ecy.2583 .Martinez-Bakker M, King AA, Rohani P (2015). “Unraveling the transmission ecology ofpolio.” PLoS Biology , (6), e1002172. doi:10.1371/journal.pbio.1002172 .McCullagh P, Nelder JA (1989). Generalized Linear Models . 2nd edition. Chapman and Hall,London.Ng B, Peshkin L, Pfeffer A (2002). “Factored particles for scalable monitoring.”
Proceedings ofthe 18th Conferenece on Uncertainty and Artificial Intelligence , pp. 370–377. .2 spatPomp : Spatiotemporal Partially Observed Markov Processes in R Park J, Ionides EL (2020). “Inference on high-dimensional implicit dynamic models usinga guided intermediate resampling filter.”
Statistics & Computing , , 1497–1522. doi:10.1007/s11222-020-09957-3 .Pitt MK, Shepard N (1999). “Filtering via simulation: Auxillary particle filters.” Jour-nal of the American Statistical Association , , 590–599. doi:10.1080/01621459.1999.10474153 .Pons-Salort M, Grassly NC (2018). “Serotype-specific immunity explains the incidence ofdiseases caused by human enteroviruses.” Science , (6404), 800–803. doi:10.1126/science.aat6777 .Ranjeva SL, Baskerville EB, Dukic V, Villa LL, Lazcano-Ponce E, Giuliano AR, Dwyer G,Cobey S (2017). “Recurring infection with ecologically distinct HPV types can explain highprevalence and diversity.” Proceedings of the National Academy of Sciences , p. 201714712. doi:10.1073/pnas.1714712114 .Rebeschini P, van Handel R (2015). “Can local particle filters beat the curse of dimensional-ity?”
The Annals of Applied Probability , (5), 2809–2866. doi:10.1214/14-AAP1061 .Roy M, Bouma MJ, Ionides EL, Dhiman RC, Pascual M (2013). “The potential eliminationof Plasmodium vivax malaria by relapse treatment: Insights from a transmission modeland surveillance data from NW India.” PLoS Neglected Tropical Diseases , , e1979. doi:10.1371/journal.pntd.0001979 .Shrestha S, Foxman B, Weinberger DM, Steiner C, Viboud C, Rohani P (2013). “Identify-ing the interaction between influenza and pneumococcal pneumonia using incidence data.” Science Translational Medicine , , 191ra84. doi:10.1126/scitranslmed.3005982 .Shrestha S, King AA, Rohani P (2011). “Statistical inference for multi-pathogen systems.” PLoS Computational Biology , , e1002135. doi:10.1371/journal.pcbi.1002135 .Sigrist F, Kunsch HR, Stahel WA (2015). “spate: An R package for spatio-temporal modelingwith a stochastic advection-diffusion process.” Journal of Statistical Software , (14), 1–23. doi:10.18637/jss.v063.i14 .Snyder C, Bengtsson T, Morzfeld M (2015). “Performance bounds for particle filtersusing the optimal proposal.” Monthly Weather Review , (11), 4750–4761. doi:10.1175/MWR-D-15-0144.1 .Tong H (1990). Non-linear Time Series: A Dynamical System Approach . Oxford SciencePubl., Oxford.Wallig M, Weston S (2019). doParallel: Foreach Parallel Adaptor for the ’parallel’ Package .R package version 1.0.15, URL https://CRAN.R-project.org/package=doParallel .Wallig M, Weston S (2020). foreach: Provides Foreach Looping Construct . R package version1.5.0, URL https://CRAN.R-project.org/package=foreach .Wesolowski A, Eagle N, Tatem AJ, Smith DL, Noor AM, Snow RW, Buckee CO (2012).“Quantifying the impact of human mobility on malaria.”
Science , (6104), 267–270. doi:10.1126/science.1223467 . idus Asfaw, Joonha Park, Allister Ho, Aaron A. King, Edward L. Ionides Proceedings of the National Academy of Sciences , (38), 11887–11892. doi:10.1073/pnas.1504964112 .Wickham H (2019). Advanced R . CRC press.Wikle CK, Zammit-Mangion A, Cressie N (2019).
Spatio-temporal Statistics with R . CRCPress.Xia Y, Bjørnstad ON, Grenfell BT (2004). “Measles metapopulation dynamics: A gravitymodel for epidemiological coupling and dynamics.”
American Naturalist , (2), 267–281. doi:10.1086/422341 . Affiliation:
Kidus Asfaw, Edward Ionides, Allister HoDepartment of StatisticsUniversity of Michigan48109 Michigan, United States of AmericaE-mail: [email protected] , [email protected] , [email protected] URL:
Aaron A. KingDepartment of Ecology & Evolutionary BiologyCenter for the Study of Complex SystemsUniversity of Michigan48109 Michigan, United States of AmericaE-mail: [email protected]
URL: http://kinglab.eeb.lsa.umich.edu/
Joonha ParkDepartment of MathematicsUniversity of Kansas66045 Kansas, United States of AmericaE-mail: [email protected]