pomp 包 COVID SEIR 模型最小二乘误差回溯

Question

我尝试为英国建模 SEIR 以评估实施的遏制措施，并在此处找到了一些带有 pomp 包的代码：https://kingaa.github.io/clim-dis/parest/parest.html 我试图将其转移到我的案例中，该案例增加了一个阶段（E）和另外三个变量。最后，我想做一个最小二乘估计来找到最佳贝塔。 Data_UK_beta0 由变量 date（从 0 到 165 的整数）和 new_cases（来自约翰霍普金斯大学数据集）组成。

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["I",,]-obs(Data_UK_pomp_beta0)
  sum(discrep^2)
}


install.packages("apricom")
library(apricom)

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


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

得到以下错误回溯（不幸的是德语，但我想消息应该很清楚:)

  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) 

我做错了什么？

Answer 1

您从apricom 导入的sse 函数与此问题无关（据我所知）。 （这也与 C(++) 代码编译无关，因此您问题中的 [compiler-errors] 标签有点误导。）

您还没有给我们提供获取您的 Data_UK_beta0 数据集的方法，所以我无法重现这一点，但我假设您实际上想要这样的东西：

beta_reg <- function (beta0) {
    params <- c(beta1=beta0, delta1=1/5.1, gamma1=1, 
                eta1=0.012649, N=67886004, S_0=67886004, E_0=5, I_0=2)
    sse_UK_beta0(params)
}
beta0 <- seq(from=1,to=40,by=1)
SSE <- sapply(beta0, beta_reg)

您还应该注意，您正在做的事情中还有许多其他潜在的棘手问题。跳出来的是，您可能会将建模中的感染流行率（I 隔间中的当前数字）与病例报告数据进行比较，这是一个（滞后的、有偏差的、不完善的）疾病发病率的度量，即每单位时间的新病例数...