Can I implement this heterogenous model using pomp?

[ibrahimStrathclyde]

I have two questions. Bear with my naivity as I am new to pomp.

1. I kept getting S4 susbsetable error when ever i ran a trajectory and wrote x["I", ,]. An example is given below

f1 <- function (beta) {
   params <- c(Beta=beta,gamma=1,N=50000,S_0=10000,I_0=10)
   sse(params)
}
beta <- seq(from=30,to=40,by=0.5)
SSE <- sapply(beta,f1)
Error in x["I", , ] : object of type 'S4' is not subsettable

2. I have a heterogeneous suscpetible model and need to simulate it stochastically. I have an R version of it as follows. Can I write this as pomp object instead?

stoch_Kseir <- function(y = init_het, parms = parms) {
  with(as.list(c(t, y, parms)), {
    x=parms[1:K]
    S <- matrix(y[1:K],K,1)
    E <- matrix(y[(K+1):(2*K)],K,1)
    I <- matrix(y[(2*K+1):(3*K)],K,1)
    R <- matrix(y[(3*K+1):(4*K)],K,1)
    H <- matrix(0,K,1)
    dN_SE <- matrix(0,K,1)
    dN_EI <- matrix(0,K,1)
    dN_IR <- matrix(0,K,1)
    lambda <- Beta/N*(sum(0.5*E+I))
    for( i in 1:K){
      dN_SE[i] <- rbinom(n = 1,round(S[i]), prob = 1 - exp(-(lambda)*(x[i]) * delta_t))
      dN_EI[i] <- rbinom(n = 1,round(E[i]), prob = 1 - exp(-mu_EI * delta_t))
      dN_IR[i] <- rbinom(n = 1,round(I[i]), prob = 1 - exp(-mu_IR * delta_t))
      S[i] <- S[i] - dN_SE[i]
      E[i] <- E[i] + dN_SE[i] - dN_EI[i]
      I[i] <- I[i] + dN_EI[i] - dN_IR[i]
      R[i] <- R[i] + dN_IR[i]
      H[i] <- H[i]+dN_EI[i]
    }
    outdf <- data.frame(S=S, E=E, I=I, R=R, H=H) # solutions as dataframe
    ## sum up the rows for S,E, I,R
    ## The susceptibles sum up for all rows and col
    Sp <- sum(outdf[1:K,1], na.rm=T)
    Ep <- sum(outdf[1:K,2], na.rm=T)
    Ip <- sum(outdf[1:K,3], na.rm=T)
    Rp <- sum(outdf[1:K,4], na.rm = T)
    Hp <- sum(outdf[1:K,5], na.rm = T)
    z  <-  data.frame(S = Sp, E = Ep, I = Ip, R = Rp, H = Hp)
    ## z stores the single compartment results
    list(
      outdf = outdf[,-which(names(outdf) == "H")],
      z = z
    )
  })
}

N <- 100000 # Total population
S <- N - 5  # Initial susceptible
E <- 4
I <- 1  # Initial infected
R <- 0  # Initial recovered
H <- 0  # Initial hospitalized
R0 <- 3 # Infection rate

mu_EI <-1/5.5 # Transmission delta=1/5.5===1/5.5*7=1.2727
mu_IR <- 1/4 # Recovery rate gamma=1/4 per day===1.75
delta_t <- 1/7 # Time step (days)

y=ini_het is the following for 10 groups of susceptibles

y
   S1    S2    S3    S4    S5    S6    S7    S8    S9   S10    E1    E2    E3    E4    E5    E6    E7    E8 
 3297  6745 10040 13209 16610 16741 13259  9977  6740  3377     0     0     0     1     2     1     0     0 
   E9   E10    I1    I2    I3    I4    I5    I6    I7    I8    I9   I10    R1    R2    R3    R4    R5    R6 
    0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0 
   R7    R8    R9   R10 
    0     0     0     0

parms <- c(x=x, K=10, N=N, Beta=Beta, mu_EI=mu_EI, mu_IR=mu_IR, delta_t=delta_t)
x
 [1] 6.145577e-05 1.967744e-03 1.500281e-02 6.418515e-02 2.044011e-01 5.650154e-01 1.169170e+00 2.105729e+00
 [9] 3.605283e+00 7.642741e+00

[kingaa]

@ibrahimStrathclyde:

With respect to your first question, you have not told us what sse is. Without this, we have no way of knowing what has gone wrong. Generally speaking, the error you see is generated whenever you try to access something that is not a vector, array, or list as if it were. That is, using the square brackets is called "subsetting" and you cannot do it unless the thing you are subsetting has entries like a vector, array, or list. [In your case, an array, because sse seems to be trying to access x as if it had 3 dimensions.] Because we know nothing about sse, however, we do not have a clue as to what x is at all. Presumably, it is internal to sse.... By the way, perhaps FAQ 1.4 is a propos here?

As for your second question, you ask us whether pomp might be suitable for your needs, but you do not explain your needs, beyond showing us a snippet of your code (presumably this function is a part of some computation you would like to perform?). It would be much more straightforward were you to explain your desires first. For example:

1. What model is this code meant to implement?
2. What are the questions that you are hoping to answer?
3. What are the data that you hope will have the answers to these questions?

Even without these things, however, it is possible to say that the answer to your specific question, "Can I write this as a pomp object instead?" is very likely to be "Yes, of course." It appears that your function is some kind of stochastic simulator, and these are what pomp is built for.

Happy to discuss further, but please have a look at the FAQ and, in particular FAQ 1.1 which has the advice:

Include in your request—at a minimum—the code that you ran that produced the error, together with a transcript of the R session that gave the errors. Better still, construct a minimal example that will reproduce the error: this allows for the most efficient solution of problems.

[ibrahimStrathclyde]

Thank you very much for your responds. Pardon my errors for not giving accurate infor for the first questions. Here is the detail

1. The code in Q1 is actually your code from some of the courses you taught previous year. I am just trying to understand the pomp and use to run the example codes from pomp tutorials.

pomp(
  data=subset(niamey,community=="A",select=-community),
  times="biweek",t0=0,
  skeleton=vectorfield(
    Csnippet("
      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("Beta","gamma","N","S_0","I_0")
) -> niameyA

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

f1 <- function (beta) {
  params <- c(Beta=beta,gamma=1,N=50000,S_0=10000,I_0=10)
  sse(params)
}

beta <- seq(from=30,to=40,by=0.5)
SSE <- sapply(beta,f1)
beta.hat <- beta[which.min(SSE)]

The error is when the sapply function acesses the x which is a vector of trajectory. I will read the FAQ as you suggested.

2. As for the second question. You are right I am simulating stochastically an SEIR model equation 1 to 4 of the paper (https://doi.org/10.1007/s00285-022-01771-x) using the idea you introduced in one of the seminars for SIS2020 or some other lecture on github.

Montalban(2022)Herd immunity under individual variation and reinfection.pdf

I am doing a research (PhD project) on how variation in susceptibility supresses epidemic outbreaks. As a toy exercise I want to simulate the model in equations 1-4 of the paper (https://doi.org/10.1007/s00285-022-01771-x) and then compare it with the reduced version of the model which is equations 8-11. Equations 1-4 modelled susceptibilities as gamma distribution with mean 1 and certain coefficient of variation. What I want to do is use typical discretized gamma discretized distibutions with say 10 group and a given theoritical cofficient of variation and then use the model 8-11 to estimate the v which corresponds to the value from the gamma discritization function. I will attach R files to show you what I am doing. I actually wanted to use pomp from the onset, but dont have enough skills to implement it in pomp. I am looking forward to hearing from you soon sir. I benefited alot from the vast resources you generously share.
Rcodes.docx

[kingaa]

Thanks for clarifying @ibrahimStrathclyde. I thought that code looked familiar, although I think it has been years. Were you in one of my courses? What is SIS2020, by the way?

Since you were so good as to send the code, it is very easy to see where the problem lies. Have a look at the documentation for trajectory, in particular, at the format option. Then you will understand why the following amendment solves the problem.

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

As for the possible value of pomp in your research, it is fair to say that you could do much or all that you envision using the package. You will need to catch up with all of the many changes that have happened in pomp since the last time you looked. I recommend that you begin with the Getting Started vignette and/or the Simulation-based Inference short course.

P.S. A word of advice: don't use MS Word to edit computer code. It's not designed for that purpose. Use a text editor.

[ibrahimStrathclyde]

Thank you once more for your response and suggestions @kingaa. By SIS2020, I am actually referring to the workshop you gave in 2020 (https://kingaa.github.io/sbied/ I watched the videos) on simulation based inference- and I am not one of the students in the workshop as I only started my PhD in 2022. I got motivated to write the code I shared ( in MS word because when i drag R script and drop it shows that it is not supported) from the those lectures. What I need to further clarify is that I actually understood how to set up a basic stochastic simulator from the course for a simple SEIR as there are four state variables. But for the code I provided there are 4 times K state variables which are later aggregated to form SEIR. I earnestly seek your help in giving me a hint on how to go about it. I will keep reading the documentation but I still need some help in this. hank you.

[kingaa]

https://kingaa.github.io/pomp/FAQ.html#vector-in-csnippet

[ibrahimStrathclyde]

That was helpful. Many thanks.

[ibrahimStrathclyde]

Good day @kingaa
Sequel to your suggestions, I have the following code:

 rSim <- Csnippet("
  double *X = &X1;
  int i;
  double S[K], E[K], I[K], R[K];
  double lambda, dN_SE, dN_EI, dN_IR;
  double x[K] = {0.02020271, 0.06187540, 0.10536052,
 0.22314355, 0.44628710, 0.73396918, 1.13943428,
 1.83258146, 2.81341072, 4.50671367};

  // Calculate effective transmission rate
  double Beta = R0 / (rho / delta + 1 / gamma);

  // Calculate lambda
  lambda = Beta / N * (sum(rho * E + I));

  // Update state variables
  for (i = 0; i < K; i++) {
    dN_SE = rbinom(S[i], 1 - exp(-lambda * x[i] * delta_t));
    dN_EI = rbinom(E[i], 1 - exp(-mu_EI * delta_t));
    dN_IR = rbinom(I[i], 1 - exp(-mu_IR * delta_t));

    S[i] -= dN_SE;
    E[i] += dN_SE - dN_EI;
    I[i] += dN_EI - dN_IR;
    R[i] += dN_IR;
  }

  // Update the state vector
  for (i = 0; i < K; i++) {
    X[i] = S[i];
    X[K + i] = E[i];
    X[2 * K + i] = I[i];
    X[3 * K + i] = R[i];
  }
")
rInit <- Csnippet("
  double *X = &X1;
  int i;
  double q[K] = {0.03333333, 0.06666667, 0.10000000,
 0.13333333, 0.16666667, 0.16666667, 0.13333333, 0.10000000,
 0.06666667, 0.03333333};

  double init_s[K], init_e[K], init_i[K], init_r[K];

  // Initialize state variables
  for (i = 0; i < K; i++) {
    init_s[i] = rmultinom(S0, q);
    init_e[i] = rmultinom(E0, q);
    init_i[i] = rmultinom(I0, q);
    init_r[i] = 0.0;
  }

  // Update the state vector
  for (i = 0; i < K; i++) {
    X[i] = init_s[i];
    X[K + i] = init_e[i];
    X[2 * K + i] = init_i[i];
    X[3 * K + i] = init_r[i];
  }
")

simulate(
  times=1:20,t0=0,
  rprocess=euler(rSim,delta.t=1),
  rinit=rInit,
  statenames=sprintf("X%d",1:K),
  paramnames=c("K", "R0", "rho", "delta", "gamma", "N"),
  params=c(K=10, R0=2,rho=0.5, delta=1/5.5, gamma=1/4,N=10000),
  format="data.frame"
)|>
  pivot_longer(-c(.id,time)) |>
  ggplot(aes(x=time,y=value,group=name,color=name))+
  geom_line()+
  theme_bw()

And the following error arised:
Error: in ‘simulate’: argument "func" is missing, with no default.
Please what func is required?

[kingaa]

I do not get the error you report when I attempt to run your codes. Please try again in a fresh session.

I do note that your call to rmultinom does not agree with the API for that function.

[ibrahimStrathclyde]

Thank you @kingaa for your kind respond as usual. I used a different computer and yes you are right. There are still issues with the rmultnomial usage. I will work on it and get back.