如何求解arg必须为NULL或deSolve包的ode函数中的字符向量错误

代码包含几个错误:

• ode呼叫中的额外名字必须被命名为
• initial_state必须作为矩阵传递
• 一个矩阵的名字不能用with
• 初始状态的数量必须等于导数的数量

另一个问题是，没有必要使用循环.R是一种向量化语言，因此状态可以直接用矩阵形式表达，这使得代码更紧凑，速度更快.

下面是代码的中间版本，只解决错误，但仍然包含循环.

然而，我也强烈建议以矩阵形式重新拟订该模式.可以在这里找到一个例子:https://github.com/tpetzoldt/covid/blob/master/models/sir_metapopulation.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 vectors like '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)