Package 'pomp'

Description

The pomp package provides facilities for inference on time series data using partially-observed Markov process (POMP) models. These models are also known as state-space models, hidden Markov models, or nonlinear stochastic dynamical systems. One can use pomp to fit nonlinear, non-Gaussian dynamic models to time-series data. The package is both a set of tools for data analysis and a platform upon which statistical inference methods for POMP models can be implemented.

Data analysis using pomp

pomp provides algorithms for:

1. Simulation of stochastic dynamical systems; see simulate.

2. Particle filtering (AKA sequential Monte Carlo or sequential importance sampling); see pfilter and wpfilter.

3. The iterated filtering methods of Ionides et al. (2006, 2011, 2015); see mif2.

4. The nonlinear forecasting algorithm of Kendall et al. (2005); see nlf.

5. The particle MCMC approach of Andrieu et al. (2010); see pmcmc.

6. The probe-matching method of Kendall et al. (1999, 2005); see probe_match.

7. Synthetic likelihood a la Wood (2010); see probe.

8. A spectral probe-matching method (Reuman et al. 2006, 2008); see spect_match.

9. Approximate Bayesian computation (Toni et al. 2009); see abc.

10. The approximate Bayesian sequential Monte Carlo scheme of Liu & West (2001); see bsmc2.

11. Ensemble and ensemble adjusted Kalman filters; see kalman.

12. Simple trajectory matching; see traj_match.

The package also provides various tools for plotting and extracting information on models and data.

Structure of the package

pomp algorithms are arranged into several levels. At the top level, estimation algorithms estimate model parameters and return information needed for other aspects of inference. Elementary algorithms perform common operations on POMP models, including simulation, filtering, and application of diagnostic probes; these functions may be useful in inference, but they do not themselves perform estimation. At the lowest level, workhorse functions provide the interface to basic POMP model components. Beyond these, pomp provides a variety of auxiliary functions for manipulating and extracting information from ‘pomp’ objects, producing diagnostic plots, facilitating reproducible computations, and so on.

Implementing a model

The basic structure at the heart of the package is the ‘pomp object’. This is a container holding a time series of data (possibly multivariate) and a model. The model is specified by specifying some or all of its basic model components. One does this using the basic component arguments to the pomp constructor. One can also add, modify, or delete basic model components “on the fly” in any pomp function that accepts them.

Documentation and examples

The package contains a number of examples. Some of these are included in the help pages. In addition, several pre-built POMP models are included with the package. Tutorials and other documentation, including a package FAQ, are available from the package website.

Useful links

• pomp homepage: https://kingaa.github.io/pomp/

• Report bugs to: https://github.com/kingaa/pomp/issues

• Frequently asked questions: https://kingaa.github.io/pomp/FAQ.html

• User guides and tutorials: https://kingaa.github.io/pomp/docs.html

• pomp news: https://kingaa.github.io/pomp/blog.html

Citing pomp

Execute citation("pomp") to view the correct citation for publications.

Author(s)

Aaron A. King

References

A. A. King, D. Nguyen, and E. L. Ionides. Statistical inference for partially observed Markov processes via the R package pomp. Journal of Statistical Software 69(12), 1–43, 2016. doi:10.18637/jss.v069.i12. An updated version of this paper is available on the pomp package website.

See the package website for more references, including many publications that use pomp.

See Also

Useful links:

• https://kingaa.github.io/pomp/

• Report bugs at https://github.com/kingaa/pomp/issues/

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

More on pomp workhorse functions: dinit(), dmeasure(), dprior(), dprocess(), emeasure(), flow(), partrans(), rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

More on pomp estimation algorithms: abc(), bsmc2(), estimation_algorithms, mif2(), nlf, pmcmc(), probe_match, spect_match

More on pomp elementary algorithms: elementary_algorithms, kalman, pfilter(), probe(), simulate(), spect(), trajectory(), wpfilter()

Description

The approximate Bayesian computation (ABC) algorithm for estimating the parameters of a partially-observed Markov process.

Usage

## S4 method for signature 'data.frame'
abc(
  data,
  Nabc = 1,
  proposal,
  scale,
  epsilon,
  probes,
  params,
  rinit,
  rprocess,
  rmeasure,
  dprior,
  ...,
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'pomp'
abc(
  data,
  Nabc = 1,
  proposal,
  scale,
  epsilon,
  probes,
  ...,
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'probed_pomp'
abc(data, probes, ..., verbose = getOption("verbose", FALSE))

## S4 method for signature 'abcd_pomp'
abc(
  data,
  Nabc,
  proposal,
  scale,
  epsilon,
  probes,
  ...,
  verbose = getOption("verbose", FALSE)
)

Arguments

Running ABC

abc returns an object of class ‘abcd_pomp’. One or more ‘abcd_pomp’ objects can be joined to form an ‘abcList’ object.

Re-running ABC iterations

To re-run a sequence of ABC iterations, one can use the abc method on a ‘abcd_pomp’ object. By default, the same parameters used for the original ABC run are re-used (except for verbose, the default of which is shown above). If one does specify additional arguments, these will override the defaults.

Continuing ABC iterations

One can continue a series of ABC iterations from where one left off using the continue method. A call to abc to perform Nabc=m iterations followed by a call to continue to perform Nabc=n iterations will produce precisely the same effect as a single call to abc to perform Nabc=m+n iterations. By default, all the algorithmic parameters are the same as used in the original call to abc. Additional arguments will override the defaults.

Methods

The following can be applied to the output of an abc operation:

abc

repeats the calculation, beginning with the last state

continue

continues the abc calculation

plot

produces a series of diagnostic plots

traces

produces an mcmc object, to which the various coda convergence diagnostics can be applied

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

Author(s)

Edward L. Ionides, Aaron A. King

References

J.-M. Marin, P. Pudlo, C. P. Robert, and R. J. Ryder. Approximate Bayesian computational methods. Statistics and Computing 22, 1167–1180, 2012. doi:10.1007/s11222-011-9288-2.

T. Toni and M. P. H. Stumpf. Simulation-based model selection for dynamical systems in systems and population biology. Bioinformatics 26, 104–110, 2010. doi:10.1093/bioinformatics/btp619.

T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. H. Stumpf. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface 6, 187–202, 2009. doi:10.1098/rsif.2008.0172.

See Also

More on methods based on summary statistics: basic_probes, nlf, probe(), probe_match, spect(), spect_match

More on pomp estimation algorithms: bsmc2(), estimation_algorithms, mif2(), nlf, pmcmc(), pomp-package, probe_match, spect_match

More on Markov chain Monte Carlo methods: pmcmc(), proposals

More on Bayesian methods: bsmc2(), dprior(), pmcmc(), prior_spec, rprior()

Description

Latent state variables that accumulate quantities through time.

Details

In formulating models, one sometimes wishes to define a state variable that will accumulate some quantity over the interval between successive observations. pomp provides a facility to make such features more convenient. Specifically, if $a$ is a state-variable named in the pomp's accumvars argument, then for each interval $(t_k,t_{k+1})$, $k=0,1,2,\dots$, $a$ will be set to zero at prior to any rprocess computation over that interval. Deterministic trajectory computation is handled slightly differently: see flow. For examples, see sir and the tutorials on the package website.

See Also

sir

More on implementing POMP models: Csnippet, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Examples

## A simple SIR model.

  ewmeas |>
    subset(time < 1952) |>
    pomp(
      times="time",t0=1948,
      rprocess=euler(
        Csnippet("
        int nrate = 6;
        double rate[nrate];     // transition rates
        double trans[nrate];    // transition numbers
        double dW;

        // gamma noise, mean=dt, variance=(sigma^2 dt)
        dW = rgammawn(sigma,dt);

        // compute the transition rates
        rate[0] = mu*pop;       // birth into susceptible class
        rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection
        rate[2] = mu;           // death from susceptible class
        rate[3] = gamma;        // recovery
        rate[4] = mu;           // death from infectious class
        rate[5] = mu;           // death from recovered class

        // compute the transition numbers
        trans[0] = rpois(rate[0]*dt);   // births are Poisson
        reulermultinom(2,S,&rate[1],dt,&trans[1]);
        reulermultinom(2,I,&rate[3],dt,&trans[3]);
        reulermultinom(1,R,&rate[5],dt,&trans[5]);

        // balance the equations
        S += trans[0]-trans[1]-trans[2];
        I += trans[1]-trans[3]-trans[4];
        R += trans[3]-trans[5];
      "),
      delta.t=1/52/20
      ),
      rinit=Csnippet("
        double m = pop/(S_0+I_0+R_0);
        S = nearbyint(m*S_0);
        I = nearbyint(m*I_0);
        R = nearbyint(m*R_0);
    "),
    paramnames=c("mu","pop","iota","gamma","Beta","sigma",
      "S_0","I_0","R_0"),
    statenames=c("S","I","R"),
    params=c(mu=1/50,iota=10,pop=50e6,gamma=26,Beta=400,sigma=0.1,
      S_0=0.07,I_0=0.001,R_0=0.93)
    ) -> ew1

  ew1 |>
    simulate() |>
    plot(variables=c("S","I","R"))

  ## A simple SIR model that tracks cumulative incidence.

  ew1 |>
    pomp(
      rprocess=euler(
        Csnippet("
        int nrate = 6;
        double rate[nrate];     // transition rates
        double trans[nrate];    // transition numbers
        double dW;

        // gamma noise, mean=dt, variance=(sigma^2 dt)
        dW = rgammawn(sigma,dt);

        // compute the transition rates
        rate[0] = mu*pop;       // birth into susceptible class
        rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection
        rate[2] = mu;           // death from susceptible class
        rate[3] = gamma;        // recovery
        rate[4] = mu;           // death from infectious class
        rate[5] = mu;           // death from recovered class

        // compute the transition numbers
        trans[0] = rpois(rate[0]*dt);   // births are Poisson
        reulermultinom(2,S,&rate[1],dt,&trans[1]);
        reulermultinom(2,I,&rate[3],dt,&trans[3]);
        reulermultinom(1,R,&rate[5],dt,&trans[5]);

        // balance the equations
        S += trans[0]-trans[1]-trans[2];
        I += trans[1]-trans[3]-trans[4];
        R += trans[3]-trans[5];
        H += trans[3];          // cumulative incidence
      "),
      delta.t=1/52/20
      ),
      rmeasure=Csnippet("
        double mean = H*rho;
        double size = 1/tau;
        reports = rnbinom_mu(size,mean);
    "),
    rinit=Csnippet("
        double m = pop/(S_0+I_0+R_0);
        S = nearbyint(m*S_0);
        I = nearbyint(m*I_0);
        R = nearbyint(m*R_0);
        H = 0;
    "),
    paramnames=c("mu","pop","iota","gamma","Beta","sigma","tau","rho",
      "S_0","I_0","R_0"),
    statenames=c("S","I","R","H"),
    params=c(mu=1/50,iota=10,pop=50e6,gamma=26,
      Beta=400,sigma=0.1,tau=0.001,rho=0.6,
      S_0=0.07,I_0=0.001,R_0=0.93)
    ) -> ew2

  ew2 |>
    simulate() |>
    plot()

  ## A simple SIR model that tracks weekly incidence.

  ew2 |>
    pomp(accumvars="H") -> ew3

  ew3 |>
    simulate() |>
    plot()

Description

Mathematically, the parts of a POMP model include the latent-state process transition distribution, the measurement-process distribution, the initial-state distribution, and possibly a prior parameter distribution. Algorithmically, each of these corresponds to at least two distinct operations. In particular, for each of the above parts, one sometimes needs to make a random draw from the distribution and sometimes to evaluate the density function. Accordingly, for each such component, there are two basic model components, one prefixed by a ‘r’, the other by a ‘d’, following the usual R convention.

Details

In addition to the parts listed above, pomp includes two additional basic model components: the deterministic skeleton, and parameter transformations that can be used to map the parameter space onto a Euclidean space for estimation purposes. There are also basic model components for computing the mean and variance of the measurement process conditional on the latent-state process.

There are thus altogether twelve basic model components:

1. rprocess, which samples from the latent-state transition distribution,

2. dprocess, which evaluates the latent-state transition density,

3. rmeasure, which samples from the measurement distribution,

4. emeasure, which computes the conditional expectation of the measurements, given the latent states,

5. vmeasure, which computes the conditional covariance matrix of the measurements, given the latent states,

6. dmeasure, which evaluates the measurement density,

7. rprior, which samples from the prior distribution,

8. dprior, which evaluates the prior density,

9. rinit, which samples from the initial-state distribution,

10. dinit, which evaluates the initial-state density,

11. skeleton, which evaluates the deterministic skeleton,

12. partrans, which evaluates the forward or inverse parameter transformations.

Each of these can be set or modified in the pomp constructor function or in any of the pomp elementary algorithms or estimation algorithms using an argument that matches the basic model component. A basic model component can be unset by passing NULL in the same way.

Help pages detailing each basic model component are provided.

See Also

workhorse functions, elementary algorithms, estimation algorithms.

More on implementing POMP models: Csnippet, accumvars, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

Several simple and configurable probes are provided with in the package. These can be used directly and as templates for custom probes.

Usage

probe_mean(var, trim = 0, transform = identity, na.rm = TRUE)

probe_median(var, na.rm = TRUE)

probe_var(var, transform = identity, na.rm = TRUE)

probe_sd(var, transform = identity, na.rm = TRUE)

probe_period(var, kernel.width, transform = identity)

probe_quantile(var, probs, ...)

probe_acf(
  var,
  lags,
  type = c("covariance", "correlation"),
  transform = identity
)

probe_ccf(
  vars,
  lags,
  type = c("covariance", "correlation"),
  transform = identity
)

probe_marginal(var, ref, order = 3, diff = 1, transform = identity)

probe_nlar(var, lags, powers, transform = identity)

Arguments

Value

A call to any one of these functions returns a probe function, suitable for use in probe or probe_objfun. That is, the function returned by each of these takes a data array (such as comes from a call to obs) as input and returns a single numerical value.

Author(s)

Daniel C. Reuman, Aaron A. King

References

B.E. Kendall, C.J. Briggs, W.W. Murdoch, P. Turchin, S.P. Ellner, E. McCauley, R.M. Nisbet, and S.N. Wood. Why do populations cycle? A synthesis of statistical and mechanistic modeling approaches. Ecology 80, 1789–1805, 1999. doi:10.2307/176658.

S. N. Wood Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102–1104, 2010. doi:10.1038/nature09319.

See Also

More on methods based on summary statistics: abc(), nlf, probe(), probe_match, spect(), spect_match

Description

Density and random generation for the Beta-binomial distribution with parameters size, mu, and theta.

Usage

rbetabinom(n = 1, size, prob, theta)

dbetabinom(x, size, prob, theta, log = FALSE)

Arguments

Details

A variable $X$ is Beta-binomially distributed if $X\sim{\mathrm{Binomial}(n,P)}$ where $P\sim{\mathrm{Beta}(\mu,\theta)}$. Using the standard $(a,b)$ parameterization, $a=\mu\,\theta$ and $b=(1-\mu)\,\theta$.

Value

C API

An interface for C codes using these functions is provided by the package. Visit the package homepage to view the pomp C API document.

See Also

More on implementing POMP models: Csnippet, accumvars, basic_components, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

blowflies is a data frame containing the data from several of Nicholson's classic experiments with the Australian sheep blowfly, Lucilia cuprina.

Usage

blowflies1(
  P = 3.2838,
  delta = 0.16073,
  N0 = 679.94,
  sigma.P = 1.3512,
  sigma.d = 0.74677,
  sigma.y = 0.026649
)

blowflies2(
  P = 2.7319,
  delta = 0.17377,
  N0 = 800.31,
  sigma.P = 1.442,
  sigma.d = 0.76033,
  sigma.y = 0.010846
)

Arguments

Details

blowflies1() and blowflies2() construct ‘pomp’ objects encoding stochastic delay-difference equation models. The data for these come from "population I", a control culture. The experiment is described on pp. 163–4 of Nicholson (1957). Unlimited quantities of larval food were provided; the adult food supply (ground liver) was constant at 0.4g per day. The data were taken from the table provided by Brillinger et al. (1980).

The models are discrete delay equations:

$R(t+1) \sim \mathrm{Poisson}(P N(t-\tau) \exp{(-N(t-\tau)/N_{0})} e(t+1) {\Delta}t)$

$S(t+1) \sim \mathrm{Binomial}(N(t),\exp{(-\delta \epsilon(t+1) {\Delta}t)})$

$N(t) = R(t)+S(t)$

where $e(t)$ and $\epsilon(t)$ are Gamma-distributed i.i.d. random variables with mean 1 and variances ${\sigma_P^2}/{{\Delta}t}$, ${\sigma_d^2}/{{\Delta}t}$, respectively. blowflies1 has a timestep (${\Delta}t$) of 1 day; blowflies2 has a timestep of 2 days. The process model in blowflies1 thus corresponds exactly to that studied by Wood (2010). The measurement model in both cases is taken to be

$y(t) \sim \mathrm{NegBin}(N(t),1/\sigma_y^2)$

i.e., the observations are assumed to be negative-binomially distributed with mean $N(t)$ and variance $N(t)+(\sigma_y N(t))^2$.

Default parameter values are the MLEs as estimated by Ionides (2011).

Value

blowflies1 and blowflies2 return ‘pomp’ objects containing the actual data and two variants of the model.

References

A.J. Nicholson. The self-adjustment of populations to change. Cold Spring Harbor Symposia on Quantitative Biology 22, 153–173, 1957. doi:10.1101/SQB.1957.022.01.017.

Y. Xia and H. Tong. Feature matching in time series modeling. Statistical Science 26, 21–46, 2011. doi:10.1214/10-sts345.

E.L. Ionides. Discussion of “Feature matching in time series modeling” by Y. Xia and H. Tong. Statistical Science 26, 49–52, 2011. doi:10.1214/11-sts345c.

S. N. Wood Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102–1104, 2010. doi:10.1038/nature09319.

W.S.C. Gurney, S.P. Blythe, and R.M. Nisbet. Nicholson's blowflies revisited. Nature 287, 17–21, 1980. doi:10.1038/287017a0.

D.R. Brillinger, J. Guckenheimer, P. Guttorp, and G. Oster. Empirical modelling of population time series: The case of age and density dependent rates. In: G. Oster (ed.), Some Questions in Mathematical Biology vol. 13, pp. 65–90, American Mathematical Society, Providence, 1980. doi:10.1007/978-1-4614-1344-8_19.

See Also

More examples provided with pomp: childhood_disease_data, compartmental_models, dacca(), ebola, gompertz(), ou2(), pomp_examples, ricker(), rw2(), verhulst()

More data sets provided with pomp: bsflu, childhood_disease_data, dacca(), ebola, parus

Examples

plot(blowflies1())
plot(blowflies2())

Description

An outbreak of influenza in an all-boys boarding school.

Details

Data are recorded from a 1978 flu outbreak in a closed population. The variable ‘B’ refers to boys confined to bed on the corresponding day and ‘C’ to boys in convalescence, i.e., not yet allowed back to class. In total, 763 boys were at risk of infection and, over the course of the outbreak, 512 boys spent between 3 and 7 days away from class (either in bed or convalescent). The index case was a boy who arrived at school from holiday six days before the next case.

References

Anonymous. Influenza in a boarding school. British Medical Journal 1, 587, 1978.

See Also

compartmental models

More data sets provided with pomp: blowflies, childhood_disease_data, dacca(), ebola, parus

Examples

if (require(tidyr) && require(ggplot2)) {

  bsflu |>
    gather(variable,value,-date,-day) |>
    ggplot(aes(x=date,y=value,color=variable))+
    geom_line()+
    labs(y="number of boys",title="boarding school flu outbreak")+
    theme_bw()

}

Description

Modified version of the Liu & West (2001) algorithm.

Usage

## S4 method for signature 'data.frame'
bsmc2(
  data,
  Np,
  smooth = 0.1,
  params,
  rprior,
  rinit,
  rprocess,
  dmeasure,
  partrans,
  ...,
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'pomp'
bsmc2(data, Np, smooth = 0.1, ..., verbose = getOption("verbose", FALSE))

Arguments

length(time(object,t0=TRUE))

Details

bsmc2 uses a version of the original algorithm (Liu & West 2001), but discards the auxiliary particle filter. The modification appears to give superior performance for the same amount of effort.

Samples from the prior distribution are drawn using the rprior component. This is allowed to depend on elements of params, i.e., some of the elements of params can be treated as “hyperparameters”. Np draws are made from the prior distribution.

Value

An object of class ‘bsmcd_pomp’. The following methods are avaiable:

plot

produces diagnostic plots

as.data.frame

puts the prior and posterior samples into a data frame

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

Author(s)

Michael Lavine, Matthew Ferrari, Aaron A. King, Edward L. Ionides

References

J. Liu and M. West. Combining parameter and state estimation in simulation-based filtering. In A. Doucet, N. de Freitas, and N. J. Gordon, (eds.), Sequential Monte Carlo Methods in Practice, pp. 197–224. Springer, New York, 2001. doi:10.1007/978-1-4757-3437-9_10.

See Also

More on Bayesian methods: abc(), dprior(), pmcmc(), prior_spec, rprior()

More on full-information (i.e., likelihood-based) methods: mif2(), pfilter(), pmcmc(), wpfilter()

More on sequential Monte Carlo methods: cond_logLik(), eff_sample_size(), filter_mean(), filter_traj(), kalman, mif2(), pfilter(), pmcmc(), pred_mean(), pred_var(), saved_states(), wpfilter()

More on pomp estimation algorithms: abc(), estimation_algorithms, mif2(), nlf, pmcmc(), pomp-package, probe_match, spect_match

Description

These functions generate B-spline basis functions. bspline_basis gives a basis of spline functions. periodic_bspline_basis gives a basis of periodic spline functions.

Usage

bspline_basis(x, nbasis, degree = 3, deriv = 0, names = NULL, rg = range(x))

periodic_bspline_basis(
  x,
  nbasis,
  degree = 3,
  period = 1,
  deriv = 0,
  names = NULL
)

Arguments

Value

If deriv>0, the derivative of that order of each of the corresponding spline basis functions are returned.

C API

Access to the underlying C routines is available: see the pomp C API document. for definition and documentation of the C API.

Author(s)

Aaron A. King

See Also

More on interpolation: covariates, lookup()

Examples

x <- seq(0,2,by=0.01)
y <- bspline_basis(x,degree=3,nbasis=9,names="basis")
matplot(x,y,type='l',ylim=c(0,1.1))
lines(x,apply(y,1,sum),lwd=2)

x <- seq(-1,2,by=0.01)
y <- periodic_bspline_basis(x,nbasis=5,names="spline%d")
matplot(x,y,type='l')

Description

LondonYorke is a data frame containing the monthly number of reported cases of chickenpox, measles, and mumps from two American cities (Baltimore and New York) in the mid-20th century (1928–1972).

ewmeas and ewcitmeas are data frames containing weekly reported cases of measles in England and Wales. ewmeas records the total measles reports for the whole country, 1948–1966. One questionable data point has been replaced with an NA. ewcitmeas records the incidence in seven English cities 1948–1987. These data were kindly provided by Ben Bolker, who writes: “Most of these data have been manually entered from published records by various people, and are prone to errors at several levels. All data are provided as is; use at your own risk.”

References

W. P. London and J. A. Yorke, Recurrent outbreaks of measles, chickenpox and mumps: I. Seasonal variation in contact rates. American Journal of Epidemiology 98, 453–468, 1973. doi:10.1093/oxfordjournals.aje.a121575.

See Also

compartmental models, bsflu

More data sets provided with pomp: blowflies, bsflu, dacca(), ebola, parus

More examples provided with pomp: blowflies, compartmental_models, dacca(), ebola, gompertz(), ou2(), pomp_examples, ricker(), rw2(), verhulst()

Examples

plot(cases~time,data=LondonYorke,subset=disease=="measles",type='n',main="measles",bty='l')
lines(cases~time,data=LondonYorke,subset=disease=="measles"&town=="Baltimore",col="red")
lines(cases~time,data=LondonYorke,subset=disease=="measles"&town=="New York",col="blue")
legend("topright",legend=c("Baltimore","New York"),lty=1,col=c("red","blue"),bty='n')

plot(
     cases~time,
     data=LondonYorke,
     subset=disease=="chickenpox"&town=="New York",
     type='l',col="blue",main="chickenpox, New York",
     bty='l'
    )

plot(
     cases~time,
     data=LondonYorke,
     subset=disease=="mumps"&town=="New York",
     type='l',col="blue",main="mumps, New York",
     bty='l'
    )

plot(reports~time,data=ewmeas,type='l')

plot(reports~date,data=ewcitmeas,subset=city=="Liverpool",type='l')

Description

Extract, set, or modify the estimated parameters from a fitted model.

Usage

## S4 method for signature 'listie'
coef(object, ...)

## S4 method for signature 'pomp'
coef(object, pars, transform = FALSE, ...)

## S4 replacement method for signature 'pomp'
coef(object, pars, transform = FALSE, ...) <- value

## S4 method for signature 'objfun'
coef(object, ...)

## S4 replacement method for signature 'objfun'
coef(object, pars, transform = FALSE, ...) <- value

Arguments

Details

coef allows one to extract the parameters from a fitted model.

coef(object,transform=TRUE) returns the parameters transformed onto the estimation scale.

coef(object) <- value sets or alters the coefficients of a ‘pomp’ object.

coef(object,transform=TRUE) <- value assumes that value is on the estimation scale, and applies the “from estimation scale” parameter transformation from object before altering the coefficients.

See Also

Other extraction methods: cond_logLik(), covmat(), eff_sample_size(), filter_mean(), filter_traj(), forecast(), logLik, obs(), pred_mean(), pred_var(), saved_states(), spy(), states(), summary(), time(), timezero(), traces()

Description

Simple SIR-type models implemented in various ways.

Usage

sir(
  gamma = 26,
  mu = 0.02,
  iota = 0.01,
  beta1 = 400,
  beta2 = 480,
  beta3 = 320,
  beta_sd = 0.001,
  rho = 0.6,
  k = 0.1,
  pop = 2100000,
  S_0 = 26/400,
  I_0 = 0.001,
  R_0 = 1 - S_0 - I_0,
  t0 = 0,
  times = seq(from = t0 + 1/52, to = t0 + 4, by = 1/52),
  seed = 329343545,
  delta.t = 1/52/20
)

sir2(
  gamma = 24,
  mu = 1/70,
  iota = 0.1,
  beta1 = 330,
  beta2 = 410,
  beta3 = 490,
  rho = 0.1,
  k = 0.1,
  pop = 1e+06,
  S_0 = 0.05,
  I_0 = 1e-04,
  R_0 = 1 - S_0 - I_0,
  t0 = 0,
  times = seq(from = t0 + 1/12, to = t0 + 10, by = 1/12),
  seed = 1772464524
)

Arguments

Details

sir() producees a ‘pomp’ object encoding a simple seasonal SIR model with simulated data. Simulation is performed using an Euler multinomial approximation.

sir2() has the same model implemented using Gillespie's algorithm.

In both cases the measurement model is negative binomial: reports is distributed as a negative binomial random variable with mean equal to rho*cases and size equal to 1/k.

This and similar examples are discussed and constructed in tutorials available on the package website.

Value

These functions return ‘pomp’ objects containing simulated data.

See Also

More examples provided with pomp: blowflies, childhood_disease_data, dacca(), ebola, gompertz(), ou2(), pomp_examples, ricker(), rw2(), verhulst()

Examples

po <- sir()
  plot(po)
  coef(po)
  
  po <- sir2()
  plot(po)
  plot(simulate(window(po,end=3)))
  coef(po)
  
  po |> as.data.frame() |> head()

Description

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

Usage

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

concat(...)

Arguments

Details

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

Examples

gompertz(sigma=2,tau=1) -> g
Np <- c(low=100,med=1000,high=10000)
lapply(
  Np,
  \(np) pfilter(g,Np=np)
) |>
  concat() -> pfs

pfs
coef(pfs)
logLik(pfs)
eff_sample_size(pfs)
cond_logLik(pfs)

pfs |> plot()

Description

The estimated conditional log likelihood from a fitted model.

Usage

## S4 method for signature 'kalmand_pomp'
cond_logLik(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'pfilterd_pomp'
cond_logLik(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'wpfilterd_pomp'
cond_logLik(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'bsmcd_pomp'
cond_logLik(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'pfilterList'
cond_logLik(object, ..., format = c("numeric", "data.frame"))

Arguments

Details

The conditional likelihood is defined to be the value of the density of

$Y(t_k) | Y(t_1),\dots,Y(t_{k-1})$

evaluated at $Y(t_k) = y^*_k$. Here, $Y(t_k)$ is the observable process, and $y^*_k$ the data, at time $t_k$.

Thus the conditional log likelihood at time $t_k$ is

$\ell_k(\theta) = \log f[Y(t_k)=y^*_k \vert Y(t_1)=y^*_1, \dots, Y(t_{k-1})=y^*_{k-1}],$

where $f$ is the probability density above.

Value

The numerical value of the conditional log likelihood. Note that some methods compute not the log likelihood itself but instead a related quantity. To keep the code simple, the cond_logLik function is nevertheless used to extract this quantity.

When object is of class ‘bsmcd_pomp’ (i.e., the result of a bsmc2 computation), cond_logLik returns the conditional log “evidence” (see bsmc2).

See Also

More on sequential Monte Carlo methods: bsmc2(), eff_sample_size(), filter_mean(), filter_traj(), kalman, mif2(), pfilter(), pmcmc(), pred_mean(), pred_var(), saved_states(), wpfilter()

Other extraction methods: coef(), covmat(), eff_sample_size(), filter_mean(), filter_traj(), forecast(), logLik, obs(), pred_mean(), pred_var(), saved_states(), spy(), states(), summary(), time(), timezero(), traces()

Description

Continue an iterative computation where it left off.

Usage

## S4 method for signature 'abcd_pomp'
continue(object, Nabc = 1, ...)

## S4 method for signature 'pmcmcd_pomp'
continue(object, Nmcmc = 1, ...)

## S4 method for signature 'mif2d_pomp'
continue(object, Nmif = 1, ...)

Arguments

See Also

mif2 pmcmc abc

Description

Incorporating time-varying covariates using lookup tables.

Usage

## S4 method for signature 'numeric'
covariate_table(..., order = c("linear", "constant"), times)

## S4 method for signature 'character'
covariate_table(..., order = c("linear", "constant"), times)

repair_lookup_table(table, t, order)

Arguments

Details

If the ‘pomp’ object contains covariates (specified via the covar argument), then interpolated values of the covariates will be available to each of the model components whenever it is called. In particular, variables with names as they appear in the covar covariate table will be available to any C snippet. When a basic component is defined using an R function, that function will be called with an extra argument, covars, which will be a named numeric vector containing the interpolated values from the covariate table.

An exception to this rule is the prior (rprior and dprior): covariate-dependent priors are not allowed. Nor are parameter transformations permitted to depend upon covariates.

repair_lookup_table applies lookup at the provided values of t and returns the resulting lookup table. If order is unsupplied, the interpolation-order from table is preserved. repair_lookup_table should be considered experimental: its interface may change without notice.

Value

covariate_table returns a lookup table suitable for inclusion of covariates in a ‘pomp’ object. Specifically, this is an object of class ‘covartable’.

repair_lookup_table returns a lookup table with entries at the provided values of t.

Extrapolation

If t is outside the range of the lookup table, the values will be extrapolated, and a warning will be issued. The type of extrapolation performed will be constant or linear according to the order flag used when creating the table.

See Also

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

More on interpolation: bsplines, lookup()

Description

A helper function to extract a covariance matrix.

Usage

## S4 method for signature 'pmcmcd_pomp'
covmat(object, start = 1, thin = 1, expand = 2.38, ...)

## S4 method for signature 'pmcmcList'
covmat(object, start = 1, thin = 1, expand = 2.38, ...)

## S4 method for signature 'abcd_pomp'
covmat(object, start = 1, thin = 1, expand = 2.38, ...)

## S4 method for signature 'abcList'
covmat(object, start = 1, thin = 1, expand = 2.38, ...)

## S4 method for signature 'probed_pomp'
covmat(object, ...)

Arguments

Value

When object is the result of a pmcmc or abc computation, covmat(object) gives the covariance matrix of the chains. This can be useful, for example, in tuning the proposal distribution.

When object is a ‘probed_pomp’ object (i.e., the result of a probe computation), covmat(object) returns the covariance matrix of the probes, as applied to simulated data.

See Also

MCMC proposals.

Other extraction methods: coef(), cond_logLik(), eff_sample_size(), filter_mean(), filter_traj(), forecast(), logLik, obs(), pred_mean(), pred_var(), saved_states(), spy(), states(), summary(), time(), timezero(), traces()

Description

Accelerating computations through inline snippets of C code

Usage

Csnippet(text)

Arguments

Details

pomp provides a facility whereby users can define their model's components using inline C code. C snippets are written to a C file, by default located in the R session's temporary directory, which is then compiled (via R CMD SHLIB) into a dynamically loadable shared object file. This is then loaded as needed.

Note to Windows and Mac users

By default, your R installation may not support R CMD SHLIB. The package website contains installation instructions that explain how to enable this powerful feature of R.

General rules for writing C snippets

In writing a C snippet one must bear in mind both the goal of the snippet, i.e., what computation it is intended to perform, and the context in which it will be executed. These are explained here in the form of general rules. Additional specific rules apply according to the function of the particular C snippet. Illustrative examples are given in the tutorials on the package website.

1. C snippets must be valid C. They will embedded verbatim in a template file which will then be compiled by a call to R CMD SHLIB. If the resulting file does not compile, an error message will be generated. Compiler messages will be displayed, but no attempt will be made by pomp to interpret them. Typically, compilation errors are due to either invalid C syntax or undeclared variables.

2. State variables, parameters, observables, and covariates must be left undeclared within the snippet. State variables and parameters are declared via the statenames or paramnames arguments to pomp, respectively. Compiler errors that complain about undeclared state variables or parameters are usually due to failure to declare these in statenames or paramnames, as appropriate.

3. A C snippet can declare local variables. Be careful not to use names that match those of state variables, observables, or parameters. One must never declare state variables, observables, covariates, or parameters within a C snippet.

4. Names of observables must match the names given given in the data. They must be referred to in measurement model C snippets (rmeasure and dmeasure) by those names.

5. If the ‘pomp’ object contains a table of covariates (see above), then the variables in the covariate table will be available, by their names, in the context within which the C snippet is executed.

6. Because the dot ‘.’ has syntactic meaning in C, R variables with names containing dots (‘.’) are replaced in the C codes by variable names in which all dots have been replaced by underscores (‘_’).

7. The headers ‘R.h’ and ‘Rmath.h’, provided with R, will be included in the generated C file, making all of the R C API available for use in the C snippet. This makes a great many useful functions available, including all of R's statistical distribution functions.

8. The header ‘pomp.h’, provided with pomp, will also be included, making all of the pomp C API available for use in every C snippet.

9. Snippets of C code passed to the globals argument of pomp will be included at the head of the generated C file. This can be used to declare global variables, define useful functions, and include arbitrary header files.

Linking to precompiled libraries

It is straightforward to link C snippets with precompiled C libraries. To do so, one must make sure the library's header files are included; the globals argument can be used for this purpose. The shlib.args argument can then be used to specify additional arguments to be passed to R CMD SHLIB. FAQ 3.7 gives an example.

C snippets are salted

To prevent collisions in parallel computations, a ‘pomp’ object built using C snippets is “salted” with the current time and a random number. A result is that two ‘pomp’ objects, built on identical codes and data, will not be identical as R objects, though they will be functionally identical in every respect.

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

spy

More on implementing POMP models: accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

dacca constructs a ‘pomp’ object containing census and cholera mortality data from the Dacca district of the former British province of Bengal over the years 1891 to 1940 together with a stochastic differential equation transmission model. The model is that of King et al. (2008). The parameters are the MLE for the SIRS model with seasonal reservoir.

Usage

dacca(
  gamma = 20.8,
  eps = 19.1,
  rho = 0,
  delta = 0.02,
  deltaI = 0.06,
  clin = 1,
  alpha = 1,
  beta_trend = -0.00498,
  logbeta = c(0.747, 6.38, -3.44, 4.23, 3.33, 4.55),
  logomega = log(c(0.184, 0.0786, 0.0584, 0.00917, 0.000208, 0.0124)),
  sd_beta = 3.13,
  tau = 0.23,
  S_0 = 0.621,
  I_0 = 0.378,
  Y_0 = 0,
  R1_0 = 0.000843,
  R2_0 = 0.000972,
  R3_0 = 1.16e-07
)

Arguments

Details

Data are provided courtesy of Dr. Menno J. Bouma, London School of Tropical Medicine and Hygiene.

Value

dacca returns a ‘pomp’ object containing the model, data, and MLE parameters, as estimated by King et al. (2008).

References

A.A. King, E.L. Ionides, M. Pascual, and M.J. Bouma. Inapparent infections and cholera dynamics. Nature 454, 877-880, 2008. doi:10.1038/nature07084.

See Also

More examples provided with pomp: blowflies, childhood_disease_data, compartmental_models, ebola, gompertz(), ou2(), pomp_examples, ricker(), rw2(), verhulst()

More data sets provided with pomp: blowflies, bsflu, childhood_disease_data, ebola, parus

Examples

# takes too long for R CMD check
  po <- dacca()
  plot(po)
  ## MLE:
  coef(po)
  plot(simulate(po))

Description

These functions are useful for generating designs for the exploration of parameter space.

profile_design generates a data-frame where each row can be used as the starting point for a profile likelihood calculation.

runif_design generates a design based on random samples from a multivariate uniform distribution.

slice_design generates points along slices through a specified point.

sobol_design generates a Latin hypercube design based on the Sobol' low-discrepancy sequence.

Usage

profile_design(
  ...,
  lower,
  upper,
  nprof,
  type = c("runif", "sobol"),
  stringsAsFactors = getOption("stringsAsFactors", FALSE)
)

runif_design(lower = numeric(0), upper = numeric(0), nseq)

slice_design(center, ...)

sobol_design(lower = numeric(0), upper = numeric(0), nseq)

Arguments

Details

The Sobol' sequence generation is performed using codes from the NLopt library by S. Johnson.

Value

profile_design returns a data frame with nprof points per profile point.

runif_design returns a data frame with nseq rows and one column for each variable named in lower and upper.

slice_design returns a data frame with one row per point. The ‘slice’ variable indicates which slice the point belongs to.

sobol_design returns a data frame with nseq rows and one column for each variable named in lower and upper.

Author(s)

Aaron A. King

References

S. Kucherenko and Y. Sytsko. Application of deterministic low-discrepancy sequences in global optimization. Computational Optimization and Applications 30, 297–318, 2005. doi:10.1007/s10589-005-4615-1.

S.G. Johnson. The NLopt nonlinear-optimization package. https://github.com/stevengj/nlopt/.

P. Bratley and B.L. Fox. Algorithm 659 Implementing Sobol's quasirandom sequence generator. ACM Transactions on Mathematical Software 14, 88–100, 1988. doi:10.1145/42288.214372.

S. Joe and F.Y. Kuo. Remark on algorithm 659: Implementing Sobol' quasirandom sequence generator. ACM Transactions on Mathematical Software 29, 49–57, 2003. doi:10.1145/641876.641879.

Examples

## Sobol' low-discrepancy design
plot(sobol_design(lower=c(a=0,b=100),upper=c(b=200,a=1),nseq=100))

## Uniform random design
plot(runif_design(lower=c(a=0,b=100),upper=c(b=200,a=1),100))

## A one-parameter profile design:
x <- profile_design(p=1:10,lower=c(a=0,b=0),upper=c(a=1,b=5),nprof=20)
dim(x)
plot(x)

## A two-parameter profile design:
x <- profile_design(p=1:10,q=3:5,lower=c(a=0,b=0),upper=c(b=5,a=1),nprof=200)
dim(x)
plot(x)

## A two-parameter profile design with random points:
x <- profile_design(p=1:10,q=3:5,lower=c(a=0,b=0),upper=c(b=5,a=1),nprof=200,type="runif")
dim(x)
plot(x)

## A single 11-point slice through the point c(A=3,B=8,C=0) along the B direction.
x <- slice_design(center=c(A=3,B=8,C=0),B=seq(0,10,by=1))
dim(x)
plot(x)

## Two slices through the same point along the A and C directions.
x <- slice_design(c(A=3,B=8,C=0),A=seq(0,5,by=1),C=seq(0,5,length=11))
dim(x)
plot(x)

Description

Evaluates the initial-state density.

Usage

## S4 method for signature 'pomp'
dinit(
  object,
  params = coef(object),
  t0 = timezero(object),
  x,
  log = FALSE,
  ...
)

Arguments

Value

dinit returns a 1-D numerical array containing the likelihoods (or log likelihoods if log=TRUE). By default, t0 is the initial time defined when the ‘pomp’ object ws constructed.

See Also

Specification of the initial-state distribution: dinit_spec

More on pomp workhorse functions: dmeasure(), dprior(), dprocess(), emeasure(), flow(), partrans(), pomp-package, rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

Description

Specification of the initial-state distribution density evaluator, dinit.

Details

To fully specify the unobserved Markov state process, one must give its distribution at the zero-time (t0). One specifies how to evaluate the log probability density function for this distribution using the dinit argument. As usual, this can be provided either as a C snippet or as an R function. In the former case, bear in mind that:

1. The goal of a this snippet is computation of a log likelihood, to be put into a variable named loglik.

2. In addition to the state variables, parameters, and covariates (if any), the variable t, containing the zero-time, will be defined in the context in which the snippet is executed.

General rules for writing C snippets can be found here.

If an R function is to be used, pass

   dinit = f

to pomp, where f is a function with arguments that can include the time t, any or all of the model state variables, parameters, and covariates. As usual, f may take additional arguments, provided these are passed along with it in the call to pomp. f must return a single numeric value, the log likelihood.

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

dinit

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

dmeasure evaluates the probability density of observations given states.

Usage

## S4 method for signature 'pomp'
dmeasure(
  object,
  y = obs(object),
  x = states(object),
  times = time(object),
  params = coef(object),
  ...,
  log = FALSE
)

Arguments

Value

dmeasure returns a matrix of dimensions nreps x ntimes. If d is the returned matrix, d[j,k] is the likelihood (or log likelihood if log = TRUE) of the observation y[,k] at time times[k] given the state x[,j,k].

See Also

Specification of the measurement density evaluator: dmeasure_spec

More on pomp workhorse functions: dinit(), dprior(), dprocess(), emeasure(), flow(), partrans(), pomp-package, rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

Description

Specification of the measurement model density function, dmeasure.

Details

The measurement model is the link between the data and the unobserved state process. It can be specified either by using one or both of the rmeasure and dmeasure arguments.

Suppose you have a procedure to compute the probability density of an observation given the value of the latent state variables. Then you can furnish

   dmeasure = f

to pomp algorithms, where f is a C snippet or R function that implements your procedure.

Using a C snippet is much preferred, due to its much greater computational efficiency. See Csnippet for general rules on writing C snippets. The goal of a dmeasure C snippet is to fill the variable lik with the either the probability density or the log probability density, depending on the value of the variable give_log.

In writing a dmeasure C snippet, observe that:

1. In addition to the states, parameters, covariates (if any), and observables, the variable t, containing the time of the observation will be defined in the context in which the snippet is executed.

2. Moreover, the Boolean variable give_log will be defined.

3. The goal of a dmeasure C snippet is to set the value of the lik variable to the likelihood of the data given the state, if give_log == 0. If give_log == 1, lik should be set to the log likelihood.

If dmeasure is to be provided instead as an R function, this is accomplished by supplying

  dmeasure = f

to pomp, where f is a function. The arguments of f should be chosen from among the observables, state variables, parameters, covariates, and time. It must also have the arguments ..., and log. It can take additional arguments via the userdata facility. f must return a single numeric value, the probability density (or log probability density if log = TRUE) of y given x at time t.

Important note

It is a common error to fail to account for both log = TRUE and log = FALSE when writing the dmeasure C snippet or function.

Default behavior

If dmeasure is left unspecified, calls to dmeasure will return missing values (NA).

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

dmeasure

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Examples

## We start with the pre-built Ricker example:
  ricker() -> po

  ## To change the measurement model density, dmeasure,
  ## we use the 'dmeasure' argument in any 'pomp'
  ## elementary or estimation function.
  ## Here, we pass the dmeasure specification to 'pfilter'
  ## as an R function.

  po |>
    pfilter(
      dmeasure=function (y, N, phi, ..., log) {
        dpois(y,lambda=phi*N,log=log)
      },
      Np=100
    ) -> pf

  ## We can also pass it as a C snippet:

  po |>
    pfilter(
      dmeasure=Csnippet("lik = dpois(y,phi*N,give_log);"),
      paramnames="phi",
      statenames="N",
      Np=100
    ) -> pf

Description

Evaluates the prior probability density.

Usage

## S4 method for signature 'pomp'
dprior(object, params = coef(object), ..., log = FALSE)

Arguments

Value

The required density (or log density), as a numeric vector.

See Also

Specification of the prior density evaluator: prior_spec

More on pomp workhorse functions: dinit(), dmeasure(), dprocess(), emeasure(), flow(), partrans(), pomp-package, rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

More on Bayesian methods: abc(), bsmc2(), pmcmc(), prior_spec, rprior()

Description

Evaluates the probability density of a sequence of consecutive state transitions.

Usage

## S4 method for signature 'pomp'
dprocess(
  object,
  x = states(object),
  times = time(object),
  params = coef(object),
  ...,
  log = FALSE
)

Arguments

Value

dprocess returns a matrix of dimensions nrep x ntimes-1. If d is the returned matrix, d[j,k] is the likelihood (or the log likelihood if log=TRUE) of the transition from state x[,j,k-1] at time times[k-1] to state x[,j,k] at time times[k].

See Also

Specification of the process-model density evaluator: dprocess_spec

More on pomp workhorse functions: dinit(), dmeasure(), dprior(), emeasure(), flow(), partrans(), pomp-package, rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

Description

Specification of the latent state process density function, dprocess.

Details

Suppose you have a procedure that allows you to compute the probability density of an arbitrary transition from state $x_1$ at time $t_1$ to state $x_2$ at time $t_2>t_1$ under the assumption that the state remains unchanged between $t_1$ and $t_2$. Then you can furnish

    dprocess = f

to pomp, where f is a C snippet or R function that implements your procedure. Specifically, f should compute the log probability density.

Using a C snippet is much preferred, due to its much greater computational efficiency. See Csnippet for general rules on writing C snippets. The goal of a dprocess C snippet is to fill the variable loglik with the log probability density. In the context of such a C snippet, the parameters, and covariates will be defined, as will the times t_1 and t_2. The state variables at time t_1 will have their usual name (see statenames) with a “_1” appended. Likewise, the state variables at time t_2 will have a “_2” appended.

If f is given as an R function, it should take as arguments any or all of the state variables, parameter, covariates, and time. The state-variable and time arguments will have suffices “_1” and “_2” appended. Thus for example, if var is a state variable, when f is called, var_1 will value of state variable var at time t_1, var_2 will have the value of var at time t_2. f should return the log likelihood of a transition from x1 at time t1 to x2 at time t2, assuming that no intervening transitions have occurred.

To see examples, consult the demos and the tutorials on the package website.

Note

It is not typically necessary (or even feasible) to define dprocess. In fact, no current pomp inference algorithm makes use of dprocess. This functionality is provided only to support future algorithm development.

Default behavior

By default, dprocess returns missing values (NA).

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

dprocess

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

Data and models for the 2014–2016 outbreak of Ebola virus disease in West Africa.

Usage

ebolaModel(
  country = c("GIN", "LBR", "SLE"),
  data = NULL,
  timestep = 1/8,
  nstageE = 3L,
  R0 = 1.4,
  rho = 0.2,
  cfr = 0.7,
  k = 0,
  index_case = 10,
  incubation_period = 11.4,
  infectious_period = 7
)

Arguments

Details

The data include monthly case counts and death reports derived from WHO situation reports, as reported by the U.S. CDC. The models are described in King et al. (2015).

The data-cleaning script is included in the R source code file ‘ebola.R’.

Model structure

The default incubation period is supposed to be Gamma distributed with shape parameter nstageE and mean 11.4 days and the case-fatality ratio ('cfr') is taken to be 0.7 (cf. WHO Ebola Response Team 2014). The discrete-time formula is used to calculate the corresponding alpha (cf. He et al. 2010).

The observation model is a hierarchical model for cases and deaths:

$p(R_t, D_t| C_t) = p(R_t | C_t) p(D_t | C_t, R_t).$

Here, $p(R_t | C_t)$ is negative binomial with mean $\rho C_t$ and dispersion parameter $1/k$; $p(D_t | C_t, R_t)$ is binomial with size $R_t$ and probability equal to the case fatality rate cfr.

References

A.A. King, M. Domenech de Cellès, F.M.G. Magpantay, and P. Rohani. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to Ebola. Proceedings of the Royal Society of London, Series B 282, 20150347, 2015. doi:10.1098/rspb.2015.0347.

WHO Ebola Response Team. Ebola virus disease in West Africa—the first 9 months of the epidemic and forward projections. New England Journal of Medicine 371, 1481–1495, 2014. doi:10.1056/NEJMoa1411100.

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

See Also

More data sets provided with pomp: blowflies, bsflu, childhood_disease_data, dacca(), parus

More examples provided with pomp: blowflies, childhood_disease_data, compartmental_models, dacca(), gompertz(), ou2(), pomp_examples, ricker(), rw2(), verhulst()

Examples

# takes too long for R CMD check
  if (require(ggplot2) && require(tidyr)) {
    
    ebolaWA2014 |>
      pivot_longer(c(cases,deaths)) |>
      ggplot(aes(x=date,y=value,group=name,color=name))+
      geom_line()+
      facet_grid(country~.,scales="free_y")+
      theme_bw()+
      theme(axis.text=element_text(angle=-90))
    
  }
  
  plot(ebolaModel(country="SLE"))
  plot(ebolaModel(country="GIN"))
  plot(ebolaModel(country="LBR"))

Description

Estimate the effective sample size of a Monte Carlo computation.

Usage

## S4 method for signature 'bsmcd_pomp'
eff_sample_size(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'pfilterd_pomp'
eff_sample_size(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'wpfilterd_pomp'
eff_sample_size(object, ..., format = c("numeric", "data.frame"))

## S4 method for signature 'pfilterList'
eff_sample_size(object, ..., format = c("numeric", "data.frame"))

Arguments

Details

Effective sample size is computed as

$\left(\sum_i\!w_{it}^2\right)^{-1},$

where $w_{it}$ is the normalized weight of particle $i$ at time $t$.

See Also

More on sequential Monte Carlo methods: bsmc2(), cond_logLik(), filter_mean(), filter_traj(), kalman, mif2(), pfilter(), pmcmc(), pred_mean(), pred_var(), saved_states(), wpfilter()

Other extraction methods: coef(), cond_logLik(), covmat(), filter_mean(), filter_traj(), forecast(), logLik, obs(), pred_mean(), pred_var(), saved_states(), spy(), states(), summary(), time(), timezero(), traces()

Description

In pomp, elementary algorithms perform POMP model operations. These operations do not themselves estimate parameters, though they may be instrumental in inference methods.

Details

There are six elementary algorithms in pomp:

• simulate which simulates from the joint distribution of latent and observed variables,

• pfilter, which performs a simple particle filter operation,

• wpfilter, which performs a weighted particle filter operation,

• probe, which computes a suite of user-specified summary statistics on actual and simulated data,

• spect, which performs a power-spectral density function computation on actual and simulated data,

• trajectory, which iterates or integrates the deterministic skeleton according to whether the latter is a (discrete-time) map or a (continuous-time) vectorfield.

Help pages detailing each elementary algorithm component are provided.

See Also

basic model components, workhorse functions, estimation algorithms.

More on pomp elementary algorithms: kalman, pfilter(), pomp-package, probe(), simulate(), spect(), trajectory(), wpfilter()

Description

Return the expected value of the observed variables, given values of the latent states and the parameters.

Usage

## S4 method for signature 'pomp'
emeasure(
  object,
  x = states(object),
  times = time(object),
  params = coef(object),
  ...
)

Arguments

Value

emeasure returns a rank-3 array of dimensions nobs x nrep x ntimes, where nobs is the number of observed variables.

See Also

Specification of the measurement-model expectation: emeasure_spec

More on pomp workhorse functions: dinit(), dmeasure(), dprior(), dprocess(), flow(), partrans(), pomp-package, rinit(), rmeasure(), rprior(), rprocess(), skeleton(), vmeasure(), workhorses

Description

Specification of the measurement-model conditional expectation, emeasure.

Details

The measurement model is the link between the data and the unobserved state process. Some algorithms require the conditional expectation of the measurement model, given the latent state and parameters. This is supplied using the emeasure argument.

Suppose you have a procedure to compute this conditional expectation, given the value of the latent state variables. Then you can furnish

  emeasure = f

to pomp algorithms, where f is a C snippet or R function that implements your procedure.

Using a C snippet is much preferred, due to its much greater computational efficiency. See Csnippet for general rules on writing C snippets.

In writing an emeasure C snippet, bear in mind that:

1. The goal of such a snippet is to fill variables named E_y with the conditional expectations of observables y. Accordingly, there should be one assignment of E_y for each observable y.

2. In addition to the states, parameters, and covariates (if any), the variable t, containing the time of the observation, will be defined in the context in which the snippet is executed.

The demos and the tutorials on the package website give examples.

It is also possible, though less efficient, to specify emeasure using an R function. In this case, specify the measurement model expectation by furnishing

  emeasure = f

to pomp, where f is an R function. The arguments of f should be chosen from among the state variables, parameters, covariates, and time. It must also have the argument .... f must return a named numeric vector of length equal to the number of observable variables. The names should match those of the observable variables.

Default behavior

The default emeasure is undefined. It will yield missing values (NA).

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

emeasure

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Description

pomp currently implements the following algorithms for estimating model parameters:

• iterated filtering (IF2)

• particle Markov chain Monte Carlo (PMCMC)

• approximate Bayesian computation (ABC)

• probe-matching via synthetic likelihood

• nonlinear forecasting

• power-spectrum matching

• Liu-West Bayesian sequential Monte Carlo

• Ensemble and ensemble-adjusted Kalman filters

Details

Help pages detailing each estimation algorithm are provided.

See Also

basic model components, workhorse functions, elementary algorithms.

More on pomp estimation algorithms: abc(), bsmc2(), mif2(), nlf, pmcmc(), pomp-package, probe_match, spect_match

Description

pomp provides a number of probability distributions that have proved useful in modeling partially observed Markov processes. These include the Euler-multinomial family of distributions and the the Gamma white-noise processes.

Usage

reulermultinom(n = 1, size, rate, dt)

deulermultinom(x, size, rate, dt, log = FALSE)

rgammawn(n = 1, sigma, dt)

Arguments

Details

If $N$ individuals face constant hazards of death in $K$ ways at rates $r_1, r_2, \dots, r_K$, then in an interval of duration $\Delta{t}$, the number of individuals remaining alive and dying in each way is multinomially distributed:

$(\Delta{n_0}, \Delta{n_1}, \dots, \Delta{n_K}) \sim \mathrm{Multinomial}(N;p_0,p_1,\dots,p_K),$

where $\Delta{n_0}=N-\sum_{k=1}^K \Delta{n_k}$ is the number of individuals remaining alive and $\Delta{n_k}$ is the number of individuals dying in way $k$ over the interval. Here, the probability of remaining alive is

$p_0=\exp(-\sum_k r_k \Delta{t})$

and the probability of dying in way $k$ is

$p_k=\frac{r_k}{\sum_j r_j} (1-p_0).$

In this case, we say that

$(\Delta{n_1},\dots,\Delta{n_K})\sim\mathrm{Eulermultinom}(N,r,\Delta t),$

where $r=(r_1,\dots,r_K)$. Draw $m$ random samples from this distribution by doing

    dn <- reulermultinom(n=m,size=N,rate=r,dt=dt),

where r is the vector of rates. Evaluate the probability that $x=(x_1,\dots,x_K)$ are the numbers of individuals who have died in each of the $K$ ways over the interval $\Delta t=$dt, by doing

    deulermultinom(x=x,size=N,rate=r,dt=dt).

Bretó & Ionides (2011) discuss how an infinitesimally overdispersed death process can be constructed by compounding a multinomial process with a Gamma white noise process. The Euler approximation of the resulting process can be obtained as follows. Let the increments of the equidispersed process be given by

    reulermultinom(size=N,rate=r,dt=dt).

In this expression, replace the rate $r$ with $r\,{\Delta{W}}/{\Delta t}$, where $\Delta{W} \sim \mathrm{Gamma}(\Delta{t}/\sigma^2,\sigma^2)$ is the increment of an integrated Gamma white noise process with intensity $\sigma$. That is, $\Delta{W}$ has mean $\Delta{t}$ and variance $\sigma^2 \Delta{t}$. The resulting process is overdispersed and converges (as $\Delta{t}$ goes to zero) to a well-defined process. The following lines of code accomplish this:

    dW <- rgammawn(sigma=sigma,dt=dt)

    dn <- reulermultinom(size=N,rate=r,dt=dW)

or

    dn <- reulermultinom(size=N,rate=r*dW/dt,dt=dt).

He et al. (2010) use such overdispersed death processes in modeling measles and the "Simulation-based Inference" course discusses the value of allowing for overdispersion more generally.

For all of the functions described here, access to the underlying C routines is available: see below.

Value

C API

An interface for C codes using these functions is provided by the package. Visit the package homepage to view the pomp C API document.

Author(s)

Aaron A. King

References

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

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

See Also

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Examples

print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1))
deulermultinom(x=dn,size=100,rate=c(1,2,3),dt=0.1)
## an Euler-multinomial with overdispersed transitions:
dt <- 0.1
dW <- rgammawn(sigma=0.1,dt=dt)
print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=dW))

Description

The mean of the filtering distribution

Usage

## S4 method for signature 'kalmand_pomp'
filter_mean(object, vars, ..., format = c("array", "data.frame"))

## S4 method for signature 'pfilterd_pomp'
filter_mean(object, vars, ..., format = c("array", "data.frame"))

Arguments

Details

The filtering distribution is that of

$X(t_k) \vert Y(t_1)=y^*_1,\dots,Y(t_k)=y^*_k,$

where $X(t_k)$, $Y(t_k)$ are the latent state and observable processes, respectively, and $y^*_t$ is the data, at time $t_k$.

The filtering mean is therefore the expectation of this distribution

$E[X(t_k) \vert Y(t_1)=y^*_1,\dots,Y(t_k)=y^*_k].$

See Also

More on sequential Monte Carlo methods: bsmc2(), cond_logLik(), eff_sample_size(), filter_traj(), kalman, mif2(), pfilter(), pmcmc(), pred_mean(), pred_var(), saved_states(), wpfilter()

Other extraction methods: coef(), cond_logLik(), covmat(), eff_sample_size(), filter_traj(), forecast(), logLik, obs(), pred_mean(), pred_var(), saved_states(), spy(), states(), summary(), time(), timezero(), traces()

Description