2

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

简单模型

 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想要在特定时间步更改单个参数和使用强制函数时使用的语句，但我认为它们不适用于这里。

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

4

3 回答 3

1

这是一种可能的方法，使用DailyTemperatureas 强制，然后使用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)

温度相关参数

参数作为列表

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

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))
于 2021-12-20T16:31:54.937 回答
0

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

有几种方法可用于二维（时间和温度）中的插值或索引。我的首选方法是创建 2D 模型，然后使用 2D 插值方法。如果参数表面是光滑的而不是像给定示例中那样随机，则此方法效果最佳。但是，为了简单起见，也可以使用舍入和表格查找。如果值不是整数，则由于舍入效应和有限的精度（IEEE 数字格式），精确比较通常不起作用，因此可以使用 match 代替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)
于 2021-12-21T13:11:23.130 回答
0

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

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

于 2021-12-21T12:19:16.793 回答