在 R 中:FME/deSolve - SIR 拟合(时变参数)

In R: FME/ deSolve - SIR fitting (time varying parameters)

我正在尝试做的事情:我有一个简单的 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,  S0 = 9990,  I0 = 10, R0 = 0, 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 函数什么都不做,它只returns每个时间步的 0.8(即使我取一个非常高的值,如 800,它表示这已经是最合适的)。我的猜测是同一变量 (beta) 的随时间变化的值,我必须像文档中那样以另一种方式处理这个问题。

非常感谢任何帮助。

我不认为按时间步长估计 beta 是个好主意。这是问题固有的,而不是 deSolveFME 的错误。如果要使用动态模型来估计与时间相关的参数,我会建议使用节点较少的合适函数,例如时间相关的线性、二次或样条曲线,例如 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, method="constant", rule=2)
  outsir <- sir_1(f_beta=f_beta,  S0 = 9990,  I0 = 10, R0 = 0, 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, 20), # 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, p = p, 
              lower=rep(0.2, 20), upper=rep(5, 20), # box constraints
              method="L-BFGS-B")

# goes in a surprisingly good direction
Fit$par

f_beta <- approxfun(x=times, y=Fit$par, method="constant", rule=2)
out2 <- sir_1(f_beta=f_beta,  S0 = 9990,  I0 = 10, R0 = 0, 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, S0, I0, R0, times) {
  # the differential equations
  sir_equations <- function(time, variables, 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, dI, dR), beta = beta))
    })
  }

  initial_values <- c(S = S0, I = I0, R = R0)
  # pass constant parameter and parameter function together as a list
  parameters <- list(
    f_beta = f_beta,
    gamma  = gamma
  )
  ode(initial_values, times, sir_equations, parameters)
}

times <- seq(0, 19)
datareal <- data.frame(
  time = times,
  I    = c(10, 32, 120, 230, 480, 567, 1040, 1743, 2300,
           2619, 3542, 4039, 4231, 6378,
           5356, 4987, 3421, 2789, 1789, 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, "gamma", exclude = TRUE)
  f_beta  <- splinefun(x = t_beta, y = y_beta)
  outsir  <- sir_1(f_beta = f_beta, gamma = gamma,
                   S0 = 9990,  I0 = 10, R0 = 0, times = times)
  modCost(model = outsir, obs = datareal)
}

## model calibration, see ?modFit
Fit <- modFit(f = sir_cost, p = p,
              # 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, n_beta)),
              # 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, "gamma", exclude = TRUE)

out2 <- sir_1(f_beta = splinefun(x = t_beta, y = y_beta), gamma,
              S0 = 9990,  I0 = 10, R0 = 0, 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"))