问题描述
我正在尝试使用 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,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,n_agegroups-1))))
})
}
shift <- function(t,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,ylab="S")
注意:这只是一个技术演示,不是一个有效的阶段结构SIR模型!