错误: 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