解决deSolve包ode函数报错：‘arg’必须为NULL或字符向量

我正在构建流行病学元种群SEIR模型，用于模拟多互联种群中传染病的传播动态，编写了Metapopulation_SEIR函数实现该模型。运行deSolve包的ode函数求解模型时触发如下错误：

Error in match.arg(method) : 'arg' must be NULL or a character vector

显式指定parms参数后仍出现相同错误，相关代码如下：

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)

尝试如下修改后仍报错：

model <- ode(initial_state, times, Metapopulation_SEIR, parms=params, connectivity_matrix)

错误原因与修复方案 #

• 核心问题：deSolve的ode函数仅支持通过parms参数传递额外模型参数，不能直接传入多个独立参数。你传入的connectivity_matrix被误识别为ode函数的method参数（该参数要求是字符向量或NULL），从而触发类型错误。

• 修复步骤：

  1. 调整模型函数的参数列表，只保留time, current_state, parms三个标准参数。
  2. 将connectivity_matrix整合到parms中，改用列表存储参数（因为向量无法容纳矩阵）。
  3. 在模型函数内部，从parms中提取connectivity_matrix及其他参数。

• 修正后的完整代码：

library(deSolve)
library(ggplot2)

# 修改后的模型函数，仅保留三个标准参数
Metapopulation_SEIR <- function(time, current_state, parms){
  with(as.list(c(current_state, parms)),{
    # 提取连通性矩阵
    connectivity_matrix <- parms$connectivity_matrix
    # 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
num_subpopulations <- 3

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)

# 将参数和连通性矩阵整合为列表
params <- list(
  beta = 0.5, 
  sigma = 0.25, 
  gamma = 0.2, 
  mu = 0.001,
  connectivity_matrix = matrix(0.1, nrow = num_subpopulations, ncol = num_subpopulations)
)
# 设置对角线值
diag(params$connectivity_matrix) <- 1

times <- seq(0, 365, by = 1)

# 正确调用ode函数，仅传递parms参数
model <- ode(initial_state, times, Metapopulation_SEIR, parms = params)

# 可添加绘图代码验证结果
head(model)

内容的提问来源于stack exchange，提问作者user2300042