以列表作为参数的 ODE 的灵敏度分析 - 结果给出标准偏差为 0

Question

注意：最初的问题是“使用包含列表的参数对 ODE 进行灵敏度分析”，因为 sensRange 函数由于参数中传递的列表而出错。问题随着列表参数的固定而演变，但出现了一个不同的问题，其中灵敏度分析显示奇怪的结果，所有运行的标准偏差为 0。

我有一个模型，使用 deSolve 包在 R 中模拟化学物质随时间的浓度。我使用的参数是化学特性（分布容积、半衰期等）以及随时间变化的重量。权重在列表中给出，由于权重随时间增加，该列表迭代 ODE 中的每个时间步长，并且化学性质以数字形式给出。

我想对这个模型进行敏感性分析，只测试化学性质。我之前没有做过敏感性分析，但我正在尝试使用 sensRange() 来遵循示例。但是，似乎 sensRange() 函数不允许将其中一个参数作为列表给出。我收到错误：

Error in yRef[, ivar] : invalid subscript type 'list'

我的模型和全局敏感性分析代码是这样设置的：

library(FME)
library(deSolve)

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)

params.sens <- c(vd = 0.2,
                 pm = 0.05,
                 halflife = 5,
                 dose = 0.001,
                 c.weight = c.weight)

solve.sens <- function(pars) {
  sens.model <- function(times, state, parameters) {
    with(as.list(c(state, parameters)), {
      
      if (times <= 5) {
        volume <- (((0.2 * times) * c.weight[times+1]) * 30)
        transferred <- pm * volume
      } else if (times > 5 & times <= 14) {
        volume <- ((0.1 * c.weight[times+1]) * 30)
        transferred <- pm * volume
      } else {
        transferred <- 0
      }
      
      intake.c <- (dose * c.weight[times+1])
      elimination.c <- concentration * vd * c.weight[times+1] * log(2) / (halflife * 12)
      
      #total
      concentration <- (intake.c + transferred - elimination.c) / (vd * c.weight[times+1])
      
      list(c(concentration))
      
    })
  }
  
  state <- c(concentration = 0.5)
  months <- seq(0, 144, 1)
  return(as.data.frame(ode(y = state, times = months, func = sens.model, parms = params.sens)))
  
}

out <- solve.sens(params.sens)

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")

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

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

我不知道如何前进，有没有人有任何提示或看到我的错误在哪里？？

编辑： 遵循 Soetart 和 Herman，2009 年的细菌生长模型，纠正我的 =/<- 错误，并将评论中的参数值添加到模型中。现在它运行没有错误，但是摘要显示了所有相同的值（平均值、最小值、最大值和所有分位数）所以我假设它不是运行正确

               x      Mean Sd       Min       Max       q05       q25       q50       q75       q95
concentration0 0  0.500000  0  0.500000  0.500000  0.500000  0.500000  0.500000  0.500000  0.500000
concentration1 1  1.246348  0  1.246348  1.246348  1.246348  1.246348  1.246348  1.246348  1.246348
concentration2 2  3.475493  0  3.475493  3.475493  3.475493  3.475493  3.475493  3.475493  3.475493
concentration3 3  7.170403  0  7.170403  7.170403  7.170403  7.170403  7.170403  7.170403  7.170403
concentration4 4 12.314242  0 12.314242 12.314242 12.314242 12.314242 12.314242 12.314242 12.314242
concentration5 5 18.890365  0 18.890365 18.890365 18.890365 18.890365 18.890365 18.890365 18.890365

Answer 1

原来post混淆了=和<-来创建命名vector，所以推荐如下代码片段：

params.sens <- c(vd = p.vd,
              pm = p.pm,
              halflife = p.hl,
              dose = concentrations$dose,
              c.weight = variables$c.weight)

经过原 poster 的编辑，此答案大部分已过时并被删除。然而，在中发现 <-（赋值）和 =（参数匹配；此处创建命名向量）之间的区别尚未完全清楚。

这是一个显示差异的示例。为避免混淆，运行它在新的 R 会话中或删除工作 space 之前：

#rm(list=ls()) # uncomment this to clear the work space
x <- c(a <- 2, b <- 3)
y <- c(d = 2, e = 3)

x
y
ls()

我们看到 y 是一个命名向量，而 x 没有命名元素。相反，a 和 b 是在用户工作 space:

中作为变量“即时”创建的

> x
[1] 2 3
> y
d e 
2 3 
> ls()
[1] "a" "b" "x" "y"

Answer 2

由于修复了分配错误，脚本现在可以正常运行并且可以重现行为。

现在我们可以看到传递给 solve.sens 的参数没有传递给 ode 函数。

解决方法是将 ode 调用中的 parms.sens 替换为 pars，并传递给调用函数：

return(as.data.frame(ode(y = state, times = months, func = sens.model, parms = pars)))

然后

plot(summ.sens, xlab = "months", ylab = "concentration")

结果：

以列表作为参数的 ODE 的灵敏度分析 - 结果给出标准偏差为 0

Sensitivity analysis for ODE with list as parameters - result gives standard deviation of 0

r

odesensitivity

desolve