Agent based models in Mata: Modelling aggregate processes, like the spread of a
> disease
using the abm package
Maarten Buis
University of Konstanz
maarten.buis@uni.kn
-------------------------------------------------------------------------------
index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------
What is an Agent Based Model?
An Agent Based Model is a simulation in which agents, that each follow
simple rules, interact with one another and thus produce a often
surprising outcome at the macro level.
The purpose of an ABM is to explore mechanisms through which actions of
the individual agents add up to a macro outcome, by varying the rules
that agents have to follow or varying with whom the agent can interact.
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------
Example: The spread of a disease
In this model each agent can be in one of three states:
Susceptible: the agent is healthy, but can get the disease.
Infectious: the agent has the disease and can spread it to others.
Recovered: the agent has had the disease and is now immune.
In the first round one or more agents get infected.
In all subsequent rounds:
- All infectious agents have a chance to become recovered, and
- all recovered agents get a chance to become susceptible.
- All agents that are still infectious contact a random number of
other agents
- Each of these contact has a chance of resulting in an infection if
someone is susceptible.
. set seed 123456789
. set scheme s1color
. clear mata
. drop _all
. run sir_main.do
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: model = sir()
: model.N(2000)
: model.tdim(150)
: model.outbreak(5)
: model.mcontacts(4)
: model.transmissibility(0.045)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(sir1, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
.
. drop _all
The number 0.045 for transmissibility was chosen to get a >> R0of
2.5, which in line with COVID-19.
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------
Adding a network
The model thus far does not need a simulation; one can compute them with
a set of differential equations.
In Stata you can do that with epimodels by Sergiy Radyakin and Paolo
Verme
However, the advantage of Agent Based Models is that it is much easier to
expand the model.
For example, what if we don't believe that people contact each other
randomly but instead it is a network that determines who contacts who.
A common model that works fairly well for a lot of social networks is a
small world network:
We assume all agents are on a circle, and that each agent is
connected to the 4 closes agents.
However, a small number of these connections are replaced by a
connection to a random agent
. clear mata
. drop _all
. set seed 28863
. run nw_main.do
. mata:
------------------------------------------------- mata (type end to exit) -----
: model = sir_nw()
: model.N(100)
: model.tdim(10)
: model.outbreak(2)
: model.degree(4)
: model.pr_long(.05)
: model.transmissibility(0.1)
: model.mindur(5)
: model.meandur(6)
: model.pr_loss(0.00)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..........
: model.export_nw()
: end
-------------------------------------------------------------------------------
.
. gen x_ego = cos((ego-1)/`n'*2*_pi)*(1+mod(ego,2)*.2)
. gen y_ego = sin((ego-1)/`n'*2*_pi)*(1+mod(ego,2)*.2)
. gen x_alter = cos((alter-1)/`n'*2*_pi)*(1+mod(alter,2)*.2)
(1 missing value generated)
. gen y_alter = sin((alter-1)/`n'*2*_pi)*(1+mod(alter,2)*.2)
(1 missing value generated)
.
.
. bys ego (alter) : gen first = _n == 1
.
. forvalues t = 1/10{
2. twoway pcspike y_ego x_ego y_alter x_alter, ///
> lcolor(gs8) || ///
> scatter y_ego x_ego if first & state`t'==1 , ///
> mcolor(yellow) msymbol(O) || ///
> scatter y_ego x_ego if first & state`t'==2 , ///
> mcolor(red) msymbol(O) || ///
> scatter y_ego x_ego if first & state`t'==3 , ///
> mcolor(black) msymbol(Oh) ///
> aspect(1) xscale(off) yscale(off) ///
> legend(order(2 "susceptible" ///
> 3 "infected" ///
> 4 "recovered") rows(1)) ///
> name(t`t', replace) title(Time `t')
3. }
. drop _all
-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------
Multiple towns on a grid
The basic SIR model could represent how a disease spreads in a specific
town
we could place multiple SIR models on a grid (like a chessboard) and see
how the disease spreads across this space
In the first round a number of people in a single town is infected
In all subsequent rounds
the disease spreads within each town as before, but then
some agents visit neighbouring town and potentially spread the
disease there
. set seed 123456789
. clear mata
. drop _all
. run grid_main.do
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: model = sir_grid()
: model.rdim(10)
: model.cdim(10)
: model.prop_close(0.15)
: model.prop_far(0.01)
:
: model.N(100)
: model.tdim(50)
: model.outbreak(5)
: model.mcontacts(4)
: model.transmissibility(0.045)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0)
:
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
: model.export_sir()
:
: end
-------------------------------------------------------------------------------
.
. replace r = -r
(5,000 real changes made)
. foreach t of numlist 1 5(5) 50 {
2. heatplot I r c if t==`t', ///
> discrete title(Time = `t') ///
> name(grid`t', replace) ///
> cuts(0(0.05)0.7) ///
> colors(plasma) ///
> yscale(off range(-10.5 -.5)) ///
> xscale(off range(.5 10.5)) ///
> ylab(none) xlab(none) ///
> plotregion(margin(zero)) ///
> aspect(1)
3. }
-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------
The ABM package
The Agent Based Models implemented in Mata as a >> class
They use a collection of helper classes available in the ABM package.
The main classes availbe in the abm package are:
o abm_pop can store a list of agents and their properties. It is
particulary helpful when at each point in time most agents'
properties are stable, but some change.
o abm_nw for storing and managing a network
o abm_grid for storing and managing a grid (like a chessboard) on
which agents can "live" and interact with one another.
o abm_util useful functions that can be imported into the model class
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------
Basic SIR model
Lets take a look at the code for the basic SIR model:
version 16.1
mata:
mata set matastrict on
class sir extends abm_util {
class abm_pop scalar agents
real scalar tdim
real scalar outbreak
real scalar removed
real scalar mcontacts
real scalar N
real scalar transmissibility
real scalar mindur
real scalar meandur
real scalar pr_heal
real scalar pr_loss
// sir_chks.do
void isid()
void istime()
// sir_set_pars.do
transmorphic N()
transmorphic tdim()
transmorphic outbreak()
transmorphic removed()
transmorphic mcontacts()
transmorphic transmissibility()
transmorphic mindur()
transmorphic meandur()
transmorphic pr_loss()
// sir_sim.do
void setup()
real matrix meet()
real scalar infect()
void progress()
void step()
void run()
// sir_export.do
void export_sir()
void export_r()
}
end
do sir_chks.do
do sir_set_pars.do
do sir_sim.do
do sir_export.do
exit
The model is implemented as a class sir
The sir class inherrits (extends) from abm_util
This gives us access to several argument checking functions and the
dots() function
The population is stored is stored in an instance of the class abm_pop.
Whatever is stored in abm_pop is assumed to persist until you store
something else. So you only have to store the disease state of an
agent when it changes. This can save a lot of reading and writing and
make the simulation run quicker.
There are various checks perfomed and we can see those
include sir_locals.do
mata:
void sir::isid(real scalar id)
{
is_int_inrange(id,1,N)
}
void sir::istime(real scalar t)
{
is_int_inrange(t,0,tdim)
}
end
The functions that set the parameters are shown
include sir_locals.do
mata:
transmorphic sir::N(| real scalar val)
{
if (args() == 1) {
is_posint(val)
N = val
}
else {
return(N)
}
}
transmorphic sir::tdim(| real scalar val)
{
if (args() == 1) {
is_posint(val)
tdim= val
}
else {
return(tdim)
}
}
transmorphic sir::outbreak(| real scalar val)
{
if (args() == 1) {
is_posint(val)
outbreak= val
}
else {
return(outbreak)
}
}
transmorphic sir::removed(| real scalar val)
{
if (args() == 1) {
is_posint(val)
removed= val
}
else {
return(removed)
}
}
transmorphic sir::mcontacts(| real scalar val)
{
if ( args()==1 ) {
if ( val < 0) _error("arguments must be positive")
mcontacts = val
}
else {
return(mcontacts)
}
}
transmorphic sir::transmissibility(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
transmissibility = val
}
else {
return(transmissibility)
}
}
transmorphic sir::mindur(| real scalar val)
{
if ( args()==1 ) {
is_posint(val)
mindur = val
}
else {
return(mindur)
}
}
transmorphic sir::meandur(| real scalar val)
{
if ( args()==1 ) {
if (val < 0) _error("argument must be positive")
meandur = val
}
else {
return(meandur)
}
}
transmorphic sir::pr_loss(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
pr_loss = val
}
else {
return(pr_loss)
}
}
end
The main simulation functions are
include sir_locals.do
mata:
void sir::run()
{
real scalar i, t
setup()
dots(0)
dots(1)
for (t=2; t<=tdim; t++) {
dots(t)
for(i=1; i<=N ; i++) {
step(i,t)
}
}
}
void sir::setup()
{
real vector id, key
real scalar i
if (N==.) _error("N has not been set")
if (tdim==.) _error("tdim has not been set")
if (outbreak==.) _error("outbreak has not been set")
if (outbreak > N) _error("the initial outbreak is larger than population")
if (removed==.) removed = 0
if (removed + outbreak > N) _error("intial outbreak + removed agents is larger than the population")
if (mcontacts==.) _error("mcontacts has not been set")
if (mcontacts > N) _error("average number of mcontacts is larger than population")
if (transmissibility==.) _error("transmissibility has not been set")
if (mindur==.) _error("mindur has not been specified")
if (meandur==.) _error("meandur has not been specified")
if (meandur <= mindur) _error("meandur must be larger than mindur")
if (pr_loss==.) pr_loss = 0
agents.abm_version("0.1.5")
agents.N(N)
agents.k(3)
agents.setup()
id = jumble((1..N)')
for (i = 1; i<= N; i++) {
key = id[i],`state',1
agents.put(key,(i<=outbreak ? `infectious' : (i <= outbreak+removed ? `removed' : `susceptible')))
key[2] = `exposes'
if (i <= outbreak) agents.put(key, J(2,0,.) )
key[2] = `dur'
agents.put(key,1)
}
// if meandur - mindur < 1 then pr_heal > 1, so there is no randomness to the duration
// and the fixed duration of the disease is mindur
pr_heal = 1/(meandur-mindur)
}
void sir::step(real scalar id, real scalar t)
{
real matrix exposed
real scalar i, k
real rowvector key
isid(id)
istime(t)
key = id, `dur',t
agents.put(key, agents.get(key, "last")+1)
progress(id,t)
key = id, `state',t
if (agents.get(key, "last") == `infectious' ) {
exposed = meet(id, t)
for(i=1; i<= cols(exposed); i++) {
exposed[2,i] = infect(exposed[1,i],t)
}
key[2] = `exposes'
agents.put(key,exposed)
}
}
void sir::progress(real scalar id, real scalar t)
{
real rowvector key1, key2, key3
isid(id)
istime(t)
key1 = id, `state',t
key2 = id, `dur',t
if (agents.get(key1, "last") == `infectious' & agents.get(key2, "last") > mindur) {
if (runiform(1,1) < pr_heal) {
agents.put(key1, `removed')
agents.put(key2,1)
key3 = id, `exposes',t
agents.put(key3,NULL)
}
}
if (agents.get(key1, "last") == `removed' & runiform(1,1) < pr_loss) {
agents.put(key1, `susceptible')
agents.put(key2,1)
}
}
real matrix sir::meet(real scalar id, real scalar t)
{
real scalar k
real colvector res
real rowvector key
isid(id)
k = rpoisson(1,1,mcontacts)
if (k == 0) {
return(J(2,0,.))
}
else{
res = ceil(runiform(1,k):*(N-1))
res = res + (res:>=id)
res = res \ J(1,k,.)
return (res)
}
}
real scalar sir::infect(real scalar id, real scalar t)
{
real vector key
real scalar infected
isid(id)
istime(t)
infected = 0
key = id,`state', t
if (agents.get(key, "last")==`susceptible') {
if (runiform(1,1) < transmissibility) {
agents.put(key, `infectious')
key[2] = `dur'
agents.put(key, 1)
key[2] = `exposes'
agents.put(key,J(2,0,.))
infected = 1
}
}
return(infected)
}
end
The "variables" in agents are actually numbered and not named. To make it
easier to read the code I use local macros defined in
// address in abm_pop class agents
local state = 1
local dur = 2
local exposes = 3
// states
local susceptible = 1
local infectious = 2
local removed = 3
The functions that allow you to export the results to Stata are
include sir_locals.do
mata:
void sir::export_sir(| string rowvector names)
{
real matrix res
pointer matrix p
string rowvector varnames
real scalar xid, n, i, j
if (args() == 1) {
if (cols(names)!= 4) _error("4 variable names need to be specified")
varnames = names
}
else {
varnames = "t", "S", "I", "R"
}
res = J(N,tdim,.)
p = agents.extract(`state', tdim)
for (i=1; i<=N; i++) {
for(j=1; j<=tdim; j++) {
res[i,j] = *p[i,j]
}
}
res = (1..tdim)', mean(res:==`susceptible')',
mean(res:==`infectious')',
mean(res:==`removed')'
xid = st_addvar("float", varnames)
n = st_nobs()
n = rows(res) - n
if (n > 0) {
st_addobs(n)
}
st_store(.,xid, res)
}
void sir::export_r()
{
real matrix res
pointer matrix p_exposes, p_state
real scalar xid, n, i, j, t
string rowvector varnames
real colvector n_ill
t = .
res = (1..tdim)',J(tdim,1,0)
n_ill = J(tdim,1,0)
p_exposes = agents.extract(`exposes',tdim)
p_state = agents.extract(`state',tdim)
for (i=1; i<=N; i++) {
if (*p_state[i,1]==`infectious') {
t=1
}
else {
t=.
}
for(j=2;j<=tdim;j++) {
if (*p_state[i,j-1]==`susceptible' &
*p_state[i,j] ==`infectious') {
t=j
n_ill[t] = n_ill[t] + 1
}
if (*p_state[i,j-1]==`infectious' &
(*p_state[i,j] ==`removed' | *p_state[i,j] ==`susceptible')) t = .
if (t!=.) {
res[t,2] = res[t,2] + rowsum((*p_exposes[i,j])[2,.])
}
}
}
res[.,2] = res[.,2]:/n_ill
varnames = "t", "reprod"
xid = st_addvar("float", varnames)
n = st_nobs()
n = rows(res) - n
if (n > 0) {
st_addobs(n)
}
st_store(.,xid, res)
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------
SIR model in a network
To add the network we make use the abm_nw class.
This takes care of creating the network and keeps track of who can
contact who.
With that the changes required are fairly small
The class definition is
version 16.1
mata:
mata set matastrict on
class sir_nw extends abm_util {
class abm_pop scalar agents
class abm_nw scalar nw // <-- new
real scalar tdim
real scalar outbreak
real scalar removed
// real scalar mcontacts // <-- no longer necessary
real scalar N
real scalar transmissibility
real scalar mindur
real scalar meandur
real scalar pr_heal
real scalar pr_loss
real scalar degree // <-- new
real scalar pr_long // <-- new
// sir_chks.do
void isid()
void istime()
// sir_set_pars.do
transmorphic N()
transmorphic tdim()
transmorphic outbreak()
transmorphic removed()
// transmorphic mcontacts() // <-- no longer necessary
transmorphic transmissibility()
transmorphic mindur()
transmorphic meandur()
transmorphic pr_loss()
transmorphic degree() // <-- new
transmorphic pr_long() // <-- new
// sir_sim.do
void setup()
// real matrix meet() // <-- no longer necessary
real scalar infect()
void progress()
void step()
void run()
// sir_export.do
void export_sir()
void export_r()
void export_nw() // <-- new
}
end
do nw_chks.do
do nw_set_pars.do
do nw_sim.do
do nw_export.do
exit
Checks and setting parameters are pretty much the same
The main simulation functions are
include nw_locals.do
mata:
void sir_nw::run()
{
real scalar i, t
setup()
dots(0)
dots(1)
for (t=2; t<=tdim; t++) {
dots(t)
for(i=1; i<=N ; i++) {
step(i,t)
}
}
}
void sir_nw::setup()
{
real vector id, key
real scalar i
if (N==.) _error("N has not been set")
if (tdim==.) _error("tdim has not been set")
if (outbreak==.) _error("outbreak has not been set")
if (outbreak > N) _error("the initial outbreak is larger than population")
if (removed==.) removed = 0
if (removed + outbreak > N) _error("intial outbreak + removed agents is larger than the population")
if (transmissibility==.) _error("transmissibility has not been set")
if (mindur==.) _error("mindur has not been specified")
if (meandur==.) _error("meandur has not been specified")
if (meandur <= mindur) _error("meandur must be larger than mindur")
if (pr_loss==.) pr_loss = 0
if (degree==.) _error("degree has not been specified") // <-- new
if (pr_long==.) _error("pr_long has not been specified") // <-- new
nw.abm_version("0.1.5") // <-- new
nw.clear() // <-- new
nw.N_nodes(0,N) // <-- new
nw.directed(0) // <-- new
nw.tdim(0) // <-- new
nw.sw(degree,pr_long) // <-- new
nw.setup() // <-- new
agents.abm_version("0.1.5")
agents.N(N)
agents.k(4)
agents.setup()
id = jumble((1..N)')
for (i = 1; i<= N; i++) {
key = id[i],`state',1
agents.put(key,(i<=outbreak ? `infectious' : (i <= outbreak+removed ? `removed' : `susceptible')))
key[2] = `exposes'
if (i <= outbreak) agents.put(key, J(2,0,.) )
key[2] = `dur'
agents.put(key,1)
}
// if meandur - mindur < 1 then pr_heal > 1, so there is no randomness to the duration
// and the fixed duration of the disease is mindur
pr_heal = 1/(meandur-mindur)
}
void sir_nw::step(real scalar id, real scalar t)
{
real matrix exposed
real scalar i, k
real rowvector key
isid(id)
istime(t)
key = id, `dur',t
agents.put(key, agents.get(key, "last")+1)
progress(id,t)
key = id, `state',t
if (agents.get(key, "last") == `infectious' ) {
exposed = nw.neighbours(id) // <-- new
exposed = exposed \ J(1,cols(exposed),.) // <-- new
for(i=1; i<= cols(exposed); i++) {
exposed[2,i] = infect(exposed[1,i],t)
}
key[2] = `exposes'
agents.put(key,exposed)
}
}
void sir_nw::progress(real scalar id, real scalar t)
{
real rowvector key1, key2, key3
isid(id)
istime(t)
key1 = id, `state',t
key2 = id, `dur',t
if (agents.get(key1, "last") == `infectious' & agents.get(key2, "last") > mindur) {
if (runiform(1,1) < pr_heal) {
agents.put(key1, `removed')
agents.put(key2,1)
key3 = id, `exposes',t
//agents.put(key3,NULL)
}
}
if (agents.get(key1, "last") == `removed' & runiform(1,1) < pr_loss) {
agents.put(key1, `susceptible')
agents.put(key2,1)
}
}
real scalar sir_nw::infect(real scalar id, real scalar t)
{
real vector key
real scalar infected
isid(id)
istime(t)
infected = 0
key = id,`state', t
if (agents.get(key, "last")==`susceptible') {
if (runiform(1,1) < transmissibility) {
agents.put(key, `infectious')
key[2] = `dur'
agents.put(key, 1)
key[2] = `exposes'
agents.put(key,J(2,0,.))
infected = 1
}
}
return(infected)
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------
SIR model on a grid
The sir class remains largely unchanged
The main trick here is that one can create a 10x10 matrix of instances of
a class, by specifying town=sir(10,10)
The abm_grid class is mainly used to find neighbouring towns.
version 16.1
mata:
mata set matastrict on
class sir_grid extends abm_util { // <-- new
class sir matrix town
class abm_grid scalar grid
real scalar rdim
real scalar cdim
real scalar prop_close
real scalar prop_medium
real scalar prop_far
real scalar tdim
real scalar outbreak
real scalar removed
real scalar mcontacts
real scalar N
real scalar transmissibility
real scalar mindur
real scalar meandur
real scalar pr_heal
real scalar pr_loss
// grid_set_pars.do
transmorphic N()
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic outbreak()
transmorphic removed()
transmorphic mcontacts()
transmorphic transmissibility()
transmorphic mindur()
transmorphic meandur()
transmorphic pr_loss()
transmorphic prop_close()
transmorphic prop_medium()
transmorphic prop_far()
// grid_chks.do
void is_valid_address()
// grid_sim.do
void setup()
void run()
void travel()
void travel_contact()
// grid_export
void export_sir()
}
class sir extends abm_util {
class abm_pop scalar agents
real scalar tdim
real scalar outbreak
real scalar removed
real scalar mcontacts
real scalar N
real scalar transmissibility
real scalar mindur
real scalar meandur
real scalar pr_heal
real scalar pr_loss
real scalar vil_id
// grid_chks.do
void isid()
void istime()
// gird_sir_set_pars.do
transmorphic N()
transmorphic tdim()
transmorphic outbreak()
transmorphic removed()
transmorphic mcontacts()
transmorphic transmissibility()
transmorphic mindur()
transmorphic meandur()
transmorphic pr_loss()
transmorphic vil_id() // <-- new
// grid_sir_sim.do
void setup()
real matrix meet()
real scalar infect()
void progress()
void step()
void run()
real scalar state() // <-- new
// grid_export.do
real matrix export_sir()
}
end
do grid_chks.do
do grid_set_pars.do
do grid_sir_set_pars.do
do grid_sir_sim.do
do grid_sim.do
do grid_export.do
exit
include grid_locals.do
mata:
transmorphic sir_grid::prop_close(| real scalar val)
{
if (args() == 1) {
is_pr(val)
prop_close = val
}
else {
return(prop_close)
}
}
transmorphic sir_grid::prop_medium(| real scalar val)
{
if (args() == 1) {
is_pr(val)
prop_medium = val
}
else {
return(prop_medium)
}
}
transmorphic sir_grid::prop_far(| real scalar val)
{
if (args() == 1) {
is_pr(val)
prop_far = val
}
else {
return(prop_far)
}
}
transmorphic sir_grid::rdim(| real scalar val)
{
if (args() == 1) {
is_posint(val)
rdim = val
}
else {
return(rdim)
}
}
transmorphic sir_grid::cdim(| real scalar val)
{
if (args() == 1) {
is_posint(val)
cdim = val
}
else {
return(cdim)
}
}
transmorphic sir_grid::N(| real matrix val)
{
if (args() == 1) {
is_posint(val)
N = val
}
else {
return(N)
}
}
transmorphic sir_grid::tdim(| real scalar val)
{
if (args() == 1) {
is_posint(val)
tdim= val
}
else {
return(tdim)
}
}
transmorphic sir_grid::outbreak(| real scalar val)
{
if (args() == 1) {
is_posint(val)
outbreak= val
}
else {
return(outbreak)
}
}
transmorphic sir_grid::removed(| real matrix val)
{
if (args() == 1) {
is_posint(val)
removed= val
}
else {
return(removed)
}
}
transmorphic sir_grid::mcontacts(| real matrix val)
{
if ( args()==1 ) {
if ( val < 0) _error("arguments must be positive")
mcontacts = val
}
else {
return(mcontacts)
}
}
transmorphic sir_grid::transmissibility(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
transmissibility = val
}
else {
return(transmissibility)
}
}
transmorphic sir_grid::mindur(| real scalar val)
{
if ( args()==1 ) {
is_posint(val)
mindur = val
}
else {
return(mindur)
}
}
transmorphic sir_grid::meandur(| real scalar val)
{
if ( args()==1 ) {
if (val < 0) _error("argument must be positive")
meandur = val
}
else {
return(meandur)
}
}
transmorphic sir_grid::pr_loss(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
pr_loss = val
}
else {
return(pr_loss)
}
}
end
include grid_locals.do
mata:
transmorphic sir::N(| real scalar val)
{
if (args() == 1) {
is_posint(val)
N = val
}
else {
return(N)
}
}
transmorphic sir::tdim(| real scalar val)
{
if (args() == 1) {
is_posint(val)
tdim= val
}
else {
return(tdim)
}
}
transmorphic sir::outbreak(| real scalar val)
{
if (args() == 1) {
is_posint(val, "zero_ok")
outbreak= val
}
else {
return(outbreak)
}
}
transmorphic sir::removed(| real scalar val)
{
if (args() == 1) {
is_posint(val, "zero_ok")
removed= val
}
else {
return(removed)
}
}
transmorphic sir::mcontacts(| real scalar val)
{
if ( args()==1 ) {
if ( val < 0) _error("arguments must be positive")
mcontacts = val
}
else {
return(mcontacts)
}
}
transmorphic sir::transmissibility(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
transmissibility = val
}
else {
return(transmissibility)
}
}
transmorphic sir::mindur(| real scalar val)
{
if ( args()==1 ) {
is_posint(val)
mindur = val
}
else {
return(mindur)
}
}
transmorphic sir::meandur(| real scalar val)
{
if ( args()==1 ) {
if (val < 0) _error("argument must be positive")
meandur = val
}
else {
return(meandur)
}
}
transmorphic sir::pr_loss(| real scalar val)
{
if ( args()==1 ) {
is_pr(val)
pr_loss = val
}
else {
return(pr_loss)
}
}
transmorphic sir::vil_id(| real scalar val)
{
if (args()==1) {
is_posint(val)
vil_id = val
}
else {
return(vil_id)
}
}
end
include grid_locals.do
mata:
void sir_grid::is_valid_address(real rowvector address)
{
if (cols(address) != 2) _error("an address ")
is_int_inrange(address[1],1,rdim)
is_int_inrange(address[2],1,cdim)
}
void sir::isid(real scalar id, real scalar from)
{
is_int_inrange(id,from,N)
}
void sir::istime(real scalar t)
{
is_int_inrange(t,0,tdim)
}
end
include grid_locals.do
mata:
void sir::setup()
{
real vector id, key
real scalar i
if (N==.) _error("N has not been set")
if (tdim==.) _error("tdim has not been set")
if (outbreak==.) _error("outbreak has not been set")
if (outbreak > N) _error("the initial outbreak is larger than population")
if (removed==.) removed = 0
if (removed + outbreak > N) _error("intial outbreak + removed agents is larger than the population")
if (mcontacts==.) _error("mcontacts has not been set")
if (mcontacts > N) _error("average number of mcontacts is larger than population")
if (transmissibility==.) _error("transmissibility has not been set")
if (mindur==.) _error("mindur has not been specified")
if (meandur==.) _error("meandur has not been specified")
if (meandur <= mindur) _error("meandur must be larger than mindur")
if (pr_loss==.) pr_loss = 0
if (vil_id==.) _error("village id has not been specified") // <-- new
agents.abm_version("0.1.5")
agents.N(N)
agents.k(3)
agents.setup()
id = jumble((1..N)')
for (i = 1; i<= N; i++) {
key = id[i],`state',1
agents.put(key,(i<=outbreak ? `infectious' : (i <= outbreak+removed ? `removed' : `susceptible')))
key[2] = `exposes'
if (i <= outbreak) agents.put(key, J(2,0,.) )
key[2] = `dur'
agents.put(key,1)
}
// if meandur - mindur < 1 then pr_heal > 1, so there is no randomness to the duration
// and the fixed duration of the disease is mindur
pr_heal = 1/(meandur-mindur)
}
void sir::run(real scalar t)
{
real scalar i
for (i=1; i<=N; i++) {
step(i,t)
}
}
void sir::step(real scalar id, real scalar t)
{
real matrix exposed
real scalar i
real rowvector key
isid(id, `internal') // <- internal to distinguish from inside and outside town
istime(t)
key = id, `dur',t
agents.put(key, agents.get(key, "last")+1)
progress(id,t)
key = id, `state',t
if (agents.get(key, "last") == `infectious' ) {
exposed = meet(id, vil_id) // <-- vil_id to keep track of town
for(i=1; i<= cols(exposed); i++) {
exposed[3,i] = infect(exposed[1,i],t) // <- exposure now needs three rows: town, id, infection
}
key[2] = `exposes'
agents.put(key,exposed)
}
}
void sir::progress(real scalar id, real scalar t)
{
real rowvector key1, key2, key3
isid(id, `internal')
istime(t)
key1 = id, `state',t
key2 = id, `dur',t
if (agents.get(key1, "last") == `infectious' & agents.get(key2, "last") > mindur) {
if (runiform(1,1) < pr_heal) {
agents.put(key1, `removed')
agents.put(key2,1)
key3 = id, `exposes',t
agents.put(key3,NULL)
}
}
if (agents.get(key1, "last") == `removed' & runiform(1,1) < pr_loss) {
agents.put(key1, `susceptible')
agents.put(key2,1)
}
}
real matrix sir::meet(real scalar id, real scalar from)
{
real scalar k
real colvector res
id = (from==vil_id ? id : 0)
isid(id, from==vil_id)
k = rpoisson(1,1,mcontacts)
if (k == 0) {
return(J(3,0,.))
}
else{
res = ceil(runiform(1,k):*(N-1))
res = res + (res:>=id)
res = res \ J(1,k,from) \ J(1,k,.)
return (res)
}
}
real scalar sir::infect(real scalar id, real scalar t)
{
real vector key
real scalar infected
isid(id, `internal')
istime(t)
infected = 0
key = id,`state', t
if (agents.get(key, "last")==`susceptible') {
if (runiform(1,1) < transmissibility) {
agents.put(key, `infectious')
key[2] = `dur'
agents.put(key, 1)
key[2] = `exposes'
agents.put(key,J(3,0,.))
infected = 1
}
}
return(infected)
}
real scalar sir::state(real scalar id, real scalar t) // <-- new
{
real vector key
key = id, `state', t
return(agents.get(key, "last"))
}
end
include grid_locals.do
mata:
void sir_grid::run()
{
real scalar t, r, c
setup()
dots(0)
dots(1)
for (t=2; t<=tdim; t++) {
for (r=1; r<=rdim; r++) {
for (c=1; c<=cdim; c++) {
town[r,c].run(t)
}
}
for (r=1; r<=rdim; r++) {
for (c=1; c<=cdim; c++) {
travel(r,c,t)
}
}
dots(t)
}
}
void sir_grid::setup()
{
real scalar r, c, j, initial
if (N==J(0,0,.)) _error("N has not been set")
if (tdim==.) _error("tdim has not been set")
if (outbreak==.) _error("outbreak has not been set")
if (removed==J(0,0,.)) removed = 0
if (mcontacts==J(0,0,.)) _error("mcontacts has not been set")
if (transmissibility==.) _error("transmissibility has not been set")
if (mindur==.) _error("mindur has not been specified")
if (meandur==.) _error("meandur has not been specified")
if (pr_loss==.) pr_loss = 0
if (prop_close==.) _error("N_close has not been specified")
if (prop_medium==.) prop_medium = 0
if (prop_far==.) prop_far = 0
if (rows(N) == 1 & cols(N) == 1){
N = J(rdim, cdim, N)
}
else if (rows(N)!=rdim | cols(N)!= cdim) {
_error("N can either be a scalar or a rdim x cdim matrix")
}
if (rows(removed) == 1 & cols(removed) == 1 ) {
removed = J(rdim,cdim,removed)
}
else if (rows(removed)!= rdim | cols(removed)!= cdim) {
_error("removed can either be a scalar or a rdim x cdim matrix")
}
if (rows(mcontacts) == 1 & cols(mcontacts)==1) {
mcontacts = J(rdim,cdim,mcontacts)
}
else if (rows(mcontacts)!=rdim | cols(mcontacts)!= cdim) {
_error("mcontacts can either be a scalar or rdim x cdim matrix")
}
grid.abm_version("0.1.5")
grid.rdim(rdim)
grid.cdim(cdim)
grid.neumann(0)
grid.torus(0)
grid.setup()
town = sir(rdim,cdim)
initial = ceil(runiform(1,1)*rdim*cdim)
j= 0
for(r=1; r<=rdim; r++) {
for(c=1; c<=cdim; c++) {
j = ++j
town[r,c].N(N[r,c])
town[r,c].tdim(tdim)
town[r,c].outbreak((j==initial)*outbreak)
town[r,c].removed(removed[r,c])
town[r,c].mcontacts(mcontacts[r,c])
town[r,c].transmissibility(transmissibility)
town[r,c].mindur(mindur)
town[r,c].meandur(meandur)
town[r,c].pr_loss(pr_loss)
town[r,c].vil_id(j)
town[r,c].setup()
}
}
}
void sir_grid::travel(real scalar r, real scalar c, real scalar t)
{
travel_contact(r,c,t,"close")
travel_contact(r,c,t,"medium")
travel_contact(r,c,t,"far")
}
void sir_grid::travel_contact(real scalar r, real scalar c, real scalar t, string scalar where)
{
real scalar i, j, k, tr, tc, N_travel, target
real vector traveling, key
real matrix neigh, exposed
if (where == "close") {
N_travel = ceil(prop_close*N[r,c])
traveling = ceil(runiform(1,N_travel)*N[r,c])
neigh = grid.find_ring(r,c,1)
}
else if (where == "medium") {
N_travel = ceil(prop_medium*N[r,c])
traveling = ceil(runiform(1,N_travel)*N[r,c])
neigh = grid.find_ring(r,c,2)
}
else {
N_travel = ceil(prop_far*N[r,c])
traveling = ceil(runiform(1,N_travel)*N[r,c])
neigh = grid.schedule()
}
k = rows(neigh)
for(i=1; i<=N_travel; i++) {
if (town[r,c].state(traveling[i],t)== `infectious') {
target = ceil(runiform(1,1)*k)
tr = neigh[target,1]
tc = neigh[target,2]
exposed = town[tr, tc].meet(traveling[i],town[r,c].vil_id)
for(j=1; j<= cols(exposed); j++) {
exposed[3,j] = town[tr,tc].infect(exposed[1,j],t)
}
key = traveling[i],`exposes', t
exposed = town[r,c].agents.get(key), exposed
town[r,c].agents.put(key,exposed)
}
}
}
end
// address in abm_pop class agents
local state = 1
local dur = 2
local exposes = 3
// states
local susceptible = 1
local infectious = 2
local removed = 3
// from
local internal = 1
local external = 0
include grid_locals.do
mata:
real matrix sir::export_sir()
{
real matrix res
pointer matrix p
real scalar i, j
res = J(N,tdim,.)
p = agents.extract(`state', tdim)
for (i=1; i<=N; i++) {
for(j=1; j<=tdim; j++) {
res[i,j] = *p[i,j]
}
}
res = (1..tdim)', mean(res:==`susceptible')',
mean(res:==`infectious')',
mean(res:==`removed')'
return(res)
}
void sir_grid::export_sir()
{
real scalar r, c, i, xid, n
real matrix res
string rowvector varnames
res = J(rdim*cdim*tdim, 6,.)
i = 0
for(r=1; r<=rdim; r++) {
for(c=1; c<=cdim; c++) {
res[|1+tdim*i++ ,1 \ (i+1)*tdim,6|] = J(tdim,1,(r,c)),town[r,c].export_sir()
}
}
varnames = "r", "c", "t", "S", "I", "R"
xid = st_addvar("float", varnames)
n = st_nobs()
n = rows(res) - n
if (n > 0) {
st_addobs(n)
}
st_store(.,xid, res)
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Final remarks
-------------------------------------------------------------------------------
Conclusion
Implementing an ABM is doable in Mata
We have seen four "helper classes" from abm:
abm_pop for storing a population of agents
abm_nw for managing a network
abm_grid for managing a grid, like a chessboard.
abm_util for importing useful functions in the model class
-------------------------------------------------------------------------------
<< index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
digression
-------------------------------------------------------------------------------
R0
14 days is a number that has floated around a lot for how long someone is
infectious when having COVID-19.
On average 4 contacts is just a convenient number
The transmissibility is then chosen to make the R0 close to 2.5, another
number that has floated around a lot for COVID-19.
If someone is 14 days infectious, and contacts 4 persons each day,
then there are 14*4=56 opportunities to infect someone.
These should result in 2.5 infections, so the chance of a contact
leading to an infection is 2.5/56 = 0.045
-------------------------------------------------------------------------------
<< index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ancillary
-------------------------------------------------------------------------------
Other Scenarios
What if the immunity against COVID-19 obtained from getting the disease
does not persist?
We can change to model to see how that impacts the spread of the disease
. clear mata
. drop _all
. run sir_main.do
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: model = sir()
: model.N(2000)
: model.tdim(250)
: model.outbreak(5)
: model.mcontacts(4)
: model.transmissibility(0.045)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0.02)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
.................................................. 200
.................................................. 250
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(anc1, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
. drop _all
Similarly, we can study the potential impact of policies encouraging
social distancing or the wearing of masks by changing the
transmissibility
. mata:
------------------------------------------------- mata (type end to exit) -----
: model.pr_loss(0)
: model.tdim(150)
: model.transmissibility(0.03)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(anc2, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
. drop _all
or the potential impact of a lockdown, which would influence the number
of people someone is in contact with
. mata:
------------------------------------------------- mata (type end to exit) -----
: model.transmissibility(0.045)
: model.mcontacts(3)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(anc3, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
. drop _all
or the impact of vaccinating a quarter of the population
. mata:
------------------------------------------------- mata (type end to exit) -----
: model.mcontacts(4)
: model.removed(500)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(anc4, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
. drop _all
-------------------------------------------------------------------------------
index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ancillary
-------------------------------------------------------------------------------
Quantifying uncertainty
ABMs are simulations, so the outcome could easily differ from run to run
One could display the outcome from multiple runs
Either because this gives a more complete description of an ABM
or because this variation is interesting in and of itself
. set seed 123456789
. clear mata
. drop _all
. run sir_main.do
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: model = sir()
: model.N(2000)
: model.tdim(150)
: model.outbreak(5)
: model.mcontacts(4)
: model.transmissibility(0.045)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0)
:
: base = "t", "S", "I", "R"
: for(i=1; i<=10; i++){
> model.run()
> names = J(1,4,"")
> for(j=1; j<=4; j++) {
> names[j] = base[j] + strofreal(i)
> }
> model.export_sir(names)
> }
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: end
-------------------------------------------------------------------------------
.
. twoway line S* t1, lcolor(green ..) || ///
> line I* t1 , lcolor(orange ..) || ///
> line R* t1, lcolor(blue ..) ///
> legend(order(1 "S" 11 "I" 21 "R") cols(3)) ///
> ylab(,angle(0)) name(anc5, replace)
-------------------------------------------------------------------------------
index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ancillary
-------------------------------------------------------------------------------
Other scenarios
We can try to quantify the impact of reducing the number of long distance
by comparing the following models:
. clear mata
. set seed 28863
. run nw_main.do
. mata:
------------------------------------------------- mata (type end to exit) -----
: model = sir_nw()
: model.N(2000)
: model.tdim(150)
: model.outbreak(10)
: model.degree(10)
: model.pr_long(.2)
: model.transmissibility(.018)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0.00)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(sir3, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
. drop _all
.
. mata:
------------------------------------------------- mata (type end to exit) -----
: model.pr_long(.05)
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
: model.export_sir()
: end
-------------------------------------------------------------------------------
.
. twoway line S I R t, name(sir4, replace) ///
> ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..)
.
. drop _all
-------------------------------------------------------------------------------
index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ancillary
-------------------------------------------------------------------------------
Other scenarios
In the real world population is not uniformly distributed over space
We could try to capture that by including a "big city" in our model.
This can be done by increasing the population in one of the towns.
. set seed 123456789
. clear mata
. drop _all
. run grid_main.do
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: model = sir_grid()
: model.rdim(10)
: model.cdim(10)
: model.prop_close(0.15)
: model.prop_far(0.01)
:
: N = J(10,10,100)
: N[5,5] =2000
: model.N(N)
: model.tdim(100)
: model.outbreak(5)
: model.mcontacts(4)
: model.transmissibility(0.045)
: model.mindur(10)
: model.meandur(14)
: model.pr_loss(0)
:
: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
: model.export_sir()
:
: end
-------------------------------------------------------------------------------
.
. replace r = -r
(10,000 real changes made)
. foreach t of numlist 1 5(5) 100 {
2. heatplot I r c if t==`t', ///
> discrete title(Time = `t') ///
> name(ancgrid`t', replace) ///
> cuts(0(0.05)0.7) ///
> colors(plasma) ///
> yscale(off range(-10.5 -.5)) ///
> xscale(off range(.5 10.5)) ///
> ylab(none) xlab(none) ///
> plotregion(margin(zero)) ///
> aspect(1)
3. }
-------------------------------------------------------------------------------
digression
-------------------------------------------------------------------------------
Class
We have a population of agents to keep track of
We need to repeatedly apply rules to them
Applying rules can be done by a function, but that function needs access
to settings and the population.
This is a good case for a class.
A class can be thought of as a tray with named slots for functions and
data
The functions are local to the class
Each function has access to all material stored in the class
. clear mata
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: class arithmatic {
>
> real scalar a
> real scalar b
>
> real scalar sum()
> real scalar prod()
> void input()
> void run()
> }
:
: void arithmatic::input(real scalar aval, real scalar bval)
> {
> a = aval
> b = bval
> }
:
: real scalar arithmatic::sum()
> {
> return(a + b)
> }
:
: real scalar arithmatic::prod()
> {
> return(a * b)
> }
:
: void arithmatic::run()
> {
> sum()
> prod()
> }
:
: // how the user interacts with this program
: math = arithmatic()
: math.input(2,4)
: math.run()
6
8
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
ancillary
-------------------------------------------------------------------------------
Alternative platforms
There are various dedicated environments for creating ABMs.
netlogo (netlogo) http://ccl.northwestern.edu/netlogo/
repast (java, python, C++, and others) https://repast.github.io/
mason (java) https://cs.gmu.edu/~eclab/projects/mason/
mesa (python) https://mesa.readthedocs.io/
agents.jl (julia) https://juliadynamics.github.io/Agents.jl/stable/
People who primarily create ABMs are not going to move away from those,
and rightly so.
The toolkits discussed in this presentation are aimed at people who
primarily use Stata and occationally want to make an ABM.
-------------------------------------------------------------------------------
index
-------------------------------------------------------------------------------
