kingaa / pomp

R package for statistical inference using partially observed Markov processes
https://kingaa.github.io/pomp
GNU General Public License v3.0
114 stars 27 forks source link

SEIR model pomp package #130

Closed Wurzelqueen closed 3 years ago

Wurzelqueen commented 3 years ago

I try to model a SEIR for UK to evaluate the implemented containment measures and found some code of @kingaa with the pomp package here: https://kingaa.github.io/clim-dis/parest/parest.html I tried to transfer this to my case which adds one stage (E) and three more variables. In the end i want to do a least squared estimation to find the optimal beta hat.

Data_UK_pomp_beta0 <- pomp(
  data= Data_UK_beta0,
  times ="date", t0=0,
  skeleton = vectorfield(
    Csnippet("
      DS=-beta1*S*I/N;
      DE= beta1*S*I/N-delta1*E;
      DI=delta1*E-(gamma1+eta1)*I;
      DR=gamma1*I;")),
  rinit = Csnippet("
      S=S_0;
      E=E_0;
      I=I_0;
      R=N-S_0-E_0-I_0;"),
  statenames = c("S","E","I","R"),
  paramnames = c("beta1","delta1","gamma1","eta1","N","S_0","E_0","I_0")) 

sse_UK_beta0 <- function (params) {
  x <- trajectory(Data_UK_pomp_beta0,params=params)
  discrep <- x["E",,]-obs(Data_UK_pomp_beta0)
  sum(discrep^2)
}

beta_reg <- function (beta0) {
  params <- c(beta1=beta0, delta1=1/5.1, gamma1=1, eta1=0.012649, N=67886004, S_0=67880000, E_0=5, I_0=2)
  sse(params)
}

beta0 <- seq(from=1,to=40,by=1)
SSE <- sapply(beta0, beta_reg)

As result i got this traceback:

 Fehler in h(simpleError(msg, call)) : 
  Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'as.matrix': Argument "dataset" fehlt (ohne Standardwert) 
7.
h(simpleError(msg, call)) 
6.
.handleSimpleError(function (cond) 
.Internal(C_tryCatchHelper(addr, 1L, cond)), "Argument \"dataset\" fehlt (ohne Standardwert)", 
    base::quote(as.matrix(dataset))) 
5.
as.matrix(dataset) 
4.
sse(params) 
3.
FUN(X[[i]], ...) 
2.
lapply(X = X, FUN = FUN, ...) 
1.
sapply(beta0, beta_reg)

What did i do wrong here?

Thanks in advance!

kingaa commented 3 years ago

Why not use the trajectory matching functionality built into pomp?

whzhuscu commented 1 year ago

Nice to meet you! I wonder if you can estimate parameters of SEIR?