Fit deterministic model using pomp

[Fuhan-Yang]

Thanks for developing such a great package.

I'm new to pomp. I'm trying to fit a deterministic version (the skeleton after removing the process stochasticity) of my stochastic model to check the differences between these two versions. I wonder if iterated filtering (IF) can fit the deterministic model as well (it should be, while all the particles are the same..right?), but it can be computationally intense to fit the deterministic model using IF.

My question is, if the goal is to compare the estimates between deterministic and stochastic version of a model, do you have any advice? Thanks!

[kingaa]

Why not use trajectory matching?

[Fuhan-Yang]

My bad. Thank you!

[Fuhan-Yang]

Hi Dr.King, the deterministic trajectory always returns NA when I tried to use trajectory(), can you please take look at my code?
Using measles as an example:

library(pomp)
library(tidyverse)

read_csv(
  "https://kingaa.github.io/sbied/stochsim/Measles_Consett_1948.csv"
) |>
  select(week,reports=cases) -> meas

sir_rinit <- Csnippet("
  S = nearbyint(eta*N);
  I = 1;
  R = nearbyint((1-eta)*N);
  H = 0;
  ")

sir_dmeas <- Csnippet("
  lik = dnbinom_mu(reports,k,rho*H,give_log);
  ")

meas_ske <- Csnippet("
   double dt = 1.0;
   DS = S - S*(1-exp(-Beta*I/N*dt));
   DI = I + S*(1-exp(-Beta*I/N*dt)) - I*(1-exp(-mu_IR*dt));
   DR = R + I*(1-exp(-mu_IR*dt));
")

meas |>
  trajectory(
    times = "week", t0 = 0,
    params=c(
      Beta=7.5,mu_IR=0.5,rho=0.5,k=10,eta=0.03,N=38000
    ),
    rinit=sir_rinit,
    dmeasure=sir_dmeas,
    skeleton = map(meas_ske,delta.t = 1),
    accumvars="H",
    statenames=c("S","I","R","H"),
    paramnames=c("Beta","mu_IR","N","eta","rho","k")
  ) -> measSIR

Thanks for the help!

[kingaa]

I saw two errors in your code. First, you loaded tidyverse after pomp, thus masking the map function. To avoid this problem, you can either load tidyverse before pomp, or explicitly specify the call to pomp::map().

The second error was that, in your skeleton component, you were not specifying how the H variable maps under the dynamics.

Compare the following.

library(pomp)
library(tidyverse)

read_csv(
  "https://kingaa.github.io/sbied/stochsim/Measles_Consett_1948.csv"
) |>
  select(week,reports=cases) -> meas

sir_rinit <- Csnippet("
  S = nearbyint(eta*N);
  I = 1;
  R = nearbyint((1-eta)*N);
  H = 0;
  ")

sir_dmeas <- Csnippet("
  lik = dnbinom_mu(reports,k,rho*H,give_log);
  ")

meas_ske <- Csnippet("
   double dt = 1.0;
   DS = S - S*(1-exp(-Beta*I/N*dt));
   DI = I + S*(1-exp(-Beta*I/N*dt)) - I*(1-exp(-mu_IR*dt));
   DR = R + I*(1-exp(-mu_IR*dt));
   DH = H + S*(1-exp(-Beta*I/N*dt));
")

meas |>
  trajectory(
    times = "week", t0 = 0,
    params=c(
      Beta=7.5,mu_IR=0.5,rho=0.5,k=10,eta=0.08,N=38000
    ),
    rinit=sir_rinit,
    dmeasure=sir_dmeas,
    skeleton = pomp::map(meas_ske,delta.t = 1),
    accumvars="H",
    statenames=c("S","I","R","H"),
    paramnames=c("Beta","mu_IR","N","eta","rho","k")
  ) -> measSIR

measSIR |> plot()

[Fuhan-Yang]

Thank you so much for the instant answer.

I tried the code above but an error returned as "Error: ‘trajectory’ is undefined for ‘object’ of class ‘tbl_df’‘tbl’‘data.frame’." This is resolved if a pomp model is used as input. While the trajectory still returned NA. Below is the code.

Not sure if this is relevant to the R or pomp version on my laptop. My pomp is Version 3.4 and R is 4.0.3. Thanks in advance!

[kingaa]

Try upgrading to recent versions. FAQ 1.1.

[Fuhan-Yang]

Thank you so much. The issue is resolved under R version 4.2.3 and pomp version 4.7.