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

This document has its origins in the Getting Started vignette.

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')`

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 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*.

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\).

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,

- \(N_t\) is the true population density at time \(t\),
- \(y_t\) is the number of individuals sampled at time \(t\),
- 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),
- the parameter \(\phi\) is proportional to our sampling effort, and also has peculiar units.

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)`

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')
```