含列表参数的ODE模型敏感性分析求助：从参数传入报错到结果标准差恒为0的问题排查

Let's break down your issues step by step and fix the code:

1. Initial Error: Invalid Subscript Type 'list' #

The first error happened because you included the c.weight list in params.sens, but sensRange is designed to work with numeric parameters only. Since you don't need to run sensitivity analysis on the weight list, we should separate it from the parameters we want to test.

2. Critical Bug: No Variation in Results (Std Dev = 0) #

The root cause here is that your solve.sens function accepts a pars argument but never uses it—you hardcoded params.sens in the ode call. This means every run of the model used the exact same initial parameters, so all results were identical.

Fixed Code #

library(FME)
library(deSolve)

# Ensure c.weight has enough elements for time steps 0-144 (145 total points)
# Added extra values to avoid index out-of-bounds errors - adjust to match your actual data
c.weight <- c(3.5, 4, 5, 5, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 38, 38, 38, 39, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 52, 54, 55, rep(55, 20))

# Only include numeric parameters for sensitivity analysis
params.sens <- c(vd = 0.2, pm = 0.05, halflife = 5, dose = 0.001)

solve.sens <- function(pars) {
  sens.model <- function(times, state, parameters) {
    with(as.list(c(state, parameters)), {
      # Handle time indexing safely to avoid out-of-bounds errors
      time_idx <- floor(times) + 1
      if (time_idx > length(c.weight)) time_idx <- length(c.weight)
      current_weight <- c.weight[time_idx]
      
      # Calculate volume and transferred amount
      if (times <= 5) {
        volume <- (((0.2 * times) * current_weight) * 30)
        transferred <- pm * volume
      } else if (times > 5 & times <= 14) {
        volume <- ((0.1 * current_weight) * 30)
        transferred <- pm * volume
      } else {
        transferred <- 0
      }
      
      # Calculate rates for ODE
      intake.c <- dose * current_weight
      elimination.c <- concentration * vd * current_weight * log(2) / (halflife * 12)
      # ODE requires derivative of concentration, not direct assignment
      dconcentration <- (intake.c + transferred - elimination.c) / (vd * current_weight)
      
      list(c(dconcentration))
    })
  }
  
  state <- c(concentration = 0.5)
  months <- seq(0, 144, 1)
  # Use the passed pars argument instead of hardcoded params.sens!
  return(as.data.frame(ode(y = state, times = months, func = sens.model, parms = pars)))
}

# Test single run
out <- solve.sens(params.sens)
head(out)

# Define parameter ranges for sensitivity analysis
parRanges <- data.frame(min = c(0.02, 0.1, 2.1), max = c(0.09, 0.2, 5.5))
rownames(parRanges) <- c("pm", "vd", "halflife")

# Run sensitivity analysis
sens <- sensRange(func = solve.sens, parms = params.sens, dist = "latin", sensvar = "concentration", parRange = parRanges, num = 50)

# Check results
head(summary(sens))
summ.sens <- summary(sens)
plot(summ.sens, xlab = "months", ylab = "concentration")

Key Modifications Explained #

• Separated weight list: Removed c.weight from params.sens to avoid list-related errors in sensRange.
• Fixed parameter passing: Changed the ode call to use the pars argument passed to solve.sens—this ensures each sensitivity run uses the varied parameters.
• Safe indexing: Added checks to prevent c.weight index out-of-bounds errors, since your time steps (0-144) require 145 weight values.
• Corrected ODE structure: Modified the model to return the derivative of concentration (dconcentration) instead of directly assigning concentration—this follows the required format for deSolve ODE functions.

After these changes, your sensitivity analysis should produce varied results with non-zero standard deviations. Adjust the c.weight values to match your actual weight data as needed.

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