Solution 2

1 Solution: Simulating competing hazards

The model we want to specify in this etivity has 3 compartments: \(I\) (infected), \(R\) (recovered) and \(M\) (dead). Like in the first etivity, infected people can recover at a rate \(\gamma\). Additionally they can also now die (transition to the \(M\) compartment) at a rate \(\mu\).

The differential equations for this model look like this: \[\begin{align} \frac{dI}{dt} & = -\gamma I -\mu I \\ \frac{dR}{dt} & = \gamma I \\ \frac{dM}{dt} & = \mu I \end{align}\]

This corresponds to the following model diagram:

As you can see, the equation describing the rate of change in the recovered (\(R\)) compartment (second line) is not affected by this addition. However, we need a new equation describing the rate of change in the deceased (\(M\)) compartment (third line). People move into this compartment from the infected compartment (\(I\)) at a rate \(\mu\) - this transition also needs to be added in the rate of change in the infected compartment \(I\) (first line).

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

Defining model input and timesteps

initial_state_values <- c(I = 1000000,   # 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.1,       # the recovery rate gamma is 0.1 days^-1
                mu = 0.2)          # the mortality rate mu is 0.2 days^-1

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

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

1.0.1 After 4 weeks, 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.2) is higher than the recovery rate (0.1), so people move more quickly from \(I\) to \(M\) than from \(I\) to \(R\).

Load the deSolve package and use it to solve the differential equations:

library(deSolve)
library(reshape2)
library(ggplot2)

output <- 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(output), 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

# ggplot automatically adds the new compartment as a separate colour/group 
# because it is contained within the output dataframe

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

# Printing the full model output:
output

##    time            I         R        M
## 1     0 1000000.0000      0.00      0.0
## 2     1  740818.2240  86393.93 172787.9
## 3     2  548811.6238 150396.13 300792.3
## 4     3  406569.6441 197810.12 395620.2
## 5     4  301194.1295 232935.29 465870.6
## 6     5  223130.0332 258956.66 517913.3
## 7     6  165298.7445 278233.75 556467.5
## 8     7  122456.2989 292514.57 585029.1
## 9     8   90717.8284 303094.06 606188.1
## 10    9   67205.4086 310931.53 621863.1
## 11   10   49786.9751 316737.67 633475.3
## 12   11   36883.0919 321038.97 642077.9
## 13   12   27323.6579 324225.45 648450.9
## 14   13   20241.8599 326586.05 653172.1
## 15   14   14995.5341 328334.82 656669.6
## 16   15   11108.9626 329630.35 659260.7
## 17   16    8229.7197 330590.09 661180.2
## 18   17    6096.7248 331301.09 662602.2
## 19   18    4516.5637 331827.81 663655.6
## 20   19    3345.9518 332218.02 664436.0
## 21   20    2478.7415 332507.09 665014.2
## 22   21    1836.2963 332721.23 665442.5
## 23   22    1360.3615 332879.88 665759.8
## 24   23    1007.7803 332997.41 665994.8
## 25   24     746.5819 333084.47 666168.9
## 26   25     553.0813 333148.97 666297.9
## 27   26     409.7326 333196.76 666393.5
## 28   27     303.5373 333232.15 666464.3
## 29   28     224.8659 333258.38 666516.8

# Printing the output at timestep/day 28:
output[output$time == 28,]

##    time        I        R        M
## 29   28 224.8659 333258.4 666516.8

# Divide the number of people who died over the 4 week period by the 
# number of people initially infected:
output[29,"M"]/1000000

## [1] 0.6665168

# Answer: proportion = 0.6665

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?

The formula is: \[\begin{align} CFR = \frac{\mu}{\mu+\gamma} \end{align}\]

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

##        mu 
## 0.6666667

# Answer: CFR = 0.6667
# This is approximately the same as the case fatality rate
# observed in the model output. The agreement would be exact
# if we allowed the simulation to run for longer.

1.0.4 Which value of \(\mu\) do you need to get a case fatality rate of 50% assuming \(\gamma\) stays fixed?

Rearranging the CFR equation above to solve for \(\mu\) gives:

\[\begin{align} CFR = \frac{\mu}{\mu+\gamma} \\ \mu = CFR(\mu+\gamma) \\ \mu = \mu CFR + \gamma CFR \\ \mu - \mu CFR = \gamma CFR \\ \mu(1-CFR) = \gamma CFR \\ \mu = \frac{\gamma CFR}{1-CFR} \end{align}\]

Calculating \(\mu\) with CFR = 0.50 and \(\gamma\) = 0.1:

\[\begin{align} \mu = \frac{0.1 * 0.50}{1-0.50} \\ \mu = 0.1 \end{align}\]

This makes sense! If \(\mu\) and \(\gamma\) are equal, then they represent two competing hazards that are also equal. Thus, half of people die and half recover, so the case fatality rate is 50%.

Double-check this by modifying the code to simulate the model using these parameters:

# The only thing we are changing is mu, so we only need to update
# the parameters vector and solve the model with this
parameters <- c(gamma = 0.1, mu = 0.1)

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

# Calculate the case fatality rate:
output[29,"M"]/1000000

## [1] 0.4981511

# This indeed gives us approximately 0.50

Again we can plot the output and compare it with what we got with our previous parameters:

output_long <- melt(as.data.frame(output), 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              

As you notice, we only see red and blue lines here representing the number of people in the \(I\) and \(M\) compartment, despite having plotted all 3 compartments. This is just because, with \(\mu\) = \(\gamma\) and \(R(0) = M(0)\) (the initial number recovered and deceased), the number of recovered and deceased people is identical at each timestep so the lines completely overlap.