R 中的 deSolve：返回导数数量错误

Question

我目前正在尝试创建一个函数，使用 deSolve 函数求解隔间模型。根据温度和降水，参数会随着时间的推移而变化，因此我使用 approxfun() 来插值参数估计值。我的虚拟数据集中有 4 个变量，当前有 122 个时间步。当我尝试运行该函数时，出现此错误：

Error in checkFunc(Func2, times, y, rho) :  The number of derivatives returned by func() (488) must equal the length of the initial conditions vector (4)

我猜这与插值参数有关，因为显然 4 个变量 * 122 个时间步长是 488，但我真的很难找出如何解决这个问题。

我将在下面粘贴我的代码！非常感谢任何帮助。

#### Fake data set up ##########################################################

fake.dates <- seq(as.Date("2020/1/1"), as.Date("2020/5/1"), 'days')
# 122 dates 

fake.temp <- sample(-3:30, size = 122, replace = T)
all.precip <- seq(0, 30, by = 0.1)
fake.precip <- sample(all.precip, size = 122, replace = T)

fake.climate <- data.frame(fake.dates, fake.temp, fake.precip)

ost.sim.constanteggs <- function (E0, L0, L3f0, L3p0, start, end, temp, precip) {
  
  require(deSolve)
  
  ## Create a sequence of dates between start and end
  date.range <- seq(as.Date(start), as.Date(end), 'days')
  
  ## Get input time points (numeric)
  global.t <- seq(1, length(date.range))
  
  ## Parameters 
  dev = pmin(1, pmax(0, -0.07258 + (0.00976 * temp)))
  mig.1 = pmax(0, 0.0066 * precip - 0.0382)
  mig.2 = pmax(0, exp(-5.48240 + (0.45392 * temp) - (0.00540 * temp^2)))
  mu.1 = pmin(1, exp(-4.38278 - (0.10640 * temp) + (0.00540 * temp^2)))
  mu.2 = pmin(1, exp(-4.38278 - (0.10640 * temp) + (0.00540 * temp^2)))
  mu.3 = pmin(1, exp((-4.864 * temp) + (0.0048 * temp^2) + (0.00008 * temp^3)))
  mu.4 = pmin(1, exp((-5.743 * temp) + (0.0068 * temp^2) + (0.00003 * temp^3)))
  mu.5 = pmin(1, 10 * exp(-6.388 - (0.2681 * temp) + (0.01633 * temp^2) - (0.00016 * temp^3)))
  
  ## Interpolate the rates
  devrate <- approxfun(dev, method = 'linear', rule = 2)
  mig1rate <- approxfun(mig.1, method = 'linear', rule = 2)
  mig2rate <- approxfun(mig.2, method = 'linear', rule = 2)
  mu1rate <- approxfun(mu.1, method = 'linear', rule = 2)
  mu2rate <- approxfun(mu.2, method = 'linear', rule = 2)
  mu3rate <- approxfun(mu.3, method = 'linear', rule = 2)
  mu4rate <- approxfun(mu.4, method = 'linear', rule = 2)
  mu5rate <- approxfun(mu.5, method = 'linear', rule = 2)
  
  print("Parameters loaded successfully")
  
  # Model function 
  para.dyn <- function(times, para.init, para.par) {
    
    with(as.list(c(para.init, para.par)), {
      
      # Parameters
      dev1 <- devrate(times)
      mig1 <- mig1rate(times)
      mig2 <- mig2rate(times)
      mu1 <- mu1rate(times)
      mu2 <- mu2rate(times)
      mu3 <- mu3rate(times)
      mu4 <- mu4rate(times)
      mu5 <- mu5rate(times)
      
      # Differential equations 
      dE <- -(mu.1 + (2 * dev)) * E + 100
      dL <- -(mu.2 + (2 * dev)) * L + (2 * dev) * E
      dL3f <- -(mu.3 + mig.1) * L3f + (2 * dev) * L
      dL3p <- -mu.4 * (L3p * (1 - mig.2)) - mu.5 * (mig.2 * L3p) + mig.1 * L3f
      
      res <- c(dE, dL, dL3f, dL3p)
      list(res)
    })
  }
  
  print('Model function loaded successfully')
  
  para.init <- c(E = E0, L = L0, L3f = L3f0, L3p = L3p0)
  print('Inital conditions for state variables loaded successfully')
  
  para.output <- lsoda(y = para.init, 
                       times = global.t, 
                       func = para.dyn,
                       parms = NULL)
  
  return(para.output)
}


ss = "2020/1/1"
ee = "2020/5/1"

test <- ost.sim.constanteggs(E0 = 100, L0 = 0, L3f0 = 0, L3p0 = 0, start = ss, end = ee, temp = fake.climate$fake.temp, precip = fake.climate$fake.precip)

我期望得到一个输出，它是每个时间点每个变量（E、L、L3f、L3p）值的数据帧。

Answer 1

方程中的参数名称输入错误，指的是全局向量，而不是插值。更正后的版本应该如下所示：

      # Differential equations
      dE <- -(mu1 + (2 * dev1)) * E + 100
      dL <- -(mu2 + (2 * dev1)) * L + (2 * dev1) * E
      dL3f <- -(mu3 + mig1) * L3f + (2 * dev1) * L
      dL3p <- -mu4 * (L3p * (1 - mig2)) - mu5 * (mig2 * L3p) + mig1 * L3f