在 R 中，SIR 模型，不是常量参数

Question

在此页面上显示了 R 中的 SIR 模型，https://rstudio-pubs-static.s3.amazonaws.com/382648_93783f69a2fd4df98ade8751c21abbad.html，它的解决方案以及 $\beta$ 和 $\gamma$ 参数的优化也被执行。（见下文）

在此代码中，假定 $\beta$ 和 $\gamma$ 在整个时间内保持不变。我想要的是有一个随时间变化的 beta，它不需要每天都改变，我们有 14 天的数据，如果它在 7 天后改变就足够了，即我们有 $\beta_1$几天[0:6] 和 $\beta_2$ 几天 [7:13] ，然后对两者执行如下优化算法，即最后我想接收一个矢量以获得 (\beta_1 的最佳值， \beta_2, \gamma) 而 gamma 一直保持不变。是否可以修改给定的代码？如果是，有人可以帮助如何修改它以接收所需的输出。

   day  cases
   0    1
   1    6
   2    26
   3    73
   4    222
   5    293
   6    258
   7    236
   8    191
   9    124
   10   69
  11    26
  12    11
  13    4

 #here beta is assumed to be constant
sir_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
dS <- -beta * I * S
dI <-  beta * I * S - gamma * I
dR <-  gamma * I
return(list(c(dS, dI, dR)))
})
}

parameters_values <- c(
beta  = 0.004, # infectious contact rate (/person/day)
gamma = 0.5    # recovery rate (/day)
)

initial_values <- c(
S = 999,  # number of susceptibles at time = 0
I =   1,  # number of infectious at time = 0
R =   0   # number of recovered (and immune) at time = 0
)

time_values <- seq(0, 10) # days


sir_values_1 <- ode(
y = initial_values,
times = time_values,
func = sir_equations,
parms = parameters_values 
)

sir_values_1

sir_values_1 <- as.data.frame(sir_values_1)
sir_values_1

sir_1 <- function(beta, gamma, S0, I0, R0, times) {
require(deSolve) # for the "ode" function

# the differential equations:
sir_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
  dS <- -beta * I * S
  dI <-  beta * I * S - gamma * I
  dR <-  gamma * I
  return(list(c(dS, dI, dR)))
  })
  }

# the parameters values:
parameters_values <- c(beta  = beta, gamma = gamma)

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

# solving
out <- ode(initial_values, times, sir_equations, parameters_values)

# returning the output:
as.data.frame(out)
}
sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = seq(0, 10))

flu <- read.table("https://uc8f29367cc06ca2f989ead2cd8e.dl.dropboxusercontent.com/cd/0/inline/BNzBF_deK5fmfGXWCB9a5YO95JkiLNFRc2Jq1w-qGNqQMXxnpn-yL-cAVoE1JQG7D4Od_SkG8YVKesqBr7wMoQHHSTNbHU_hhyahK7up0EDEft-u7Vf4xZJvu4cTNuUjXFb-QaHlOfBPnFhKspeb7RbO/file", header = TRUE)

predictions <- sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = flu$day)
predictions

model_fit <- function(beta, gamma, data, N = 763, ...) {
I0 <- data$cases[1] # initial number of infected (from data)
times <- data$day   # time points (from data)
# model's predictions:
predictions <- sir_1(beta = beta, gamma = gamma,   # parameters
                   S0 = N - I0, I0 = I0, R0 = 0, # variables' intial values
                   times = times)                # time points
# plotting the observed prevalences:
 with(data, plot(day, cases, ...))
# adding the model-predicted prevalence:
with(predictions, lines(time, I, col = "red"))
}

predictions <- sir_1(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = flu$day)
predictions

ss <- function(beta, gamma, data = flu, N = 763) {
I0 <- data$cases[1]
times <- data$day
 predictions <- sir_1(beta = beta, gamma = gamma,   # parameters
                   S0 = N - I0, I0 = I0, R0 = 0, # variables' intial values
                   times = times)                # time points
sum((predictions$I[-1] - data$cases[-1])^2)
}
ss(beta = 0.004, gamma = 0.5)

beta_val <- seq(from = 0.0016, to = 0.004, le = 100)
ss_val <- sapply(beta_val, ss, gamma = 0.5)

min_ss_val <- min(ss_val)
min_ss_val

beta_hat <- beta_val[ss_val == min_ss_val]
beta_hat

 plot(beta_val, ss_val, type = "l", lwd = 2,
 xlab = expression(paste("infectious contact rate ", beta)),
 ylab = "sum of squares")
# adding the minimal value of the sum of squares:
abline(h = min_ss_val, lty = 2, col = "grey")
# adding the estimate of beta:
abline(v = beta_hat, lty = 2, col = "grey")

ss(beta = 0.004, gamma = 0.5)

ss2 <- function(x) {
ss(beta = x[1], gamma = x[2])
}
ss2(c(0.004, 0.5))



starting_param_val <- c(0.004, 0.5)
ss_optim <- optim(starting_param_val, ss2)

Answer 1

这当然是可能的。您所需要的只是梯度函数中的 if 语句：

beta <- if (time<6) beta1 else beta2

或

beta <- ifelse(time<6, beta1, beta2))

并确保您的参数向量包括 beta1 和 beta2。

在 R 中，SIR 模型，不是常量参数

In R, SIR model,not constant parameters

r

variables

optimization

parameters

ode