Preliminaries

## load pomp package
library("pomp")
## load Tristan da Cunha data
data(fluTdc1971)

This loads the pomp package, and the Tristan da Cunha data from the fitR package.

The deterministic skeleton

## define deterministic skeleton
seitlDeterSkel <- Csnippet("
    double trans[5];

    double beta = R0 / D_inf;
    double epsilon = 1 / D_lat;
    double nu = 1 / D_inf;
    double tau = 1 / D_imm;

    double N = S + E + I + T + L;

    trans[0] = beta * I / N * S;
    trans[1] = epsilon * E;
    trans[2] = nu * I;
    trans[3] = alpha * tau * T;
    trans[4] = (1 - alpha) * tau * T;

    DS = -trans[0] + trans[4];
    DE = trans[0] - trans[1];
    DI = trans[1] - trans[2];
    DT = trans[2] - trans[3] - trans[4];
    DL = trans[3];
    DInc = trans[1];
")

This defines the so-called deterministic skeleton. In pomp, the stochastic and deterministic versions of the same model can be specified in the same pomp object, and then used by different fitting methods. For example, one can get a first approximate fit using a fast method and the deterministic model before performing more computationally extensive methods on the stochastic models around those solutions. The deterministic skeleton defines how the model would behave if there were no noise in the system.

The model code is implemented as C code which is saved as a character string and converted into R code using Csnippet.

The first line

    double trans[5];

defines an array trans (of doubles, or floating point numbers), which will hold 5 elements. These are defined below according to the 5 transitions in the model.

The following section

    double beta = R0 / D_inf;
    double epsilon = 1 / D_lat;
    double nu = 1 / D_inf;
    double tau = 1 / D_imm;

    double N = S + E + I + T + L;

defines some auxiliary variables, to be used below when defining the transitions.

The following section

    trans[0] = beta * I / N * S;
    trans[1] = epsilon * E;
    trans[2] = nu * I;
    trans[3] = alpha * tau * T;
    trans[4] = (1 - alpha) * tau * T;

defines the strength of the 5 transitions in the model: infection, incubation, recovery, temporary immunity and long-term immunity.

Lastly, it is defined how the transitions change the compartments:

    DS = -trans[0] + trans[4];
    DE = trans[0] - trans[1];
    DI = trans[1] - trans[2];
    DT = trans[2] - trans[3] - trans[4];
    DL = trans[3];
    DInc = trans[1];

The stochastic model

## define stochastic model, for use with euler, see ?euler
seitlStochSim <- Csnippet("
    double rate[5];
    double dN[5];

    double beta = R0 / D_inf;
    double epsilon = 1 / D_lat;
    double nu = 1 / D_inf;
    double tau = 1 / D_imm;

    double N = S + E + I + T + L;

    rate[0] = beta * I / N;
    rate[1] = epsilon;
    rate[2] = nu;
    rate[3] = alpha * tau;
    rate[4] = (1 - alpha) * tau;

    reulermultinom(1, S, &rate[0], dt, &dN[0]);
    reulermultinom(1, E, &rate[1], dt, &dN[1]);
    reulermultinom(1, I, &rate[2], dt, &dN[2]);
    reulermultinom(2, T, &rate[3], dt, &dN[3]);

    S += -dN[0] + dN[4];
    E += dN[0] - dN[1];
    I += dN[1] - dN[2];
    T += dN[2] - dN[3] - dN[4];
    L += dN[3];
    Inc += dN[1];
")

This defines the stochastic model. This, again, is implemented as C code which is saved as a character string and later converted into R code using Csnippet. It defines transitions which are then used by the pomp function euler, which in turn uses the Euler-Maruyama approximation to the full stochastic system. This calculates the number of events happening in a given fixed time step by sampling from a multinomial distribution with the probabilities of each event given by the product of its rate and the time step. This happens in the following four lines of code:

    reulermultinom(1, S, &rate[0], dt, &dN[0]);
    reulermultinom(1, E, &rate[1], dt, &dN[1]);
    reulermultinom(1, I, &rate[2], dt, &dN[2]);
    reulermultinom(2, T, &rate[3], dt, &dN[3]);

Instead of calculating the rates of transition as in the deterministic skeleton above, reulermultinom draws random numbers from the compartment specified in the second argument to get the number of individuals who transition from that compartment to another one. These numbers are stored in the array N, and the numbers in each compartment adjusted in the following segment of code:

    S += -dN[0] + dN[4];
    E += dN[0] - dN[1];
    I += dN[1] - dN[2];
    T += dN[2] - dN[3] - dN[4];
    L += dN[3];
    Inc += dN[1];

The number of susceptibles decreases by the number of individuals randomly drawn from S at rate rate[0] (infection rate) and increases by the number of individuals randomly drawn from T at rate rate[4] (rate of loss of temporary immunity), and so on.

Random point observations

## define sampling random point observations
seitlGenObsPoint <- Csnippet("
    obs = rpois(rho * Inc);
")

(you may notice that the function in the example looks slightly different, to prevent numerical problems later)

This, again, defines a character string with C code, to be used with Csnippet. It defines a function to randomly generate point observations from model incidence, like our function rPointObs earlier. In this case, it uses a random draw from a Poisson distribution with mean rho*Inc, and stores it in the variable obs.

Point observation probability density

## define point observation probability density
seitlPointLike <- Csnippet("
    lik = dpois(obs, rho * Inc, give_log);
")

(you may notice that the function in the example looks slightly different, to prevent numerical problems later)

This, again, defines a character string with C code converted ot R using Csnippet. It defines a function to evaluate the probability density of point observations of the model incidence, like our function dPointObs earlier. In this case, it evaluates the probability following a Poisson probability distribution with mean rho * Inc and stores it in the variable lik. It can do both the logged and natural likelihood; which of the two is used will be given by give_log.

Prior density

seitlPrior <- Csnippet("
  lik = dunif(R0, 1, 50, 1) +
          dunif(D_lat, 0, 10, 1) +
          dunif(D_inf, 0, 15, 1) +
          dunif(D_imm, 0, 50, 1) +
          dunif(alpha, 0, 1, 1) +
          dunif(rho, 0, 1, 1);

  lik = give_log ? lik : exp(lik);
")

This, again, defines a character string with C code converted to R using Csnippet. It defines a function to evaluate the prior probability density of a set of parameters. The first segment of code,

lik <- dunif(R0, 1, 50, 1) +
  dunif(D_lat, 0, 10, 1) +
  dunif(D_inf, 0, 15, 1) +
  dunif(D_imm, 0, 50, 1) +
  dunif(alpha, 0, 1, 1) +
  dunif(rho, 0, 1, 1)

calculates the logged probability density according to uniform distributions with different bounds (e.g., 1 and 50 for R0), and takes the sum. The last arguments of 1 in dunif indicate that we want the logarithms of the probability densities. The last line,

lik <- give_log ? lik:exp(lik)

returns the calculated log-likelihood if give_log is 1, and the exponential of that number (that is, the natural likelihood), if it is 0. The value of give_log is set by whichever function calls the snippet.

Constructing the pomp object

seitlPomp <- pomp(
  data = fluTdc1971[, c("time", "obs")],
  skeleton = vectorfield(seitlDeterSkel),
  rprocess = euler(step.fun = seitlStochSim, delta.t = 0.1),
  rmeasure = seitlGenObsPoint,
  dmeasure = seitlPointLike,
  dprior = seitlPrior,
  partrans = parameter_trans(
    log = c("R0", "D_inf", "D_lat", "D_imm", "alpha", "rho")
  ),
  times = "time",
  t0 = 1,
  accumvars = "Inc",
  paramnames = c("R0", "D_inf", "D_lat", "D_imm", "alpha", "rho"),
  statenames = c("S", "E", "I", "T", "L", "Inc"),
  obsnames = c("obs")
)

This constructs the pomp object seitlPomp from the C functions defined above. The individual options are:

data <- fluTdc1971[, c("time", "obs")]

This specifies the data set. We select the two relevant columns, time and obs.

skeleton <- vectorfield(seitlDeterSkel)

These define the deterministic skeleton, as given by the C code above. If you are interested in an explanation of the vectorfield, see ?pomp.

rprocess <- euler(step.fun = seitlStochSim, delta.t = 0.1
)

This defines the stochastic (random) process, or the function that is used to sample trajectories. In this case, we use the Euler-Maruyama method (see above), with the step function given by our C code above, and a time step of delta.t = 0.1.

rmeasure <- Csnippet(seitlGenObsPoint)
dmeasure <- Csnippet(seitlPointLike)

These define the functions to randomly sample form the measurement process (equivalent to rPointObs earlier in the course) or evaluate the probability density at a measurement point (equivalent to dPointObs earlier in the course).

times <- "time"
t0 <- 1

These indicate that times in the data set (fluTdc1971) are given by the “time” column, and that the simulations are to start with the initial state at time 1.

paramnames <- c("R0", "D_inf", "D_lat", "D_imm", "alpha", "rho")
statenames <- c("S", "E", "I", "T", "L", "Inc")
obsnames <- c("obs")

These define the names of the parameters, states, and the observation variable.

Return to the pomp session.

 

This web site and the material contained in it were originally created in support of an annual short course on Model Fitting and Inference for Infectious Disease Dynamics at the London School of Hygiene & Tropical Medicine. All material is under a MIT license. Please report any issues or suggestions for improvement on the corresponding GitHub issue tracker. We are always keen to hear about any uses of the material here, so please do get in touch using the Discussion board if you have any questions or ideas, or if you find the material here useful or use it in your own teaching.