Agent based models in Mata: Modelling aggregate processes, like the spread of a > disease

Maarten Buis University of Konstanz
maarten.buis@uni.kn -------------------------------------------------------------------------------

index >>

-------------------------------------------------------------------------------

Table of Contents
    What is an Agent Based Model?

        What is an Agent Based Model?

        Example: The spread of a disease

            R0

        Other Scenarios(ancillary)

        Quantifying uncertainty(ancillary)

        Adding a network

        Other scenarios(ancillary)
 
    How to implement an ABM in Mata

        Basic SIR model

            Class

        SIR model in a network
 
    Final remarks

        Conclusion

        Example of an ABM on a grid: Schelling's segregation model
       (ancillary)

        Alternative platforms
 
 
-------------------------------------------------------------------------------

<< 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
. 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 "infected" /// > 3 "susceptible" /// > 4 "recovered") rows(1)) /// > name(t`t', replace) title(Time `t') 3. }
. drop _all
    Close knit networks tend to slow down the spread of a disease, while the
    long distance ties tend

So that is one way in which advise against unnecesary travel or "social buble" rules in countries like Belgium works
-------------------------------------------------------------------------------

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

run abm_pop.mata

mata:

mata set matastrict on

class sir {
	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                          posint()
	void                          pr()
	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()
	void                          dots()
	
	// 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
The population is stored is stored in another class abm_pop, which is available from https://github.com/maartenteaches/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::posint(real scalar val)
{
	if (val <= 0 | val != ceil(val)) {
		_error("argument must be a positive integer")
	} 
}

void sir::pr(real scalar val)
{
	if (val <0 | val > 1) {
		_error("argument must be between 0 and 1")
	}
}

void sir::isid(real scalar id)
{
	if (id < 0 | id > N | id!=ceil(id)) {
		_error("argument is not a valid id")
	}
	
}

void sir::istime(real scalar t)
{
	if (t<0 | t > tdim | t!=ceil(t)) {
		_error("argument is not a valid time")
	}
}
end

The functions that set the parameters are shown
include sir_locals.do

mata:

transmorphic sir::N(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		N = val
	}
	else {
		return(N)
	}
}

transmorphic sir::tdim(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		tdim= val
	}
	else {
		return(tdim)
	}
}


transmorphic sir::outbreak(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		outbreak= val
	}
	else {
		return(outbreak)
	}
}

transmorphic sir::removed(| real scalar val)
{
	if (args() == 1) {
		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 ) {
		pr(val)
		transmissibility = val
	}
	else {
		return(transmissibility)
	}
}

transmorphic sir::mindur(| real scalar val)
{
	if ( args()==1 ) {
		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 ) {
		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.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)
}



void sir::dots(real scalar t) 
{
	if (t==0) {
		printf("----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5\n")
	}
	else if (mod(t, 50) == 0) {
		printf(". %9.0f\n", t)
	}
	else {
		printf(".")
	}
	displayflush()
}
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, which is available from https://github.com/maartenteaches/abm_nw.
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

run abm_pop.mata
run abm_nw.mata                                          // <-- new

mata:

mata set matastrict on

class sir_nw {
	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                          posint()
	void                          pr()
	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()
	void                          dots()
	
	// 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.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.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)
}

void sir_nw::dots(real scalar t) 
{
	if (t==0) {
		printf("----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5\n")
	}
	else if (mod(t, 50) == 0) {
		printf(". %9.0f\n", t)
	}
	else {
		printf(".")
	}
	displayflush()
}
end
 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Final remarks
-------------------------------------------------------------------------------

Conclusion
Implementing an ABM is doable in Mata
We have seen two "helper classes":
abm_pop for storing a population of agents abm_nw for managing a network
There is also the abm_grid class to manage a grid, like a chessboard, which is available from https://github.com/maartenteaches/abm_grid.
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Final remarks
-------------------------------------------------------------------------------

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/
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

-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
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
Recently, there have an increasing number of news stories about that immunity against COVID-19 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
 
 
-------------------------------------------------------------------------------

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(anc4, 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

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

Example of an ABM on a grid: Schelling's segregation model
In this model there are blue and red agents, that live in cells on a square grid. Each cell can contain only one agent, but not all cells are occupied. An agent is happy if the proportion of agents of the same colour in her neighbourhood is above a certain cut-off value. If she is unhappy she will move to the nearest available cell that does satisfy her needs.
The point of this simulation that even with a low cut-off value, like 30%, you will get fairly seggregated neighbourhoods.

. clear mata
. drop _all
. set seed 1234567
. . run schelling_main.mata
. . mata: ------------------------------------------------- mata (type end to exit) ----- : model = schelling()
: model.rdim(10)
: model.cdim(10)
: model.tdim(5)
: model.nred(45)
: model.nblue(45)
: model.limit(.3)
: model.run()
: model.summtable()
Average proportion of neighbours with same color and average proportion of happy agents
----------------------------------- Iteration prop same prop happy ----------------------------------- 0 .543 .867 1 .688 .944 2 .725 .978 3 .738 .978 4 .768 .989 5 .768 .989 -----------------------------------
: model.extract(("r","c","t", "id", "red", "happy"))
: end -------------------------------------------------------------------------------
. . replace r = -r (540 real changes made)
. . forvalues t = 0/5 { 2. twoway scatter r c if red == 0 & t == `t' & happy == 1, /// > msymbol(O) mcolor(blue) msize(*2) /// > mlabel(id) mlabpos(0) || /// > scatter r c if red == 1 & t == `t' & happy == 1, /// > msymbol(O) mcolor(red) msize(*2) /// > mlabel(id) mlabpos(0) || /// > scatter r c if red == 0 & t == `t' & happy == 0, /// > msymbol(Oh) mcolor(blue) msize(*2) /// > mlabel(id) mlabpos(0) || /// > scatter r c if red == 1 & t == `t' & happy == 0, /// > msymbol(Oh) mcolor(red) msize(*2) /// > mlabel(id) mlabpos(0) /// > yscale(off range(-10.5 -.5)) /// > xscale(off range(.5 10.5)) /// > ylab(none) xlab(none) /// > plotregion(margin(zero)) /// > aspect(1) /// > xline(.5(1)9.5) /// > yline(-.5(-1)-9.5) /// > legend(off) /// > name(gr`t', replace) title(Iteration `t') 3. }
 
 
-------------------------------------------------------------------------------

index

-------------------------------------------------------------------------------