错误: in ‘mif2’: ‘dmeasure’ with log=TRUE returns illegal value

[xiaozuo986]

Thank you for your time at first!
I'm trying to estimate parameters in a epidemic model via the mif2 algorithm. Unfortunately there is often an error message that the dmeasure returns non-finite values during the mif2 process. I found that one of the state variables “cases” is negative by using the Rprintf statement, and the parameter “rho” related to this state is also negative, but I have set the parameter “rho” to be finite (0.0.1). I don't know why this happens.The code and data are:

simul_pf <- Csnippet("
                     // compute transmission rate
                     double betaIN = exp(b1*season1+b2*season2+b3*season3+b4*season4+b5*season5+b6*season6+b4*season4*bH*RH);

                     // gamma white noise
                     double dW = rgammawn(sigPRO,dt);

                     // force of infection
                     double foi = (betaOUT+((I1+q0*I2)/pop)*betaIN)*(dW/dt);

                     double dBS1 = (delta*pop+dpopdt)*dt;
                     double dS2S1 = muS2S1*S2*dt;
                     double dS1E = F*S1*dt;

                     double dEI1 = muEI1*E*dt;
                     double dI1S2 = muI1S2*I1*dt; // ask about the fIR
                     double dS2I2 = F*S2*dt;
                     double dI2S2 = muI2S2*I2*dt;
                     double dS1D = delta*S1*dt;
                     double dED = delta*E*dt;
                     double dI1D = delta*I1*dt;
                     double dS2D = delta*S2*dt;
                     double dI2D = delta*I2*dt;
                     double dKK = ((foi-K)/(tau/2.0))*dt;
                     double dFF = ((K-F)/(tau/2.0))*dt;

                     // compute equations
                     S1 += dBS1 + dS2S1 + dS1E - dS1D;
                     E += dS1E - dEI1 - dED;
                     I1 += dEI1 - dI1S2 - dI1D;
                     S2 += dI1S2 - dS2S1 - dS2I2 +dI2S2 - dS2D;
                     I2 += dS2I2 - dI2S2 - dI2D;
                     K += dKK;
                     F += dFF;
                     cases += rho*dEI1;
                     W += (dW-dt)/sigPRO;
                     ")

############ rmeas #################
rmeas_pf <- Csnippet("
                     double size = 1.0/sigOBS/sigOBS;
                     PF = rnbinom_mu(size,cases);
                     ")

############ dmeas #################
dmeas_pf <- Csnippet("
                     double size = 1.0/sigOBS/sigOBS;
                     lik = dnbinom_mu(PF,size,cases,1);
                     if (!give_log) lik = exp(lik);
                     if (!R_FINITE(lik)) Rprintf(\"%lg %lg %lg \\n\",PF,size,cases,lik);
                     ")

############ fromEst #################
fromEst_pf <- Csnippet("
                       TsigOBS = expit(sigOBS);
                       TsigPRO = expit(sigPRO);
                       TmuS2S1 = exp(muS2S1);
                       TmuEI1 = exp(muEI1);
                       TmuI1S2 = exp(muI1S2);
                       TmuI2S2 = exp(muI2S2);
                       TbetaOUT = exp(betaOUT);
                       Trho = expit(rho);
                       Ttau = expit(tau);
                       Tq0 = expit(q0);
                       TS1_0 = expit(S1_0);
                       TE_0 = expit(E_0);
                       TI1_0 = expit(I1_0);
                       TS2_0 = expit(S2_0);
                       TI2_0 = expit(I2_0);
                       TK_0 = expit(K_0);
                       TF_0 = expit(F_0);
                       double sum = TS1_0 + TE_0 + TI1_0 + TS2_0 + TI2_0;
                       TS1_0 /= sum;
                       TE_0 /= sum;
                       TI1_0 /= sum;
                       TS2_0 /= sum;
                       TI2_0 /= sum;
                       ")

############ fromEst #################

toEst_pf <- Csnippet("
                     TsigOBS = logit(sigOBS);
                     TsigPRO = logit(sigPRO);
                     TmuS2S1 = log(muS2S1);
                     TmuEI1 = log(muEI1);
                     TmuI1S2 = log(muI1S2);
                     TmuI2S2 = log(muI2S2);
                     TbetaOUT = log(betaOUT);
                     Trho = logit(rho);
                     Ttau = logit(tau);
                     Tq0 = logit(q0);
                     TS1_0 = logit(S1_0);
                     TE_0 = logit(E_0);
                     TI1_0 = logit(I1_0);
                     TS2_0 = logit(S2_0);
                     TI2_0 = logit(I2_0);
                     TK_0 = logit(K_0);
                     TF_0 = logit(F_0);

                     double sum = S1_0 + E_0 + I1_0 + S2_0 + I2_0;
                     TS1_0 /= log(S1_0/sum);
                     TE_0 /= log(E_0/sum);
                     TI1_0 /= log(I1_0/sum);
                     TS2_0 /= log(S2_0/sum);
                     TI2_0 /= log(I2_0/sum);
                     ")

############ initlz #################

initlz_pf <- Csnippet("
                      double m = pop/(S1_0 + E_0 + I1_0 + S2_0 + I2_0);

                      S1 = nearbyint(m*S1_0);
                      E = nearbyint(m*E_0);
                      I1 = nearbyint(m*I1_0);
                      S2 = nearbyint(m*S2_0);
                      I2 = nearbyint(m*I2_0);

                      K = K_0;
                      F = F_0;
                      cases = 0;
                      W = 0;
                      ")
dat <- read.csv("~/all_ahmedabad.csv")
now.num <- 3
z <- read.csv("~/paramGrid_surat.csv")
z <- subset(z, rho > 0.01)
param_pf <- as.numeric(z[now.num,])
names(param_pf) <- colnames(z)
par_names <- c("sigOBS","sigPRO","muS2S1","muEI1","muI1S2","muI2S2","betaOUT","delta","rho","tau","q0")
vp_names <- c("S1_0","E_0","I1_0","I2_0","S2_0","K_0","F_0")
sp_names <- c("b1","b2","b3","b4","b5","b6","bH")

pomp(
  data=subset(dat,select=c("time","PF")),
  params=param_pf,
  times="time",
  t0=with(dat,2*time[1]-time[2]),
  covar=covariate_table(time=dat$time,
    temp=dat$temp,
    RH=dat$RH,
    pop=dat$pop,
    dpopdt=dat$dpopdt,
    season1=dat$season1,
    season2=dat$season2,
    season3=dat$season3,
    season4=dat$season4,
    season5=dat$season5,
    season6=dat$season6,
    times ="time"
  ),
  rprocess = euler(step.fun = simul_pf, delta.t=1/365),
  rmeasure = rmeas_pf,
  dmeasure = dmeas_pf,
  rinit=initlz_pf,
  parameter_trans(toEst = toEst_pf,fromEst = fromEst_pf),
  accumvars = c("W","cases"),
  statenames = c("cases","S1","E","I1","S2","I2","K","F","W"),
  paramnames = c(par_names,vp_names,sp_names)
) -> po

############## data ########################
y <- read.csv("~/paramGrid_surat.csv")
param <- as.numeric(y[now.num,])
param.names <- colnames(y)
names(param) <- param.names

############## simulate #####################
rdd <- rw_sd(
  sigOBS=0.03,
  sigPRO=0.03,
  muS2S1=0.03,
  muEI1=0.03,
  muI2S2=0.03,
  muI2S2=0.03,
  delta=0.03,
  rho=ivp(0.03),
  tau=0.03,
  betaOUT=ivp(0.00001),
  S1_0=ivp(0.03),
  E_0=ivp(0.03),
  I1_0=ivp(0.03),
  S2_0=ivp(0.03),
  I2_0=ivp(0.03),
  K_0=ivp(0.03),
  F_0=ivp(0.03),
  b1=0.03,
  b2=0.03,
  b3=0.03,
  b4=0.03,
  b5=0.03,
  b6=0.03,
  bH=0.03
)

  seed <- ceiling(runif(1,min=1,max=2^30))
  set.seed(seed)
  param <- as.numeric(y[1,])
  names(param) <- colnames(y)

  tryCatch(
    mif2(
      po,
      Np=1000,
      Nmif=50,
      cooling.type="geometric",
      cooling.fraction.50=0.5,
      transform=TRUE,
      start=param,
      rw.sd=rdd,
      pred.mean=TRUE,
      filter.mean=TRUE,
      max.fail=500
    ),
    error = function(e) e
  ) -> mifout

[xiaozuo986]

This is the data used：
all_ahmedabad.csv
paramGrid_surat.csv

[kingaa]

Your transformation functions are not inverse to one another. For example, if you have

T_sigPRO = logit(sigPRO);

in the snippet that implements transformation to the estimation scale, then you would want

sigPRO = expit(T_sigPRO);

in the snippet that transforms back onto the model scale. See the manual.

P.S. I'm surprised the code that you've furnished runs at all. You must be using a very old version of the package. Why is that?

[xiaozuo986]

Thank you very much for your reply!

The version of "pomp" that I am using is 5.6. Unfortunately, even after making the changes to the code as you suggested, I am still not getting the correct results. The modified code is as follows:

############ fromEst #################
fromEst_pf <- Csnippet("
 sigOBS = expit(T_sigOBS );
 sigPRO = expit(T_sigPRO );
 muS2S1 = exp(T_muS2S1 );
 muEI1 = exp(T_muEI1 );
 muI1S2 = exp(T_muI1S2 );
 muI2S2 = exp(T_muI2S2 );
 betaOUT = exp(T_betaOUT );
 rho = expit(T_rho );
 tau = expit(T_tau );
 q0 = expit(T_q0 );
 S1_0 = expit(T_S1_0 );
 E_0 = expit(T_E_0 );
 I1_0 = expit(T_I1_0 );
 S2_0 = expit(T_S2_0 );
 I2_0 = expit(T_I2_0 );
 K_0 = expit(T_K_0 );
 F_0 = expit(T_F_0 );
 double sum = T_S1_0 + T_E_0 + T_I1_0 + T_S2_0 + T_I2_0;
 S1_0 = T_S1_0/sum;
 E_0 = T_E_0/sum;
 I1_0 = T_I1_0/sum;
 S2_0 = T_S2_0/sum;
 I2_0 = T_I2_0/sum;
 ")

############ fromEst #################

toEst_pf <- Csnippet("
 T_sigOBS = logit(sigOBS);
 T_sigPRO = logit(sigPRO);
 T_muS2S1 = log(muS2S1);
 T_muEI1 = log(muEI1);
 T_muI1S2 = log(muI1S2);
 T_muI2S2 = log(muI2S2);
 T_betaOUT = log(betaOUT);
 T_rho = logit(rho);
 T_tau = logit(tau);
 T_q0 = logit(q0);
 T_S1_0 = logit(S1_0);
 T_E_0 = logit(E_0);
 T_I1_0 = logit(I1_0);
 T_S2_0 = logit(S2_0);
 T_I2_0 = logit(I2_0);
 T_K_0 = logit(K_0);
 T_F_0 = logit(F_0);
double sum = S1_0 + E_0 + I1_0 + S2_0 + I2_0;
                     T_S1_0 = log(S1_0/sum);
                     T_E_0 = log(E_0/sum);
                     T_I1_0 = log(I1_0/sum);
                     T_S2_0 = log(S2_0/sum);
                     T_I2_0 = log(I2_0/sum);
                     ")

The error message is:

NOTE: The provided objects ‘transform’,‘start’,‘pred.mean’,‘filter.mean’,‘max.fail’ are available for use by POMP basic components.
70 159.591 -15087.7 
error: in ‘mif2’: ‘dmeasure’ with log=TRUE returns illegal value.
Log likelihood, data, states, and parameters are:
   time:         1997
 loglik:          NaN
     PF:           70
  cases:     -15087.7
     S1:       979665
      E:      15483.3
     I1:  1.03117e+06
     S2:  1.55717e+06
     I2:      5349.31
      K:     0.255111
      F:     0.307045
      W:   -0.0593478
 sigOBS:    0.0791581
 sigPRO:     0.220959
 muS2S1:      2.00935
  muEI1:       56.747
 muI1S2:      8.26055
 muI2S2:      86.2213
betaOUT:    0.0960373
    rho:  -0.00846253
    tau:     0.679695
     bH:    -0.838488
     q0:     0.845467
     b1:     -2.85551
     b2:      4.44979
     b3:     -3.16204
     b4:     0.944407
     b5:     0.521347
     b6:     -3.78877
  delta:    0.0333043
   S1_0:    0.0868384
    E_0:     0.213538
   I1_0: -0.000980317
   S2_0:     0.111821
   I2_0:     0.012409
    K_0:     0.292909
    F_0:     0.319137

[kingaa]

You were not passing the parameter transformations properly: like all arguments to pomp() (and other pomp functions), the arguments should be given by name. This explains why the nonsense code you had in your original post did not cause errors: it was never called. You chose to ignore the various warnings that pomp() gave you about mysterious objects being made available for use by basic components. In general, it is never safe to ignore a warning unless you understand the warning and know that it is safe to ignore it.

Also, the code you were running looks like it has been lifted with minimal changes from code designed to work under a very old version of the package. In particular, the transform, max.fail, start, pred.mean, filter.mean arguments are all meaningless now: some of the warnings that you felt happy to ignore were informing you of this fact. In a forthcoming release, these warnings will become errors. Had you been using the most recent version of pomp (as FAQ 1.1 instructs you to do before submitting a help request), you would have been informed of this. Do read the FAQ. Also, do install the latest version of pomp. In general, frequently updating your library of R packages is a standard practice that you would do well to adopt.

I've had a pretty thorough look at your code, which you haven't tried to explain, but it seems familiar enough to me that I can guess something about what you are trying to do. I've made quite a number of corrections, some of which are based on this kind of informed speculation. So read the following carefully. In particular, I don't think your rho parameter is an initial-value parameter and I don't think you want F_0 or K_0 to be negative.

You'll also note that I've changed the style of the code: good style makes code easier to read and also less prone to errors.

Running this under pomp v 5.11, I can now see that your model assumptions are problematic. In particular, for allowable values of the parameters, you can have negative values of the state variables, which are incompatible with the data. This happens because there is nothing stopping you removing more individuals from a compartment than are in the compartment.

Have a look at FAQ 3.6. It is not necessary to use the Euler-multinomial approximation, but you need to do so something to make sure that your compartment occupancy variables do not go negative.

library(tidyverse)
library(pomp)

read_csv("all_ahmedabad.csv") -> dat

read_csv("paramGrid_surat.csv") |>
  filter(rho > 0.01) |>
  slice(3) -> param_pf

dat |>
  select(-PF) |>
  covariate_table(times="time") |>
  repair_lookup_table(
    t=with(dat,c(2*time[1]-time[2],time))
  ) -> covars_pf

dat |>
  select(time,PF) |>
  pomp(
    times="time",
    t0=with(dat,2*time[1]-time[2]),
    params=param_pf,
    covar=covars_pf,
    rprocess = euler(
      step.fun=Csnippet(r"{
        // transmission rate
        double betaIN = exp(b1*season1+b2*season2+b3*season3+b4*season4+b5*season5+b6*season6+b4*season4*bH*RH);
        // force of infection
        double dW = rgammawn(sigPRO,dt);
        double foi = (betaOUT+((I1+q0*I2)/pop)*betaIN)*(dW/dt);
        // individual terms
        double dBS1 = (delta*pop+dpopdt)*dt;
        double dS2S1 = muS2S1*S2*dt;
        double dS1E = F*S1*dt;
        double dEI1 = muEI1*E*dt;
        double dI1S2 = muI1S2*I1*dt; // ask about the fIR
        double dS2I2 = F*S2*dt;
        double dI2S2 = muI2S2*I2*dt;
        double dS1D = delta*S1*dt;
        double dED = delta*E*dt;
        double dI1D = delta*I1*dt;
        double dS2D = delta*S2*dt;
        double dI2D = delta*I2*dt;
        double dKK = ((foi-K)/(tau/2.0))*dt;
        double dFF = ((K-F)/(tau/2.0))*dt;
        // balance equations
        S1 += dBS1 + dS2S1 + dS1E - dS1D;
        E += dS1E - dEI1 - dED;
        I1 += dEI1 - dI1S2 - dI1D;
        S2 += dI1S2 - dS2S1 - dS2I2 + dI2S2 - dS2D;
        I2 += dS2I2 - dI2S2 - dI2D;
        if (S1 < 0 || E < 0 || I1 < 0 || S2 < 0 || I2 < 0)
          Rprintf("%lg %lg %lg %lg %lg\n",S1,E,I1,S2,I2);
        K += dKK;
        F += dFF;
        cases += rho*dEI1;
        W += (dW-dt)/sigPRO;}"
      ),
      delta.t=1/365
    ),
    rmeasure = Csnippet(r"{
      double size = 1.0/sigOBS/sigOBS;
      PF = rnbinom_mu(size,cases);}"
    ),
    dmeasure = Csnippet(r"{
      double size = 1.0/sigOBS/sigOBS;
      lik = dnbinom_mu(PF,size,cases,give_log);}"
    ),
    rinit=Csnippet(r"{
      double m = pop/(S1_0 + E_0 + I1_0 + S2_0 + I2_0);
      S1 = nearbyint(m*S1_0);
      E = nearbyint(m*E_0);
      I1 = nearbyint(m*I1_0);
      S2 = nearbyint(m*S2_0);
      I2 = nearbyint(m*I2_0);
      K = K_0;
      F = F_0;
      cases = 0;
      W = 0;}"
    ),
    partrans=parameter_trans(
      logit=c("sigOBS","sigPRO","rho","tau","q0"),
      log=c("muS2S1","muEI1","muI1S2","muI2S2","betaOUT","K_0","F_0"),
      barycentric=c("S1_0","E_0","I1_0","S2_0","I2_0")
    ),
    accumvars = c("W","cases"),
    statenames = c("cases","S1","E","I1","S2","I2","K","F","W"),
    paramnames = c(
      c("sigOBS","sigPRO","muS2S1","muEI1","muI1S2","muI2S2","betaOUT","delta","rho","tau","q0"),
      c("S1_0","E_0","I1_0","I2_0","S2_0","K_0","F_0"),
      c("b1","b2","b3","b4","b5","b6","bH")
    )
  ) -> po

set.seed(974650257)

po |>
  pfilter(Np=1000) |>
  replicate(n=5) |>
  concat() -> pf

plot(pf)

pf |>
  logLik() |>
  logmeanexp(se=TRUE,ess=TRUE)

po |>
  mif2(
    Np=1000,
    Nmif=50,
    cooling.fraction.50=0.5,
    rw.sd=rw_sd(
      sigOBS=0.03,sigPRO=0.03,
      muS2S1=0.03,muEI1=0.03,muI2S2=0.03,muI2S2=0.03,
      delta=0.03,
      rho=0.03,tau=0.03,
      S1_0=ivp(0.03),E_0=ivp(0.03),I1_0=ivp(0.03),S2_0=ivp(0.03),I2_0=ivp(0.03),
      K_0=ivp(0.03),F_0=ivp(0.03),
      b1=0.03,b2=0.03,b3=0.03,b4=0.03,b5=0.03,b6=0.03,bH=0.03
    )
  ) -> mifout

[kingaa]

As a simple example of what I'm talking about, here's something I knocked up:

library(pomp)

simulate(
  t0=-1/52,
  times=seq(0,10,by=1/52),
  params=c(
    gamma = 26, mu = 0.02, iota = 0.1,
    Beta = 400, Beta_sd = 0.01,
    rho = 0.6, k = 100,
    pop = 1e6,
    S_0 = 25/400, I_0 = 0.001, R_0 = 375/400
  ),
  globals = Csnippet(r"{
    void eeulermultinom
    (int m, double size, const double *rate, double dt, double *trans)
    {
      double lambda = 0.0;
      int j, k;
      for (k = 0; k < m; k++) lambda += rate[k];
      if (lambda > 0.0) {
        size = size*(1-exp(-lambda*dt));
        for (k = 0; k < m; k++)
          trans[k] = size*rate[k]/lambda;
      } else {
        for (k = 0; k < m; k++) trans[k] = 0.0;
      }
    }}"
  ),
  rprocess = euler(
    step.fun=Csnippet(r"{
      double dW = rgammawn(Beta_sd,dt);
      double rate[6];
      double trans[6];
      rate[0] = mu*pop;
      rate[1] = Beta*(I+iota)/pop*dW/dt;
      rate[2] = mu;
      rate[3] = gamma;
      rate[4] = mu;
      rate[5] = mu;
      trans[0] = rate[0]*dt;
      eeulermultinom(2,S,&rate[1],dt,&trans[1]);
      eeulermultinom(2,I,&rate[3],dt,&trans[3]);
      eeulermultinom(1,R,&rate[5],dt,&trans[5]);
      // balance equations
      S += trans[0] - trans[1] - trans[2];
      I += trans[1] - trans[3] - trans[4];
      R += trans[3] - trans[5];
      C += trans[3];
      W += (dW-dt)/Beta_sd;}"
    ),
    delta.t=1/365
  ),
  rmeasure = Csnippet(r"{
      reports = rnbinom_mu(k,rho*C);}"
  ),
  rinit=Csnippet(r"{
      double m = pop/(S_0+I_0+R_0);
      S = m*S_0;
      I = m*I_0;
      R = m*R_0;
      C = 0;
      W = 0;}"
  ),
  partrans=parameter_trans(
    logit=c("rho"),
    log=c("gamma","mu","k","Beta","Beta_sd","iota"),
    barycentric=c("S_0","I_0","R_0")
  ),
  obsnames = "reports",
  accumvars = c("W","C"),
  statenames = c("S","I","R","C","W"),
  paramnames = c(
    "pop","rho","k","gamma","mu","Beta","Beta_sd","iota",
    "S_0","I_0","R_0"
  )
) |>
  plot()

[kingaa]

This eeulermultinom function is available as part of the regular pomp C API as of the new release 5.11.1.0. As a result, in that version, you can drop the globals argument in the above.

[xiaozuo986]

The code is now running smoothly without any errors. I modified the "slice" function in your code to precisely extract the 10th row, which was a crucial step for my specific requirements.There is still an extensive amount for me to master and understand about the "pomp" package. The entire process of exchanging ideas with you has been an absolute delight and a highly rewarding experience.

Thank you so very much for your extremely detailed and comprehensive explanation. It has been of immense help to me!

[kingaa]

Thank you for providing your code and making it possible for me to help. Importantly, your contribution stimulated improvements to the package.

Keep in mind that just because your code runs without errors does not mean that it is correct. If simply changing the parameter set eliminated the errors, there is reason to doubt the correctness of the inferences. Beyond the syntax errors and use of defunct options, the serious issue we uncovered was the fact that, in certain parameter regimes the compartment occupancies become negative.

Finally, take my advice: keep your R packages updated and read the FAQ before posting.

[xiaozuo986]

I apologize for not explaining clearly at the beginning. The code and data utilized in this project are actually derived from a specific paper (https://doi.org/10.1038/s41467-022-28145-7). If you want to know more details of the model you can check this article.
I have taken your advice and upgraded the R package to the latest version. Thanks again for your patience!