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

Question

我想做的是：我有一个简单的SIR模型，具有时变的传输速率β，我已经在R中实现了这一点(感谢@tpetzoldt)。我们有一个N=10000种群，伽马也是固定的。

  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函数不执行nothing，它只返回每个时间步骤的0.8 (即使我的值非常高，比如800，它说这已经是最适合的了)。我的猜测是，对于相同变量(beta)的时变值，我必须以另一种方式来处理，就像在documentation.

中那样

任何帮助都是非常感谢的。

Answer 1

我不认为估计每个时间步骤的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, 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) 

Answer 2

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

下面的示例演示如何为此目的使用样条，但当然也可以(也可能更可取)使用具有某种机械意义的函数。

作为一个副作用，它也有可能识别伽马，而不是事先固定它。尽管如此，这仍然是一个技术演示，但我没有提出科学问题，一个时间依赖的测试版是否会有任何意义。

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"))

​

​