NA in traj_objfun()

[Fuhan-Yang]

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!!

[kingaa]

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!

[Fuhan-Yang]

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.