On this page a SIR model in R is shown, https://rstudio-pubs-static.s3.amazonaws.com/382648_93783f69a2fd4df98ade8751c21abbad.html, the solution of it and the optimization of the $\beta$ and $\gamma$ parameter is also executed. (see below)
In this code both $\beta$ and $\gamma$ are assumed to be constant over the whole time. What I want is to to have a time varying beta, it does not need to change each day, we have fourteen days of data, it would suffice if it would change after seven days, i.e we have $\beta_1$ for days[0:6] and $\beta_2$ for days[7:13] and then do the optimization algorithm like below for both, i.e. in the end I want to receive a vector for the optimal values of (\beta_1, \beta_2, \gamma) whereas gamma stayed constant the whole time. Would it be possible with a modification of the code given? If yes could someone help how to modify it to receive the desired output.
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) 
没有评论:
发表评论