Help needed: New to pomp

[ibrahimStrathclyde]

Dear Professor King,

I am emailing to ask some help as regards your popular package POMP. I am reading the documentation online and trying to use some examples to write my codes using the pomp package. But I as a novice I have spent hours trying to figure out the errors. I don’t seem able to fix it. Below is the code:

library(pomp)
open.sir.ode <- Csnippet("
  double Beta =R0*(1/(rho/delta+1/gamma))
  DS= -beta*(rho*E+I)*(S)^(1+v^2) ;  
  DE= beta*(rho*E+I)*(S)^(1+v^2)-delta*E;
  DI=delta*E-gamma*I;
  DR=gamma*I;")

init1 <- Csnippet("
  double epsilon
  S=1-epsilon;
  E=epsilon*0.8;
  I=epsilon*0.2;
  R= 0;")

pomp(data=data.frame(time=1:46,data=0),
  times="time",t0=0,
  skeleton=vectorfield(open.sir.ode),
  rinit=init1,
  params = c(gamma=0.25,rho=0.5,Beta=0.9,delta=0.333,
    epsilon=0.0001,phi=0009,v=1),
  statenames=c("S","E","I","R"),
  paramnames=c("gamma","rho","Beta","delta","epsilon", "phi","v")
) -> open.sir

params1<-c(gamma,rho,Beta,delta,epsilon,phi,v)

x <-trajectory(open.sir,params=params1,format="data.frame")

library(ggplot2)
ggplot(data=x,mapping=aes(x=time,y=I))+geom_line()

[kingaa]

Here are some things I notice.

I do not see how R0 supposed to be set in your code.

What is beta and is it supposed to be different from Beta?

To raise a number to a power in C, one uses the pow() function.

In C, you must end a complete line with a semicolon.

You have declared epsilon and Beta to be both parameters and local variables.

If you are new to C, you should look at some introductory tutorials about the language. There are many good ones.

Compare the following:

library(ggplot2)
library(pomp)

open.sir.ode <- Csnippet("
  DS = -Beta*(rho*E+I)*pow(S,1+v*v);  
  DE =  Beta*(rho*E+I)*pow(S,1+v*v)-delta*E;
  DI = delta*E-gamma*I;
  DR = gamma*I;")

init1 <- Csnippet("
  S = 1-epsilon;
  E = epsilon*0.8;
  I = epsilon*0.2;
  R = 0;")

trajectory(
  times=1:46,
  t0=0,
  skeleton=vectorfield(open.sir.ode),
  rinit=init1,
  params = c(gamma=0.25,rho=0.5,Beta=0.9,delta=0.333,
    epsilon=0.0001,phi=0009,v=1),
  statenames=c("S","E","I","R"),
  paramnames=c("gamma","rho","Beta","delta","epsilon","phi","v"),
  format="data.frame"
) -> x

x |> 
  ggplot(aes(x=time,y=I))+
  geom_line()

[ibrahimStrathclyde]

Thank you Prof sir. It worked. I used the R0 to calculate Beta in the main code. And the beta and Beta are the same actually. confused with the example in the online code i was using.I will study C before going back to the POMP.
Thanks
Ibrahim.

[kingaa]

No need to master all of C programming. I think with just a slightly better idea of C syntax you will be able to write C snippets that do what you want!

[kingaa]

Also, when you say "it worked", I guess that you mean it ran without errors. Was it doing what you had in mind? That is another question!

[ibrahimStrathclyde]

Yes sir I mean it ran without errors. I am trying to do my first model fitting using the package. So basically I am just trying to reproduce the examples you gave in the Niamey pandemic with a different model I am working with in my PG program - before fitting my actual data. I am doing sort of "window shopping" with the package trying to learn to do something different. And thanks for the clarification with the C programming.
Ibrahim