NA in traj_objfun() #190

Fuhan-Yang · 2023-03-29T05:15:13Z

Fuhan-Yang
Mar 29, 2023

I'm using traj_objfun() to fit the deterministic version of my stochastic model. While I found the results are not consistent, especially, traj_objfun() mostly returns NA, given I've checked my deterministic model and ensured it is valid. Can you please give me advice on this? Below is my model:

I'm building a SIRS model. Below is my code in pomp:

library(tidyverse)
library(pomp)
# simulated data #
dt = data.frame(daynum = c(3653:(3653 + 3652)),
                ili = round(rbinom(3653,size = 120,
                                   prob = 0.5 + rnorm(3652,sd = 0.01)) + 2 * sin(2*pi/365*c(1:3653))))

### initial condition ###
sirs_init = Csnippet("
  S = N-300;
  I = 300;
  R = 0;
")

### process model ###
sirs_binom_step = Csnippet("

  double dN_SI = rbinom(S,beta*I/N*dt);;
  double dN_IR = rbinom(I,mu_IR*dt);
  double dN_R = rbinom(R,mu_R*dt);
    
  // flow
  S = S - dN_SI + dN_R;
  I = I - dN_IR + dN_SI;
  R = R - dN_R + dN_IR;
 
  // delta-J
  J = dN_SI;
")

### measure model ###
inc_dmeas = Csnippet("
  lik = dbinom(ili, J, rho, give_log);
")

inc_rmeas <- Csnippet("
  ili = rbinom(J, rho);
")

##### plug in pomp ####
sirs_sto_annual_beta = 
  pomp(data = dt,
       times = 'daynum',
       t0 = 1,
       rinit = sirs_init,
       rprocess = euler(sirs_binom_step,delta.t = 1),
       rmeasure = inc_rmeas,
       dmeasure = inc_dmeas,
       obsnames = 'ili',
       statenames = c('S','I','R','J'),
       paramnames = c('beta','mu_IR','mu_R','N','rho'),
       accumvars = 'J')

### skeleton ###
sirs_annual_beta_skeleton = Csnippet("
  
  double dt = 1.0;
 
  // flow
  DS = S - S*beta*I/N*dt + R*mu_R*dt;
  DI = I - I*mu_IR*dt + beta*S*I/N*dt;
  DR = R - R*mu_R*dt + I*mu_IR*dt;
   
  // delta-J 
  DJ = S*beta*I/N*dt;
   
")

## build traj_func ##
traj_func = 
  traj_objfun(sirs_sto_annual_beta,
              est = c('beta','rho','mu_R'),
              params = c(beta = 0.2,
                         rho = 0.01,mu_R = 4/150,
                         mu_IR = 1/5,N = 1e7),
              rinit = sirs_init,
              skeleton = pomp::map(sirs_annual_beta_skeleton,delta.t = 1),
              paramnames = c('beta','mu_IR','mu_R','N','rho'),
              statenames = c('S','I','R','J'))

traj_func(c(beta = 0.2,mu_R = 4/150,rho = 0.01))

I also used traj_func() to do a global search on likelihood, it returns NA most of the time and few valid numbers on discrete parameter combinations.

Can you please tell me what might go wrong? Thank you so much!!

Answered by kingaa

Mar 29, 2023

I notice that you round J in ili_nll_objfun, but not in traj_func. Notice that dbinom insists (correctly) on size being an integer:

dbinom(50,size=c(100,100.1),prob=0.5)
[1] 0.07958924        NaN
Warning message:
In dbinom(50, size = c(100, 100.1), prob = 0.5) : NaNs produced

The usual rounding function in C is nearbyint. Cf. the following:

### measure model ###
inc_dmeas <- Csnippet("
  lik = dbinom(ili, nearbyint(J), rho, give_log);
")

## build traj_func using simulated data ##
data.frame(
  daynum = seq(from=3653,to=3653 + 3652)
) |>
  mutate(
    ili = round(
      rbinom(
        n = length(daynum),
        size = 120,
        prob = 0.5 + rnorm(n=length(daynum),sd = 0.01)
      ) …

View full answer

kingaa · 2023-03-29T19:10:04Z

kingaa
Mar 29, 2023
Maintainer

I notice that you round J in ili_nll_objfun, but not in traj_func. Notice that dbinom insists (correctly) on size being an integer:

dbinom(50,size=c(100,100.1),prob=0.5)
[1] 0.07958924        NaN
Warning message:
In dbinom(50, size = c(100, 100.1), prob = 0.5) : NaNs produced

The usual rounding function in C is nearbyint. Cf. the following:

### measure model ###
inc_dmeas <- Csnippet("
  lik = dbinom(ili, nearbyint(J), rho, give_log);
")

## build traj_func using simulated data ##
data.frame(
  daynum = seq(from=3653,to=3653 + 3652)
) |>
  mutate(
    ili = round(
      rbinom(
        n = length(daynum),
        size = 120,
        prob = 0.5 + rnorm(n=length(daynum),sd = 0.01)
      ) +
        2 * sin(2*pi/365*daynum)
    )
  ) |>
  traj_objfun(
    times = "daynum",
    t0 = 1,
    rinit = sirs_init,
    skeleton = pomp::map(sirs_annual_beta_skeleton, delta.t = 1),
    dmeasure = inc_dmeas,
    accumvars = "J",
    params = c(
      beta = 0.2,seas = 0.02,
      rho = 0.01,mu_R = 4/150,
      mu_IR = 1/5,N = 1e7
    ),
    est = c("beta","seas","rho","mu_R"),
    paramnames = c("beta","seas","mu_IR","mu_R","N","rho"),
    statenames = c("S","I","R1","R2","R3","R4","J")
  ) -> traj_func

traj_func(c(beta = 0.2,seas = 0.02,mu_R = 4/150,rho = 0.01))
[1] 14292.71

I also notice that you have extraneous code in your example. Please reduce your examples to their essentials. Not only is this a courtesy to those you are asking for help, and to the readers of this discussion, but it may help you find the error yourself!

1 reply

Fuhan-Yang Mar 31, 2023
Author

Thank you so much for your time checking my code. Sorry about the lengthy question. The discussion is updated as only the core issue is included.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

NA in traj_objfun() #190

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 1 reply

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

NA in traj_objfun() #190

Fuhan-Yang Mar 29, 2023

Replies: 1 comment · 1 reply

kingaa Mar 29, 2023 Maintainer

Fuhan-Yang Mar 31, 2023 Author

Fuhan-Yang
Mar 29, 2023

Replies: 1 comment 1 reply

kingaa
Mar 29, 2023
Maintainer

Fuhan-Yang Mar 31, 2023
Author