COVID-19 modelling - for fun

library(deSolve)   # package to solve the model
library(reshape2)  # package to change the shape of the model output
library(ggplot2)   # package for plotting

1

initial_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

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

3 Next will be competeing model: Competing Hazards in Compartmental Models

initial_state_values <- c(I = 100000000,   # at the start, there are 10^6 
                                         # infected people
                          R = 0,         # no one has recovered
                          M = 0)         # no one has died 

parameters <- c(gamma = 0.071,       # the recovery rate gamma is 14 days^-1
                mu = 0.038)          # the mortality rate mu is 

times <- seq(from = 0, to = 100, by = 1) # model the course of the
                                         # infection over 100 days

# remember that the rates are in units of days^-1

The model function describing this set of differential equations looks like this:

# As always, the model function takes as input arguments
# (in the following order): time, state and parameters

cohort_model <- function(time, state, parameters) {  

    with(as.list(c(state, parameters)), {  # tell R to unpack the state
                                           # variable names (I, R and M)
                                           # and model parameters 
                                           # (gamma and mu) from the inputs
                                           # "state" and "parameters"
        
    # The differential equations
      dI <- -gamma * I - mu * I  # included -mu * I
      dR <- gamma * I            # no change
      dM <- mu * I               # a new line to describe the rate of 
                                 # change in the M compartment
      
    # Return the number of people in each compartment at each
    # timestep (in the same order as the input state variables)
    return(list(c(dI, dR, dM))) 
    })
  
}

1.0.1 After 100 days, do you expect more people to have recovered or more people to have died, and why?

We expect more people to die than to recover because the mortality rate (0.038) is higher than the recovery rate (0.071), so people move more quickly from I to M than from I to R.

output1 <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = cohort_model,
                            parms = parameters))

Plotting the output over time indeed shows that more people have indeed died than recovered:

output_long <- melt(as.data.frame(output1), id = "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(colour = "Compartment")        # add legend title

1.0.2 Based on the model output, what proportion of the initially infected cohort died before recovering over the 4 week period?

output1

##     time            I        R        M
## 1      0 1.000000e+08        0        0
## 2      1 8.967304e+07  6726734  3600224
## 3      2 8.041254e+07 12758801  6828654
## 4      3 7.210837e+07 18167940  9723686
## 5      4 6.466177e+07 23018479 12319749
## 6      5 5.798418e+07 27368104 14647718
## 7      6 5.199618e+07 31268546 16735278
## 8      7 4.662655e+07 34766191 18607257
## 9      8 4.181145e+07 37902635 20285917
## 10     9 3.749360e+07 40715180 21791223
## 11    10 3.362165e+07 43237275 23141077
## 12    11 3.014955e+07 45498914 24351531
## 13    12 2.703602e+07 47526994 25436983
## 14    13 2.424402e+07 49345636 26410340
## 15    14 2.174035e+07 50976467 27283179
## 16    15 1.949524e+07 52438883 28065881
## 17    16 1.748197e+07 53750275 28767753
## 18    17 1.567662e+07 54926241 29397143
## 19    18 1.405770e+07 55980766 29961537
## 20    19 1.260597e+07 56926390 30467645
## 21    20 1.130415e+07 57774359 30921488
## 22    21 1.013678e+07 58534760 31328463
## 23    22 9.089957e+06 59216634 31693410
## 24    23 8.151241e+06 59828091 32020668
## 25    24 7.309465e+06 60376403 32314131
## 26    25 6.554620e+06 60868092 32577288
## 27    26 5.877727e+06 61309004 32813269
## 28    27 5.270737e+06 61704383 33024881
## 29    28 4.726430e+06 62058931 33214639
## 30    29 4.238333e+06 62376865 33384801
## 31    30 3.800642e+06 62661967 33537391
## 32    31 3.408152e+06 62917626 33674222
## 33    32 3.056193e+06 63146883 33796923
## 34    33 2.740581e+06 63352465 33906953
## 35    34 2.457563e+06 63536817 34005620
## 36    35 2.203771e+06 63702131 34094098
## 37    36 1.976189e+06 63850372 34173439
## 38    37 1.772109e+06 63983305 34244586
## 39    38 1.589104e+06 64102510 34308386
## 40    39 1.424998e+06 64209405 34365597
## 41    40 1.277839e+06 64305261 34416900
## 42    41 1.145877e+06 64391218 34462905
## 43    42 1.027543e+06 64468298 34504159
## 44    43 9.214287e+05 64537418 34541153
## 45    44 8.262731e+05 64599400 34574327
## 46    45 7.409442e+05 64654981 34604074
## 47    46 6.644272e+05 64704823 34630750
## 48    47 5.958121e+05 64749517 34654671
## 49    48 5.342828e+05 64789596 34676122
## 50    49 4.791077e+05 64825535 34695357
## 51    50 4.296304e+05 64857764 34712606
## 52    51 3.852627e+05 64886664 34728074
## 53    52 3.454767e+05 64912579 34741944
## 54    53 3.097995e+05 64935819 34754382
## 55    54 2.778066e+05 64956658 34765535
## 56    55 2.491177e+05 64975345 34775537
## 57    56 2.233914e+05 64992103 34784506
## 58    57 2.003218e+05 65007130 34792548
## 59    58 1.796347e+05 65020605 34799760
## 60    59 1.610839e+05 65032688 34806228
## 61    60 1.444488e+05 65043524 34812027
## 62    61 1.295317e+05 65053241 34817227
## 63    62 1.161550e+05 65061954 34821891
## 64    63 1.041597e+05 65069768 34826073
## 65    64 9.340317e+04 65076774 34829823
## 66    65 8.375746e+04 65083057 34833185
## 67    66 7.510787e+04 65088691 34836201
## 68    67 6.735151e+04 65093744 34838905
## 69    68 6.039615e+04 65098274 34841330
## 70    69 5.415906e+04 65102337 34843504
## 71    70 4.856608e+04 65105980 34845454
## 72    71 4.355068e+04 65109247 34847203
## 73    72 3.905322e+04 65112176 34848770
## 74    73 3.502021e+04 65114803 34850176
## 75    74 3.140369e+04 65117159 34851437
## 76    75 2.816064e+04 65119272 34852568
## 77    76 2.525250e+04 65121166 34853582
## 78    77 2.264469e+04 65122864 34854491
## 79    78 2.030618e+04 65124388 34855306
## 80    79 1.820917e+04 65125754 34856037
## 81    80 1.632872e+04 65126979 34856693
## 82    81 1.464246e+04 65128077 34857281
## 83    82 1.313034e+04 65129062 34857808
## 84    83 1.177437e+04 65129945 34858280
## 85    84 1.055844e+04 65130737 34858704
## 86    85 9.468071e+03 65131447 34859085
## 87    86 8.490308e+03 65132084 34859425
## 88    87 7.613517e+03 65132655 34859731
## 89    88 6.827272e+03 65133168 34860005
## 90    89 6.122223e+03 65133627 34860251
## 91    90 5.489983e+03 65134039 34860471
## 92    91 4.923035e+03 65134408 34860669
## 93    92 4.414635e+03 65134739 34860846
## 94    93 3.958738e+03 65135036 34861005
## 95    94 3.549920e+03 65135302 34861148
## 96    95 3.183322e+03 65135541 34861276
## 97    96 2.854581e+03 65135755 34861390
## 98    97 2.559790e+03 65135947 34861493
## 99    98 2.295441e+03 65136119 34861585
## 100   99 2.058392e+03 65136274 34861668
## 101  100 1.845823e+03 65136412 34861742

output1[output1$time == 100,]

##     time        I        R        M
## 101  100 1845.823 65136412 34861742

output1[102,"M"]/100000000

## [1] NA

1.0.3 Now use the competing hazards formula given in the video lecture to calculate the case fatality rate. Does this agree with your answer to the previous question?

# Calculate this by referring to the values in the parameters vector:
parameters["mu"]/(parameters["mu"]+parameters["gamma"])

##        mu 
## 0.3486239