参数值作为另一个向量的函数。解解答案

【问题标题】：Parameter values as a function of another vector. deSolve参数值作为另一个向量的函数。解解
【发布时间】：2022-01-22 05:14:05
【问题描述】：

我希望建立一个人口动态模型，其中每个参数值对应于当天的温度。例如

简单模型

 library(deSolve)
set.seed(1)

pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)



lv_model <- function(pars, times = seq(0, 50, by = 1)) {
  # initial state 
  state <- c(x = 1, y = 2)
  # derivative
  deriv <- function(t, state, pars) {
    with(as.list(c(state, pars)), {
      d_x <- alpha * x - beta * x * y
      d_y <- delta * beta * x * y - gamma * y
      return(list(c(x = d_x, y = d_y)))
    })
  }
  # solve
  ode(y = state, times = times, func = deriv, parms = pars)
}
lv_results <- lv_model(pars = pars, times = seq(0, 50, by = 1))

我现在想使用一系列每日温度 DailyTemperature<-floor(runif(50,0,40))

并使参数值成为温度的函数

TraitTemperature<-seq(1,40,1)

#trait responses to temperature
alpha<- abs(rnorm(40,mean = 0.5,sd=1))
beta<- abs(rnorm(40,mean = 0.2,sd=0.5))
delta<-abs(rnorm(40,mean=1,sd=2))
gamma<- seq(0.025,1,0.025)
parameters<-as.data.frame(cbind(TraitTemperature,alpha,beta,delta,gamma))

这样对于迭代的每个时间步，它都会查看每日温度，然后在参数数据框中找到相应的温度值。

回顾档案，我看到 if/else 语句用于在特定时间步更改单个参数和使用强制函数，但我认为它们不适用于这里。

我希望这是有道理的，我对如何使它起作用的想法很感兴趣。到目前为止，我还尝试使用 for loop 遍历每日温度列表，然后使用 match 函数来识别值，但这并没有利用每日时间步长。

【问题讨论】：

对deSolve 没有太多经验，但我确实使用迭代方法进行了很多此类动态建模。因此，解决此问题的另一种方法可能是将您的 DE 转换为一种格式，其中 y 在时间 t 的值是时间 t-1 的函数。然后在循环中迭代函数。如果速度是个问题，最好在 Rcpp 中进行此迭代，因为 R 在这种情况下会变得有点慢。
如果我理解正确的话，这就是我们所说的强制。您可以在 deSolve 帮助页面 ?forcings 或例如以下页面中找到更多相关信息：tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html
有几种方法可以做到这一点。一种想法是根据温度为参数创建 4 个信号，但如果信号的索引（例如温度）与时间向量完全对应，也可以通过索引访问来实现（见下文）。另一种方法是使用 simecol 包中的approxTime1，它能够一次返回整个参数值向量。最后，它也可以通过回调来完成，其中parms 是一个进行任意插值的函数。

标签： r for-loop iteration ode desolve

【解决方案1】：

这是一种可能的方法，使用DailyTemperature 作为强制，然后使用parameters 数据框作为查找表。这里不需要match，只要温度是整数，数据框中的温度对应数据框的行索引即可。

我想补充一点，原则上，不连续的强迫会使模型变慢，因为 ODE 根据定义是连续的。幸运的是，由于求解器非常强大，它应该适用于实际应用：

library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {

  pars <- parameters[DailyTemperature[floor(t + 1)], 2:5]
  #print(pars)
  
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y
    list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
  })
}


state <- c(x = 1, y = 2)
times = seq(0, 50, by = 1)

# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
  TraitTemperature = seq(1,40,1),
  alpha = abs(rnorm(40,mean = 0.5,sd=1)),
  beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(40,mean=1,sd=2)),
  gamma = seq(0.025,1,0.025)
)


DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)

参数列表

在上面的示例中，从ode 传递的pars 变量只是被从全局变量parameters 和DailyTemperature 派生的pars 覆盖。这可行，但也可以考虑将两者都作为列表传递给deriv 函数。

deriv <- function(t, state, p) {
  
  parameters <- p[[1]]
  DailyTemperature <- p[[2]]

  parms <- parameters[DailyTemperature[floor(t + 1)], 2:5]
  # ...
}

然后：

out <- ode(y = state, times = times, func = deriv,
  parms = list(parameters, DailyTemperature))

【讨论】：

非常感谢！我将以 0.1 的间隔使用温度值，所以我猜我会在代码的第 4 行使用match 函数？ ``` pars
您还可以使用 approxfun 的插值来获取给定温度的参数行，或者更快地获取直接计算指数的公式。
所以我已经了解了 approxfun 我遇到的问题是使其成为时间和 DailyTemperature 的函数。例如alphainput<- approxfun(alpha,rule = 2) 将遍历 alpha 向量，但不一定针对 DailyTemperature。我对索引也不是很好。有什么想法吗？非常感谢您的帮助

【解决方案2】：

这似乎适用于使用当前可重现代码进行索引：

set.seed(1)
deriv <- function(t, state, pars) {
  pars<- parameters[match(parameters$TraitTemperature[parameters[2:5]],DailyTemperature),]
  #print(pars)
  
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y
    list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
  })
}

state <- c(x = 1, y = 2)
times = seq(0, 50, by = 1)

# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
  TraitTemperature = seq(1,40,1),
  alpha = abs(rnorm(40,mean = 0.5,sd=1)),
  beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(40,mean=1,sd=2)),
  gamma = seq(0.025,1,0.025)
)

DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)

但是如果我增加值的分辨率，例如


# Parameter datasets 
parameters <- data.frame(
  TraitTemperature = seq(0.1,40,0.1),
  alpha = abs(rnorm(400,mean = 0.5,sd=1)),
  beta = abs(rnorm(400,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(400,mean=1,sd=2)),
  gamma = seq(0.0025,1,0.0025)
)

# random daily temperature dataset 
DailyTemperature <- round(runif(51, 0, 40),1) # one more because start zero

然后我在某些时间步得到 NA

【讨论】：

这里最好再问一个问题。尽管如此，我仍然可以为这个问题写第二个答案，但为了避免让稍后可能阅读此内容的 SO 用户感到困惑，我建议您提取具体问题，然后再问一次。然后，我们可以将 Q+A 移至新线程，从此处删除过时的代码。

【解决方案3】：

OP 扩展了她/他的问题，因此它可能是启动新线程的首选方式，但为了提供快速反馈，我尝试给出另一个答案。

在二维（时间和温度）中存在多种插值或索引方法。我的首选方法是创建 2D 模型，然后使用 2D 插值方法。如果参数表面是平滑的，而不是像给定示例中那样随机，则此方法效果最佳。但是，为了简单起见，也可以使用舍入和表格查找。如果值不是整数，由于舍入效应和有限的精度（IEEE 数字格式），精确比较通常不起作用，因此可以使用which.min 代替匹配：

DailyTemperature <- round(runif(51, 0, 40), 1)

TraitTemperature <- seq(0, 40, by=0.1)
N <- length(TraitTemperature)

parameters <- data.frame(
  TraitTemperature = TraitTemperature,
  alpha = abs(rnorm(N, mean = 0.5, sd=1)),
  beta = abs(rnorm(N, mean = 0.2, sd=0.5)),
  delta = abs(rnorm(N, mean=1,sd=2)),
  gamma = seq(0.025, 1, length.out=N)
)


t <- 17
actualTemp <- DailyTemperature[floor(t+1)]
actualTemp
pars <- parameters[which.min(abs(actualTemp - parameters$TraitTemperature)), 1:5]

head(pars)

【讨论】：