THIS IS A TEST OF THE REPORTING AND PUBLISHING CAPABILITIES OF THE reticulate PACKAGE LINKING PYTHON TO RMARKDOWN. THIS REPORT IS NOT INTENDED TO PROVIDE ANY INSIGHT INTO THE CURRENT COVID-19 PANDEMIC IN MIAMI-DADE COUNTY. INSPECT AT YOUR OWN RISK.

Setup

We need to load some R packages before we can perform our computing.

library(reticulate)
library(tidyverse)

The shell script to install the packages we need is in src/install_python_packages_20200403.sh. This creates a conda environment called env_sflCHIME and activates it. At current, I have been unable to run these commands on Windows.

On Mac, I have a Conda environment created for this project. In the terminal, I type which python and it tells me

which python
/Users/gabrielodom/anaconda3/envs/env_sflCHIME/bin/python

If it says

/usr/bin/python

(or some other place), then something is wrong.

I specifically want this python environment, so I “point” reticulate to the particular environment that we installed with Anaconda, replete with the version of Python we selected and all associated libraries. (NOTE: this line must be executed before any other reticulate functions or Python calls, otherwise it will have no effect.)

use_python(
    "/Users/gabrielodom/anaconda3/envs/env_sflCHIME/bin/python",
    required = TRUE
)

We also have some Python libraries to load.

from functools import reduce
from typing import Tuple, Dict, Any
import pandas as pd
import streamlit as st
import numpy as np
import altair as alt

IT WORKS!!!!!! (This seriously took a full weekend to figure out, and then I called in the big guns: Tim N. and Athina H. Note that this is not compatible with environments created via Anaconda, at least not by default.)


Parameters

Now, we want to set our initial parameters. We organise them the same way as on the CHARM app.

Hospital Parameters

# Regional Population
S_default = 2761581.0

# Hospital Market Share
Penn_market_share = 0.15

# Currently Hospitalised Patients (at current hospital?)
known_cases = 400.0
current_hosp = known_cases

# Currently Hospitalised Patients
known_infections = 1000.0
initial_infections = known_infections

Spread and Contact Parameters

# Doubling Time
doubling_time = 6.0

# Social Distancing (% reduction in social contact)
relative_contact_rate = 0.0

Severity Rates

# Hospitalisation (% of total infections)
hosp_rate = 0.05

# ICU (% of total infections)
icu_rate = 0.02

# Ventilated (% of total infections)
vent_rate = 0.01

Severity Times

# Infectious Days
recovery_days = 14.0

# Average hospital length of stay (days)
hosp_los = 7.0

# Average ICU length of stay (days)
icu_los = 9.0

# Average length of stay on a ventilator (days)
vent_los = 10.0

Directly Estimated Parameters

Given the parameters above, we calculate the following disease spread parameters.

SIR Initial Values

total_infections = current_hosp / Penn_market_share / hosp_rate
print(total_infections)
## 53333.333333333336
detection_prob = initial_infections / total_infections
print(detection_prob)
## 0.01875
Infct = initial_infections / detection_prob
Susc, Infct, Recov = S_default - Infct, Infct, 0.0
print(Susc)
## 2708247.6666666665
print(Infct)
# print(Recov)
# type(Recov)
## 53333.333333333336
intrinsic_growth_rate = 2 ** (1 / doubling_time) - 1
print(intrinsic_growth_rate)
## 0.12246204830937302

NOTE: the total infections and initial infections are equal at \(t = 0\).

SIR Transition Rates

# type(Recov)
# print(Recov)

gamma = 1.0 / recovery_days
# mean recovery rate, gamma, (in 1/days).
print(gamma)
## 0.07142857142857142
beta = (intrinsic_growth_rate + gamma) / Susc * (1.0 - relative_contact_rate)
# Contact rate, beta
# {rate based on doubling time} / {initial Susc}
print(beta)
## 7.15926472029741e-08
r_t = beta / gamma * Susc 
# Current Reproduction Number
# r_t is r_0 after distancing
print(r_t)
## 2.7144686763312222
r_naught = r_t / (1.0 - relative_contact_rate)
# Initial Reproduction Number
print(r_naught)
## 2.7144686763312222
doubling_time_t = 1.0 / np.log2(beta * Susc - gamma + 1.0) 
# Current Doubling Time
# doubling time after distancing
print(doubling_time_t)

# type(Recov)
## 5.999999999999998

As an aside, if I want to access any of these values in R (for instance, \(R_0\)), I can call them as py$r_naught, which yields 2.7144687.


Model Specification

We can now specify the SIR model at outset and over time.

# type(Recov)

# SIR at Outset
def sir(y, beta, gamma, N):
    S, I, R = y
    Sn = (-beta * S * I) + S
    In = (beta * S * I - gamma * I) + I
    Rn = gamma * I + R
    if Sn < 0:
        Sn = 0.0
    if In < 0:
        In = 0.0
    if Rn < 0:
        Rn = 0.0

    scale = N / (Sn + In + Rn)
    return Sn * scale, In * scale, Rn * scale
    
# type(Recov)

# Dynamic SIR 
def sim_sir(S, I, R, beta, gamma, n_days, beta_decay = None):
    N = S + I + R
    s, i, r = [S], [I], [R]
    for day in range(n_days):
        y = S, I, R
        S, I, R = sir(y, beta, gamma, N)
        if beta_decay:
            beta = beta * (1 - beta_decay)
        s.append(S)
        i.append(I)
        r.append(R)
            
    s, i, r = np.array(s), np.array(i), np.array(r)
    return s, i, r

# type(Recov)

# Test Function
def test_fun(x):
  y = x + 2
  return y
  
# type(Recov)

Run the model, assuming no decay in contact rate over the next 30 days.

n_days = 100
beta_decay = 0.0

test_fun(3)
# type(Susc)
# type(Infct)
# type(Recov)
## 5
y = Susc, Infct, Recov
# print(y)
# type(y)
N = Susc + Infct + Recov
print(N)
## 2761581.0
type(N)
## <class 'float'>
sir(y, beta, gamma, N)
## (2697906.8336139764, 59864.6425764999, 3809.5238095238096)
s, i, r = sim_sir(
  Susc, Infct, Recov, beta, gamma, n_days, beta_decay = beta_decay
)

ERROR: Error in py_call_impl(callable, dots$args, dots$keywords) : TypeError: unsupported operand type(s) for +: 'float' and 'type'
SOLUTION: I can’t have a Python object named R.

We now inspect the results of this model

res_df <- tibble(
    days  = seq_len(py$n_days + 1),
    Susc  = py$s,
    Infct = py$i,
    Recov = py$r
)

ggplot(data = res_df) +
    aes(x = days) +
    geom_line(aes(y = Susc), colour = "black") + 
    geom_line(aes(y = Infct), colour = "red") + 
    geom_line(aes(y = Recov), colour = "green")