R - deSolve 包(ode 函数):根据时间改变 SIR 模型中的参数矩阵 [英] R - deSolve package (ode function): change a matrix of parameters in SIR model according to time
问题描述
我正在尝试使用 deSolve 包中的函数 ode 来模拟病毒在人群中的传播.我的模型的基础是一个 SIR 模型，我在这里发布了一个更简单的模型演示，它只包含三个状态S(敏感)、I(感染)和 R(恢复).每个州在我的代码中由 m*n 矩阵表示，因为我的人口中有m个年龄组和n个亚群.>
问题是:在模拟期间，会有几次疫苗接种活动将处于状态S的人转移到状态I.每个疫苗接种活动都以开始日期、结束日期、其覆盖率和持续时间为特征.我想要做的是，一旦 时间 t 落入一项疫苗接种活动的开始日期 和结束日期 的区间，代码会计算有效疫苗接种率(也是一个m*n矩阵，基于覆盖率和持续时间)并乘以S(m*n 矩阵)，得到一个由人过渡到状态 I 的矩阵.现在，我正在使用 if() 来确定 时间 t 是否在 开始日期 和 结束日期 之间>:
#初始化有效接种率矩阵irrate_matrix = 矩阵(数据 = 代表(0，m*n)，nrow = m，ncol = n)for (i in 1:length(tbegin)){如果 (t>=tbegin[i] & t<=tend[i]){for (j in 1:n){irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]}}}这里，irrate_matrix是m*n有效接种率矩阵，m = 2是年龄组数，n = 2 是子种群的数量，tbegin = c(5, 20, 35) 是开始3 次疫苗接种活动的日期，tend = c(8, 23, 38) 是 3 次疫苗接种活动的结束日期，covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) 是每个亚群的每种疫苗接种的覆盖率(例如，covir[1] = 0.35 是亚群1第一次接种疫苗的覆盖率，而covir[4] = 0.225是覆盖率subpopulation2) 的第一次接种，duration = c(4, 4, 4) 是每次接种的持续时间(以天为单位).
在计算irrate_matrix后，我将其带入导数，因此我有:
dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)dR = as.matrix(gammaS*I) - as.matrix(mu*R)我想从第 0 天到第 50 天以 1 天为步长进行模拟，因此:
times = seq(0, 50, 1)我的代码的当前问题是:每次 时间 t 接近 tbegin[i] 或tend[i]，模拟变得慢得多，因为它在这个时间点比在任何其他时间点迭代更多的轮次.例如，一旦 时间 t 达到 tbegin[1] = 5，模型会在时间点 5 迭代多轮.我附上了打印这些迭代的屏幕截图(screenshot1 和 screenshot2).我发现这就是为什么我的更大模型现在需要很长时间运行的原因.
我曾尝试使用事件"tpetzoldt 在这个问题中提到的 deSolve 函数 stackoverflow:将参数的值更改为时间的函数.然而，我发现改变一个参数矩阵并在每次有疫苗接种活动时改变它对我来说很不方便.
我正在寻找关于以下方面的解决方案:
如何在有疫苗接种活动时将我的
irrate_matrix更改为非零矩阵，并在没有疫苗接种时将其设为零矩阵?(每次接种都要计算)同时，如何避免在任何
tbegin[i]或tend[i]处迭代多轮，从而使代码运行得更快?(我认为我不应该使用if()但我不知道我应该如何处理我的情况)如果我需要使用强制"或事件"功能，你能否也告诉我如何有多个强制"/事件"?在模型中?现在，我有一个事件"在我更大的模型中用于将病毒引入人群，如:
virusevents = data.frame(var = "I1", time = 2, value = 1, method = "add")欢迎任何好主意，非常感谢直接提供一些代码！提前致谢！
作为参考，我在这里发布了整个演示:
library(deSolve)####################################(1) 定义sir 函数####################################sir_basic <- 函数 (t, x, params){ # 检索初始状态S = 矩阵(数据 = x[(0*m*n+1):(1*m*n)], nrow = m, ncol = n)I = 矩阵(数据 = x[(1*m*n+1):(2*m*n)], nrow = m, ncol = n)R = 矩阵(数据 = x[(2*m*n+1):(3*m*n)], nrow = m, ncol = n)with(as.list(params), {N = as.matrix(S + I + R)# 打印当前迭代打印(粘贴0(时间总人口"，t，是"，总和(N)))# 通过检查时间 t 计算 irrate_matrixirrate_matrix = 矩阵(数据 = 代表(0，m*n)，nrow = m，ncol = n)for (i in 1:length(tbegin)){如果 (t>=tbegin[i] & t<=tend[i]){for (j in 1:n){irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]}}}#衍生品dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)dR = as.matrix(gammaS*I) - as.matrix(mu*R)衍生物＜- c(dS, dI, dR)列表(衍生物)})}###################################(2) 表征参数####################################m = 2 # 年龄组数n = 2 # 子种群数tbegin = c(5, 20, 35) # 开始日期趋势 = c(8, 23, 38) # 结束日期持续时间 = c(4, 4, 4) # 持续时间covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) # 覆盖率b = 0.0006 # 每日出生率mu = 0.0006 # 每日死亡率gammaS = 0.05 # 从 I 到 R 的转换率参数 = c(m = m, n = n,tbegin = tbegin，tend = 趋向，持续时间 = 持续时间，covir = covir，b = b, mu = mu, gammaS = gammaS)#######################################(3) 初始状态########################################初始化 = c(S = c(20000, 40000, 10000, 20000),I = rep(0, m*n),R = rep(0, m*n))#######################################(4) 运行模拟########################################时间 = seq(0, 50, 1)traj <- ode(func = sir_basic,y = 初始化，参数 = 参数，次 = 次)情节(轨迹)元素明智的操作对于矩阵和向量是相同的，所以 as.matrix 转换是多余的，因为没有 true 使用矩阵乘法.与 rep 相同:无论如何都会回收零.
实际上，CPU 时间已经减少到 50%.相比之下，使用带有 approxTime 的外部强制代替内部的 if 和 for 会使模型变慢(未显示).
简化代码
sir_basic2 <- 函数 (t, x, params){ # 检索初始状态S = x[(0*m*n+1):(1*m*n)]我 = x[(1*m*n+1):(2*m*n)]R = x[(2*m*n+1):(3*m*n)]with(as.list(params), {N = S + I + R# 打印当前迭代#print(paste0(时间总人口", t, 是", sum(N)))# 通过检查时间 t 计算 irrate_matrixirrate_matrix = 矩阵(数据 = 0，nrow = m，ncol = n)for (i in 1:length(tbegin)){如果 (t >= tbegin[i] & t <= 趋向 [i]){for (j in 1:n){irrate_matrix[, j] = -log(1-covir[(j-1) * length(tbegin)+i])/duration[i]}}}#衍生品dS = b*N - irrate_matrix*S - mu*SdI = irrate_matrix*S - gammaS*I - mu*IdR = gammaS*I - mu*R列表(c(dS，dI，dR))})}基准
每个模型版本运行 10 次.模型 sir_basic 是原始实现，其中 print 行被禁用以进行公平比较.
system.time(对于(我在 1:10)traj <- ode(func = sir_basic,y = 初始化，参数 = 参数，次 = 次))系统时间(对于(我在 1:10)traj2 <- ode(func = sir_basic2,y = 初始化，参数 = 参数，次 = 次))情节(traj，traj2)摘要(traj - traj2)当我使用 method="adams" 而不是默认的 lsoda 求解器时，我观察到了另一个显着的加速，但这对于您的完整模型可能会有所不同.>
I am trying to simulate the transmission of viruses in a population using the function ode from the deSolve package. The basic of my model is a SIR model and I posted a much simpler demo of my model here, which consists of only three states S(susceptible), I(infectious) and R(recovered). Each state is represented by a m*n matrix in my code, since I have m age groups and n subpopulations in my population.
The problem is: during the simulation period, there will be several vaccination activities that transfer people in state S to state I. Each vaccination activity is characterized by a begin date, an end date, its coverage rate and duration. What I want to do is once the time t falls into the interval of begin date and end date of one vaccination activity, the code calculates the effective vaccination rate (also a m*n matrix, based on coverage rate and duration) and times it with S (m*n matrix), to get a matrix of people transited to state I. Right now, I am using if() to decide if time t is between a begin date and a end date:
#initialize the matrix of effective vaccination rate
irrate_matrix = matrix(data = rep(0, m*n), nrow = m, ncol = n)
for (i in 1:length(tbegin)){
if (t>=tbegin[i] & t<=tend[i]){
for (j in 1:n){
irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]
}
}
}
Here, irrate_matrix is the m*n effective vaccination rate matrix, m = 2 is the number of age groups, n = 2 is the number of subpopulations, tbegin = c(5, 20, 35) is the begin date of 3 vaccination activities, tend = c(8, 23, 38) is the end date of 3 vaccination activities, covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) is the coverage rate of each vaccination for each subpopulation (e.g., covir[1] = 0.35 is the coverage rate of the first vaccination for subpopulation1, while covir[4] = 0.225 is the coverage rate of the first vaccination for subpopulation2) and duration = c(4, 4, 4) is the duration of each vaccination (in days).
After calculating irrate_matrix, I take it into derivatives and therefore I have:
dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)
dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)
dR = as.matrix(gammaS*I) - as.matrix(mu*R)
I want to do a simulation from day 0 to day 50, by 1-day step, thus:
times = seq(0, 50, 1)
The current issue with my code is: every time the time t comes to a time point close to a tbegin[i] or tend[i], the simulation becomes much slower since it iterates at this time point for much more rounds than at any other time point. For example, once the time t comes to tbegin[1] = 5, the model iterates at time point 5 for many rounds. I attached screenshots from printing out those iterations (screenshot1 and screenshot2). I find this is why my bigger model takes a very long running time now.
I have tried using the "events" function of deSolve mentioned by tpetzoldt in this question stackoverflow: change the value of a parameter as a function of time. However, I found it's inconvenient for me to change a matrix of parameters and change it every time there is a vaccination activity.
I am looking for solutions regarding:
How to change my
irrate_matrixto non-zero matrix when there is a vaccination activity and let it be zero matrix when there is no vaccination? (it has to be calculated for each vaccination)At the same time, how to make the code run faster by avoiding iterating at any
tbegin[i]ortend[i]for many rounds? (I think I should not useif()but I do not know what I should do with my case)If I need to use "forcing" or "events" function, could you please also tell me how to have multiple "forcing"/"events" in the model? Right now, I have had an "events" used in my bigger model to introduce a virus to the population, as:
virusevents = data.frame(var = "I1", time = 2, value = 1, method = "add")
Any good idea is welcome and directly providing some codes is much appreciated! Thank you in advance!
For reference, I post the whole demo here:
library(deSolve)
##################################
###(1) define the sir function####
##################################
sir_basic <- function (t, x, params)
{ # retrieve initial states
S = matrix(data = x[(0*m*n+1):(1*m*n)], nrow = m, ncol = n)
I = matrix(data = x[(1*m*n+1):(2*m*n)], nrow = m, ncol = n)
R = matrix(data = x[(2*m*n+1):(3*m*n)], nrow = m, ncol = n)
with(as.list(params), {
N = as.matrix(S + I + R)
# print out current iteration
print(paste0("Total population at time ", t, " is ", sum(N)))
# calculate irrate_matrix by checking time t
irrate_matrix = matrix(data = rep(0, m*n), nrow = m, ncol = n)
for (i in 1:length(tbegin)){
if (t>=tbegin[i] & t<=tend[i]){
for (j in 1:n){
irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]
}
}
}
# derivatives
dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)
dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)
dR = as.matrix(gammaS*I) - as.matrix(mu*R)
derivatives <- c(dS, dI, dR)
list(derivatives)
})
}
##################################
###(2) characterize parameters####
##################################
m = 2 # the number of age groups
n = 2 # the number of sub-populations
tbegin = c(5, 20, 35) # begin dates
tend = c(8, 23, 38) # end dates
duration = c(4, 4, 4) # duration
covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) # coverage rates
b = 0.0006 # daily birth rate
mu = 0.0006 # daily death rate
gammaS = 0.05 # transition rate from I to R
parameters = c(m = m, n = n,
tbegin = tbegin, tend = tend, duration = duration, covir = covir,
b = b, mu = mu, gammaS = gammaS)
##################################
#######(3) initial states ########
##################################
inits = c(
S = c(20000, 40000, 10000, 20000),
I = rep(0, m*n),
R = rep(0, m*n)
)
##################################
#######(4) run simulations########
##################################
times = seq(0, 50, 1)
traj <- ode(func = sir_basic,
y = inits,
parms = parameters,
times = times)
plot(traj)
Element wise operations are the same for matrices and vectors, so the as.matrix conversions are redundant, as no true matrix multiplication is used. Same with the rep: the zero is recycled anyway.
In effect, CPU time reduces already to 50%. In contrast, use of an external forcing with approxTime instead of the inner if and for made the model slower (not shown).
Simplified code
sir_basic2 <- function (t, x, params)
{ # retrieve initial states
S = x[(0*m*n+1):(1*m*n)]
I = x[(1*m*n+1):(2*m*n)]
R = x[(2*m*n+1):(3*m*n)]
with(as.list(params), {
N = S + I + R
# print out current iteration
#print(paste0("Total population at time ", t, " is ", sum(N)))
# calculate irrate_matrix by checking time t
irrate_matrix = matrix(data = 0, nrow = m, ncol = n)
for (i in 1:length(tbegin)){
if (t >= tbegin[i] & t <= tend[i]){
for (j in 1:n){
irrate_matrix[, j] = -log(1-covir[(j-1) * length(tbegin)+i])/duration[i]
}
}
}
# derivatives
dS = b*N - irrate_matrix*S - mu*S
dI = irrate_matrix*S - gammaS*I - mu*I
dR = gammaS*I - mu*R
list(c(dS, dI, dR))
})
}
Benchmark
Each model version is run 10 times. Model sir_basic is the original implementation, where print line was disabled for a fair comparison.
system.time(
for(i in 1:10)
traj <- ode(func = sir_basic,
y = inits,
parms = parameters,
times = times)
)
system.time(
for(i in 1:10)
traj2 <- ode(func = sir_basic2,
y = inits,
parms = parameters,
times = times)
)
plot(traj, traj2)
summary(traj - traj2)
I observed another considerable speedup, when I use method="adams" instead of the default lsoda solver, but this may differ for your full model.
这篇关于R - deSolve 包(ode 函数):根据时间改变 SIR 模型中的参数矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
