traj_objfun question

[Steven-Lau-lib]

I am using traj_objfun for the preliminary estimation of parameter values, but I encountered an issue during the process. I used dmeasure to calculate the likelihood, but the parameters corresponding to the minimum value did not fit the data well.
However, when I estimate the parameters using the least squares method, the estimated parameters fit the data well, but I noticed that their likelihood is not the best.

What could be the reason that my model does not estimate the parameters well using the likelihood-based method?
What I'm building is a model with age structure. This is the code for each part of the pomp.

rinit_fn <- function (
  Ntotal,
  prop_Total1,
  prop_Total2,
  prop_Total3,
  prop_Total4,
  prop_Total5,
  prop_Total6,
  prop_Total7,
  prop_Total8,
  prop_Total9,
  prop_Total10,
  prop_Total11,
  prop_Total12,
  prop_Total13,
  prop_Total14,
  prop_Total15,
  prop_Total16,
  prop_Total17,
  prop_Total18,
  prop_Total19,
  prop_Total20,
  prop_Total21,
  prop_Total22,
  ...
) {
  Num_Ag <- 22
  prop_Total <- c(
    prop_Total1,prop_Total2,prop_Total3,prop_Total4,prop_Total5,
    prop_Total6,prop_Total7,prop_Total8,prop_Total9,prop_Total10,
    prop_Total11,prop_Total12,prop_Total13,prop_Total14,prop_Total15,
    prop_Total16,prop_Total17,prop_Total18,prop_Total19,prop_Total20,
    prop_Total21,prop_Total22
  )

  state <- c(
    S1 = round(Ntotal * prop_Total * 0.2),
    E1 = round(Ntotal * prop_Total * 0.1),
    I1 = round(Ntotal * prop_Total * 0.1),
    S2 = round(Ntotal * prop_Total * 0.2),
    E2 = round(Ntotal * prop_Total * 0.1),
    I2 = round(Ntotal * prop_Total * 0.1),
    R = round(Ntotal * prop_Total * 0.2),
    Vwp = numeric(Num_Ag),
    Vap = numeric(Num_Ag),
    SV1 = numeric(Num_Ag),
    SV2 = numeric(Num_Ag),
    C1 = numeric(Num_Ag),
    C2 = numeric(Num_Ag)
  )

  names(state) <- c(
    paste0(
      rep(
        c(
          "S1","E1","I1",
          "S2","E2","I2",
          "R",
          "Vwp","Vap",
          "SV1","SV2",
          "C1","C2"
        ),
        each = Num_Ag
      ),
      1:Num_Ag
    )
  )
  state
}

dobs <- Csnippet("
  int i;
  double mu;
  double rho[22];
  const double *C1_vec = &C11;
  const double *C2_vec = &C21;
  double f = 0.0;

  for (i = 0; i < 22; i++) {
    if (i < 3) {
      rho[i] = rho1;
    } else if (i < 8) {
      rho[i] = rho1 * Ratio_2_1;
    } else if (i < 10) {
      rho[i] = rho1 * Ratio_3_1;
    } else if (i < 22) {
      rho[i] = rho1 * Ratio_4_1;
    }
   }

  // If all report variables are available
  if (ISNA(report0_) && ISNA(report1_2) && ISNA(report3_4) && ISNA(report5_9) && ISNA(report10_14) && ISNA(report15_24) && ISNA(report25_)){
    f += (give_log) ? 0 : 1;
   } else {
    // Likelihood for report0_
    mu = 0.0;
    for (i = 0; i < 3; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report0_), 1.0 / k, mu, 1);

    // Likelihood for report1_2
    mu = 0.0;
    for (i = 3; i < 5; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report1_2), 1.0 / k, mu, 1);

    // Likelihood for report3_4
    mu = 0.0;
    for (i = 5; i < 8; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report3_4), 1.0 / k, mu, 1);

    // Likelihood for report5_9
    mu = 0.0;
    for (i = 8; i < 9; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report5_9), 1.0 / k, mu, 1);

    // Likelihood for report10_14
    mu = 0.0;
    for (i = 9; i < 10; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report10_14), 1.0 / k, mu, 1);

    // Likelihood for report15_24
    mu = 0.0;
    for (i = 10; i < 12; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report15_24), 1.0 / k, mu, 1);

    // Likelihood for report25_
    mu = 0.0;
    for (i = 12; i < 22; i++) {
      mu += (C1_vec[i] + C2_vec[i] * Ratio_sec) * rho[i];
    }
    f += dnbinom_mu(nearbyint(report25_), 1.0 / k, mu, 1);

   }
  lik = (give_log) ? f : exp(f);
")


skeleton_fn <- Csnippet("
  int i, j;
  int nages = 22, nages_cmat = 16; // No of age classes in the model, and in the contact matrix (without infant classes)

  double plambda[nages]; // Force of infection in each age group

  const double *S1_vec = &S11;
  const double *E1_vec = &E11;
  const double *I1_vec = &I11;
  const double *S2_vec = &S21;
  const double *E2_vec = &E21;
  const double *I2_vec = &I21;
  const double *R_vec = &R1;
  const double *Vwp_vec = &Vwp1;
  const double *Vap_vec = &Vap1;
  const double *SV1_vec = &SV11;
  const double *SV2_vec = &SV21;

  const double *CR= &C_1_1;
  const double *pmu= &mu1;
  const double *pdelta= &pdelta1;

  double pop = 0.0;
  // Total population in each group
  double N_vec[nages]; // Total population in each age group

  for (i = 0; i < nages; i++) {
    N_vec[i] = S1_vec[i] + E1_vec[i] + I1_vec[i] + S2_vec[i] +
      E2_vec[i] + I2_vec[i] + R_vec[i] + Vwp_vec[i] + Vap_vec[i] +
      SV1_vec[i] + SV2_vec[i];
    pop += N_vec[i];
   }

  // Create augmented contact matrix with seasonality factor
  // Necessary because more age classes in the model than in the contact data
  // Matrix of contact rates with seasonal forcing
  double Cseas[nages_cmat][nages_cmat];
  for (i = 0; i < nages_cmat; i++)
    for (j = 0; j < nages_cmat; j++)
      Cseas[i][j] = CR[nages_cmat*i+j];

  // Augmented contact matrix with 3 extra rows for infant classes
  double Cseas_aug[nages][nages_cmat];
  // Add two rows in the augmented contact matrix
  for (i = 0; i < nages; i++) {
    for (j = 0; j < nages_cmat; j++) {
      if (i <= 3) {
        Cseas_aug[i][j] = Cseas[0][j];
      } else  if (i < 8) {
        Cseas_aug[i][j] = Cseas[1][j];
      } else {
        Cseas_aug[i][j] = Cseas[i - 6][j];
      }
    }
   }

  // No infected and population size without infant age classes
  double I1_tilde[nages_cmat], I2_tilde[nages_cmat], N_tilde[nages_cmat];

  // Sum no infected in 0-2 y
  I1_tilde[0] = I1_vec[0] + I1_vec[1] + I1_vec[2] + I1_vec[3];
  // Sum no infected in 0-2 y
  I2_tilde[0] = I2_vec[0] + I2_vec[1] + I2_vec[2] + I2_vec[3];
  // Sum population in 0-2 y
  N_tilde[0] = N_vec[0] + N_vec[1] + N_vec[2] + N_vec[3];

  // Sum no infected in 2-5 y
  I1_tilde[1] = I1_vec[4] + I1_vec[5] + I1_vec[6] + I1_vec[7];
  // Sum no infected in 2-5 y
  I2_tilde[1] = I2_vec[4] + I2_vec[5] + I2_vec[6] + I2_vec[7];
  // Sum population in 2-5 y
  N_tilde[1] = N_vec[4] + N_vec[5] + N_vec[6] + N_vec[7];

  for (i = 2; i < nages_cmat; i++) {
    I1_tilde[i] = I1_vec[i + 6];
    I2_tilde[i] = I2_vec[i + 6];
    N_tilde[i] = N_vec[i + 6];
   }

  double *q = &FOI_1;
  double FOI[nages];
  for (i = 0; i < nages; i++) {
    if (i < 9) {
      FOI[i] = q[0];
    } else if (i < 11) {
      FOI[i] = q[0]*q[1];
    } else {
      FOI[i] = q[0]*q[1]*q[2];
    }
   }


  // Force of infection in each age group
  for (i = 0; i < nages; i++) {
    plambda[i] = 0.0; // Initialize values
    for (j = 0; j < nages_cmat; j++) {
      plambda[i] += FOI[i] * Cseas_aug[i][j] * (I1_tilde[j] + FOI_sec * I2_tilde[j]) / N_tilde[j];
    }
   }

  double *dS1 = &DS11;
  double *dE1 = &DE11;
  double *dI1 = &DI11;
  double *dS2 = &DS21;
  double *dE2 = &DE21;
  double *dI2 = &DI21;
  double *dR = &DR1;
  double *dVwp = &DVwp1;
  double *dVap = &DVap1;
  double *dSV1 = &DSV11;
  double *dSV2 = &DSV21;

  double *dC1 = &DC11;
  double *dC2 = &DC21;

  double pcoveragep;
  double pcoverageb;
  for (i = 0; i < nages; i++) {
    if (i == 1) {
      pcoveragep = primary_cov;
    } else {
      pcoveragep = 0;
    }
    if (i == 6  || i == 7) {
      pcoverageb = booster_cov;
    } else {
      pcoverageb = 0;
    }

    if (i == 0) {
      dS1[i] = pop*pbirth - plambda[i] * S1_vec[i] -
        pcoveragep * S1_vec[i] - pcoverageb * S1_vec[i] -
        (pdelta[i] + pmu[i])*S1_vec[i];  // Primary susceptible
      dE1[i] = plambda[i] * (S1_vec[i] + SV1_vec[i]) - paramsigma * E1_vec[i] -
        (pdelta[i] + pmu[i])*E1_vec[i];    // Primary exposed
      dI1[i] = paramsigma * E1_vec[i] - paramgamma * I1_vec[i] -
        (pdelta[i] + pmu[i])*I1_vec[i];    // Primary infected

      // Primary susceptible
      dS2[i] = - plambda[i] * S2_vec[i] + 1/paramomegaI * R_vec[i]
        + 1/(paramomegaI*Ratio_Vap_I) * Vap_vec[i]
        + 1/(paramomegaI*Ratio_Vwp_I) * Vwp_vec[i]
        - (pcoveragep+pcoverageb) * S2_vec[i]
        - (pdelta[i] + pmu[i])*S2_vec[i];
      // Primary exposed
      dE2[i] = plambda[i] * (S2_vec[i] + paramepsilonwp * Vwp_vec[i]
                             + paramepsilonap * Vap_vec[i] + SV2_vec[i])
        - paramsigma * E2_vec[i] - (pdelta[i] + pmu[i])*E2_vec[i];
      // Primary infected
      dI2[i] = paramsigma * E2_vec[i] - paramgamma * I2_vec[i] -
        (pdelta[i] + pmu[i])*I2_vec[i];
      // Recovered
      dR[i] = paramgamma * (I1_vec[i] + I2_vec[i]) -
        1/paramomegaI * R_vec[i] - (pdelta[i] + pmu[i])*R_vec[i];

      // Vaccinated
      dVwp[i] = (1-prop_ap)*pcoveragep * parameff * (S1_vec[i]+S2_vec[i])
        - 1/(paramomegaI*Ratio_Vwp_I) * Vwp_vec[i] -
        plambda[i] *paramepsilonwp * Vwp_vec[i] - (pdelta[i] + pmu[i])*Vwp_vec[i];
      // Vaccinated
      dVap[i] = prop_ap*pcoveragep * parameff * (S1_vec[i]+S2_vec[i])
        + pcoverageb * parameff * (S1_vec[i]+S2_vec[i]+SV1_vec[i]+SV2_vec[i])
        - 1/(paramomegaI*Ratio_Vap_I) * Vap_vec[i] -
        plambda[i] *paramepsilonap * Vap_vec[i]
        - (pdelta[i] + pmu[i])*Vap_vec[i];

      // Vaccine failure
      dSV1[i] = pcoveragep * (1-parameff) * S1_vec[i]
        + pcoverageb * (1-parameff) * S1_vec[i]
        - plambda[i] * SV1_vec[i]- pcoverageb * parameff * SV1_vec[i]
        - (pdelta[i] + pmu[i])*SV1_vec[i];
      // Vaccine failure
      dSV2[i] = pcoveragep * (1-parameff) * S2_vec[i]
        + pcoverageb * (1-parameff) * S2_vec[i]
        - plambda[i] * SV2_vec[i] - pcoverageb * parameff * SV2_vec[i]
        - (pdelta[i] + pmu[i])*SV2_vec[i];
    } else {
      // Primary susceptible
      dS1[i] = pdelta[i-1]*S1_vec[i-1] - plambda[i] * S1_vec[i] - pcoveragep * S1_vec[i] - pcoverageb * S1_vec[i] - (pdelta[i] + pmu[i])*S1_vec[i];
      // Primary exposed
      dE1[i] = pdelta[i-1]*E1_vec[i-1] + plambda[i] * (S1_vec[i]+SV1_vec[i])
        - paramsigma * E1_vec[i] - (pdelta[i] + pmu[i])*E1_vec[i];
      // Primary infected
      dI1[i] = pdelta[i-1]*I1_vec[i-1] + paramsigma * E1_vec[i]
        - paramgamma * I1_vec[i] - (pdelta[i] + pmu[i])*I1_vec[i];

      // Primary susceptible
      dS2[i] = pdelta[i-1]*S2_vec[i-1]  - plambda[i] * S2_vec[i]
        + 1/paramomegaI * R_vec[i] + 1/(paramomegaI*Ratio_Vap_I) * Vap_vec[i]
        + 1/(paramomegaI*Ratio_Vwp_I) * Vwp_vec[i]
        - (pcoveragep+pcoverageb) * S2_vec[i]
        - (pdelta[i] + pmu[i])*S2_vec[i];
      // Primary exposed
      dE2[i] = pdelta[i-1]*E2_vec[i-1]
        + plambda[i] * (S2_vec[i] + paramepsilonwp * Vwp_vec[i] + paramepsilonap * Vap_vec[i] + SV2_vec[i])
        - paramsigma * E2_vec[i] - (pdelta[i] + pmu[i])*E2_vec[i];
      // Primary infected
      dI2[i] = pdelta[i-1]*I2_vec[i-1] + paramsigma * E2_vec[i]
        - paramgamma * I2_vec[i] - (pdelta[i] + pmu[i])*I2_vec[i];
      // Recovered
      dR[i] = pdelta[i-1]*R_vec[i-1] + paramgamma * (I1_vec[i] + I2_vec[i])
        - 1/paramomegaI * R_vec[i] - (pdelta[i] + pmu[i])*R_vec[i];

      // Vaccinated
      dVwp[i] = pdelta[i-1]*Vwp_vec[i-1] +
        (1-prop_ap) * pcoveragep * parameff * (S1_vec[i]+S2_vec[i])
        - 1/(paramomegaI*Ratio_Vwp_I) * Vwp_vec[i] -
        plambda[i] *paramepsilonwp * Vwp_vec[i] -
        (pdelta[i] + pmu[i])*Vwp_vec[i];
      // Vaccinated
      dVap[i] = pdelta[i-1]*Vap_vec[i-1] + prop_ap*pcoveragep * parameff * (S1_vec[i]+S2_vec[i])
        + pcoverageb * parameff * (S1_vec[i]+S2_vec[i]+SV1_vec[i]+SV2_vec[i])
        - 1/(paramomegaI*Ratio_Vap_I) * Vap_vec[i]
        - plambda[i] * paramepsilonap * Vap_vec[i]
        - (pdelta[i] + pmu[i])*Vap_vec[i];

      // Vaccine failure
      dSV1[i] = pdelta[i-1]*SV1_vec[i-1] + pcoveragep * (1-parameff) * S1_vec[i] +
        pcoverageb * (1-parameff) * S1_vec[i] - plambda[i] * SV1_vec[i] -
        pcoverageb * parameff * SV1_vec[i] - (pdelta[i] + pmu[i])*SV1_vec[i];
      // Vaccine failure
      dSV2[i] = pdelta[i-1]*SV2_vec[i-1] + pcoveragep * (1-parameff) * S2_vec[i]
        + pcoverageb * (1-parameff) * S2_vec[i] - plambda[i] * SV2_vec[i]
        - pcoverageb * parameff * SV2_vec[i] - (pdelta[i] + pmu[i])*SV2_vec[i];
    }

    dC1[i] = paramsigma * E1_vec[i];   // Cumulative infections
    dC2[i] = paramsigma * E2_vec[i];   // Cumulative infections

   }
")

trans_to_fn <- Csnippet("
  T_paramomegaI = logit((paramomegaI-10)/15);
  T_FOI_1 = logit(FOI_1/0.5);
  ")

trans_from_fn <- Csnippet("
  paramomegaI = 15*expit(T_paramomegaI)+10;
  FOI_1 = 0.5*expit(T_FOI_1);
  ")

This is the data used in the model and 10,000 sets of parameters.
modeldata.zip

load("./data/modeldata.RData")

simufromyear <- 1955
simutoyear <- 2014

SEIRV <- pomp(
  cfile="SEIRV_source",
  cdir = "./",
  data = data.frame(
    time=c(0,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$time/365),
    report0_=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report0_),
    report1_2=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report1_2),
    report3_4=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report3_4),
    report5_9=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report5_9),
    report10_14=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report10_14),
    report15_24=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report15_24),
    report25_=c(NA,Reported_Case[Reported_Case$year>=simufromyear & Reported_Case$year<=simutoyear,]$report25_)
  ),
  times = "time",
  t0 =-100,
  covar = covariate_table(covars1, times="time", order="constant"),
  covarnames = c("pbirth",
    paste0("mu",1:22),
    "pseas","primary_cov","booster_cov","prop_ap",
    paste0("C_", rep(1:16, each = 16), "_", rep(1:16, times = 16))),
  rinit = rinit_fn,
  skeleton = vectorfield(skeleton_fn),
  params = paramslist,
  paramnames = names(paramslist),
  partrans=parameter_trans(
    toEst = trans_to_fn,
    fromEst = trans_from_fn,
    logit = estpars[!(estpars %in% c("paramomegaI","FOI_1"))]),
  statenames = c(
    paste0(
      rep(
        c(
          "S1","E1","I1",
          "S2","E2","I2",
          "R",
          "Vwp","Vap",
          "SV1","SV2",
          "C1","C2"
        ),
        each = 22
      ),1:22
    )
  ),
  accumvars = c(paste0("C1",1:22),paste0("C2",1:22)),
  obsnames = c("report0_","report1_2","report3_4","report5_9","report10_14","report15_24","report25_"),
  dmeasure = dobs,
  verbose = FALSE
)

This is the one I selected the best one among the 10,000 randomly generated parameter values ​​and ran the trajectory, which is the result of report1_2

Finally, this is the running environment. Thank you very much!

-----------------------------------
 contents of site Rprofile file 
-----------------------------------
# Things you might want to change

# options(papersize="a4")
# options(editor="notepad")
# options(pager="internal")

# set the default help type
# options(help_type="text")
  options(help_type="html")

# set a site library
# .Library.site <- file.path(chartr("\\", "/", R.home()), "site-library")

# set a CRAN mirror
# local({r <- getOption("repos")
#       r["CRAN"] <- "http://my.local.cran"
#       options(repos=r)})

# Give a fortune cookie, but only to interactive sessions
# (This would need the fortunes package to be installed.)
#  if (interactive()) 
#    fortunes::fortune()


-----------------------------------
 sessionInfo 
-----------------------------------
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8  LC_CTYPE=Chinese (Simplified)_China.utf8    LC_MONETARY=Chinese (Simplified)_China.utf8 LC_NUMERIC=C                                LC_TIME=Chinese (Simplified)_China.utf8    

time zone: Asia/Shanghai
tzcode source: internal

attached base packages:
[1] grid      parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] doFuture_1.0.1        circumstance_0.0.10.1 egg_0.4.5             gridExtra_2.3         ggsci_3.2.0           future.apply_1.11.2   future_1.34.0         tidyr_1.3.1           reshape_0.8.9         dplyr_1.1.4           doParallel_1.0.17    
[12] iterators_1.0.14      foreach_1.5.2         openxlsx_4.2.7.1      lubridate_1.9.3       socialmixr_0.3.2      pomp_5.11             subplex_1.9           ggplot2_3.5.1         reshape2_1.4.4        plyr_1.8.9            magrittr_2.0.3       
[23] nloptr_2.1.1         

loaded via a namespace (and not attached):
 [1] dotCall64_1.2     gtable_0.3.5      spam_2.11-0       lattice_0.22-6    vctrs_0.6.5       tools_4.4.1       generics_0.1.3    curl_5.2.3        tibble_3.2.1      fansi_1.0.6       pkgconfig_2.0.3   data.table_1.16.0 lifecycle_1.0.4   farver_2.1.2     
[15] compiler_4.4.1    stringr_1.5.1     fields_16.3       munsell_0.5.1     codetools_0.2-20  maps_3.4.2        pillar_1.9.0      wpp2017_1.2-3     parallelly_1.38.0 countrycode_1.6.0 tidyselect_1.2.1  zip_2.3.1         digest_0.6.37     mvtnorm_1.3-1    
[29] stringi_1.8.4     purrr_1.0.2       listenv_0.9.1     labeling_0.4.3    colorspace_2.1-1  cli_3.6.3         utf8_1.2.4        withr_3.0.1       scales_1.3.0      oai_0.4.0         timechange_0.3.0  httr_1.4.7        globals_0.16.3    deSolve_1.40     
[43] coda_0.19-4.1     viridisLite_0.4.2 rlang_1.1.4       Rcpp_1.0.13       glue_1.7.0        xml2_1.3.6        rstudioapi_0.16.0 jsonlite_1.8.9    R6_2.5.1

[kingaa]

These codes are very complicated, so it is hard for me to say for sure, but it doesn't seem that you've provided the codes you used for the least-squares fitting. Are you sure that the vectorfield and the integrator are identical in both calculations? You can of course write a dmeasure component that corresponds to least squares. Have you verified that doing that gives the same results?

[Steven-Lau-lib]

Thank you for your reply. Indeed, this code is quite long because there are many compartments and multiple age-related parameters.
First, I will add the random sampling method.

estpars <- c("FOI_1","FOI_2_1","FOI_3_2",
             "FOI_sec",
             "parameff",
             "paramomegaI",
             "Ratio_Vwp_I","Ratio_Vap_I",
             "rho1","Ratio_2_1","Ratio_3_1","Ratio_4_1","Ratio_sec",
             "paramepsilonwp","paramepsilonap",
             "k")

start.range <- subset(data.frame(FOI_1=c(0,0.5),
                                 FOI_2_1=c(0,1),
                                 FOI_3_2=c(0,1),
                                 FOI_sec = c(0,1),
                                 parameff = c(0.5,1),
                                 paramomegaI=c(10,25),
                                 Ratio_Vwp_I=c(0,1),
                                 Ratio_Vap_I=c(0,1),
                                 rho1=c(0,0.2),
                                 Ratio_2_1 = c(0,1),
                                 Ratio_3_1 = c(0,1),
                                 Ratio_4_1 = c(0,1),
                                 Ratio_sec = c(0,0.05),
                                 paramepsilonwp=c(0,0.01),
                                 paramepsilonap=c(0,0.01),
                                 k=c(0,0.1)),
                      select = estpars)

simu_time <- 1000/0.5
# simu_time <- 1000
start.values <- sobol_design(lower = setNames(as.numeric(start.range[1, ]), names(start.range[1, ])),
                             upper = setNames(as.numeric(start.range[2, ]), names(start.range[2, ])),
                             nseq = simu_time)

start.values <- start.values[start.values$Ratio_Vap_I<start.values$Ratio_Vwp_I,]
simu_time <- nrow(start.values)

ofun <- traj_objfun(SEIRV, est = estpars)

Next, the method of likelihood calculation is described.

lik <- numeric(simu_time)

for(i in 1:simu_time){
  lik[i] <- ofun(start.values[i,])
}

The following is the code for the least squares method. I refer to the previous course.

lsm <- numeric(simu_time)
for (i in 1:simu_time) {
  coef(SEIRV, estpars) <- as.numeric(start.values[i,])
  
  x <- trajectory(SEIRV,format="array")
  C1 <- x[paste0("C1",1:22),,]
  C1_vec <- rbind(apply(C1[1:3,],2,sum),
                  apply(C1[4:5,],2,sum),
                  apply(C1[6:8,],2,sum),
                  C1[9,],
                  C1[10,],
                  apply(C1[11:12,],2,sum),
                  apply(C1[13:22,],2,sum))
  C1_weights <- c(SEIRV@params["rho1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_2_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_2_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_3_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_3_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_4_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_4_1"]) 
  obs1_vec <- sweep(C1_vec, 1, C1_weights, `*`)  
  
  
  
  C2 <- x[paste0("C2",1:22),,]
  C2_vec <- rbind(apply(C2[1:3,],2,sum),
                  apply(C2[4:5,],2,sum),
                  apply(C2[6:8,],2,sum),
                  C2[9,],
                  C2[10,],
                  apply(C2[11:12,],2,sum),
                  apply(C2[13:22,],2,sum))
  C2_weights <- c(SEIRV@params["rho1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_2_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_2_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_3_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_3_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_4_1"],
                  SEIRV@params["rho1"]*SEIRV@params["Ratio_4_1"])*SEIRV@params["Ratio_sec"]  
  obs2_vec <- sweep(C2_vec, 1, C2_weights, `*`)  
  obs_vec <- obs1_vec + obs2_vec
  
  discrep <- obs_vec-obs(SEIRV)
  lsm[i] <- sum(discrep[,-1]^2)
}
start.values$lsm <- lsm

The initial compressed file already contains a set of calculated results.

[Steven-Lau-lib]

Hi Professor, I have been wondering if there is an error in my dmeasure function, so I used the Rprintf statement to output the process value.

    Rprintf(\"rho1 = %lg, Ratio_2_1 = %lg, Ratio_3_1 = %lg, Ratio_4_1 = %lg | \\n\", rho1,Ratio_2_1,Ratio_3_1,Ratio_4_1);

Then I realized that in ofun, transform was done by default. So I adjusted the program and it finally worked.

ofun <- traj_objfun(SEIRV, est = estpars)

lik <- numeric(simu_time)
for(i in 1:simu_time){
  coef(SEIRV, estpars) <- as.numeric(start.values[i,estpars])
  lik[i] <- ofun(coef(SEIRV, estpars, transform = TRUE))
}
start.values$lik <- lik

I hadn't noticed before that traj_objfun would perform a transform by default.

[kingaa]

Ah very good. I'm glad you found the error yourself, using a very effective approach. Yes, in general pomp uses transformations wherever it makes sense to do so, when they have been provided.

[Steven-Lau-lib]

Thank you. Pomp is very interesting, I need more time to get familiar with and learn it.

[kingaa]

If you ever feel that the documentation is not very clear, please do feel free to suggest changes, or even just to point out places where it is confusing or not sufficiently clear. This will help improve it for everyone!