r - 使用 deSolve 包中的 ODE 的现实年龄结构化模型
问题描述
我正在尝试使用 deSolve 包中的 ODE 模拟一个现实的年龄结构模型,其中所有个人都可以在时间步长结束时转移到下一个年龄组(而不是以给定的速率连续老化)。
例如,考虑具有易感 (S) 和传染 (I) 两个状态的模型,每个状态分为 4 个年龄组(S1、S2、S3、S4 和 I1、I2、I3、I4),S1 中的所有个体都应该在时间步结束时进入 S2,S2 中的应该进入 S3,依此类推。
我尝试分两步完成,第一步是求解 ODE,第二步是在时间步长结束时将个人转移到下一个年龄组,但没有成功。
以下是我的尝试之一:
library(deSolve)
times <- seq(from = 0, to = 100, by = 1)
n_agecat <- 4
#Initial number of individuals in each state
S_0 = c(999,rep(0,n_agecat-1))
I_0 = c(1,rep(0,n_agecat-1))
si_initial_state_values <- c(S = S_0,
I = I_0)
# Parameter values
si_parameters <- c(beta = 0.01) #contact rate assuming random mixing
si_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
n_agegroups <- 4
S <- state[1:n_agegroups]
I <- state[(n_agegroups+1):(2*n_agegroups)]
# Total population
N <- S+I
# Force of infection
lambda <- beta * I/N
# Solving the differential equations
dS <- -lambda * S
dI <- lambda * S
# Trying to shift all individuals into the following age group
S <- c(0,S[-n_agecat])
I <- c(0,I[-n_agecat])
return(list(c(dS, dI)))
})
}
output <- as.data.frame(ode(y = si_initial_state_values,
times = times,
func = si_model,
parms = si_parameters))
任何指导将不胜感激,在此先感谢您!
解决方案
我看过你的模型。原则上,在事件函数中实现 shift 是可行的,但主要模型仍然存在几个问题:
- 消亡:如果每个时间步移动年龄组并且第一个元素填充为零,则所有内容都在 4 个时间步内移到末尾,并且人口消亡。
- 感染:在您的情况下,感染者只能感染同一年龄组,因此您需要在计算 lambda 之前对“年龄”组进行汇总。
- 最后,什么是“年龄”组?你想知道感染后的时间吗?
总而言之,有几种选择:我个人更喜欢离散模型进行这种模拟,即差分方程、年龄结构化矩阵模型或基于个体的模型。
如果你想让它成为一个 ODE,我建议让易感者一起作为一个状态,并仅将受感染者实现为阶段结构。
这是一个简单的例子,请检查:
library(deSolve)
times <- seq(from = 0, to = 100, by = 1)
n_agegroups <- 14
n_agecat <- 14
# Initial number of individuals in each state
S_0 = c(999) # only one state
I_0 = c(1, rep(0,n_agecat-1)) # several stages
si_initial_state_values <- c(S = S_0,
I = I_0)
# Parameter values
si_parameters <- c(beta = 0.1) # set contact parameter to a higher value
si_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
S <- state[1]
I <- state[2:(n_agegroups + 1)]
# Total population
N <- S + sum(I)
# Force of infection
#lambda <- beta * I/N # old
lambda <- beta * sum(I) / N # NEW
# Solving the differential equations
dS <- -lambda * S
dI <- lambda * S
list(c(dS, c(dI, rep(0, n_agegroups-1))))
})
}
shift <- function(t, state, p) {
S <- state[1]
I <- state[2:(n_agegroups + 1)]
I <- c(0, I[-n_agecat])
c(S, I)
}
# output time steps (note: ode uses automatic simulation steps!)
times <- 1:200
# time step of events (i.e. shifting), not necessarily same as times
evt_times <- 1:200
output <- ode(y = si_initial_state_values,
times = times,
func = si_model,
parms = si_parameters,
events=list(func=shift, time=evt_times))
## default plot function
plot(output, ask=FALSE)
## plot totals
S <- output[,2]
I <- rowSums(output[, -(1:2)])
par(mfrow=c(1,2))
plot(times, S, type="l", ylim=c(0, max(S)))
lines(times, I, col="red", lwd=1)
## plot stage groups
matplot(times, output[, -(1:2)], col=rainbow(n=14), lty=1, type="l", ylab="S")
注意:这只是一个技术演示,不是一个有效的阶段结构化 SIR 模型!
推荐阅读
- c++ - 是否有计划向 std::optional 添加“期望”?
- php - 如何在忽略表情符号的同时匹配非字母数字的正则表达式特殊字符?
- java - 服务器未运行 body multipart/form-data
- sql-server - SQL Server 中的存储过程奇怪错误
- javascript - 如何在 iOS 设备上选择图像期间保持 socket.io 连接处于活动状态
- scala - 如何使用类型参数扩展案例类?
- javascript - 获取div的id
- keycloak - Keycloak - how to implement delegated administration
- c++ - When looping through all processes, the process name stays the same but the process ID changes
- android - 如何使 TextView 在具有更高屏幕比例的设备上响应?