在 R 中：FME/deSolve - SIR 拟合（时变参数）
In R: FME/ deSolve - SIR fitting (time varying parameters)
我正在尝试做的事情：我有一个简单的 SIR 模型，具有随时间变化的传输率 beta，我已经在 R 中实现了它（感谢@tpetzoldt）。我们有一个人口N=10000，gamma也是固定的。
sir_1 <- function(f_beta, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- f_beta(time)
gamma <- f_gamma(time)
with(as.list(variables), {
dS <- -beta * I * S/10000
dI <- beta * I * S/10000 - 1/5 * I
dR <- 1/5 * I
return(list(c(dS, dI, dR), beta=beta))
})
}
# time dependent parameter functions
parameters_values <- list(
f_beta = f_beta
)
# the initial values of variables
initial_values <- c(S = S0, I = I0, R = R0)
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
f_beta <- approxfun(x=times, y=seq(0.901, 0.92, by=0.001), rule=2)
out <- as.data.frame(sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times))
现在我有了一些“真实”数据，使用 FME 包我想在每个时间步获得最佳 beta 参数
datareal <- cbind(time = times, I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,
5356, 4987, 3421, 2789, 1789,1156))
sir_cost <- function (f_beta) {
outsir <- as.data.frame(sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times))
costf <- modCost(model = outsir, obs = datareal)
}
p <- rep(0.8, 20)
Fit <- modFit(f = sir_cost, p = p)
Fit
$par
[1] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
我的问题：
- 对于我在每个时间步取 0.8 的初始值，但是 Fit 函数什么都不做，它只returns每个时间步的 0.8（即使我取一个非常高的值，如 800，它表示这已经是最合适的）。我的猜测是同一变量 (beta) 的随时间变化的值，我必须像文档中那样以另一种方式处理这个问题。
非常感谢任何帮助。
我不认为按时间步长估计 beta 是个好主意。这是问题固有的，而不是 deSolve 或 FME 的错误。如果要使用动态模型来估计与时间相关的参数，我会建议使用节点较少的合适函数，例如时间相关的线性、二次或样条曲线，例如 3-5 节而不是 20 节。然后将 approxfun 替换为该函数并将其插入。模型拟合是一门艺术，因此请使用起始值和求解器。还有，读书。
注意以下只是技术演示:
library("deSolve")
library("FME")
sir_1 <- function(f_beta, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- parameters$f_beta(time)
with(as.list(variables), {
dS <- -beta * I * S/10000
dI <- beta * I * S/10000 - 1/5 * I
dR <- 1/5 * I
return(list(c(dS, dI, dR), beta=beta))
})
}
initial_values <- c(S = S0, I = I0, R = R0)
parameters <- list(f_beta=f_beta)
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
# use method "constant" to leave beta constant over time step
f_beta <- approxfun(x=times, y=seq(0.901, 0.92, by=0.001), method="constant", rule=2)
out <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
plot(out)
datareal <- cbind(time = times, I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,
5356, 4987, 3421, 2789, 1789,1156))
plot(out, obs=datareal)
sir_cost <- function (p) {
f_beta <- approxfun(x=times, y=p, method="constant", rule=2)
outsir <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
modCost(model = outsir, obs = datareal)
}
# Play with start values!!!
p <- rep(0.8, 20)
# e.g.: consider random start values
set.seed(123)
p <- runif(20, min=0.8, max=1.2)
# try other solvers, especially such with true box constraints
Fit <- modFit(f = sir_cost, p = p,
lower=rep(0.2, 20), upper=rep(5, 20), # box constraints
method="Port")
summary(Fit) # system is singular (that is what we expected)
# use another solver. Note: it takes a while
Fit <- modFit(f = sir_cost, p = p,
lower=rep(0.2, 20), upper=rep(5, 20), # box constraints
method="L-BFGS-B")
# goes in a surprisingly good direction
Fit$par
f_beta <- approxfun(x=times, y=Fit$par, method="constant", rule=2)
out2 <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
# compare with data
plot(out, out2, obs=datareal)
# but see how unstable beta is
plot(out2)
用随时间变化的参数拟合模型可能是个好主意，但如果有理由这样做，我建议限制参数的数量并使用一种平滑函数。
以下示例展示了如何为此目的使用样条曲线，但当然也可以（并且可能更可取）使用具有某种机械意义的函数。
作为副作用，还可以识别伽玛而不是先验地固定它。尽管如此，这仍然是一个技术演示，但我将科学问题悬而未决，即随时间变化的测试版是否有意义。
library("FME")
sir_1 <- function(f_beta, gamma, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- parameters$f_beta(time)
gamma <- parameters$gamma
with(as.list(variables), {
dS <- -beta * I * S / 10000
dI <- beta * I * S / 10000 - gamma * I
dR <- gamma * I
# return vector of derivatives, and beta as auxiliary variable
return(list(c(dS, dI, dR), beta = beta))
})
}
initial_values <- c(S = S0, I = I0, R = R0)
# pass constant parameter and parameter function together as a list
parameters <- list(
f_beta = f_beta,
gamma = gamma
)
ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
datareal <- data.frame(
time = times,
I = c(10, 32, 120, 230, 480, 567, 1040, 1743, 2300,
2619, 3542, 4039, 4231, 6378,
5356, 4987, 3421, 2789, 1789, 1156)
)
## define parameter as a vector: gamma and beta
t_beta <- c(0, 12, 16, 19) # consider more or less knots
n_beta <- length(t_beta)
y_beta <- rep(1, n_beta)
p <- c(gamma = 1/5, y_beta) # combine all parameters in one vector
## a small helper function for parameter selection
select <- function(p, which, exclude = FALSE) {
parnames <- names(p)
p[(which == parnames) != exclude]
}
## check the helper function
select(p, "gamma")
select(p, "gamma", excl=TRUE)
## cost function, see ?modCost help page
sir_cost <- function (p) {
gamma <- select(p, "gamma")
y_beta <- select(p, "gamma", exclude = TRUE)
f_beta <- splinefun(x = t_beta, y = y_beta)
outsir <- sir_1(f_beta = f_beta, gamma = gamma,
S0 = 9990, I0 = 10, R0 = 0, times = times)
modCost(model = outsir, obs = datareal)
}
## model calibration, see ?modFit
Fit <- modFit(f = sir_cost, p = p,
# lower bound to avoid negative values of beta
lower = c(gamma = 0, rep(0.0, n_beta)),
# note: high sensitivity wrt. upper bound
upper = c(gamma=1, rep(2.0, n_beta)),
# an algorithm that supports box constraints
method = "Port")
## all parameters were identifiable
summary(Fit)
## smaller time steps to obtain a curves
times <- seq(0, 19, 0.1)
## split components of fitted parameters
gamma <- select(Fit$par, "gamma")
y_beta <- select(Fit$par, "gamma", exclude = TRUE)
out2 <- sir_1(f_beta = splinefun(x = t_beta, y = y_beta), gamma,
S0 = 9990, I0 = 10, R0 = 0, times = times)
## show fitted curves and compare simulation with data
## see ?plot.deSolve help page
plot(out2, obs = datareal, which = c("S", "R", "I", "beta"),
las = 1, obspar = list(pch = 16, col = "red"))
我正在尝试做的事情：我有一个简单的 SIR 模型，具有随时间变化的传输率 beta，我已经在 R 中实现了它（感谢@tpetzoldt）。我们有一个人口N=10000，gamma也是固定的。
sir_1 <- function(f_beta, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- f_beta(time)
gamma <- f_gamma(time)
with(as.list(variables), {
dS <- -beta * I * S/10000
dI <- beta * I * S/10000 - 1/5 * I
dR <- 1/5 * I
return(list(c(dS, dI, dR), beta=beta))
})
}
# time dependent parameter functions
parameters_values <- list(
f_beta = f_beta
)
# the initial values of variables
initial_values <- c(S = S0, I = I0, R = R0)
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
f_beta <- approxfun(x=times, y=seq(0.901, 0.92, by=0.001), rule=2)
out <- as.data.frame(sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times))
现在我有了一些“真实”数据，使用 FME 包我想在每个时间步获得最佳 beta 参数
datareal <- cbind(time = times, I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,
5356, 4987, 3421, 2789, 1789,1156))
sir_cost <- function (f_beta) {
outsir <- as.data.frame(sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times))
costf <- modCost(model = outsir, obs = datareal)
}
p <- rep(0.8, 20)
Fit <- modFit(f = sir_cost, p = p)
Fit
$par
[1] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
我的问题：
- 对于我在每个时间步取 0.8 的初始值，但是 Fit 函数什么都不做，它只returns每个时间步的 0.8（即使我取一个非常高的值，如 800，它表示这已经是最合适的）。我的猜测是同一变量 (beta) 的随时间变化的值，我必须像文档中那样以另一种方式处理这个问题。
非常感谢任何帮助。
我不认为按时间步长估计 beta 是个好主意。这是问题固有的，而不是 deSolve 或 FME 的错误。如果要使用动态模型来估计与时间相关的参数，我会建议使用节点较少的合适函数，例如时间相关的线性、二次或样条曲线，例如 3-5 节而不是 20 节。然后将 approxfun 替换为该函数并将其插入。模型拟合是一门艺术，因此请使用起始值和求解器。还有，读书。
注意以下只是技术演示:
library("deSolve")
library("FME")
sir_1 <- function(f_beta, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- parameters$f_beta(time)
with(as.list(variables), {
dS <- -beta * I * S/10000
dI <- beta * I * S/10000 - 1/5 * I
dR <- 1/5 * I
return(list(c(dS, dI, dR), beta=beta))
})
}
initial_values <- c(S = S0, I = I0, R = R0)
parameters <- list(f_beta=f_beta)
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
# use method "constant" to leave beta constant over time step
f_beta <- approxfun(x=times, y=seq(0.901, 0.92, by=0.001), method="constant", rule=2)
out <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
plot(out)
datareal <- cbind(time = times, I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,
5356, 4987, 3421, 2789, 1789,1156))
plot(out, obs=datareal)
sir_cost <- function (p) {
f_beta <- approxfun(x=times, y=p, method="constant", rule=2)
outsir <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
modCost(model = outsir, obs = datareal)
}
# Play with start values!!!
p <- rep(0.8, 20)
# e.g.: consider random start values
set.seed(123)
p <- runif(20, min=0.8, max=1.2)
# try other solvers, especially such with true box constraints
Fit <- modFit(f = sir_cost, p = p,
lower=rep(0.2, 20), upper=rep(5, 20), # box constraints
method="Port")
summary(Fit) # system is singular (that is what we expected)
# use another solver. Note: it takes a while
Fit <- modFit(f = sir_cost, p = p,
lower=rep(0.2, 20), upper=rep(5, 20), # box constraints
method="L-BFGS-B")
# goes in a surprisingly good direction
Fit$par
f_beta <- approxfun(x=times, y=Fit$par, method="constant", rule=2)
out2 <- sir_1(f_beta=f_beta, S0 = 9990, I0 = 10, R0 = 0, times = times)
# compare with data
plot(out, out2, obs=datareal)
# but see how unstable beta is
plot(out2)
用随时间变化的参数拟合模型可能是个好主意，但如果有理由这样做，我建议限制参数的数量并使用一种平滑函数。
以下示例展示了如何为此目的使用样条曲线，但当然也可以（并且可能更可取）使用具有某种机械意义的函数。
作为副作用，还可以识别伽玛而不是先验地固定它。尽管如此，这仍然是一个技术演示，但我将科学问题悬而未决，即随时间变化的测试版是否有意义。
library("FME")
sir_1 <- function(f_beta, gamma, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- parameters$f_beta(time)
gamma <- parameters$gamma
with(as.list(variables), {
dS <- -beta * I * S / 10000
dI <- beta * I * S / 10000 - gamma * I
dR <- gamma * I
# return vector of derivatives, and beta as auxiliary variable
return(list(c(dS, dI, dR), beta = beta))
})
}
initial_values <- c(S = S0, I = I0, R = R0)
# pass constant parameter and parameter function together as a list
parameters <- list(
f_beta = f_beta,
gamma = gamma
)
ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
datareal <- data.frame(
time = times,
I = c(10, 32, 120, 230, 480, 567, 1040, 1743, 2300,
2619, 3542, 4039, 4231, 6378,
5356, 4987, 3421, 2789, 1789, 1156)
)
## define parameter as a vector: gamma and beta
t_beta <- c(0, 12, 16, 19) # consider more or less knots
n_beta <- length(t_beta)
y_beta <- rep(1, n_beta)
p <- c(gamma = 1/5, y_beta) # combine all parameters in one vector
## a small helper function for parameter selection
select <- function(p, which, exclude = FALSE) {
parnames <- names(p)
p[(which == parnames) != exclude]
}
## check the helper function
select(p, "gamma")
select(p, "gamma", excl=TRUE)
## cost function, see ?modCost help page
sir_cost <- function (p) {
gamma <- select(p, "gamma")
y_beta <- select(p, "gamma", exclude = TRUE)
f_beta <- splinefun(x = t_beta, y = y_beta)
outsir <- sir_1(f_beta = f_beta, gamma = gamma,
S0 = 9990, I0 = 10, R0 = 0, times = times)
modCost(model = outsir, obs = datareal)
}
## model calibration, see ?modFit
Fit <- modFit(f = sir_cost, p = p,
# lower bound to avoid negative values of beta
lower = c(gamma = 0, rep(0.0, n_beta)),
# note: high sensitivity wrt. upper bound
upper = c(gamma=1, rep(2.0, n_beta)),
# an algorithm that supports box constraints
method = "Port")
## all parameters were identifiable
summary(Fit)
## smaller time steps to obtain a curves
times <- seq(0, 19, 0.1)
## split components of fitted parameters
gamma <- select(Fit$par, "gamma")
y_beta <- select(Fit$par, "gamma", exclude = TRUE)
out2 <- sir_1(f_beta = splinefun(x = t_beta, y = y_beta), gamma,
S0 = 9990, I0 = 10, R0 = 0, times = times)
## show fitted curves and compare simulation with data
## see ?plot.deSolve help page
plot(out2, obs = datareal, which = c("S", "R", "I", "beta"),
las = 1, obspar = list(pch = 16, col = "red"))