我正在尝试创建一个流行病学 SEIR 模型并随着时间的推移模拟模型隔间。更具体地说，下面的代码定义了一个名为 Metapopulation_SEIR 的函数。该函数代表了一种称为 Metapopulation SEIR 模型的区室模型，该模型常用于流行病学中，用于模拟传染病在多个相互关联的人群中的传播。参见代码：
library(deSolve)
library(ggplot2)
# Metapopulation SEIR model function
Metapopulation_SEIR <- function(time, current_state, params, connectivity_matrix){
with(as.list(c(current_state, params)),{
# Calculate total population size in each subpopulation
N <- apply(current_state, 1, sum)
# Initialize rate of change vectors
dS <- rep(0, nrow(current_state))
dE <- rep(0, nrow(current_state))
dI <- rep(0, nrow(current_state))
dR <- rep(0, nrow(current_state))
dM <- rep(0, nrow(current_state))
# Loop through each subpopulation
for (i in 1:nrow(current_state)) {
# Calculate total number of individuals in subpopulation i
Ni <- N[i]
# Calculate rates of change for each compartment in subpopulation i
dSi <- -beta * S[i] * I[i] / Ni
dEi <- (beta * S[i] * I[i] / Ni) - sigma * E[i]
dIi <- sigma * E[i] - gamma * I[i] - mu * I[i]
dRi <- gamma * I[i]
dMi <- mu * I[i]
# Update rate of change vectors
dS[i] <- dSi
dE[i] <- dEi
dI[i] <- dIi
dR[i] <- dRi
dM[i] <- dMi
# Calculate movement of individuals between subpopulations
for (j in 1:nrow(current_state)) {
dS[i] <- dS[i] + connectivity_matrix[i, j] * (beta * S[j] * I[j] / Ni)
dE[i] <- dE[i] + connectivity_matrix[i, j] * (sigma * E[j])
dI[i] <- dI[i] + connectivity_matrix[i, j] * (gamma * I[j])
dR[i] <- dR[i] + connectivity_matrix[i, j] * (gamma * R[j])
dM[i] <- dM[i] + connectivity_matrix[i, j] * (mu * I[j])
}
}
# Return a list of rates of change for each compartment
return(list(c(dS, dE, dI, dR, dM)))
})
}
# Example inputs
# Number of subpopulations
num_subpopulations <- 3
# Initial compartment counts for each subpopulation
initial_state <- matrix(c(
S1 = 900, E1 = 10, I1 = 5, R1 = 85,
S2 = 950, E2 = 5, I2 = 2, R2 = 43,
S3 = 850, E3 = 15, I3 = 7, R3 = 100
), ncol = 4, byrow = TRUE)
# Parameters
params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)
# Connectivity matrix (specifies the movement of individuals between subpopulations)
# Example: Fully connected network where individuals can move between all subpopulations
connectivity_matrix <- matrix(0.1, nrow = num_subpopulations, ncol = num_subpopulations)
diag(connectivity_matrix) <- 1 # Individuals stay within their own subpopulation
# Time points
times <- seq(0, 365, by = 1)
# Solve the model
model <- ode(initial_state, times, Metapopulation_SEIR, params, connectivity_matrix)
最后一行出现以下错误：
Error in match.arg(method) : 'arg' must be NULL or a character vector
我尝试将论证命名如下：
model <- ode(initial_state, times, Metapopulation_SEIR, parms=params, connectivity_matrix)
但仍然遇到同样的错误。 感谢您的帮助！
代码包含多个错误：
odeinitial_statewith进一步的问题是，没有必要使用循环。 R是一种矢量化语言，因此状态可以直接以矩阵形式表示，这使得代码更紧凑，速度更快。
这里是代码的中间版本，解决了错误，但我也强烈建议以矩阵形式重新制定模型。
library(deSolve)
# Metapopulation SEIR model function
Metapopulation_SEIR <- function(time, current_state, params, connectivity_matrix, n_sub){
# tp: access to names in state **matrix** not possible with "with",
# this works only for parms
#with(as.list(c(current_state, params)),{
with(as.list(params),{
# tp: reformat state vector as matrix
dim(current_state) <- c(n_sub, 5)
# tp: extract states
S <- current_state[,1]
E <- current_state[,2]
I <- current_state[,3]
R <- current_state[,4]
M <- current_state[,5]
# Calculate total population size in each subpopulation
N <- apply(current_state, 1, sum)
# Initialize rate of change vectors
dS <- rep(0, nrow(current_state))
dE <- rep(0, nrow(current_state))
dI <- rep(0, nrow(current_state))
dR <- rep(0, nrow(current_state))
dM <- rep(0, nrow(current_state))
# tp: loop is not necessary, R is a vectorized language
# Loop through each subpopulation
for (i in 1:nrow(current_state)) {
# Calculate total number of individuals in subpopulation i
Ni <- N[i]
# Calculate rates of change for each compartment in subpopulation i
dSi <- -beta * S[i] * I[i] / Ni
dEi <- (beta * S[i] * I[i] / Ni) - sigma * E[i]
dIi <- sigma * E[i] - gamma * I[i] - mu * I[i]
dRi <- gamma * I[i]
dMi <- mu * I[i]
# Update rate of change vectors
dS[i] <- dSi
dE[i] <- dEi
dI[i] <- dIi
dR[i] <- dRi
dM[i] <- dMi
# Calculate movement of individuals between subpopulations
for (j in 1:nrow(current_state)) {
dS[i] <- dS[i] + connectivity_matrix[i, j] * (beta * S[j] * I[j] / Ni)
dE[i] <- dE[i] + connectivity_matrix[i, j] * (sigma * E[j])
dI[i] <- dI[i] + connectivity_matrix[i, j] * (gamma * I[j])
dR[i] <- dR[i] + connectivity_matrix[i, j] * (gamma * R[j])
dM[i] <- dM[i] + connectivity_matrix[i, j] * (mu * I[j])
}
}
# Return a list of rates of change for each compartment
return(list(c(dS, dE, dI, dR, dM)))
})
}
# Example inputs
# Number of subpopulations
num_subpopulations <- 3
# Initial compartment counts for each subpopulation
# tp: add 5th column for M states
initial_state <- matrix(c(
S1 = 900, E1 = 10, I1 = 5, R1 = 85, M1=0,
S2 = 950, E2 = 5, I2 = 2, R2 = 43, M2=0,
S3 = 850, E3 = 15, I3 = 7, R3 = 100, M3=0
), ncol = 5, byrow = TRUE)
# tp: reformat initial_state into vector
initial_state <- as.vector(initial_state)
# Parameters
params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)
# Connectivity matrix (specifies the movement of individuals between subpopulations)
# Example: Fully connected network where individuals can move between all subpopulations
connectivity_matrix <- matrix(0.1, nrow = num_subpopulations, ncol = num_subpopulations)
diag(connectivity_matrix) <- 1 # Individuals stay within their own subpopulation
# Time points
times <- seq(0, 365, by = 1)
# Solve the model
# tp: add n_sub argument
model <- ode(initial_state, times, Metapopulation_SEIR, params, "lsoda",
connectivity_matrix = connectivity_matrix, n_sub=num_subpopulations)
# tp: alternative: works too, but extra arguments must always be named
model <- ode(initial_state, times, Metapopulation_SEIR, params,
connectivity_matrix = connectivity_matrix, n_sub=num_subpopulations)