在 R：FME/ deSolve - SIR 拟合时变参数

如何解决在 R：FME/ deSolve - SIR 拟合时变参数

我想要做的：我有一个简单的 SIR 模型，具有随时间变化的传输速率 beta，我已经在 R 中实现了它（感谢 @tpetzoldt）。我们有 N=10000 的人口，gamma 也是固定的。

  sir_1 <- function(f_beta,S0,I0,R0,times) {

  # the differential equations
  sir_equations <- function(time,variables,parameters) {
  beta  <- f_beta(time)
  gamma <- f_gamma(time)
  with(as.list(variables),{
  dS <- -beta * I * S/10000
  dI <-  beta * I * S/10000 - 1/5 * I
  dR <-  1/5 * I
  return(list(c(dS,dI,dR),beta=beta))
  })
  }

# time dependent parameter functions
parameters_values <- list(
f_beta  = f_beta
)

# the initial values of variables
initial_values <- c(S = S0,I = I0,R = R0)



out <- ode(initial_values,times,sir_equations,parameters)
}

times <- seq(0,19)


f_beta <- approxfun(x=times,y=seq(0.901,0.92,by=0.001),rule=2)


out <- as.data.frame(sir_1(f_beta=f_beta,S0 = 9990,I0 = 10,R0 = 0,times = times))

现在我有了一些“真实”数据，我想通过 FME 包在每个时间步获得最佳的 beta 参数

 datareal <- cbind(time = times,I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,5356,4987,3421,2789,1789,1156))



 sir_cost <- function (f_beta) {
 outsir <- as.data.frame(sir_1(f_beta=f_beta,times = times))
 costf <- modCost(model = outsir,obs = datareal)
 }


p <- rep(0.8,20)
Fit <- modFit(f = sir_cost,p = p)

Fit
$par
[1] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

我的问题：

1. 对于我在每个时间步长取 0.8 的初始值，但是 Fit 函数什么都不做，它只为每个时间步长返回 0.8（即使我取了一个非常高的值，比如 800，它说这已经是最合适的了）。我的猜测是对于同一个变量 (beta) 的时变值，我必须像文档中那样以另一种方式处理这个问题。

非常感谢任何帮助。

解决方法

我认为按时间步长估算 beta 不是一个好主意。这是问题所固有的，而不是 deSolve 或 FME 的错。如果使用动态模型来估计时间相关参数，我会建议使用具有较少结点的合适函数，例如时间相关的线性、二次或样条曲线，例如 3-5 节而不是 20 节。然后用该函数替换 approxfun 并将其插入。模型拟合是一门艺术，因此请使用起始值和求解器。还有，看书。

请注意，以下只是一个技术演示：

library("deSolve")
library("FME")

sir_1 <- function(f_beta,S0,I0,R0,times) {
  # the differential equations
  sir_equations <- function(time,variables,parameters) {
    beta  <- parameters$f_beta(time)
    with(as.list(variables),{
      dS <- -beta * I * S/10000
      dI <-  beta * I * S/10000 - 1/5 * I
      dR <-  1/5 * I
      return(list(c(dS,dI,dR),beta=beta))
    })
  }

  initial_values <- c(S = S0,I = I0,R = R0)
  parameters <- list(f_beta=f_beta)
  out <- ode(initial_values,times,sir_equations,parameters)
}

times <- seq(0,19)
# use method "constant" to leave beta constant over time step
f_beta <- approxfun(x=times,y=seq(0.901,0.92,by=0.001),method="constant",rule=2)
out <- sir_1(f_beta=f_beta,S0 = 9990,I0 = 10,R0 = 0,times = times)
plot(out)

datareal <- cbind(time = times,I=c(10,32,120,230,480,567,1040,1743,2300,2619,3542,4039,4231,6378,5356,4987,3421,2789,1789,1156))
plot(out,obs=datareal)

sir_cost <- function (p) {
  f_beta <- approxfun(x=times,y=p,rule=2)
  outsir <- sir_1(f_beta=f_beta,times = times)
  modCost(model = outsir,obs = datareal)
}

# Play with start values!!!
p <- rep(0.8,20)

# e.g.: consider random start values
set.seed(123)
p <- runif(20,min=0.8,max=1.2)


# try other solvers,especially such with true box constraints
Fit <- modFit(f = sir_cost,p = p,lower=rep(0.2,20),upper=rep(5,# box constraints
              method="Port")
summary(Fit) # system is singular (that is what we expected)

# use another solver. Note: it takes a while
Fit <- modFit(f = sir_cost,# box constraints
              method="L-BFGS-B")

# goes in a surprisingly good direction
Fit$par

f_beta <- approxfun(x=times,y=Fit$par,rule=2)
out2 <- sir_1(f_beta=f_beta,times = times)

# compare with data
plot(out,out2,obs=datareal)

# but see how unstable beta is
plot(out2) 

,

拟合具有时间相关参数的模型可能是一个好主意，但如果有理由这样做，我建议限制参数的数量并使用一种平滑函数。

以下示例显示了如何为此目的使用样条曲线，但当然也可以（并且可能更可取）使用具有某种机械意义的函数。

作为副作用，还可以识别伽马而不是先验地修复它。尽管如此，这仍然是一个技术演示，但我将科学问题悬而未决，即依赖于时间的测试版是否有意义。

library("FME")

sir_1 <- function(f_beta,gamma,parameters) {
    beta  <- parameters$f_beta(time)
    gamma <- parameters$gamma
    with(as.list(variables),{
      dS <- -beta * I * S / 10000
      dI <-  beta * I * S / 10000 - gamma * I
      dR <-  gamma * I
      # return vector of derivatives,and beta as auxiliary variable
      return(list(c(dS,beta = beta))
    })
  }

  initial_values <- c(S = S0,R = R0)
  # pass constant parameter and parameter function together as a list
  parameters <- list(
    f_beta = f_beta,gamma  = gamma
  )
  ode(initial_values,19)
datareal <- data.frame(
  time = times,I    = c(10,1156)
)

## define parameter as a vector: gamma and beta
t_beta <- c(0,12,16,19)      #  consider more or less knots
n_beta <- length(t_beta)
y_beta <- rep(1,n_beta)
p      <- c(gamma = 1/5,y_beta) # combine all parameters in one vector

## a small helper function for parameter selection
select <- function(p,which,exclude = FALSE) {
  parnames <- names(p)
  p[(which == parnames) != exclude]
}

## check the helper function
select(p,"gamma")
select(p,"gamma",excl=TRUE)

## cost function,see ?modCost help page
sir_cost <- function (p) {
  gamma   <- select(p,"gamma")
  y_beta  <- select(p,exclude = TRUE)
  f_beta  <- splinefun(x = t_beta,y = y_beta)
  outsir  <- sir_1(f_beta = f_beta,gamma = gamma,obs = datareal)
}

## model calibration,see ?modFit
Fit <- modFit(f = sir_cost,# lower bound to avoid negative values of beta
              lower = c(gamma = 0,rep(0.0,n_beta)),# note: high sensitivity wrt. upper bound
              upper = c(gamma=1,rep(2.0,# an algorithm that supports box constraints
              method = "Port")

## all parameters were identifiable
summary(Fit)

## smaller time steps to obtain a curves
times <- seq(0,19,0.1)

## split components of fitted parameters
gamma  <- select(Fit$par,"gamma")
y_beta <- select(Fit$par,exclude = TRUE)

out2 <- sir_1(f_beta = splinefun(x = t_beta,y = y_beta),times = times)

## show fitted curves and compare simulation with data
## see ?plot.deSolve help page
plot(out2,obs = datareal,which = c("S","R","I","beta"),las = 1,obspar = list(pch = 16,col = "red"))