我正在尝试使用 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))
任何指导将不胜感激,在此先感谢您!