Validation and Inference of Agent Based Models

In the above model there exists a set of parameters which achieve equilibrium in the sheep and wolf populations. This is often found through trial and error. Another, more sophisticated method involves looping through many thousands of simulations, changing the set of parameters each time. The sheep and wolf populations are analysed in each run. The output is some metric describing how close the populations are to sustaining over time. This is demonstrated below. In each run, the parameters are adjusted, the model is run for a set number of time ticks, and the final numbers of both wolf and sheep compared to the initial populations. The data can be analysed with appropriate statistical models to identify stationarity. One such method is to view the two populations as two time series. The oscillation of the species populations can be viewed as seasonality. Given enough simulation time, the difference in population numbers from the start to the end of the time series is compared to determine if the system has reached equilibrium.

The BehaviorSpace tool in ABM software NetLogo allows a simulation to be run multiple times with changing parameters. A simple demonstration of this tool for analysis is when the ‘energy from grass’ parameter is changed. This parameter determines the energy that sheep gather from eating grass. Because sheep lose energy from movement, a low amount of energy gained from grass is likely to result in a larger decrease in the sheep population than if a large amount of energy was gained from eating grass. The Wolf-Sheep predation ABM allows this theory to be tested in a simulation based manner, followed by the application of robust statistical analysis methods from the output data to answer questions about the model. Keeping all other parameters fixed, the grass energy parameter was updated. The simulation was run for 1000 ticks and the population of the two species recorded at each tick.

When the grass energy is low, both species die out extremely quickly. This is due to the sheep not gaining enough energy to move, subsequently dying out. The wolves, now with no species to predate, also die out. At a grass energy of four, an equilibrium appears to be reached. The sheep can gain enough energy to continually move, which sustains the populations of the wolves predating them. If the grass energy is changed to 7, the sheep population varies far more. The population increases sharply due to the high energy levels of the sheep (who then reproduce more frequently), which results in a high level of competition for the available grass. Due to this competition, the population drops off to an extremely low level at $\sim$340 ticks. This results in the wolf population dropping to zero from the lack of prey. With no predators remaining, the sheep population increases to a high level and remains at $\sim$500 for the remainder of the simulation. In the last run the grass energy is set to 10. This has a similar effect as the previous parameter, although the sheep population reaches a higher level before remaining constant.

While this ‘brute-force’ approach is often used to analyse complex systems for the purposes of parameter inference, it is computationally inefficient. When looping through samples of each parameter, much of the simulation time is spent producing data far from the desired outcome. For a set of parameters $\theta$ of dimension $P$, each sampled at N equal intervals $a\times\frac{max(\theta)}{N},a=\{1,...N\}$, the required number of simulations quickly becomes prohibitively large. Consider the case of the wolf-sheep predation model with 7 parameters. If we would like to sample each parameter at 100 even intervals, the number of simulations to cover all possible combinations becomes:

Clearly this number of simulations is unreasonable. This can be shortened by reducing the number of sampled intervals for each parameter or disregarding some combinations of variables, but the issue remains of a large amount of computation time being wasted sampling parameters which do not produce simulated data close to the desired outcome of the ABM.

While ABM is a ’bottom up’ approach focusing on the aggregate behaviour of agents, Equation Based Modelling (EBM) is a top down approach. EBM is based on an interrelation between a series of partial differential equations (Parunak, Savit & Riolo, 1998). The variability in the system is examined over time and/or space. Validation is enabled by comparing the output of the model to the real world system behaviour. In comparison, ABM allows validation of both the individual and aggregate system behaviour. EBM is suited for modelling physical processes where the behaviour of the components is not as important as that of the overall system.

Generally, EBM’s are a set of differential equations and parameters. Extending the Lotka–Volterra model, time can be eliminated from each equation to produce a single equation, and then solved for a constant quantity $V$:

This model will give the number of two populations $x$ and $y$ given the set of parameters $\alpha$, $\beta$, $\gamma$ and $\delta$. A major limitation of this in ecological models is that neither population can reach zero - they always recover after falling to a low level. However, the more flexible structure of an ABM allows either population to go extinct.

While Agent Based Modelling is well suited for many classes of complex systems, it is better suited for some contexts than others. In problems with a large number of homogeneous agents, it is often the case that more accurate results can be obtained in a shorter amount of time using an aggregate solution. Tracking the behaviour of individual agents is not always necessary. For example, if we are concerned with the amount of water in a container, it is more useful to track the water level directly and model it using a set of differential equations. Parameters in the equations such as temperature can be adjusted and the effects observed. Where there are only a few interacting agents in a system, a similar approach can be taken by deriving a model with differential equations. One such example is in the behaviour of a tennis ball falling onto different surfaces. Here, the physical behaviours of the model are understood with a rigorous mathematical structure such as kinematic equations. This can be used to developed a set of differential equations for effectively describing the behaviour of the ball given different materials and surfaces.

Agent Based Modelling is well suited for systems with many heterogeneous agents interacting in an environment. In many populations, individual agents possess unique properties which are not shared across the entire population. Their behaviours are guided by both the environment, the history of their actions and their current properties. ABM is well suited for a system where it is useful to track individual agents and their actions. In this manner it is much more powerful than using equation based modelling. Where the behaviour of agents relies on their history of interactions, the agents can be modelled as intelligent beings with the ability to learn, reason about situations, and deal with uncertainty. An example of this is in the behaviour of pedestrians in a transportation network. Humans are intelligent beings, with each individual’s decision making process different from the rest. The reason for one individual choosing to travel a certain route is likely different to another individual, and the long term behaviours of each pedestrian differ substantially. While it would be impractical to model such a process using EBM, pedestrian modelling is well suited to an ABM approach, allowing for individual actions carried out by several heterogeneous agents.

Recent advances in the field of artificial intelligence have allowed for the construction of agents in an ABM, such that they are able to learn from their past and take an action based on a complex set of variables. One such influential paradigm to enable this is BDI - Belief-Desire-Intention (Singh, Padgham & Logan, 2016). BDI recognises “that an agent can be identified as having: a set of beliefs about its environment and about itself; a set of desires which are computational states which it wants to maintain, and a set of intentions which are computational states which the agent is trying to achieve” (O’Hare, Jennings, 1996). More recent developments have expanded on this philosophy, going beyond simple reactive behaviour to embedding cognitive complexity (Abrams, 2013). While many ABMs are not suited to this more complex agent intelligence, the recent advances in modelling highly complex systems allows for the embedding of more advanced levels of agent intelligence to achieve a more realistic simulation.

Although ABM is much less mathematically oriented than EBM, it is important to establish a mathematical framework for representing the system and the agent properties. There exist a number of mathematical formalisations for Agent Based Models. A series of mathematical objects can be developed to represent a theoretical abstraction of agent based simulations. These objects capture the features of simulations and allow for the development of mathematical theory to explain the output of a simulation (Laubenbacher et al, 2007). An ABM can be described as a Markov Chain where the state of the system at time $t$ is given by the state at all nodes in the system: $X_{t}={x_{i,t}},\;i=1,...,N,t=1,...,T$. The state at a given node $i$ for time $t+1$ is given by

That is, the state at time $t+1$ in the system is reliant only on the previous state. This is an important feature of an ABM, allowing for the model at each tick of the simulation to make agent decisions and update the environment based only on the current state. Storing every prior state and set of decisions, while updating the current state using this prior data, becomes costly in terms of computation and memory as a system increases in complexity.

Perhaps the most common framework for describing agent based simulations is Cellular Automata (CA). They are typically defined over a regular grid such as the two dimensional $\mathbb{Z}^{2}$. Each grid point $(i,j)$ represents a site or node. Each node has a state $x_{i,j}(t)$, where $t$ denotes a time step. The neighborhood N for each site consists of a set of nodes that can influence the future state of the node $x_{i,j}$. Based on the current state $x_{i,j}(t)$ and the current state of all nodes in the neighborhood $N$, a function $f(i,j)$ computes the next state $x_{i,j}(t+1)$ of the site at $(i,j)$. Laubenbacher et al (2007) defines this as:

where $\bar{x}_{i,j}$ denotes the tuples consisting of all the states $x_{i^{\prime},j^{\prime}}(t)$ with $(i^{\prime},j^{\prime})\in N$.

Another modelling framework is that of finite dynamical systems. This framework represents an ABM as a time-discrete dynamical system of a finite state set. It can be seen that the representation of an ABM in this way is mathematically rich enough to allow the derivation of formal results.

As the foundation of Agent Based Modelling lies in simulation based on observed data in addition to data output, statistics plays a major role at nearly every stage of the modelling process. Perhaps the most important step in ABM based research lies in the analysis of data from the simulation to derive parameters for understanding the data generating process. Other analysis may include descriptive, predictive, or inferential statistics. These models can include, but are not limited to; spatial, network, survival, or time series analysis. Most ABM frameworks provide a method to import and export data into more robust statistical packages such as R or Python. These packages provide advanced modelling tools which allow the researcher to extract insights from the ABM.

To date, there is no standard format for the documentation of an ABM. Limited replicability is an issue in scientific research using ABMs. To maximise the reproducibility of an ABM, a format for their design should be adopted which is popular in the scientific community. By standardising the description of a model in this way, it can be more easily reproduced by another party or revisited by the original developer. One such format is ODD - Overview, Design, Details, first proposed by Grimm et al (2006) as a standard protocol for describing agent based models. It was subsequently reviewed and updated in 2010 (Grimm et al, 2010). ODD describes a set of key aspects of the model:

1. 1.

   Purpose
   What is the purpose of the model? What are the desired outcomes?

2. 2.

   Entities, State Variables and Scales
   What entities exist in the model? Agents, environment, collectives, spatial units. How are the entities characterised? What properties do they hold? What are the spatial and temporal resolutions of the model?

3. 3.

   Process overview and scheduling
   What actions are assigned to each entity, and in which order? How are the state variables updated? How is time modelled - discretely, continuously, or as a continuum where both can occur?

4. 4.

   Design Concepts
   Emergence, Adaptation, Fitness, Prediction, Sensing, Interaction, Stochasticity, Collectives, Observation.

5. 5.

   Initialisation
   What is the initial state of the model world? How many entities are created and what are their default properties?

6. 6.

   Input Data
   Is external data used as input to represent the change in a process over time?

ABMs developed in the course of this dissertation will follow the ODD protocol. For each, thorough documentation is provided to enable reproduction of the model by researchers.

Several ABM toolkits are available to enable the development of agent based model simulations. Perhaps the most popular of these is NetLogo. NetLogo was developed in 1999 by Uri Wilesnky. It uses a philosophy of ‘low threshold, no ceiling’, with a goal of making the software easy to use for beginners while remaining flexible enough to develop advanced models. It contains a number of example models for fields such as network analysis, social science, and biology. All ABM models in this report are developed in NetLogo. The RNetLogo and pyNetLogo libraries in R and Python respectively allow for programmatic, headless control of a NetLogo model. Statistical models were developed in both R and Python, with the simulation element communicating with the NetLogo model for sending commands and receiving output.

Bayesian inference is a method of statistical inference grounded in Bayes Rule. It states that the probability of some event A, given that event B has occurred, is given by:

In the context of Bayesian parameter estimation, we aim to find the posterior distribution of the parameters given the data:

where $p(\theta|x)$ is the posterior distribution,
$p(x|\theta)$ is the likelihood,
$p(\theta)$ is the prior, and
$p(x)$ is the probability of the data. Often the denominator $p(x)$ involves computing an integral to normalise the posterior to a probability distribution:

In some cases, such as models where a conjugate prior is available, this is straightforward to calculate numerically. If no conjugate prior is available, there are a wide range of computational methods for sampling from the posterior. One of these is Markov chain Monte Carlo (MCMC). MCMC constructs a Markov chain with the target posterior dsitibution as its equilibrium distribution. By sampling from this chain for a long period of time, a sample from the posterior is obtained. A common algorithm for this is known as Metropolis Hastings, developed by Metropolis, Rosenbluth, Rosenbluth, Teller and Teller (1953). Formally, the algorithm uses a proposal distribution $g(x^{{}^{\prime}}|x)$ as the conditional probability of accepting a state $x^{{}^{\prime}}$ given $x$, and an acceptance ratio $A(x^{{}^{\prime}}|x)$ as the probability of accepting the proposed state $x^{{}^{\prime}}$.

MCMC can be used for parameter inference in any class of statistical model where the likelihood $g(x^{{}^{\prime}}|x)$ can be evaluated. Consider a simple generalised linear model of the form $\pi=\frac{1}{1+\exp[-(\beta_{0}+\beta_{1}x_{i})]}$. Here the parameters are the intercept $\beta_{0}$ and a single coefficient $\beta_{1}$. The likelihood follows a Binomial distribution of the form

The use of this likelihood in the calculation of the acceptance probability will construct a Markov chain that converges on the posterior distributions for $\beta_{0}$ and $\beta_{1}$.

For a large class of statistical models, the likelihood is either prohibitively expensive to compute or intractable to calculate. This arises often in complex systems where stochastic processes are involved and as a result, methods such as MCMC cannot be used. An alternative method of approximating the posterior distribution is required for these problems, known as likelihood free inference.

The likelihood in a statistical model is a function of the parameters given the data, $f(\theta|x)$. In simple models following a known distribution this is easily calculated. However, in complex models the integral often has no closed form solution allowing it to be evaluated in a finite number of operations. An example of this is the g and k distribution.

The g-and-k distribution is known as a quantile distribution, defined in terms of its inverse cumulative distribution function. It has five parameters - $\theta=\{a,b,c,g,k\}$, and its inverse cumulative distribution function is given by (Drovandi, Pettitt, 2011):

For a quantile distribution such as this, the likelihood cannot be calculated analytically, only numerically. In practice, simulating from the distribution to compute the likelihood is much more straightforward than the numerical calculation. We draw from the distribution by simulating values for substitution into the above function. This idea is the basis of likelihood free methods - the acceptance or rejection of samples in a simulation based model. There is ongoing development in these methods due to the advent of big data and the use of simulation models to evaluate parameters. Two such popular methods are known as Synthetic Likelihood and Approximate Bayesian Computation.

ABC aims to approximate a posterior distribution without the existence of the likelihood function $P(y|\theta)$. In essence this is achieved by simulating data with a set of prior parameters $\theta^{{}^{\prime}}$. The simulated data is then compared with the observed data in some capacity - either in totality or via a set of summary statistics. If the simulated data is deemed close enough to the observed data (using a defined distance measure), the set of sampled prior parameters belongs to the approximated posterior distribution $P(\theta|d(y^{*},y)\leq\epsilon)$. Formally this can be denoted as

where $\mathbbm{1}_{\Delta(x,x^{*})\leq\epsilon}$ denotes the set of sampled points where the distance measure meets the epsilon threshold.

While ABC algorithms are well grounded in mathematical theory, they make a number of assumptions and approximations. These need to be carefully assessed before considering the use of ABC for a given model. As $\epsilon\rightarrow\infty$ and $N\rightarrow\infty$, $\theta_{sim}=\theta$. Practically, we need to set $\epsilon\neq 0$ to approximate the posterior in a reasonable amount of time. With a high dimensional parameter space, the probability of simulating prior parameters equal to that of the posterior is extremely small. As $\epsilon$ increases, the posterior distribution shifts more towards the prior. A key component of ABC is choosing $\epsilon$ such that the posterior is calculated both accurately and within a reasonable amount of time. One such approach is to run the algorithm a large number of times and consider the top $k\%$ of the sampled parameters as the acceptance area. Another is to run the algorithm multiple times until a desired $k\%$ of $N$ is attained.

A second potential issue is in the use of summary statistics. In most cases the complete dataset is high in dimensionality. When comparing the observed and simulated data, most samples will be rejected and the algorithm will need to be run for an impractical amount of time. If the data is reduced to a set of summary statistics such that $dim(s)<<dim(y^{*})$, this curse of dimensionality is avoided. However, a poor choice of summary statistics will result in a biased posterior that is not truly representative. Additionally, the credible intervals will be inflated due to a loss of information (“Approximate Bayesian computation,” n.d.). This can be resolved by choosing a sufficient statistic, which completely explains the original dataset using a smaller subset of statistics.

The simplest ABC algorithm, rejection sampling, is a straightforward procedure of simulating data with some given set of parameters and evaluating the data against a set of observed data. The algorithm proceeds as follows:

Direct evaluation of the likelihood is not carried out in the above algorithm. For almost any model where data is observed, data can be simulated, and there exists a set of parameter values to be inferred, rejection ABC can be used as a straightforward approach to approximate the posterior distribution. This is achieved by noting that the acceptance probability of a given parameter is proportional to the probability of generating the observed data $y_{obs}$, under the model $p(y|\theta)$ for a fixed parameter vector $\theta$. For a fixed $\theta$, if we generate a dataset from the model $y\sim p(y|\theta)$, then the probability of generating our observed dataset exactly, so that $y=y_{obs}$, is equal to $p(y_{obs}|\theta)$ (Sisson, 2018).

A limitation of the rejection algorithm is a low acceptance rate when the prior is very different to the posterior. Beaumont et al (2009) proposed a particle filtering based method known as ABC-SMC. This algorithm alleviates the issue of low acceptance rates by avoiding low rejection rates of the sampling region and gradually evolving towards the target posterior through population filters. In high dimensional data from a complex system this results in a large gain in efficiency. The algorithm proceeds as follows:

$K_{\theta,t}$ in step 14 of the algorithm is chosen as a conditional density that serves as a transition kernel to “move” sampled parameters and then appropriately weight accepted values. In contexts of real-valued parameters, for example, $K_{\theta,t}(\theta|\theta^{*})$ might be taken as a multivariate normal or $t$ density centred at or near $\theta^{*}$, and whose scales may decrease as $t$ increases. (Bonassi & West, 2015). Kernel choice is discussed further in Section 5.5.1.

The algorithm constructs $T$ intermediary distributions with increasingly small tolerance schedules. The tolerances $\epsilon$ are chosen such that gradual evolution towards the target posterior is achieved. The algorithm continues until $N$ particles have been accepted in population $T$. The final set of particles is an approximation of the posterior distribution. The inclusion of a weight calculation for each particle in a population enables the sampling from population $t-1$ with probabilities equal to the normalised weights. These weights are calculated such that particles further from the population mean are sampled and refined.

The ABC-SMC algorithm addresses the main drawback of rejection ABC in the inefficiency of sampling parameters with a high distance from the posterior, and hence rejecting a large number of samples. The repeated sampling from posterior approximations results in a distribution closer to the posterior distribution (Beaumont, 2010). One drawback exists in the decreasing epsilon tolerance. If the quantile of the distances from which to reduce epsilon is appropriately chosen, the number of accepted samples in each population will not decrease significantly as $\epsilon$ decreases. However, setting this quantile too high will reduce $\epsilon$ such that a larger proportion of samples are rejected, eliminating a key advantage of the SMC approach. This can be addressed by either calculating $\epsilon$ at each population using an adaptive method (described in Section 5.5.3), or including conditional density estimation for the final population sample.

An extension to the original SMC algorithm was proposed by Bonassi & West (2015). A joint kernel $K_{t}(x,\theta|x^{*},\theta^{*})$ is used on the joint distribution of accepted values $(x,\theta)$ to raise the importance of proposals where the simulated data $x_{i}$ is close to $x_{obs}$. This addition is known as ABC-SMC with Adaptive Weights. The algorithm proceeds as follows:

The added step of updating the data based weights adds computation time, although this is usually negligible in comparison to the simulation time (Bonassi & West, 2015). The addition of adaptive weights increases the proportion of accepted samples in comparison to Algorithm 2, and hence decreases the overall number of simulations. As the simulations tend to dominate computation time, ABC-SMC with Adaptive Weights tends to be more computationally efficient than ABC-SMC.

There are three methods of controlling the balance between computational speed and posterior accuracy in the ABC-SMC algorithm - the choice of epsilon in each population, the size of the populations, and the number of populations. If the model complexity is small or the computation time is not a concern, the epsilon array can be set to small values, with a large number of populations and a large number of particles in each population. For a complex system with multiple parameters and high dimensional data, the computation time of the algorithm increases exponentially. The two principal reasons for this are 1). The time to run the simulation with the prior parameters and 2). Looping until the epsilon threshold is reached for each sample of the current population. Adjusting the two settings mentioned above will help to reduce the computation time.

An additional measure is that of parallelisation. In the advent of big data there now exist several methods for conducting complex analysis on cluster computing systems. A cluster system comprises of two or more computers working together to perform tasks (ESDS, 2014). They are particularly suited to complex computational tasks which can be divided into many smaller tasks for each computer (node). For the analysis of large datasets it is important to design statistical algorithms which are well suited to this infrastructure. In some cases this is straightforward. In the rejection ABC algorithm, the number of desired loops can be divided between each node and the results merged at the end of the computation. With more complex ABC algorithms with dependencies amongst steps this is less clear. In each population of ABC-SMC there is a dependence on the distance, weight and parameter values from the previous population. Simulations within the current population $t$ can be executed in parallel, and the results sent back to the master node for evaluation until the population threshold $N$ is reached. As the simulation time dominates when using ABC to infer for an ABM, this method of parallelisation can result in a significant decrease in computation time.

Two ABC models were developed to infer the mean and standard deviation of a normal distribution - Rejection ABC and ABC-SMC with Adaptive Weights. The true mean was taken to be approximately 4. A Normal prior $\theta^{1}$ was set for the mean and a Gamma prior $\theta^{2}$ for the standard deviation.

Rejection ABC with a epsilon of $\epsilon=0.5$ and $N=30,000$ results in a high rejection rate of 96.5%. This can be adjusted through a more informative prior, although in many cases this is not available. Another adjustment is in epsilon. If this is increased, more samples will be accepted with the tradeoff of a less accurate posterior. The acceptance rate of 3.5% appears to give an accurate posterior.

The ABC-SMC model has a longer running time, but reduces the variance of both posterior estimates $\mu$ and $\sigma$. The acceptance rate at the final population is 50%, which is significantly higher than the rejection model. The increase in computation time is due to the multiple populations and increased time spent computing weights. For a complex model, the reduced time spent simulating data will outweigh this. The posterior mean estimates are slightly closer to the true mean.

A challenge in designing an ABM to represent a complex system is representing the system in an accurate manner. Where there is observed data for a system, an ABM can be said to be validated against such a system if its output (simulated data) is close to the system’s observed data. ABC is a suitable approach for this validation. If an ABM is designed with control from a set of parameters, ABC can be used to run the simulation with a set of prior parameters. If the ABC is modelled correctly, it will converge on a set of parameters which will produce simulated data close to the observed data from the system that the ABM is modelled from. In comparison to the ‘brute force’ method described in Section 3.1, significantly less time is spent simulating data using parameters far from the true parameters of the system. Further models will focus on the use of ABC algorithms to infer parameters for an Agent Based Model. While the final pedestrian model is relatively simple in nature, it serves as a basis for the addition of larger amounts of data, a higher dimension of parameters, and a more complex ABM environment.