How to add the exp(cubic spline) function to the deterministic model？

[LBJlove1]

Dear Professor, I'm a big fan of pomp. I want to add the time-varying function beta=exp(cubic spline) to a simple deterministic SIR model, I tried many times but failed, I would like to ask for your help. Can you give me some advice or give me a simple example?

##########################
library(ggplot2)
library(tidyverse)
library(pomp)
library(plyr)
library(readxl)
datanba <- read_xlsx("C:/Users/1/Desktop/data.xlsx")
pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
      double beta;
      double AA[3];
      AA = bspline_basis(times,nbasis=4,degree=3)
      beta = exp(b0*AA[0]+b1*AA[1]+b2*AA[2]);
      
      DS = -beta*S*I/N;
      DI = beta*S*I/N-gamma*I;
      DR = gamma*I;")),
  rinit=Csnippet("
      S = S_0;
      I = I_0;
      R = N-S_0-I_0;"),
  statenames=c("S","I","R"),
  paramnames=c("b0","b1","b2","gamma","N","S_0","I_0")) -> lg

#############################
Error: error in building shared-object library from C snippets: in ‘Cbuilder’: compilation error: cannot compile shared-object library ‘C:/Users/1/AppData/Local/Temp/RtmpukQwHv/6504/pomp_3ba766659a8ad6327a277867added0ff.dll’: status = 1
compiler messages:
gcc  -I"D:/R-42~1.2/include" -DNDEBUG     -I"D:/rtools42/x86_64-w64-mingw32.static.posix/include"     -O2 -Wall  -std=gnu99 -mfpmath=sse -msse2 -mstackrealign  -c C:/Users/1/AppData/Local/Temp/RtmpukQwHv/6504/pomp_3ba766659a8ad6327a277867added0ff.c -o C:/Users/1/AppData/Local/Temp/RtmpukQwHv/6504/pomp_3ba766659a8ad6327a277867added0ff.o
C:/Users/1/AppData/Local/Temp/RtmpukQwHv/6504/pomp_3ba766659a8ad6327a277867added0ff.c: In function '__pomp_skelfn':
C:/Users/1/AppData/Local/Temp/RtmpukQwHv/6504/pomp_3ba766659a8ad6327a277867added0ff.c:62:12: warning: implicit declaration of function 'bspline_basis' [-Wimplicit-function-declaration]
   62 |       AA = bspline_basis(times,nbasis=4,degree=3);
      |            ^~~~~~~~~~~~~

[kingaa]

That's something that is quite easily accomplished. As the help page (?bspline_basis) points out, there is access to B-splines in the pomp C API. Have you had a look there?

Keep in mind that C has syntax that is different from that of R. In particular, your R-style call to bspline_basis has two problems:

1. There is no bspline_basis function in the pomp C API.
2. Arguments in C are passed by position only, never by name.

Try your hand at using the bspline_eval function and let me know if you run into trouble.

An alternative approach, which may prove faster, is to define the B-spline basis as a matrix using a call to pomp::bspline_basis, then to pass this as a covariate to whichever pomp elementary or estimation functions you are using. The basis functions you define in this way will simply be available to your C snippets just as the state variables and parameters are. There are examples of this included in the package. In particular, the sir, sir2, and dacca examples all include (periodic) B-spline bases among their covariates. Look at the source code to see how this is done.

[LBJlove1]

Very emotional your reply, I tried to build the following simple model, but encountered a problem in estimating parameters with maximum likelihood, I don't know if it is a problem with my model or a problem with the optimization function optim, thank you again!

library(tidyverse)
library(pomp)

time <- 0:14
measles <- c(1,3,7,25,72,222,282,256,233,189,123,70,25,11,4)
datanba <- data.frame(time,measles)
pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
       double beta;
       double seas[9];
       int nbasis = 9;

       beta = exp(b0*seas[0]+b1*seas[1]*b2*seas[2]+
               b3*seas[3]+b4*seas[4]+b5*seas[5]+
               b6*seas[6]+b7*seas[7]+b8*seas[8]);

       DS = -beta*S*I/N;
       DI = beta*S*I/N-gamma*I;
       DR = gamma*I;")),
  rinit=Csnippet("
       S = S_0;
       I = I_0;
       R = N-S_0-I_0;"),
  covar = covariate_table(
    t=seq(from=0,to=15,by=0.1),
    seas = bspline_basis(
      x=t,
      nbasis=9,
      degree=3,
      ),
    times = "t"
  ),
  statenames=c("S","I","R"),
  paramnames=c(
    "b0","b1","b2","b3","b4","b5","b6","b7","b8",
    "gamma","N","S_0","I_0")
) -> lg

poisson.loglik <- function (params) {
  x <- trajectory(lg,params=params,format="array")
  sum(dpois(x=obs(lg),lambda=x["I",,],log=TRUE))
}

f2 <- function (par) {
  params <- c(
    b0=exp(par[1]),b1=exp(par[2]),
    b2=exp(par[3]),b3=exp(par[4]),
    b4=exp(par[5]),b5=exp(par[6]),
    b6=exp(par[7]),b7=exp(par[8]),b8=exp(par[9]),
    gamma=expit(par[10]),
    N=763,S_0=762,I_0=1
  )
  -poisson.loglik(params)
}

optim(
  fn=f2,
  par=c(log(0.6),log(0.6),log(0.6),log(0.6),
    log(0.6),log(0.6),log(0.6),log(0.6),log(0.6),
    logit(0.6))
) -> fit2

coef(lg) <- c(
  b0=exp(fit2$par[1]),b1=exp(fit2$par[2]),
  b2=exp(fit2$par[3]),b3=exp(fit2$par[4]),b4=exp(fit2$par[5]),
  b5=exp(fit2$par[6]),b6=exp(fit2$par[7]),b7=exp(fit2$par[8]),
  b8=exp(fit2$par[9]),
  gamma=exp(fit2$par[10]),N=763,S_0=762,I_0=1
)

x <- trajectory(lg,format ="data.frame")

ggplot(data=join(as.data.frame(lg),x,by='time'),
  mapping=aes(x=time))+
  geom_point(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+
  ylab("I")

[LBJlove1]

And I can't use bspline_eval function in a C fragment, only periodic_bspline_basis_eval function, I don't understand what the reason is.

[kingaa]

I'm not sure what the precise problem was, since you included no error message, but I can see that you are not accessing the seasonality covariates you defined. The array seas[9] that you declare is local (and uninitialized): using it will result in undefined behavior. See FAQ 3.1 for advice on how to include a vector of covariates. You can also, per my advice above, mimic the sir, sir2, and dacca examples.

Before posting, please read FAQ 1.1. Also, note that I have edited your code for clarity.

Finally, I see that you are constructing an objective function for trajectory matching. Why not use pomp's facility for this, traj_objfun? There is an example of its use in the Getting Started vignette and in the manual.

[LBJlove1]

Thank you very much for your guidance, but when I use the function bspline_eval, the program reports an error, but when I use the function periodic_bspline_basis_eval, the program can run and get the result, the code is as follows:

library(ggplot2)
library(tidyverse)
library(pomp)
library(plyr)
library(readxl)
time <- 0:14
measles <- c(1,3,7,25,72,222,282,256,233,189,123,70,25,11,4)
datanba <- data.frame(time,measles)
pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
       double beta3;
       double seas3[9];
       int nbasis = 9;
       bspline_eval(seas3,t,15,9,3,0,1);
       beta3 = exp(b0*seas3[0]+b1*seas3[1]*b2*seas3[2]+b3*seas3[3]+b4*seas3[4]+
                b5*seas3[5]+b6*seas3[6]+b7*seas3[7]+b8*seas3[8]);
       DS = -beta3*S*I/N;
       DI = beta3*S*I/N-gamma*I;
       DR = gamma*I;")),
    rinit=Csnippet("
       S = S_0;
       I = I_0;
       R = N-S_0-I_0;"),
     statenames=c("S","I","R"),
     paramnames=c("b0","b1","b2","b3","b4","b5","b6","b7","b8","gamma","N","S_0","I_0")) -> lg
 
logit <- function(p) 
       log(p/(1-p)) 
 expit <- function(x) 
     1/(1+exp(-x))
 
poisson.loglik <- function (params) {
       x <- trajectory(lg,params=params,format="array")
       sum(dpois(x=obs(lg),lambda=x["I",,],log=TRUE))
     }
 
   
f6 <- function (par) {
       params <- c(b0=exp(par[1]),b1=exp(par[2]),b2=exp(par[3]),b3=exp(par[4]),
                                   b4=exp(par[5]),b5=exp(par[6]),b6=exp(par[7]),
                   b7=exp(par[8]),b8=exp(par[9]),gamma=expit(par[10]),
                                   N=763,S_0=762,I_0=1)
       -poisson.loglik(params)
     }
 
optim(fn=f6,par=c(log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),logit(0.8))) -> fit2

coef(lg) <- c(b0=exp(fit2$par[1]),b1=exp(fit2$par[2]),b2=exp(fit2$par[3]),b3=exp(fit2$par[4]),b4=exp(fit2$par[5]),
              b5=exp(fit2$par[6]),b6=exp(fit2$par[7]),b7=exp(fit2$par[8]),b8=exp(fit2$par[9]),
              gamma=exp(fit2$par[10]),N=763,S_0=762,I_0=1)
 
x<-trajectory(lg,format ="data.frame")

ggplot(data=join(as.data.frame(lg),x,by='time'),
       mapping=aes(x=time))+
  geom_point(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+
  ylab("I")`

错误: in ‘trajectory’: 无法载入共享目标对象‘C:/Users/1/AppData/Local/Temp/Rtmp6BDrHb/9760/pomp_5663e0e130bb1d564c91330ce871760f.dll’：:
  LoadLibrary failure:  找不到指定的模块。

When I use function periodic_bspline_basis_eval, results as following:

`library(ggplot2)
library(tidyverse)
library(pomp)
library(plyr)
library(readxl)
time <- 0:14
measles <- c(1,3,7,25,72,222,282,256,233,189,123,70,25,11,4)
datanba <- data.frame(time,measles)
pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
       double beta3;
       double seas3[9];
       int nbasis = 9;
       periodic_bspline_basis_eval(t,6,3,9,seas3);
       beta3 = exp(b0*seas3[0]+b1*seas3[1]*b2*seas3[2]+b3*seas3[3]+b4*seas3[4]+
                b5*seas3[5]+b6*seas3[6]+b7*seas3[7]+b8*seas3[8]);
       DS = -beta3*S*I/N;
       DI = beta3*S*I/N-gamma*I;
       DR = gamma*I;")),
    rinit=Csnippet("
       S = S_0;
       I = I_0;
       R = N-S_0-I_0;"),
     statenames=c("S","I","R"),
     paramnames=c("b0","b1","b2","b3","b4","b5","b6","b7","b8","gamma","N","S_0","I_0")) -> lg
 
logit <- function(p) 
       log(p/(1-p)) 
 expit <- function(x) 
     1/(1+exp(-x))
 
poisson.loglik <- function (params) {
       x <- trajectory(lg,params=params,format="array")
       sum(dpois(x=obs(lg),lambda=x["I",,],log=TRUE))
     }
 
   
f6 <- function (par) {
       params <- c(b0=exp(par[1]),b1=exp(par[2]),b2=exp(par[3]),b3=exp(par[4]),
                                   b4=exp(par[5]),b5=exp(par[6]),b6=exp(par[7]),
                   b7=exp(par[8]),b8=exp(par[9]),gamma=expit(par[10]),
                                   N=763,S_0=762,I_0=1)
       -poisson.loglik(params)
     }
 
optim(fn=f6,par=c(log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),log(0.8),logit(0.8))) -> fit2

coef(lg) <- c(b0=exp(fit2$par[1]),b1=exp(fit2$par[2]),b2=exp(fit2$par[3]),b3=exp(fit2$par[4]),b4=exp(fit2$par[5]),
              b5=exp(fit2$par[6]),b6=exp(fit2$par[7]),b7=exp(fit2$par[8]),b8=exp(fit2$par[9]),
              gamma=exp(fit2$par[10]),N=763,S_0=762,I_0=1)
 
x<-trajectory(lg,format ="data.frame")

ggplot(data=join(as.data.frame(lg),x,by='time'),
       mapping=aes(x=time))+
  geom_point(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+
  ylab("I")`

But my goal is to use a cubic spline function to construct the infectious term beta in the model, and then calculate beta through parameter estimation. There is no need for periodic spline functions to construct. I have reviewed the file you recommended, but only have examples of periodic spline basis functions. I am not sure which function to use to find a simple cubic spline basis function. Could you please give me a simple example？

[kingaa]

Why not follow my advice (as I gave you above) and construct the bsplines using bspline_basis in R, pass these to your pomp object as covariates, and then access them in the C snippets? The only problem with the code you showed above is that you fail to access the covariates you've constructed.

If you want to access the bspline_eval function inside the C snippet, then you'll have to understand it first.

I don't know why the issue of periodic splines is entering into your thinking. I did not mention anything about them in my last reply and I understand that you want ordinary (non-periodic) splines.

[LBJlove1]

Thank you very much for your suggestion, I have successfully solved the problem. At the same time, I have the following small question：

1. I have learned the IF2 algorithm in the stochastic model and implemented it, so can a deterministic model also use this method to make estimates?
2. Can I build a deterministic model in rprocess？

Looking forward to your response

[LBJlove1]

After building an exp(spline) beta model, I had some problems with parameter estimation. The following code works and the fit works fine,

########################Data#######################
library(tidyverse)
library(plyr)
library(pomp)
library(ggsci)


time <- 1:54
reports <- c(1,1,1,1,1,1,10,23,45,70,148,155,58,178,340,297,741,521,512,520,405,321,367,390,297,
             190,92,131,73,61,62,41,35,32,46,97,34,69,43,24,31,14,26,34,71,23,14,60,15,10,8,5,3,1)
datanba <- data.frame(time,reports)

################### Construct a deterministic SIR model #####################

pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
        double beta;
         int nbasis = 9;
  
         beta = exp(b0**&seas_1+b1**&seas_2+b2**&seas_3+
                 b3**&seas_4+b4**&seas_5+b5**&seas_6+
                 b6**&seas_7+b7**&seas_8+b8**&seas_9);
 
         DS = -beta*S*I/N;
         DI = beta*S*I/N-gamma*I;
         DR = gamma*I;")),
  rinit=Csnippet("
         S = S_0;
         I = I_0;
         R = N-S_0-I_0;"),
  covar = covariate_table(
    t=seq(from=0,to=55,by=0.1),
    seas = bspline_basis(
      x=t,
      nbasis=9,
      degree=3,
    ),
    times = "t"
  ),
  statenames=c("S","I","R"),
  paramnames=c(
    "b0","b1","b2","b3","b4","b5","b6","b7","b8",
    "gamma","N","S_0","I_0")
) -> lg

###################Objective function##################

sse <- function (params) {
  x <- trajectory(lg,params=params,format="array")
  discrep <- x["I",,]-obs(lg)
  sum(discrep^2)
}

f2 <- function (par) {
  params <- c(
    b0=exp(par[1]),b1=exp(par[2]),
    b2=exp(par[3]),b3=exp(par[4]),
    b4=exp(par[5]),b5=exp(par[6]),
    b6=exp(par[7]),b7=exp(par[8]),b8=exp(par[9]),
    gamma=expit(par[10]),N=50000,S_0=50000-exp(par[11]),I_0=exp(par[11])
  )
  sse(params)
}

######################Minimizing objective function##################

optim(
  fn=f2,
  par=c(log(2),log(0.1),log(0.1),log(0.2),
        log(0.1),log(0.1),log(0.1),log(0.1),log(0.1),logit(0.1),log(2)
  ),
  method = "Nelder-Mead",
  control = list(maxit = 20000)
) -> fit2

coef(lg) <- c(
  b0=exp(fit2$par[1]),b1=exp(fit2$par[2]),
  b2=exp(fit2$par[3]),b3=exp(fit2$par[4]),b4=exp(fit2$par[5]),
  b5=exp(fit2$par[6]),b6=exp(fit2$par[7]),b7=exp(fit2$par[8]),
  b8=exp(fit2$par[9]),gamma=expit(fit2$par[10]),
  N=50000,S_0=50000-exp(fit2$par[11]),I_0=exp(fit2$par[11])
)

x <- trajectory(lg,format ="data.frame")


ggplot(data=join(as.data.frame(lg),x,by='time'),
       mapping=aes(x=time))+
  geom_point(aes(y=reports),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+
  ylab("I")+
  theme_bw()

##############Visualization##############

data=join(as.data.frame(lg),x,by='time')
model.pred <- data$I
raply(10000,rnbinom(n=length(model.pred),mu = model.pred,size=10)) ->simdat


aaply(simdat,2,quantile,probs=c(0.025,0.5,0.975))->quantiles


times <- 1:54
quantiles  <- data.frame(times,quantiles)

colnames(quantiles) <- c("times","lower","mean","higher")
data <- data.frame(quantiles,datanba$reports)


ggplot(data,aes(x=times))+
  geom_line(aes(y=mean),color=pal_jco()(1))+
  geom_ribbon(aes(ymin=lower,ymax=higher),fill=pal_jco()(1),alpha=0.3)+
  geom_point(aes(y=datanba.reports),color="red")+
  theme_bw()

coef(lg)
          b0           b1           b2           b3           b4           b5           b6           b7           b8 
8.063935e-04 4.585929e-03 4.767480e-01 8.773362e-04 8.388647e-05 3.511879e-04 1.589562e-01 8.800913e-03 1.234533e-03 
       gamma            N          S_0          I_0 
1.000000e+00 5.000000e+04 4.999312e+04 6.877639e+00 

but when I set gamma to 0.2, S to 49999 and estimate the other parameters, I get the following error. I just changed the values of two parameters, and the result was that I didn't get the correct result and got a warning message. I don't understand what is the reason for this, is it because after changing the parameters, the integrator does not work? And when using traj_objfun functions to make estimates, the following problems will also be encountered.
Please see warnings:

1: In lsoda(y, times, func, parms, ...) :
  repeated convergence test failures on a step, but integration was successful - inaccurate Jacobian matrix?
2: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go
3: in ‘trajectory’: abnormal exit from ODE integrator, istate = -5
4: In lsoda(y, times, func, parms, ...) :
  repeated convergence test failures on a step, but integration was successful - inaccurate Jacobian matrix?
5: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go

48: in ‘trajectory’: abnormal exit from ODE integrator, istate = -5
49: In lsoda(y, times, func, parms, ...) :
  repeated convergence test failures on a step, but integration was successful - inaccurate Jacobian matrix?
50: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go

########################Data#######################
library(tidyverse)
library(dplyr)
library(plyr)
library(pomp)
library(ggsci)


time <- 1:54
reports <- c(1,1,1,1,1,1,10,23,45,70,148,155,58,178,340,297,741,521,512,520,405,321,367,390,297,
             190,92,131,73,61,62,41,35,32,46,97,34,69,43,24,31,14,26,34,71,23,14,60,15,10,8,5,3,1)
datanba <- data.frame(time,reports)

################### Construct a deterministic SIR model #####################

pomp(
  data=datanba,
  times="time",t0=0,
  skeleton=vectorfield(
    Csnippet("
        double beta;
         int nbasis = 9;
  
         beta = exp(b0**&seas_1+b1**&seas_2+b2**&seas_3+
                 b3**&seas_4+b4**&seas_5+b5**&seas_6+
                 b6**&seas_7+b7**&seas_8+b8**&seas_9);
 
         DS = -beta*S*I/N;
         DI = beta*S*I/N-gamma*I;
         DR = gamma*I;")),
  rinit=Csnippet("
         S = S_0;
         I = I_0;
         R = N-S_0-I_0;"),
  covar = covariate_table(
    t=seq(from=0,to=55,by=0.1),
    seas = bspline_basis(
      x=t,
      nbasis=9,
      degree=3,
    ),
    times = "t"
  ),
  statenames=c("S","I","R"),
  paramnames=c(
    "b0","b1","b2","b3","b4","b5","b6","b7","b8",
    "gamma","N","S_0","I_0")
) -> lg

###################Objective function##################

sse <- function (params) {
  x <- trajectory(lg,params=params,format="array")
  discrep <- x["I",,]-obs(lg)
  sum(discrep^2)
}

f2 <- function (par) {
  params <- c(
    b0=exp(par[1]),b1=exp(par[2]),
    b2=exp(par[3]),b3=exp(par[4]),
    b4=exp(par[5]),b5=exp(par[6]),
    b6=exp(par[7]),b7=exp(par[8]),b8=exp(par[9]),
    gamma=0.2,N=50000,S_0=49999,I_0=1
  )
  sse(params)
}

######################Minimizing objective function##################

optim(
  fn=f2,
  par=c(log(2),log(0.1),log(0.1),log(0.2),
        log(0.1),log(0.1),log(0.1),log(0.1),log(0.1)
  ),
  method = "Nelder-Mead",
  control = list(maxit = 20000)
) -> fit2

coef(lg) <- c(
  b0=exp(fit2$par[1]),b1=exp(fit2$par[2]),
  b2=exp(fit2$par[3]),b3=exp(fit2$par[4]),b4=exp(fit2$par[5]),
  b5=exp(fit2$par[6]),b6=exp(fit2$par[7]),b7=exp(fit2$par[8]),
  b8=exp(fit2$par[9]),gamma=0.2,
  N=50000,S_0=49999,I_0=1
)

x <- trajectory(lg,format ="data.frame")

ggplot(data=join(as.data.frame(lg),x,by='time'),
       mapping=aes(x=time))+
  geom_point(aes(y=reports),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+
  ylab("I")+
  theme_bw()

[kingaa]

With respect to the questions you ask above,

1. I have learned the IF2 algorithm in the stochastic model and implemented it, so can a deterministic model also use this method to make estimates?

It might be possible, but this is not the case for which the algorithm is best suited. When the latent state process is deterministic, I recommend you use trajectory matching, via the traj_objfun method in pomp.

2. Can I build a deterministic model in rprocess？

In principle, the answer is "yes". Note that the goal of an rprocess C snippet is different from that of a skeleton C snippet. Read the manual for details.

Now for the more recent questions.

Before addressing the issue you ran into, note that you can streamline your workflow quite a bit (and eliminate sources of common errors) by using traj_objfun to construct an objective function. Compare the following to your code. In particular, note the simpler skeleton snippet, the integration of your error metric via a dmeasure snippet, and the integration of parameter transformations into the objective function.

library(tidyverse)
library(pomp)

data.frame(
  time=1:54,
  reports=c(
    1,1,1,1,1,1,10,23,45,70,148,155,58,178,340,297,741,521,512,
    520,405,321,367,390,297,190,92,131,73,61,62,41,35,32,46,97,
    34,69,43,24,31,14,26,34,71,23,14,60,15,10,8,5,3,1
  )
) |>
  traj_objfun(
    times="time",t0=0,
    skeleton=vectorfield(
      Csnippet("
      double beta;
      beta = exp(b0*seas_1+b1*seas_2+b2*seas_3+
                 b3*seas_4+b4*seas_5+b5*seas_6+
                 b6*seas_7+b7*seas_8+b8*seas_9);
      DS = -beta*S*I/N;
      DI = beta*S*I/N-gamma*I;
      DR = gamma*I;"
      )
    ),
    rinit=Csnippet("
      S = N*(1-I_0);
      I = N*I_0;
      R = 0;"
    ),
    dmeasure=Csnippet("
      double discrep = I-reports;
      if (give_log) lik = -discrep*discrep;
      else lik = exp(-discrep*discrep);"
    ),
    covar = covariate_table(
      t=seq(from=0,to=55,by=0.1),
      seas = bspline_basis(x=t,nbasis=9,degree=3),
      times = "t"
    ),
    statenames=c("S","I","R"),
    paramnames=c(
      paste0("b",0:8),
      "gamma","N","I_0"
    ),
    partrans=parameter_trans(
      log=c(paste0("b",0:8),"gamma"),
      logit=c("I_0")
    ),
    params=c(
      b0=2.0,b1=0.1,b2=0.1,
      b3=0.2,b4=0.1,b5=0.1,
      b6=0.1,b7=0.1,b8=0.1,
      gamma=0.2,N=50000,
      I_0=1/50000
    ),
    est=c(paste0("b",0:8),"gamma","I_0")
  ) -> f2

Now we fit to the data.

fit2 <- optim(
  fn=f2,
  par=coef(f2,c(paste0("b",0:8),"gamma","I_0"),transform=TRUE),
  control=list(maxit=10000)
)
fit2$convergence
f2(fit2$par)
coef(f2)

bind_cols(
  f2 |> trajectory(format="d"),
  reports = f2 |> as_pomp() |> obs() |> t()
) |>
  ggplot(aes(x=time))+
  geom_point(aes(y=reports),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+ylab("I")+
  theme_bw()

With fixed $\gamma$ and $I_0$:

f2 |>
  traj_objfun(
    params=c(
      b0=2.0,b1=0.1,b2=0.1,
      b3=0.2,b4=0.1,b5=0.1,
      b6=0.1,b7=0.1,b8=0.1,
      gamma=1,N=50000,
      I_0=1/50000
    ),
    est=c(paste0("b",0:8))
  ) -> f3

fit3 <- optim(
  fn=f3,
  par=coef(f3,paste0("b",0:8),transform=TRUE),
  control=list(maxit=10000)
)
fit3$convergence
f3(fit3$par)
coef(f3)

bind_cols(
  f3 |> trajectory(format="d"),
  reports = f3 |> as_pomp() |> obs() |> t()
) |>
  ggplot(aes(x=time))+
  geom_point(aes(y=reports),color='black')+
  geom_line(aes(y=I),color='red')+
  xlab("time")+ylab("I")+
  theme_bw()

Note that the issue you encountered no longer appears. The warnings you were seeing were generated by the ODE integrator. Presumably, the problem was due to inadmissible values in your parameters or initial conditions. These in turn were probably due to the inappropriate parameter transformations you were applying.

[LBJlove1]

Thank you very much for your reply and the code modification. After successfully estimating the parameters b0-b8, I tried to visualize beta(t), but I found that R0=beta/gamma is always less than 1, is this phenomenon reasonable?

coef(f2)
          b0           b1           b2           b3           b4           b5           b6 
1.757295e-01 1.278786e-04 8.250521e-01 2.059958e-02 1.285904e-02 2.171483e-01 6.859711e-02 
          b7           b8        gamma            N          I_0 
7.838992e-02 6.433913e-01 1.092065e+00 5.000000e+04 1.052940e-05 

[kingaa]

How are you computing $R_0$? Is there a right way? Is there a unique correct formula? Is there a right way in which $R_0<1$ and one has a deterministic outbreak? I leave these things for you to contemplate. This is not the place to get into a discussion of that.

I am closing this issue now. Feel free to reopen if you have more issues with pomp, or to move to the Discussions forum.

[LBJlove1]

I have encountered some problems again. When I use pomp to build a stochastic model to simulate the COVID-19 epidemic, no matter how I change the model parameters, the following error is reported in pfilter:

library(tidyverse)
library(pomp)

time <- 1:23
reports <- c(2,4,10,12,26,40,32,58,52,36,38,48,54,37,25,18,18,6,3,6,3,2,1)
covid_yz <- data.frame(time,reports)

covid_yz |>
  ggplot(aes(x=time,y=reports))+
  geom_line()+
  geom_point()

seir_step <- Csnippet("
  double rate[12];
  double dN[12];

  double Beta;

  Beta = exp(b0**&seas_1+b1**&seas_2+b2**&seas_3+
             b3**&seas_4+b4**&seas_5+b5**&seas_6+
             b6**&seas_7+b7**&seas_8+b8**&seas_9);

  rate[0] = Beta*(1-pd)*(1-em)*pm*(I+K1*E+K2*A)/N;
  rate[1] = (1-Beta)*q*(1-pd)*(1-em)*pm*(I+K1*E+K2*A)/N;
  rate[3] = lambda;
  rate[4] = f1*sigma;
  rate[5] = (1-f1)*sigma;
  rate[6] = gamma2;
  rate[7] = gamma1*rho;
  rate[8] = gamma1*(1-rho);
  rate[9] = b*(1-f2);
  rate[10] = b*f2;
  rate[11] = gamma3;

  reulermultinom(2,S,&rate[0],dt,&dN[0]);
  dN[0]=(1-q)*dN[0];
  dN[2]=q*dN[0];
  reulermultinom(1,Sq,&rate[3],dt,&dN[3]);
  reulermultinom(2,E,&rate[4],dt,&dN[4]);
  reulermultinom(1,A,&rate[6],dt,&dN[6]);
  reulermultinom(2,I,&rate[7],dt,&dN[7]);
  reulermultinom(2,Eq,&rate[9],dt,&dN[9]);
  reulermultinom(1,H,&rate[11],dt,&dN[11]);

  S += -dN[0]-dN[1]-dN[2]+dN[3]+dN[9];
  Sq += dN[1]-dN[3];
  E += dN[0]-dN[4]-dN[5];
  A += dN[4]-dN[6];
  I += dN[5]-dN[7]-dN[8];
  Eq += dN[2]-dN[9]-dN[10];
  H += dN[8]+dN[10]-dN[11];
  R += dN[6]+dN[7]+dN[11];
")

seir_init <- Csnippet("
  S = N-1000;
  Sq = 0;
  E = 0;
  A = 0;
  I = 1000;
  Eq = 0;
  H = 0;
  R = 0;
")


seir_dmeas <- Csnippet("
  lik = dnbinom_mu(reports,k,I,give_log);
")

seir_rmeas <- Csnippet("
  reports = rnbinom_mu(k,I);
")

covar <- covariate_table(
  t=seq(from=0,to=23,by=0.1),
  seas = bspline_basis(
    x=t,
    nbasis=9,
    degree=3,
    ),
  times = "t"
)

covid_yz |>
  pomp(
    times = "time",
    t0 = 0,
    rprocess = euler(seir_step,delta.t=1),
    rinit = seir_init,
    rmeasure = seir_rmeas,
    dmeasure = seir_dmeas,
    covar = covar,
    paramnames = c("N","b0","b1","b2","b3","b4","b5","b6","b7","b8",
      "q","pd","em","pm","K1","K2","lambda","f1","sigma","gamma2","gamma1",
      "rho","b","f2","gamma3","k"),
    statenames = c("S","Sq","E","A","I","Eq","H","R")
  ) -> covid_yz

covid_yz |>
  simulate(
    params=c(
      N=4515600,
      b0=10,b1=-0.7,b2=-0.7,b3=-1.17,b4=-1.16,
      b5=-1.48,b6=-1.733,b7=-1,b8=2.36,
      q=0.2,pd=0.2,em=0.2,pm=0.2,
      K1=0.2,K2=0.2,
      lambda=0.2,f1=0.2,sigma=0.2,
      gamma1=0.2,gamma2=0.2,gamma3=0.2,
      rho=0.02,b=0.2,f2=0.2,k=100
    ),
    nsim=200,
    format="data.frame",
    include.data=TRUE
  ) -> sims

sims |>
  ggplot(aes(x=time,y=reports,group=.id,color=.id=="data"))+
  geom_line()+
  guides(color="none")

covid_yz |>
  pfilter(
    Np=5000,
    params=c(
      N=4515600,
      b0=10,b1=-0.7,b2=-0.7,b3=-1.17,b4=-1.16,b5=-1.48,b6=-1.733,b7=-1,b8=2.36,
      q=0.2,pd=0.2,em=0.2,pm=0.2,K1=0.2,K2=0.2,
      lambda=0.2,f1=0.2,sigma=0.2,gamma2=0.2,gamma1=0.2,
      rho=0.02,b=0.2,f2=0.2,gamma3=0.2,k=100
    )
    ) -> pf

The error is as follows, but I am not sure why the value of state S is NA. Is it a problem with the model parameter settings or the transition probability between models.

错误: in ‘pfilter’: ‘dmeasure’ with log=TRUE returns illegal value.
Log likelihood, data, states, and parameters are:
   time:            2
 loglik:           NA
reports:            4
      S:           NA
     Sq:           NA
      E:           NA
      A:           NA
      I:           NA
     Eq:           NA
      H:           NA
      R:           50
      N:   4.5156e+06
     b0:           10
     b1:         -0.7
     b2:         -0.7
     b3:        -1.17
     b4:        -1.16
     b5:        -1.48
     b6:       -1.733
     b7:           -1
     b8:         2.36
      q:          0.2
     pd:          0.2
     em:          0.2
     pm:          0.2
     K1:          0.2
     K2:          0.2
 lambda:          0.2
     f1:          0.2
  sigma:          0.2
 gamma2:          0.2
 gamma1:          0.2
    rho:         0.02
      b:          0.2
     f2:          0.2
 gamma3:          0.2
      k:          100
此外: There were 50 or more warnings (use warnings() to see the first 50)

[kingaa]

Thanks for providing sufficiently so that I can easily reproduce the problem. Notice that something is already not right in the simulate call:

> sims |> filter(.id!="data") |> head()
  time .id reports     seas_1    seas_2    seas_3      seas_4 seas_5 seas_6
1    1   1     908 0.06729953 0.6074902 0.3222514 0.002958823      0      0
2    1   2     836 0.06729953 0.6074902 0.3222514 0.002958823      0      0
3    1   3     819 0.06729953 0.6074902 0.3222514 0.002958823      0      0
4    1   4     714 0.06729953 0.6074902 0.3222514 0.002958823      0      0
5    1   5     911 0.06729953 0.6074902 0.3222514 0.002958823      0      0
6    1   6     813 0.06729953 0.6074902 0.3222514 0.002958823      0      0
  seas_7 seas_8 seas_9  S Sq  E A   I Eq   H R
1      0      0      0 NA NA NA 0 833 NA 165 2
2      0      0      0 NA NA NA 0 852 NA 146 2
3      0      0      0 NA NA NA 0 813 NA 185 2
4      0      0      0 NA NA NA 0 815 NA 183 2
5      0      0      0 NA NA NA 0 806 NA 186 8
6      0      0      0 NA NA NA 0 804 NA 192 4

I looked into it a little. Examine the following code. It prints the values of the rate array. From the very beginning, rate[1] is negative, which results in NaN values.

seir_step <- Csnippet(r"{
  double rate[12];
  double dN[12];

  double Beta = exp(b0*seas_1+b1*seas_2+b2*seas_3+
                    b3*seas_4+b4*seas_5+b5*seas_6+
                    b6*seas_7+b7*seas_8+b8*seas_9);

  rate[0] = Beta*(1-pd)*(1-em)*pm*(I+K1*E+K2*A)/N;
  rate[1] = (1-Beta)*q*(1-pd)*(1-em)*pm*(I+K1*E+K2*A)/N;
  rate[2] = q*rate[0];
  rate[0] -= rate[2];
  rate[3] = lambda;
  rate[4] = f1*sigma;
  rate[5] = (1-f1)*sigma;
  rate[6] = gamma2;
  rate[7] = gamma1*rho;
  rate[8] = gamma1*(1-rho);
  rate[9] = b*(1-f2);
  rate[10] = b*f2;
  rate[11] = gamma3;

  for (int i = 0; i < 12; i++)
    Rprintf("%ld %lg \n",i,rate[i]);
  Rprintf("t=%lg\n\n",t);

  reulermultinom(3,S, &rate[0], dt,&dN[0]);
  reulermultinom(1,Sq,&rate[3], dt,&dN[3]);
  reulermultinom(2,E, &rate[4], dt,&dN[4]);
  reulermultinom(1,A, &rate[6], dt,&dN[6]);
  reulermultinom(2,I, &rate[7], dt,&dN[7]);
  reulermultinom(2,Eq,&rate[9], dt,&dN[9]);
  reulermultinom(1,H, &rate[11],dt,&dN[11]);

  S  += -dN[0]-dN[1] -dN[2]+dN[3]+dN[9];
  Sq +=  dN[1]-dN[3];
  E  +=  dN[0]-dN[4] -dN[5];
  A  +=  dN[4]-dN[6];
  I  +=  dN[5]-dN[7] -dN[8];
  Eq +=  dN[2]-dN[9] -dN[10];
  H  +=  dN[8]+dN[10]-dN[11];
  R  +=  dN[6]+dN[7] +dN[11];
}")

Observe that I also changed the code in two places to (1) eliminate a needless *& operation, which has no effect (since * and & are mutually inverse) and (2) eliminate an error that would arise as soon as you had fixed this one. The second argument to reulermultinom has to be an integer (i.e., a double precision number with no integer part). Cf. rbinom.

[kingaa]

I am closing this issue now. Feel free to reopen if more discussion is warranted.