Fit deterministic model using pomp #188

Fuhan-Yang · 2023-03-15T20:22:39Z

Fuhan-Yang
Mar 15, 2023

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!

Answered by kingaa

Mar 20, 2023

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…

View full answer

kingaa · 2023-03-15T21:48:03Z

kingaa
Mar 15, 2023
Maintainer

Why not use trajectory matching?

1 reply

Fuhan-Yang Mar 16, 2023
Author

My bad. Thank you!

Fuhan-Yang · 2023-03-20T14:37:10Z

Fuhan-Yang
Mar 20, 2023
Author

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!

0 replies

kingaa · 2023-03-20T16:43:56Z

kingaa
Mar 20, 2023
Maintainer

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

0 replies

Fuhan-Yang · 2023-03-20T20:00:51Z

Fuhan-Yang
Mar 20, 2023
Author

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!

2 replies

kingaa Mar 21, 2023
Maintainer

Try upgrading to recent versions. FAQ 1.1.

Fuhan-Yang Mar 21, 2023
Author

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fit deterministic model using pomp #188

{{title}}

Replies: 4 comments 3 replies

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Fit deterministic model using pomp #188

Fuhan-Yang Mar 15, 2023

Replies: 4 comments · 3 replies

kingaa Mar 15, 2023 Maintainer

Fuhan-Yang Mar 16, 2023 Author

Fuhan-Yang Mar 20, 2023 Author

kingaa Mar 20, 2023 Maintainer

Fuhan-Yang Mar 20, 2023 Author

kingaa Mar 21, 2023 Maintainer

Fuhan-Yang Mar 21, 2023 Author

Fuhan-Yang
Mar 15, 2023

Replies: 4 comments 3 replies

kingaa
Mar 15, 2023
Maintainer

Fuhan-Yang Mar 16, 2023
Author

Fuhan-Yang
Mar 20, 2023
Author

kingaa
Mar 20, 2023
Maintainer

Fuhan-Yang
Mar 20, 2023
Author

kingaa Mar 21, 2023
Maintainer

Fuhan-Yang Mar 21, 2023
Author