```
import Pkg; Pkg.add("Distributions")
import Pkg; Pkg.add("Plots")
import Pkg; Pkg.add("LsqFit")
import Pkg; Pkg.add("PyPlot")
```

In this post I just wanted to play with julia and R together on the topic of agent based models. Agent based models have some particular advantages in that they capture the random effects of individuals in the model (and their shifting states). The downside of agent based models is of course that you have to program all of those transitions in your model. Regardless of that, they do show a much wider range of values and allow for very fined tune interactions and effects that you might miss with a deterministic model. I’m also using this as a way of exploring julia and the agents package which is built for this purpose.

## Prep

JuliaCall is a great R package that allows you to call Julia from R.

To start, we can take the following Agent Based SIR model from https://raw.githubusercontent.com/epirecipes/sir-julia/master/script/abm/abm.jl which provides a done of great epimodels in a couple of different languages. Regardless, the set-up of this program is very well done. We have agents who can take on one of three states, S, I, R. The probability of an infectious encounter is 5% and agents have an average of 10 contacts per period. The infectious period is 9 days. In theory, if we wanted to add some additional detail we could add a quarantine state and captured the effect of having some portion of agents move to quarantine and no longer transmit. That will be for another day.

```
using Agents
using Random
using DataFrames
using Distributions
using DrWatson
using StatsPlots
using BenchmarkTools
function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end;
mutable struct Person <: AbstractAgent
::Int64
id::Symbol
statusend
function init_model(β::Float64,c::Float64,γ::Float64,N::Int64,I0::Int64)
= @dict(β,c,γ)
properties = ABM(Person; properties=properties)
model for i in 1:N
if i <= I0
= :I
s else
= :S
s end
= Person(i,s)
p = add_agent!(p,model)
p end
return model
end;
function agent_step!(agent, model)
transmit!(agent, model)
recover!(agent, model)
end;
function transmit!(agent, model)
# If I'm not susceptible, I return
!= :S && return
agent.status = rand(Poisson(model.properties[:c]))
ncontacts for i in 1:ncontacts
# Choose random individual
= random_agent(model)
alter if alter.status == :I && (rand() ≤ model.properties[:β])
# An infection occurs
= :I
agent.status break
end
end
end;
function recover!(agent, model)
!= :I && return
agent.status if rand() ≤ model.properties[:γ]
= :R
agent.status end
end;
susceptible(x) = count(i == :S for i in x)
infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x);
= 0.1
δt = 400
nsteps = nsteps*δt
tf = 0:δt:tf;
t
= 0.05
β = 10.0*δt
c = rate_to_proportion(0.11,δt);
γ
= 1000
N = 10;
I0
= init_model(β,c,γ,N,I0)
abm_model
= [(:status, f) for f in (susceptible, infected, recovered)]
to_collect = run!(abm_model, agent_step!, nsteps; adata = to_collect);
abm_data, _
:t] = t; abm_data[!,
```

Looking at a single iteration:

```
<- JuliaCall::julia_eval("abm_data")
out
%>%
out ggplot(aes(step, infected_status))+
geom_line()+
labs(
title = "Single Run of Agent Based SIR Model"
)
```

Ok, that’s good, but in reality, we want to run the same model many times in order to get a range of possible values. I’ll do this with a simple function. Ideally, we could use glue to make some of the parameters inputs to our function so that we could model the effect of changing the number of contacts from 10 per day to 5 per day which could signify physical distancing.

```
<- function(contacts = 10, contacts_after_int = 5){
run_simulation <- formatC(contacts, digits = 1)
contacts_in <- formatC(contacts_after_int, digits = 1)
contacts_after_int_in <- glue::glue("
julia_script using Agents
using Random
using DataFrames
using Distributions
using DrWatson
using StatsPlots
using BenchmarkTools
function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end;
mutable struct Person <: AbstractAgent
id::Int64
status::Symbol
end
function init_model(β::Float64,c::Float64,γ::Float64,N::Int64,I0::Int64)
properties = @dict(β,c,γ)
model = ABM(Person; properties=properties)
for i in 1:N
if i <= I0
s = :I
else
s = :S
end
p = Person(i,s)
p = add_agent!(p,model)
end
return model
end;
function agent_step!(agent, model)
transmit!(agent, model)
recover!(agent, model)
end;
function transmit!(agent, model)
# If I'm not susceptible, I return
agent.status != :S && return
ncontacts = rand(Poisson(model.properties[:c]))
for i in 1:ncontacts
# Choose random individual
alter = random_agent(model)
if alter.status == :I && (rand() ≤ model.properties[:β])
# An infection occurs
agent.status = :I
break
end
end
end;
function recover!(agent, model)
agent.status != :I && return
if rand() ≤ model.properties[:γ]
agent.status = :R
end
end;
susceptible(x) = count(i == :S for i in x)
infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x);
δt = 0.1
nsteps = 400
tf = nsteps*δt
t = 0:δt:tf;
β = 0.03
co = δt> 14.0 ? {contacts_in} : {contacts_after_int_in}
c = co*δt
γ = rate_to_proportion(0.11,δt);
N = 1000
I0 = 10;
abm_model = init_model(β,c,γ,N,I0)
to_collect = [(:status, f) for f in (susceptible, infected, recovered)]
abm_data, _ = run!(abm_model, agent_step!, nsteps; adata = to_collect);
abm_data[!,:t] = t;
")
<- tempfile(fileext = ".jl")
tmp cat(julia_script, file = tmp)
julia_source(tmp)
<- JuliaCall::julia_eval("abm_data")
out
out }
```

Now we can run it a few times with a few different contact rates. Just to make things interesting, I am adding the option to but a before and after contact rate to allow people to simulate different intervention and the effect of waiting.

We have no social distancing, moderate distancing, and severe social distancing.

```
<- map(1:25, ~run_simulation(contacts = 10,
out_fill_none contacts_after_int = 10))
<- map(1:25, ~run_simulation(contacts = 10,
out_fill_10 contacts_after_int = 7))
<- map(1:25, ~run_simulation(contacts = 10,
out_fill_05 contacts_after_int = 2))
```

And then we can visualise our outputs (of course using our curve statistics to better capture the range of possibilities).

```
<- out_fill_10 %>%
steps10 bind_rows(.id = ".draws") %>%
group_by(step) %>%
curve_interval(infected_status, .width = c(.5, .8, .95)) %>%
mutate(contact = 7)
<- out_fill_05 %>%
steps05 bind_rows(.id = ".draws") %>%
group_by(step) %>%
curve_interval(infected_status, .width = c(.5, .8, .95)) %>%
mutate(contact = 5)
<- out_fill_none %>%
step_none bind_rows(.id = ".draws") %>%
group_by(step) %>%
curve_interval(infected_status, .width = c(.5, .8, .95)) %>%
mutate(contact = 10)
%>%
steps10 bind_rows(steps05) %>%
bind_rows(step_none) %>%
ggplot(aes(x = step, y = infected_status, group = contact)) +
geom_hline(yintercept = 1, color = "gray75", linetype = "dashed") +
geom_lineribbon(aes(ymin = .lower, ymax = .upper)) +
scale_fill_brewer() +
labs(
title = "Simulated SIR Curve for Infections",
subtitle = "Using an Agent Based Model",
y = "Infected Individuals"
)
```

It’s amazing but not unsurprising to see the impact of reducing the number of contacts per day on disease transmission.

## Reuse

## Citation

```
@online{dewitt2020,
author = {Michael DeWitt and Michael DeWitt},
title = {Julia {ABM} {SIR}},
date = {2020-08-09},
url = {https://michaeldewittjr.com/programming/2020-08-09-julia-abm-sir/julia-abm-sir.html},
langid = {en}
}
```