COVID-19 modelling - for fun
BInh Thang Tran
5/22/2020
library(deSolve) # package to solve the model
library(reshape2) # package to change the shape of the model output
library(ggplot2) # package for plottinginitial_number_infected <- 100000000 # Number of population in Vietnam
initial_number_recovered <- 0 # the initial number of people in "0"
# the recovered state
recovery_rate <- 0.071 # the rate of recovery gamma,
# in units of days^-1
follow_up_duration <- 100 # the duration to run the model for,
# in units of days# The initial state values are stored as a vector
# and each value is assigned a name.
initial_state_values <- c(I = initial_number_infected,
R = initial_number_recovered)
# Parameters are also stored as a vector with assigned names and values.
parameters <- c(gamma = recovery_rate)
# In this case we only have one parameter, gamma.# The times vector creates a sequence of timepoints at which we want to
# calculate the number of people in the I and R compartment.
times <- seq(from = 0, to = follow_up_duration, by = 1) cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dI <- -gamma*I
dR <- gamma*I
return(list(c(dI, dR)))
})
}output <- as.data.frame(ode(y = initial_state_values, times = times, func = cohort_model, parms = parameters))output## time I R
## 1 0 100000000.00 0
## 2 1 93146189.21 6853811
## 3 2 86762125.67 13237874
## 4 3 80815613.74 19184386
## 5 4 75276664.45 24723336
## 6 5 70117344.29 29882656
## 7 6 65311634.18 34688366
## 8 7 60835298.34 39164702
## 9 8 56665762.09 43334238
## 10 9 52781997.97 47218002
## 11 10 49164419.69 50835580
## 12 11 45794783.38 54205217
## 13 12 42656095.57 57343904
## 14 13 39732527.49 60267473
## 15 14 37009335.23 62990665
## 16 15 34472785.41 65527215
## 17 16 32110085.92 67889914
## 18 17 29909321.39 70090679
## 19 18 27859393.09 72140607
## 20 19 25949963.00 74050037
## 21 20 24171401.63 75828598
## 22 21 22514739.49 77485261
## 23 22 20971621.85 79028378
## 24 23 19534266.56 80465733
## 25 24 18195424.89 81804575
## 26 25 16948344.90 83051655
## 27 26 15786737.40 84213263
## 28 27 14704744.29 85295256
## 29 28 13696908.94 86303091
## 30 29 12758148.71 87241851
## 31 30 11883729.34 88116271
## 32 31 11069241.01 88930759
## 33 32 10310576.18 89689424
## 34 33 9603908.79 90396091
## 35 34 8945675.06 91054325
## 36 35 8332555.41 91667445
## 37 36 7761457.83 92238542
## 38 37 7229502.19 92770498
## 39 38 6734005.79 93265994
## 40 39 6272469.78 93727530
## 41 40 5842566.57 94157433
## 42 41 5442128.11 94557872
## 43 42 5069134.94 94930865
## 44 43 4721706.02 95278294
## 45 44 4398089.23 95601911
## 46 45 4096652.51 95903347
## 47 46 3815875.70 96184124
## 48 47 3554342.80 96445657
## 49 48 3310734.87 96689265
## 50 49 3083823.36 96916177
## 51 50 2872463.95 97127536
## 52 51 2675590.70 97324409
## 53 52 2492210.78 97507789
## 54 53 2321399.37 97678601
## 55 54 2162295.05 97837705
## 56 55 2014095.43 97985905
## 57 56 1876053.14 98123947
## 58 57 1747472.01 98252528
## 59 58 1627703.59 98372296
## 60 59 1516143.86 98483856
## 61 60 1412230.23 98587770
## 62 61 1315438.64 98684561
## 63 62 1225280.97 98774719
## 64 63 1141302.53 98858697
## 65 64 1063079.81 98936920
## 66 65 990218.33 99009782
## 67 66 922350.64 99077649
## 68 67 859134.47 99140866
## 69 68 800251.02 99199749
## 70 69 745403.33 99254597
## 71 70 694314.80 99305685
## 72 71 646727.77 99353272
## 73 72 602402.28 99397598
## 74 73 561114.76 99438885
## 75 74 522657.02 99477343
## 76 75 486835.10 99513165
## 77 76 453468.34 99546532
## 78 77 422388.48 99577612
## 79 78 393438.77 99606561
## 80 79 366473.22 99633527
## 81 80 341355.84 99658644
## 82 81 317959.96 99682040
## 83 82 296167.58 99703832
## 84 83 275868.82 99724131
## 85 84 256961.29 99743039
## 86 85 239349.65 99760650
## 87 86 222945.08 99777055
## 88 87 207664.84 99792335
## 89 88 193431.89 99806568
## 90 89 180174.43 99819826
## 91 90 167825.62 99832174
## 92 91 156323.17 99843677
## 93 92 145609.07 99854391
## 94 93 135629.30 99864371
## 95 94 126333.53 99873666
## 96 95 117674.87 99882325
## 97 96 109609.65 99890390
## 98 97 102097.22 99897903
## 99 98 95099.67 99904900
## 100 99 88581.71 99911418
## 101 100 82510.49 99917490#number of people have recovered after 2 months
output[output$time == 60, c("time","R")] ## time R
## 61 60 98587770#proportion of the total popultion corespond to
output[output$time == 60,"R"]/(output[output$time == 60,"I"]+
output[output$time == 60,"R"]) # note this is only on a separate line ## [1] 0.9858777#plot
# First turn the output dataset into a long format,
# so that the number in each compartment at each timestep
# are all in the same column
output_long <- melt(as.data.frame(output), id = "time")
# Plot the number of people in each compartment over time
ggplot(data = output_long, # specify object containing data to plot
aes(x = time, y = value, colour = variable, group = variable)) + # assign columns to axes and groups
geom_line() + # represent data as lines
xlab("Time (days)")+ # add label for x axis
ylab("Number of people") + # add label for y axis
labs(title = paste("Number infected and recovered over time when gamma =", parameters["gamma"],"days^-1")) # add title
Based on the plot, at what timepoint were infected and recovered individuals equal in number?
output[output$time == 10,]## time I R
## 11 10 49164420 50835580#4.1 What is the recovery rate in 7 days?
parameters <- c(gamma = 0.142)
output <- as.data.frame(ode(y = initial_state_values, times = times, func = cohort_model, parms = parameters))
# Plotting the output
output_long <- melt(as.data.frame(output), id = "time")
# turn output dataset into long format
ggplot(data = output_long, # specify object containing data to plot
aes(x = time,
y = value,
colour = variable,
group = variable)) + # assign columns to axes and groups
geom_line() + # represent data as lines
xlab("Time (days)")+ # add label for x axis
ylab("Number of people") + # add label for y axis
labs(title = paste("Number infected and recovered over time when gamma =",
parameters["gamma"],"days^-1")) + # add title
scale_color_brewer(palette = "Set1")