The R codes for this document are provided as a script.

This document has its origins in the Getting Started vignette.

Introduction

The R package pomp provides facilities for modeling POMPs, a toolbox of statistical inference methods for analyzing data using POMPs, and a development platform for implmenting new POMP inference methods. The basic data-structure provided by pomp is the object of class pomp, alternatively known as a pomp object. It is a container that holds real or simulated data and a POMP model, possibly together with other information such as model parameters, that may be needed to do things with the model and data. A real pomp data analysis begins with constructing one or more pomp objects to hold the data and the model or models under consideration. Here, we’ll illustrate this process a dataset of Parus major abundance in Wytham Wood, near Oxford (McCleery and Perrins 1991).

Download and plot the data:

loc <- url("https://kingaa.github.io/sbied/intro/parus.csv")
dat <- read.csv(loc)
head(dat)
 year pop
 1960 148
 1961 258
 1962 185
 1963 170
 1964 267
 1965 239
plot(pop~year,data=dat,type='o')

A stochastic Ricker model

Let’s suppose we wish to investigate the extent to which the Ricker model (Ricker 1954) explains these data. First, let’s recall the details of the Ricker (1954) model.

The deterministic Ricker map

The Ricker map describes the deterministic dynamics of a simple population: \[N_{t+1} = r\,N_{t}\,\exp(-c\,N_{t}).\] Here, \(N_t\) is the population density at time \(t\) and \(r\) is a fixed value, related to the population’s intrinsic capacity to increase. \(N\) is a state variable, \(r\) and \(c\) are parameters. If we know \(r\) and the initial condition \(N_0\), the deterministic Ricker equation predicts the future population density at all times \(t=1,2,\dots\). We can view the initial condition, \(N_0\) as a special kind of parameter, an initial-value parameter.

Process noise

We can model process noise in this system by making the growth rate \(r\) into a random variable. For example, if we assume that the intrinsic growth rate is log-normally distributed, \(N\) becomes a stochastic process governed by \[N_{t+1} = r\,N_{t}\,\exp(-c\,N_{t}+\varepsilon_{t}), \qquad \varepsilon_{t}\;\sim\;\mathrm{Normal}\left(0,\sigma\right),\] where the new parameter \(\sigma\) is the standard deviation of the noise process \(\varepsilon\).

Measurement error

Let’s suppose that the Ricker model is our model for the dynamics of a real population. However, we cannot know the exact population density at any time, but only estimate it through sampling.

Let’s model measurement error by assuming the measurements, \(y_t\), are Poisson with mean \(\phi\,N_t\): \[y_{t}\;\sim\;\mathrm{Poisson}\left(\phi\,N_{t}\right)\]

In this equation,

  1. \(N_t\) is the true population density at time \(t\),
  2. \(y_t\) is the number of individuals sampled at time \(t\),
  3. the choice of units for \(N\) is peculiar and depends on the parameters (e.g., \(N=\log(r)/c\) is the equilibrium of the deterministic model),
  4. the parameter \(\phi\) is proportional to our sampling effort, and also has peculiar units.

Implementing the Ricker model in pomp

The call to construct a pomp object is, naturally enough, pomp. Documentation on this function can be had by doing ?pomp. Learn about the various things you can do once you have a pomp object by doing methods?pomp and following the links therein. Read an overview of the package as a whole with links to its main features by doing package?pomp. A complete index of the functions in pomp is returned by the command library(help=pomp). Finally, the home page for the pomp project is https://kingaa.github.io/pomp/; there you have access to the complete source code, tutorials, manuals, a news blog, a facility for reporting issues and making help requests, etc.

Now, to construct our pomp object:

library(pomp)
parus <- pomp(dat,times="year",t0=1959)

The times argument specifies that the column of dat labelled “year” gives the measurement times; t0 is the “zero-time”, the time at which the state process will be initialized. We’ve set it to one year prior to the beginning of the data. Plot the new pomp object:

plot(parus)

Adding in the process model simulator

We can add the stochastic Ricker model to parus by writing a Csnippet that simulates one realization of the stochastic process, from an arbitary time \(t\) to \(t+1\), given arbitrary states and parameters. We provide this to pomp in the form of a Csnippet, a little snippet of C code that performs the computation. The following does this.

stochStep <- Csnippet("
  e = rnorm(0,sigma);
  N = r*N*exp(-c*N+e);
")
pomp(parus,rprocess=discrete_time(step.fun=stochStep,delta.t=1),
     paramnames=c("r","c","sigma"),statenames=c("N","e")) -> parus

Note that in the above, we use the exp and rnorm functions from the R API. In general any C function provided by R is available to you. pomp also provides a number of C functions that are documented in the header file, pomp.h, that is installed with the package. See the Csnippet documentation (?Csnippet) to read more about how to write them. Note too that we use discrete.time.sim here because the model is a stochastic map. We specify that the time step of the discrete-time process is delta.t, here, 1 yr.

At this point, we have what we need to simulate the state process:

sim <- simulate(parus,params=c(N_0=1,e_0=0,r=12,c=1,sigma=0.5),
                format="data.frame")
plot(N~year,data=sim,type='o')