First off, this article is not going to show you the trajectory of the current coronavirus pandemic. I think we have enough people doing this with results all over the place.

I became interested in epidemiology as a mathematics undergraduate student in the early 1990s, and one of my first jobs was calculating confidence intervals for data in EpiInfo for a team of practicing acute epidemiologists. I dusted off the epidemiology in my head and started looking for code examples that would get me a head start on doing some analysis. As usual, I found a lot of completed analyses, a lot of self-congratulatory complicated code, but very little that was simple enough to be a building block for my knowledge. For me, this can lead to the unfortunate situation of accepting that code works without knowing everything that it is doing, and I learned long ago not to fall into this trap. It will bite you badly at some point.

I am a big fan of Brady Haran’s Numberphile YouTube channel, and he posted a nice simple video (https://youtu.be/k6nLfCbAzgo) of mathematician Ben Sparks working through how to build a SIR (Susceptible, Infected, Recovered), epidemiological model, using Geogebra.

Solving the SIR model means solving an ordinary differential equation problem, and most analytical software handles this pretty succinctly. I work for SAS, and most of the time, my clients want to use SAS to solve their analytical problems. My background and experience are a bit more varied, and I like to explore how to do things using various tools as it helps me help them incorporate their code into SAS workflows in the “use the tool you love” mentality. This exercise has reinforced to me that SAS PROC MODEL is a deeply developed tool with a rich set of solvers all in one place, backed up by SASs vast academic center of excellence.

So now would be an excellent time to take a look at the YouTube video and the Geogebra workbook: https://www.geogebra.org/m/aqpv5df7. I created this Geogebra workbook as precisely as in the video as possible, so the parameters are not formally named.









There are several ways to construct this example using SAS software. SAS PROC IML contains a rich set of matrix-based tools if you wanted to build an ODE solver anywhere from scratch to the fully parameterized solver that simply needs to be called. For this example, I decided to focus on PROC MODEL, which is an interface to a large selection of solvers that can also be simple or very complex. Lastly, for SIR models in SAS, the tried and true data step is also an option. I have provided a link to a paper at the bottom of this section that demonstrates the data step approach.

DATA d_scaffold; DO time = 1 TO 20; OUTPUT; END; RUN; PROC MODEL data ​=d_scaffold plots=all outmodel=sir; parameters transm 1.9 recov 0.18 ; dependent S_T 0.99 R_T 0 I_T 0.01 ; DERT.S_T = -transm * S_T * I_T; DERT.I_T = transm * S_T * I_T - recov * I_T; DERT.R_T = recov * I_T; SOLVE S_T I_T R_T / out=SIR_SOLVE ; RUN;QUIT; ods graphics / reset width=7in height=5in imagemap; proc sgplot data=sir_solve; series x=time y=S_T / lineattrs=(color="blue") smoothconnect; series x=time y=I_T / lineattrs=(color="red") smoothconnect; series x=time y=R_T / lineattrs=(color="green") smoothconnect; run; ods graphics / reset;

And the output from the graph of the simulated S, I and R functions:

This PharmaSug paper is a great place to start, and I gleaned much from it for this article: https://www.pharmasug.org/proceedings/china2014/SP/PharmaSUG-China-2014-SP03.pdf

If you do not have access to SAS software, you can get a free learning version here:













There are also a number of libraries that solve ordinary differential equations. I opted for the deSolve library. Notice that the R code is very similar to the PROC MODEL code in its general configuration. You set some initial values, then specify the first derivatives explicitly. The ode function then uses these to solve the system.

library(deSolve) library(ggplot2) N <- 1 Istart <- 0.01 Sstart <- N-Istart Rstart <- 0 transm <- 3.25 # beta recov <- 0.3 # gamma maxT <- 20 init <- c(S = Sstart, I = Istart, R = Rstart) parameters <- c(transm, recov) times <- seq(0, maxT, by = 1) ## Create an SIR function (adapted from deSolve Examples) sir <- function(time, state, parameters) { with(as.list(c(state, parameters)), { dS <- -transm * S * I dI <- transm * S * I - recov * I dR <- recov * I return(list(c(dS, dI, dR))) }) } # Run ODE Solver results <- ode(y = init, times = times, func = sir, parms = parameters) #plot SIR Function plotds <- data.frame(results) names(plotds) <- c('time','S','I','R') ggplot(data=plotds) + geom_line(aes(x=time,y=S),color='green') + geom_line(aes(x=time,y=I),color='darkred') + geom_line(aes(x=time,y=R),color='blue') + theme_bw()

The best way to understand the output of this simulation is to simply graph it:









The python version of this code is again very similar to both the SAS and R manifestations. I chose the odeint function from scipy.integrate as the solver, and I am sure there are many other tools out there that are accessible from python.

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt import pandas as pd N = 1 Istart = 0.01 Sstart = N - Istart Rstart = 0 transm = 3.25 recov = 0.30 maxT = 20 t = np.linspace(0, maxT, maxT) def deriv(SIR, t, N, transm, recov): S, I, R = SIR dSdt = -transm * S * I dIdt = transm * S * I - recov * I dRdt = recov * I return dSdt, dIdt, dRdt SIR0 = Sstart, Istart, Rstart ret = odeint(deriv, SIR0, t, args=(N, transm, recov)) S, I, R = ret.T plotData = pd.DataFrame(ret.T) plotData = plotData.transpose() plotData.columns = ['S','I','R'] with pd.plotting.plot_params.use('x_compat',True): plotData['S'].plot(color='g') plotData['I'].plot(color='r') plotData['R'].plot(color='b')

And we would expect the graph from Python to look the same as the others, and so it does:

Slightly Off-Topic: SAS, R and Python Integrated

If you have SAS Viya, you can also access the SAS solvers using R or Python code using the libraries provided by our GitHub repositories.

SAS has put a lot of energy into embracing what our customers want when it comes to other free and open-source analytical solutions. What this means is SAS is working to build interfaces in both directions so that an analyst can have the most amount of flexibility in how they work with code. Another example is model management. Our model management solution. SAS Open Model Manager, accepts and keeps up with code written in just about anything and streamlines the relationship between modeling and DevOps.

The times have definitely changed. I listened to an internal tech talk recently given by a SAS employee who mentioned that they had never written a data step but do their work entirely in Python, running on a SAS Viya cluster. If you haven't kept up with SAS over the last ten years, you will be very surprised by the additions to the platform. As always though, SAS remains backward compatibility as much as is possible. You can still run those programs you wrote in 1995 without worry:)

All the code in this article is available here: https://github.com/scoyote/SimpleEpiModels