pomp 包 COVID SEIR 模型最小二乘误差回溯
pomp package COVID SEIR model least squares errors traceback
我尝试为英国建立 SEIR 模型以评估实施的遏制措施,并在此处找到了一些带有 pomp 包的代码:https://kingaa.github.io/clim-dis/parest/parest.html
我试图将其转移到我的案例中,该案例增加了一个阶段 (E) 和另外三个变量。最后我想做一个最小二乘估计来找到最佳贝塔。
Data_UK_beta0 由变量日期(从 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)
我做错了什么?
您从 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
隔间中的当前数字)与 病例报告 数据,它是 发病率 的(滞后的、有偏差的、不完善的)衡量指标,即每单位时间的新病例数 ...
我尝试为英国建立 SEIR 模型以评估实施的遏制措施,并在此处找到了一些带有 pomp 包的代码:https://kingaa.github.io/clim-dis/parest/parest.html 我试图将其转移到我的案例中,该案例增加了一个阶段 (E) 和另外三个变量。最后我想做一个最小二乘估计来找到最佳贝塔。 Data_UK_beta0 由变量日期(从 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)
我做错了什么?
您从 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
隔间中的当前数字)与 病例报告 数据,它是 发病率 的(滞后的、有偏差的、不完善的)衡量指标,即每单位时间的新病例数 ...