Sensitivity analysis of agent-based models: a new protocol

Agent-based models (ABMs) are increasingly used in the management sciences. Though useful, ABMs are often critiqued: it is hard to discern why they produce the results they do and whether other assumptions would yield similar results. To help researchers address such critiques, we propose a systematic approach to conducting sensitivity analyses of ABMs. Our approach deals with a feature that can complicate sensitivity analyses: most ABMs include important non-parametric elements, while most sensitivity analysis methods are designed for parametric elements only. The approach moves from charting out the elements of an ABM through identifying the goal of the sensitivity analysis to specifying a method for the analysis. We focus on four common goals of sensitivity analysis: determining whether results are robust, which elements have the greatest impact on outcomes, how elements interact to shape outcomes, and which direction outcomes move when elements change. For the first three goals, we suggest a combination of randomized finite change indices calculation through a factorial design. For direction of change, we propose a modification of individual conditional expectation (ICE) plots to account for the stochastic nature of the ABM response. We illustrate our approach using the Garbage Can Model, a classic ABM that examines how organizations make decisions.


Introduction
Simulation modeling has become increasingly important in studying organizational behavior (Carley 2002;Harrison et al. 2007).Among the several simulation approaches available, agent-based models (ABMs) play a special role (Prietula and Carley 1994).Since early classics such as the Garbage Can Model (Cohen et al. 1972), researchers have increasingly used agent-based simulation techniques to address relevant organizational, strategic and operational questions (Prietula et al. 1998;Luo et al. 2018;Barnes et al. 2020).As Anderson (1999) emphasizes, ABMs allow researchers to examine open systems-common in management situations-whose behavior cannot be described analytically by equations derived from energy conservation principles or decision-theoretical axioms.
We can broadly define ABMs as computational models in which aggregate outcomes emerge from agents' properties, behaviors and interactions, without the imposition of any top-down constraint.This makes ABMs extremely flexible, as it is relatively easy to vary the building blocks of an ABM.Yet the flexibility of ABMs does not come for free: results are often hard to interpret (Rahmandad and Sterman 2008).Lacking a closed-form expression that links inputs to outputs, agent-based modelers often struggle to assess whether the results of their ABMs are robust and their conclusions are valid.
Addressing these issues is one purpose of sensitivity analysis.However, as our literature review shows, sensitivity analysis is often omitted or only performed in a partial, ad-hoc manner in work that employs ABMs.This is understandable as ABMs raise a number of challenges for sensitivity analysis.For instance, the flexibility of ABMs makes it possible to vary not only parameters but also agents' behavioral rules.While sensitivity analysis with respect to parameters is relatively straightforward, it is less clear how to perform sensitivity analysis on elements that do not belong to a well-defined numerical space.Yet varying non-parametric elements, not just parametric elements, is surely essential in any sensitivity analysis of an ABM.After all, ABMs are often regarded as axiomatic systems that generate outcomes which should be regarded as propositions or theorems (Gallegati and Richiardi 2009).We cannot test the robustness and stability of ABMs' theoretical findings without varying their non-parametric elements.(Henceforth, we use the term non-parametric elements to indicate elements of an ABM simulator other than parameters; thus, we do not use the term non-parametric in the sense of models that are free of assumptions about the frequency distributions of the variables being assessed, a usage found in statistics.) To strengthen future work that employs ABMs, we develop a systematic process for conducting sensitivity analysis on ABMs.Our first contribution is to propose a general conceptual structure for ABMs that identifies their "moving parts"-the elements that can be subjected to a sensitivity analysis.This conceptual structure is helpful for researchers to catalogue all the assumptions that underlie their model and to assess the breadth of their sensitivity analysis.
It is a common misconception that sensitivity analysis is just about showing that the main conclusions of a paper are robust to a range of assumptions.Robustness is only one of several goals.Sensitivity analyses can reveal which elements of a model, or combinations of elements, have the greatest impact on the results, and how strongly various elements interact with each other to influence model outcomes.In our second, more technical, contribution, we propose a design that allows the simultaneous variation of parametric and non-parametric ABM elements to address simultaneously multiple aims of a sensitivity analysis.
A further challenge of sensitivity analysis, which usually involves running many variants of a core model, is the visualization of results.Consequently, we devote special attention to graphical representation.In this context, a third contribution of our paper is a modification of the well-known individual conditional expectation (ICE) plots, a modification that accounts for the stochastic nature of ABM's response, by adding a test on the difference of mean values to validate the statistical significance of the insights regarding direction of change.We call the modification "stochastic individual conditional expectation" plots, or S-ICE (see Appendix 3 for greater details).
To illustrate our approach, we apply our sensitivity analysis process to the Garbage Can Model (GCM) in the implementation of Fioretti andLomi (2008, 2010).We choose this ABM because it is well-known in the management sciences and the software implementation is publicly available.We show that a non-parametric element is the most important driver of the results, demonstrating that a sensitivity analysis focused only on parameters misses important aspects of ABMs.We then show that interactions between model elements are relevant and offer new insights on the managerial interpretation of the findings.These insights would be missed by approaches to sensitivity analysis that vary "one parameter at a time."As we illustrate with our analysis of the GCM, sensitivity analysis can reveal new research insights and point to fruitful extensions for future modeling.

Related literature
Our paper builds on two related literatures, one on agent-based modeling and another one on sensitivity analysis.These subjects are vast and this review cannot claim to be exhaustive.In this section, we provide synthetic overviews of these fields in order to highlight the research gap that we address.Appendix 1 provides a more comprehensive review of the literature, with several additional details.

Agent-based models in the management sciences
ABMs have been used to address a number of topics in the management sciences.A first use of ABMs is for theory development.An early classic is the garbage can model (GCM).The model arose in association with the "garbage can" theory of organizational decision-making (Cohen et al. 1972).In the garbage can theory, organizations are viewed as "collections of choices looking for problems" (Cohen et al. 1972).Each opportunity for choice is like a garbage can into which problems, solutions, and decision makers have been dumped, and what emerges from the collection is "organized anarchy."To any scholar who has spent time in real organizations, the garbage can is a welcome complement to traditional models.Over the years, it was noted that the computer model presented in the original article did not reflect the corresponding verbal theory (Bendor et al. 2001).Several authors have reformulated or extended the initial simulation model, replicating and generalizing the original results (Masuch and LaPotin 1989;Fioretti andLomi 2008, 2010;Lomi and Harrison 2012;Troitzsch 2012).While a classic, the garbage can model (GCM) still remains an influential work in the management sciences (Glynn et al. 2020), with substantial spillovers into other disciplines (Simshauser 2018).Another classical group of ABMs with a similar focus is the family of NK models (Levinthal 1997;Rivkin and Siggelkow 2003;Baumann et al. 2018), which applies techniques first developed in evolutionary biology to study the interrelationship between organizational design and market selection forces.ABMs have also been used to support theoretical investigations in fields that range from innovation diffusion (Garcia and Jager 2011;Fibich and Gibori 2010), knowledge transfer (Levine and Prietula 2012), and organizational learning (e.g., Levinthal 1997) to management accounting (e.g., Wall 2016) and organizational design (e.g., Dosi et al. 2003;Clement and Puranam 2018).
A second use of ABMs is as simulation tools for reproducing the behavior of actual complex systems.Examples of this vast literature are works such as Amini et al. (2012) and Stummer et al. (2015), in which ABM is used to simulate product diffusion, Utomo et al. (2018), on modeling agri-food supply chains, Barnes et al. (2010), Ayer et al. (2019), Barnes et al. (2020) on modeling disease transmission.Often, agent-based models are also part of hybrid simulations; we refer to the reviews of Brailsford et al. (2019), Robertson (2019), Currie et al. (2020) for additional details.Further discussion can also be found in Appendix 1.

Sensitivity analysis
Broadly speaking, sensitivity analysis can be thought of as the exploration of a mathematical or numerical model.The model is typically regarded as a black box that processes a set of inputs and calculates one or more quantities of interest (outputs).Thus, the exploration is not performed by a direct inspection of the model.Instead, the properties of the model are obtained indirectly, by investigating how the output changes given variations in the inputs.
Despite much theoretical progress [see the handbook by Ghanem et al. (2016)], sensitivity analysis is often a neglected task (Saltelli et al. 2020).The recent investigation by Saltelli et al. (2019), who inspect how sensitivity analysis is performed to help scientific modeling across several disciplines, shows that either sensitivity analysis is not applied or it is applied unsystematically.Reviews focused on sensitivity analysis for ABMs reach similar conclusions: in a survey of papers published in Ecological Modelling and in the Journal of Artificial Societies and Social Simulation, Thiele et al. (2014) find that 88% and 76% of the papers published in the years 2009 and 2010, respectively, do not include a serious sensitivity analysis.Similarly, Utomo et al. (2018) survey agri-food supply chain ABMs, finding that 28% of papers do not incorporate any form of sensitivity analysis, 68% of them perform a basic analysis, and only 4% of them apply a more systematic approach.
We have also performed a closer investigation of the works that have applied some form of sensitivity analysis in agent-based modeling; it is reported in Appendix 1.The analysis shows that the majority of authors use sensitivity analysis to check if the main qualitative conclusions of their work are robust when some pre-selected parameters are set to values different from those included in their baseline scenario.To assess robustness, modelers typically vary oneparameter-at-a-time, plot the output variable against different parameter values, and show qualitatively whether conclusions are robust.Only a minority of authors consider other goals beyond robustness, vary multiple parameters at the same time, and use more sophisticated quantitative sensitivity methods.However, realistic ABMs would benefit substantially from a systematic method to identify the importance of computationally expensive assumptions.For example, if the outcome of the model was not sensitive to replacing an expensive assumption with a computationally cheaper one, this would be really useful for model building.At the same time, research on ABMs with a theoretical focus could benefit from the systematic approach to sensitivity analysis that we propose, as it would allow a researcher to develop new insights about the importance and interaction of model elements.Yet, authors are not rigorous about mapping goals to methods, and do not vary models elements other than parameters (and usually only a subset of parameters is pre-selected without a formal procedure), possibly due to the lack of a systematic approach.We aim to fill these gaps in the remainder of our work.
In this respect, our work is related to Lorscheid et al. (2012).The authors propose a systematic approach to opening the black-box of simulation models.They are interested in the whole pipeline of model analysis, including the formulation of the research question.Our work focuses on the sensitivity analysis of ABMs.Given this difference in goals, the two papers focus on different aspects.Lorscheid et al. (2012) include a pre-experimental phase to decide the number of simulations that need to be run to achieve stable results, and they also consider an iterative approach to select interesting parameter ranges.This can be seen as a premise to our approach, which addresses the goals that can be answered by sensitivity analysis; we elaborate on the nature of ABMs creating a conceptual map of the elements of an ABM that can be subjected to sensitivity analysis.We also exploit recent results showing that new sensitivity measures can be extracted from the design; we consider new visualization methods that address different sensitivity goals.

The elements of agent-based models
In this section, we propose a conceptual structure for ABMs, with the goal of facilitating sensitivity analyses.Our structure classifies the elements of an ABM into a number of sets and subsets.Admittedly, it is not the only conceptual framework that can describe an ABM, and sometimes the distinction between elements might be blurred, subjective or dependent on the specific model under scrutiny.Nonetheless, our aim is to classify the "moving parts" of an ABM that can be subjected to sensitivity analysis.We first describe the conceptual structure at a theoretical level, and then illustrate our framework using the Garbage Can Model.

Theoretical structure
Our conceptual framework is represented in the diagram of Fig. 1.
All elements of the ABM are contained in the outer rectangle.They can be classified into six sets and subsets: -Principles (A) are high-level elements that define the nature of the ABM.Principles do not include specific algorithmic implementations; rather, they are conceptual guidelines that influence the modeler in formulating specific procedures or in choosing certain parameters.This is illustrated in Fig. 1, in which principles are represented outside of the rectangle containing all practical assumptions of an ABM.We consider principles to be out of the scope of sensitivity analysis, in the sense that varying principles would lead to a different ABM rather than a sensitivity analysis of a given ABM.-Assumptions are low-level elements that define a specific implementation of an ABM.As such, we would like to include as many assumptions as possible within the scope of a sensitivity analysis.While most assumptions are instances of procedures and parameters, some are not.Such assumptions are identified with the letter B in Fig. 1.Telling principles from assumptions is not always straightforward in practice: changing some assumptions may lead to analyzing a new model.We address this issue further in the discussion section.-Parameters are a specific subset of assumptions.They are cardinal quantities that influence the evolution of the model but are determined outside of the simulation run, either at the initialization stage or through a direct intervention of the modeler.Some parameters (C) characterize the environment or define general properties of the simulation (e.g., the number of time steps), while other parameters are closely associated to agents and determine their properties (D).
Fig. 1 We identify four types of elements: principles, assumptions, parameters, and procedures.Some parameters are associated with agents, as are some procedures.The capital letters in the diagram identify relevant sets discussed in the text -Procedures are also a subset of assumptions.They can be defined as algorithmic prescriptions that determine the time evolution of the ABM.Unlike parameters, procedures are not cardinal, nor are they drawn from an easily defined set of possibilities.Procedures and parameters are sometimes intertwined, in the sense that certain procedures regulate the distribution of parameters, and certain parameters co-determine the effects of procedures.Certain procedures define the main mechanisms through which the model works or through which agents' attributes are initialized (E), while other procedures define what an agent "does", and are also commonly known as behavioral rules (F).
Changing the assumptions that characterize agents makes it possible to conduct a sensitivity analysis "at the agent level", that is to evaluate the impact of agents on simulation outcomes.1

Illustrating the structure: the garbage can model
To illustrate the conceptual structure outlined in Sect.3.1, we consider the Garbage Can Model (GCM).We first provide a summary description of the GCM and then we map out the elements of the model following the conceptual structure of Fig. 1.We focus on the Fioretti andLomi (2008, 2010) implementation, which takes the most careful steps to match the ABM to the original verbal theory of Cohen et al. (1972).In Appendix 2 we include more details on the software implementation, while we refer to Fioretti andLomi (2008, 2010) for greater details.

Description of the GCM
Fioretti and Lomi model four classes of agents: participants, who can make decisions; choice opportunities, openings for participants to make decisions; problems concerning conditions or people inside or outside the organization; and solutions in search of problems. 2The organization is depicted as cells on a grid, and over time, participants, choice opportunities, problems, and solutions move-randomly and independently-on the grid.When at least one choice opportunity, at least one participant, and at least one solution happen to be collocated in a cell, the participants can make a decision.This is easy if, by chance, no problem is in the cell.The participants just declare a decision, even though it solves no problem!Cohen et al. (1972) refer to this as a decision by oversight.Matters are more challenging when one or more problems are in the cell along with at least one choice opportunity, participant, and solution.Then participants can make a decision only if the solution is good enough to solve the problem(s).In particular, participants are characterized by a level of ability, a cardinal variable that takes a value between two extremes (minimum and maximum ability).Likewise, solutions have a certain efficiency, and problems have a certain difficulty, both taking values between a minimum and a maximum.If the sum of the abilities of the participants who are present multiplied by the efficiency of the most effective solution in the cell is greater than the sum of the difficulties of problems in the cell, then the participants can make a decision that solves the problems.Cohen et al. (1972 refer to this as a decision by resolution.A decision by resolution is the most desirable outcome, as it associates a solution to a problem. If the participants in a choice opportunity lack a solution good enough to solve their problems (because their ability is too low and/or the efficiency of the solutions is too low), they are blocked-unable to make a decision or to move on.They remain stuck until another choice opportunity, moving randomly, happens into their cell.The newly-arrived choice opportunity then grabs the most difficult problem and wanders off with it, freeing up the participants and solutions.Cohen et al. (1972) refer to this as a flight.In subsequent wandering, the now-freed participants might very well stumble onto the same vexing problem again.
A final twist in the model concerns organizational structure.The model described so far corresponds to an anarchy setting.In their ABM, Fioretti andLomi (2008, 2010) also allow a hierarchy setting.In a hierarchy, choice opportunities are always assigned a ranking from most important to least important.Participants can also be ranked relative to one another, as can problems and solutions.The rankings are then used as follows: -When hierarchy is applied to participants (that is, the participant structure is hierarchical rather than anarchic), each participant is allowed to be part of choice opportunities as important as, or less important than, herself.So the most important participant can take part in all choice opportunities, but less important participants can take part only in less important choice opportunities.-When hierarchy applies to problems (i.e., the problem structure is hierarchical), each problem can be considered in choice opportunities as important as, or less important than, itself.The most important problems can be considered in all choice opportunities, but less important problems are on the table only in less important choice opportunities.-When hierarchy applies to solutions (i.e., the solution structure is hierarchical), each solution can be considered in choice opportunities as important as, or less important than, itself.The most important solutions can enter all choice opportunities, but less important solutions are up for consideration only in less important choice opportunities.
The model allows researchers the flexibility to apply hierarchy rather than anarchy just to participants, just to problems, just to solutions, or to any combination of the three.3

Elements of the GCM
To illustrate our conceptual structure of Fig. 1, in this section we map various elements of the GCM.
-Principles (A): A key principle of the GCM is organized anarchy: decision makers work on problems they have stumbled upon, using solutions that happen to be available.This can be contrasted with game-theoretic models of industrial organization, in which agents with clear preferences seek rational solutions to well-defined problems.Another distinguishing principle of the GCM is independence between the objects of decision-making.Indeed, choice opportunities, participants (decision-makers), solutions and problems exist independently of one another.This principle distinguishes the GCM from other models where, instead, solutions exist only attached to problems.-Assumptions that are neither parameters nor procedures (B): The assumption that agents move on a regular square grid as opposed to a more sophisticated structure such as a network fits into this class.-Parameters that are not agent parameters (C): A first example of an element in this class is the number of choice opportunities, as it determines how many agents of this type exist, but not their properties.Note that the number of choice opportunities is an assumption rather than a principle: in the original implementation of the GCM (Cohen et al. 1972) there was only one choice opportunity, but some scholars have since argued that multiple choice opportunities are more in line with the verbal theory of the GCM (Bendor et al. (2001;Fioretti and Lomi (2010).Thus, the number of choice opportunities does not seem to lead to conceptually different models.Another example of an element in class C is the size of the grid on which agents wander.-Agent parameters (D): The minimum and maximum levels of ability, efficiency and difficulty determine key attributes of the agents, so we consider them as "agent parameters".-Procedures that are not agent procedures (E): As an example of a procedure in this class, we consider the rule for assigning values to ability, efficiency and difficulty.This procedure assigns these values to agents by sampling the interval delimited by the minimum and maximum values above uniformly at random.While it determines an agent's parameter, it does not define what an agent "does", and so we do not think of it as an agent procedure.-Agent procedures-or behavioral rules (F): The participant, problem and solution structures determine the access of participants, problems and solutions to choice opportunities, and so clearly influence how agents behave.Thus, we consider these procedures as behavioral rules.
archical access to choice opportunities.We prefer, and use, more intuitive terms: participant structure, problem structure, and solution structure.
Footnote 3 (continued) 4 A systematic process for sensitivity analysis: six steps In this section, we outline a process in six steps to make sensitivity analysis of ABMs systematic.We deliberately keep the discussion as general as possible, leaving all mathematical details and most references to Appendices 2 and 3. Also, for illustration purposes, we shall focus on traditional parametric elements in this section, leaving the extension to non-parametric elements to Sect. 5. We visualize the steps in Fig. 2.

Output of interest
The first step of the process is choosing the output(s) of interest.Indeed, as noted in Lee et al. (2015), usually ABMs produce a multiplicity of outputs.One or more of the outputs could be quantities of interest, as long as they are deemed relevant by the analyst or by the decision maker.Moreover, the ABM response is frequently stochastic.For instance, in the Garbage Can ABM, the number of decisions by resolution or by oversight is stochastic with respect to the simulation inputs.Formally, in stochastic simulators the output Y is the conditional distribution of the quantity of interest given inputs X.Not infrequently, analysts are interested in one or more summary statistics or functions of such distribution.For example, quantities of interest can be a moment (the expected value, the variance of Y), a quantile (the median, the 95 or 99 quantile) or the probability that Y is above a given threshold.All these quantities can be the target of a sensitivity analysis.

Goal
After the analyst has determined the quantities of interest, she decides on the goal of the sensitivity analysis (this is also known as setting in the sensitivity analysis literature).A variety of goals are possible.The analyst may be wishing to increase her own understanding of the model behavior, or might be asked to test theoretical aspects as part of a broader research investigation, or may be required to deliver robust managerial insights to a stakeholder.Broadly, one may be interested in the relative importance of alternative model elements (factor prioritization), or in whether they increase or decrease the quantity of interest (direction of change), or in how they interact (interaction quantification), or in whether conclusions drawn from the model are robust with respect to variations in the inputs (robustness analysis).

Elements
A third step is to decide which elements of the ABM to vary.This phase requires a critical review of the model, its main principles, assumptions, parameters, and procedures.For instance, varying a given element may not make managerial sense in a given application, and in this case it would be safe to exclude that element from the sensitivity analysis.The structure laid out in Sect. 3 is a useful reference to the researcher for classifying the element(s) at hand.It helps in balancing which elements are left out of the analysis (e.g., it is hazardous to focus entirely on parameters and to ignore procedures-see Sect.7 for further discussion).

Sensitivity method/design
A fourth step is to choose the most appropriate method for each goal-element combination, and the related experimental design, that is, the choice of how to sample points in the input space.The choice of the method also implies the choice of the scale, local or global.Let us start from a local sensitivity analysis, namely a sensitivity analysis around a particular point in the input space (also known as scenario).
A natural choice for Agent-Based Models is to use finite difference methods.Finite differences are given by the difference between the output at a base scenario and the output at alternative scenarios, obtained by varying one or more of the inputs.The relevant sensitivity measures are main, interaction and total effects (mathematically defined in Appendix 2).They capture, respectively, the individual influence of each input, the influence of an input when varied together with other inputs net of its individual effect, and the sum of the two.Total effects can answer a factor prioritization goal, while main and interaction effects can be used to address direction of change and interaction quantification.Also, the examination of model output across scenarios addresses robustness.
The evaluation of the model at several locations in the input space is the basis for a global sensitivity analysis.Differently from a local analysis, a global sensitivity analysis implies considering the response of the model with inputs sampled across the entire input space, or at least across the portion of the input space that the analyst deems most relevant Saltelli et al. (2008)-see Appendix 2 for mathematical aspects.This exploration is often performed through various forms of Monte Carlo experimental designs.The selection of the global sensitivity method depends on the goal.For factor prioritization, variance methods are a classic set of tools.According to these methods, the most important inputs are those that explain most output variance.Moment-independent methods are a more recent approach that exploits the full distribution of the output given the inputs.For direction of change, popular techniques are, among others, gradients computed at randomized locations and partial dependence functions.For interaction quantification, the analyst can use high-order variance-based methods: instead of looking at which individual inputs explain variance the most, these methods consider the contribution to variance of two or more inputs that are varied together.Finally, the examination of the output behavior across the Monte Carlo-generated scenarios helps the analyst in conducting a robustness analysis.

Assignment of values
Once the goal and the method/design of the analysis have been identified, the researcher selects numerical values for the parameters.The type of assignment depends on the design.For instance, if the analyst is considering a local analysis, she needs to assign reference values that form a basic scenario (usually called base case) and variation ranges/levels to the inputs to reach other scenarios (sometimes called best case or worst case).Conversely, if the analysis is global and the design is probabilistic, the analyst needs to specify values in the form of ranges (supports) and distributions for the inputs.Value assignment is a delicate step in ABMs: One may have constraints on parameter changes dictated by procedures or other model elements.That is, a choice for a procedure may limit the range of "reasonable" values that a parameter can take-Sect.7 provides further discussion.

Results communication/visualization
The last step is the choice of how to visualize the results of a sensitivity analysis.For any combination of goal/element/method, there may be one or more appropriate visualization tools.For example, for factor prioritization and interaction quantification, parametric elements and calculation of main, total and interaction effects (local analysis), bar-charts or tornado diagrams Eschenbach (1992) are suitable.If instead one is interested in direction of change, a variety of tools are available, ranging from spiderplots Eschenbach (1992) to ICE plots Goldstein et al. (2015).In this respect, we note that ICE plots have been defined for simulators whose output is deterministic.Most ABMs have a stochastic response.We then propose a modification of ICE plots that takes the stochastic nature of the response into account.Specifically, because the response of a (stochastic) ABM simulator is, in fact, a distribution, we introduce a test for the statistical significance of the difference in mean values.The test allows us to distinguish whether a small but non-null change in a mean (or conditional mean) is indeed significant or just the consequence of numerical noise.

An experimental design for including non-parametric elements
Local and global methods have been mainly conceived for inputs which are continuous or discrete real numbers.However, ABMs include also elements such as behavioral rules and interaction procedures, that do not belong to a well-defined numerical space.
While it is not always necessary to vary all the elements of an ABM simultaneously, it is important to stress that a researcher who is performing a sensitivity analysis involving only parameters is implicitly fixing a substantial portion of an agent-based simulation.Suppose one performs a global sensitivity analysis on all of the model parameters.As Fig. 1 highlights, only a fraction of the ABM would be tested (subsets C and D).Using the example of the GCM, the effect of the participant, solution and problem structures on simulation outcomes would be ignored.
To address these issues, we treat the alternative specifications of a non-parametric element as the levels of a categorical variable (see Appendix 3 for all technical details).The levels may or may not have an ordinal meaning.We then propose to perform a full factorial design that comprises all possible combinations of (a) levels of the categorical variables associated with non-parametric elements and (b) two or more discrete parameter values for each parameter (the choice of a full factorial design is not restrictive; computationally cheaper designs can be also exploited-see Sect.7).
This design makes it possible to address sensitivity analysis goals in a way that is in between a local and a global approach.In particular, it allows the researcher to compute finite difference measures (main, interaction and total effects) at several locations of the input space and to obtain global measures as appropriate functions of these local effects.
For example, as we show in Appendix 3, from the variance of the main effects it is possible to estimate the total order variance-based sensitivity indices, a recommended global importance measure (Saltelli et al. 2008).This allows the researcher to address the factor prioritization goal.The same sample allows also the calculation of moment-independent global sensitivity measures.We should emphasize that obtaining global measures as appropriate functions of local measures is by no means novel in sensitivity analysis (Morris 1991).Nonetheless, our approach extends this idea to non-parametric elements, making it possible for example to assess the contribution of a behavioral rule to output variance.
Our design also allows the researcher to address additional goals besides factor prioritization.Regarding direction of change, for a given input, the design allows one to consider all pairs of points that are only different because of that input.One can then compute the resulting Newton quotients at all such pairs of points.We propose to visualize this information in Stochastic Individual Conditional Expectation (S-ICE) plots that take into account the stochastic nature of the ABM response.We will describe these plots when presenting our example of the GCM in Sect.6 (Appendix 3 reports additional details).
Our design also allows the researcher to quantify interactions of a parameter and a procedure.To do so, the researcher evaluates the simulator response at all pairs of points that differ in two or more inputs.It is then possible to calculate all corresponding interaction effects as well as second-order Newton quotients.We plot them as histograms to facilitate visualization.Such visualization shows heterogeneity of interactions across scenarios, enriching the insights coming from reporting a single number that represents a global interaction between inputs.
Finally, regarding robustness, we show how S-ICE plots can be used to test the robustness of conjectures on the model (we refer the reader to the example of the GCM in the next section).

The garbage can model as an illustration
This section offers an illustration of our approach through its application to the GCM.

Output of interest
Recall from Sect.3.3 that decisions in the GCM can happen by oversight, flight or resolution.A notable outcome of many simulations is that decisions by resolution are far from ubiquitous.That is, decisions as envisioned in more traditional models of organizations, in which participants gather to apply a solution directly to a problem, are not at all the only thing that happens.Quite often, simulated participants make empty, symbolic decisions that solve no problems (oversight) or get bogged down in intractable issues (flight).The portion of decisions by resolution is therefore a natural quantity of interest in GCM simulations and a focus of previous literature (Takahashi 1997;Fioretti andLomi 2008, 2010).In this section, our sensitivity analysis focuses solely on the fraction of decisions by resolution for reasons of space; the same analysis can be conducted for other outputs.To deal with stochasticity, for each set of inputs we perform 100 simulations with different random seeds, and report the average fraction of decisions by resolution across these simulations.Each simulation lasts 1000 time steps; this is enough for results to be sufficiently stable, in the sense that we observe that the output of interest shows little variability across simulation runs.Concerning this issue, the recent work of Vandin et al. (2021) proposes a way to automatically select the number and duration of simulation runs, provided some requirements.

Goal
We consider obtaining insights regarding all four main goals: factor prioritization, direction of change, interaction quantification and robustness analysis.

Elements
We select some elements of the GCM to illustrate the range of types of elements shown in Fig. 1.Of the elements that we consider (see Table 1), the number of opportunities and the minima and maxima of ability, efficiency and difficulty are parametric elements (D and F in Fig. 1, respectively), while participant, solution and problem structures (E) are procedures and thus non-parametric elements.
In contrast, we decide not to vary other elements of the GCM: these are the principles (A), which by definition are not candidates for sensitivity analysis; the size of the grid (D), which comprises 195 cells (15 in one direction and 13 in the other direction), and the number of participants, solutions and problems (D), which are 25 each; and the assignment rule of ability/difficulty/efficiency (C).We also do not vary the structure of the grid-a regular square lattice (B).
Looking back at Fig. 1, we notice that we vary elements in all main sets of our classification, suggesting that our exploration of the model touches all main classes of elements.Our conceptual structure also suggests which parts of the model are left untested.Since we do not subject the assignment rule of ability/difficulty/efficiency to our sensitivity analysis, we are not varying any non-agent procedure (C in Fig. 1).

Sensitivity method/design
We use two variance-based indicators (first order and total order indexes) and two moment-independent methods, namely the and Ku importance measures Baucells and Borgonovo (2013) (see Appendix 2 for the definitions).Using alternative measures has the advantage that they look at alternative properties of the output distribution.If alternative measures concur in indicating an input as important, this indication is more robust than the one coming from a single sensitivity measure.

Assignment of values
We assign participant structure, problem structure, and solution structure (which are non-parametric elements) to categorical variables X 1 , X 2 , and X 3 .These can take two values: X i = {A, H} , with X i = A representing anarchy and X i = H represent- ing hierarchy ( i = 1, 2, 3 ).Only considering these two options is consistent with the organizational structures of the original Garbage Can Model and with the findings of Fioretti and Lomi (2008), who identified Anarchy and Hierarchy as the most relevant structures.The remaining elements that we test are parametric.According to our design, we select a few discrete values that the parameters can take.In particular, the number of choice opportunities can take a base value of 25 (as for all other classes of agents), and an incremented value of 40, which is sufficiently different to the baseline to generate a substantial change in the output.We assign three values to each of the remaining parametric variables: the two extremes and the central value We consider all possible combinations of elements, giving rise to a full factorial design with 11, 664 input configurations (scenarios).Evaluating the ABM at these scenarios, allows one to compute 5,832 main effects for each element taking two values, X 1 ∶ X 4 , and 11,664 for each element taking three values, X 5 ∶ X 10 .Also, it is possible to calculate a total of 664, 848 pairwise interaction effects.

Results communication/visualization
For result visualization we use the tools discussed in Sect. 5.

Results of the numerical experiments
Across the 11,664 scenarios, we obtain an average of 10.3% decisions by resolution, 65.4% by oversight and 24.3% by flight, corroborating the original insights of the GCM.From now on, we focus on the share of decisions by resolution.

Factor prioritization
Figure 3 displays a Pareto-chart visualizing the ranking of the elements.Consistently across the four sensitivity measures, problem structure is the most important element.The number of choice opportunities is the second most important element in affecting the share of decisions by resolution, while solution and participant structure are ranked third and have a similar importance.
Overall, assumptions regarding non-parametric elements impact the output of the GCM more than assumptions regarding parameters.One may argue that this result is driven by the interval we chose for parameters -that is, between 0 and 2 (or 0.0 and 0.2) for minimum values, X 5 , X 7 , X 9 , and between 8 and 10 (or 0.8 and 1.0) for maximum values, X 6 , X 7 , X 10 .The argument would be that choosing relatively Fig. 3 Ranking of the GCM inputs according to first-and total-order variance-based, i and Ku i momentindependent measures.Sorting is based on the total-order indexes 2 i narrow intervals in which parameters are varied could artificially result in procedures looking more important than parameters.To test the robustness of our result, we assigned alternative ranges to minimum and maximum values, X 5 ∶ X 10 .In par- ticular, we repeated the calculations letting X 5 and X 9 take values between 0 and 4 and X 6 and X 10 take values between 6 and 10 (we also let X 7 and X 8 take values between 0.0 and 0.4, and 0.6 and 1.0, respectively).These intervals are almost as wide as they can be relative to each other, as it would not make sense to have a minimum value of, say, ability, be larger than its maximum value.Re-running the analysis, the findings remained unchanged.

Direction of change
Figure 4 displays the S-ICE plots for the inputs.These graphs can be read as follows.Consider the top left panel as a reference.The horizontal axis reports participant structure ( X 1 = A is anarchy, while X 1 = H is hierarchy).The vertical axis reports the expected share of decisions by resolution conditional on the value of X 1 in all scenarios (Because there are several scenarios, it is difficult to distinguish small black dots on the vertical axis, which instead appear as a continuous vertical black line).The large black dots on each vertical line represent the average percentage of decisions by resolution across all scenarios.For instance, given a participant structure that is anarchic (A), on average about 10% of the decisions are taken by resolution.This fraction increases on average to about 11% when participant struc- ture equals H-technically, these two dots form the graph of the corresponding partial dependence function, which, in our design case, is discrete.This information is complemented by the lines joining the smaller black dots at X 1 = A, H , whose purpose is to show the change in the portion of decisions by resolution, scenario by scenario.These lines join pairs of small black dots such that an input in a given scenario is switched to an alternative value in the other scenario, with all other elements unchanged.More specifically, in the first panel, the lines join pairs of conditional expectations obtained when participant structure is A in a given scenario and switched to H in the other scenario, with all other elements fixed.As we mentioned, for each scenario we run 100 simulations to estimate the conditional mean of the share of decisions by resolution in that scenario.Hence, the value we are dealing with is a sample mean; to assess whether the observed change is different from zero, we need to perform a statistical test.We use a two-sample t-test at 5% significance level.If the null hypothesis is not rejected the corresponding line in Fig. 4 is gray; otherwise it is blue or red depending on whether the change is positive or negative.In the first panel, we observe a majority of blue lines, indicating that the switch of participant structure from anarchy to hierarchy increases the portion of decisions by resolution in the majority of scenarios.At the same time, not all lines are blue, i.e., there are scenarios in which the opposite occurs and thus the output is not monotonic in this variable.
We now comment on the other panels.The solutions structure panel (first row, second panel in Fig. 4) delivers insights similar to the participants structure panel.The problem structure panel (first row, third panel) yields instead a different message.When problem structure is A (anarchy), we expect about 11.5% decisions by resolution; this number decreases to about 9% when problem structure is hierarchi- cal.Thus, the switch of problem structure from A to H leads to a decrease (significant percentage-wise) in the average portion of decisions taken by resolution.Note that virtually all lines are red: for almost all combinations of all other inputs, the switch of problem structure from anarchy to hierarchy leads to a decrease in the portion of decisions by resolution.(See the paragraph on robustness below for additional comments on this property.)Moving to the following panel, it is possible to see that the number of choice opportunities, on average, has a positive effect on decisions by resolution; however, the plot shows two well clustered sets of positive and negative realizations.This is a sign that this element is involved in interactions with other inputs.In fact, it can be shown (see Appendix 3) that if an input is binary and the model is additive in that input, then there is only one possible slope for the one-way lines in an S-ICE plot.Thus, changes in slope denote non-additivity.Finally, panels for minimum and maximum values of ability, efficiency and difficulty ( X 5 ∶ X 10 ) exhibit a weaker but regular impact, with ability and efficiency param- eters having on average a positive effect on decisions by resolution, and difficulty parameters having on average a negative effect.Overall, the panels in Fig. 4 signal a non-monotonic behavior of the quantity of interest as a function of the inputs.Traditional local sensitivity analysis methods that vary one input at a time would thus be inadequate to study this output.

Interaction quantification
We then come to our third goal, understanding interaction effects quantitatively.Overall, there are 45 possible pairwise interactions among the 10 variables we focused on.The first question we address is whether pairwise interactions are significant in determining the model outcome.To answer this question, we fit a linear regression model with all 45 interaction terms on the input-output data (we use the subroutine fitlm.mavailable in Matlab).Results show that 43 pairwise interactions are statistically significant at the 1% threshold for decisions by resolution.Thus, interactions do matter in the GCM.Another indication in this sense comes from the sensitivity measures in Fig. 3.A visual inspection of the values of the first-and total-order variance-based sensitivity indexes shows several cases with a large discrepancy between these measures.For example, the values of the first-and totalorder indexes for number of opportunities and participant structure in the Pareto chart of Fig. 3 reveals low first-order indexes and large total-order indexes, suggesting that these inputs owe a significant part of their importance to their interactions with the remaining inputs.
We gain information about the signs of interactions through second-order finite difference interaction effects.In Fig. 5, we report the histograms of the normalized second-order Newton quotients for the five most important pairwise interactions.For each interaction, in the legend we also display the mean value, the standard deviation, and the percentage of the scenarios that lead to a positive interaction.Note that interactions are considered positive when increasing both inputs together leads to an increase in output, by more than what could be expected by raising the two inputs individually.This definition only applies to elements that have a notion of ordinality; when considering categorical variables, this interpretation does not apply.To study interactions between numerical and categorical elements, by convention we consider {A, H} as an ordered set in which A is lower.Therefore, positive interactions between a parameter and a categorical variable mean that output increases when the parameter is increased together with one out of problem/solution/participant structures switching from anarchy to hierarchy across two scenarios.
We find that the number of opportunities is involved in three of the five most important interactions.Specifically, the strongest interaction is between number of opportunities and problem structure.The interaction is positive in 98.8% of the cases and has an average effect of +2.8% on the share of decision by resolution, a large impact when considering that the average number of decision by resolutions is 10.4% .The second most important interaction is the one between participant struc- ture and problem structure.This interaction is mostly positive ( 96% of the scenar- ios).Note that, had the analyst focused only on parameters, the sensitivity analysis would not have revealed the most important interactions.

Robustness
Let us suppose for a moment that the experiments on the GCM were carried out with the goal of validating the following two conjectures: 1. C1 A hierarchical problem structure decreases efficiency in decision making.

C2 A higher number of opportunities increases efficiency in decision-making.
Here, what we mean by an increase in efficiency is an increase in the fraction of decisions taken by resolution.
We can check if these conjectures are true by considering each possible line in the S-ICE plots.Conjecture C1 is confirmed by the experiments of the ABM.Indeed, in the third panel of the top row of Fig. 4, all lines are red and none are blue.There would a few cases in which the change in problem structure from anarchy to hierarchy slightly increases the expected number of decision by resolutions.However, it turns out that none of these increases is statistically significant, and indeed these lines are colored grey in our representation.These cases are a byproduct of the stochasticity of the model and we expect them to drop to zero for a large enough number of replications.In sum, we cannot find any statistically significant case in which a hierarchical problem structure decreases the share of decisions by resolution.Consequently, there is not enough statistical evidence to reject the first conjecture.
Conversely, conjecture C2 is not robust.Indeed, the cluster of red lines in the fourth panel in Fig. 4 shows that the opposite of the conjecture occurs in 26.5% of the scenarios, and many of these cases are statistically significant.In particular, we find that this behavior is critically dependent on the interaction with problem structure.When problem structure is hierarchic, the number of opportunities has a consistently positive effect (in 100% of the scenarios); yet, the effect is mixed when problem structure is anarchic, with a positive effect in 47% of the cases and a negative effect in the remaining.

Some managerial insights
Managerial insights are obtained when a question of interest for an operations research/management science investigator is answered by interrogating the ABM at hand.To illustrate, consider a researcher who is using the GCM to understand the combination of elements necessary to achieve a higher share of decisions by resolution (this number might be considered a proxy of the effectiveness of an organization).As is observable in the problem structure panel of Fig. 4, while the impact of a hierarchic problem structure on decisions by resolution is almost always negative, sometimes this effect is milder.This is because of an interaction effect: the negative impact can be mitigated by an increased number of opportunities and anarchic participant and solution structures (Fig. 5).Thus, the combination that leads to the highest share of decisions by resolution is the one with hierarchic participants and solutions structures, anarchic decisions structure, and a low number of choice opportunities.Note that, had the researcher performed a simple series of one-at-atime sensitivities (equivalent to considering only individual effects), the researcher would have naively picked the structure in which problems and solutions freely float in the organization as optimal.
More generally, interaction analysis reveals a further insight: the strongest interaction effect is the one between number of opportunities and problem structure, which can be seen as complementary variables given the sign of their interaction.Suppose that an organization has a hierarchic problem structure and a large number of choice opportunities.Because the interaction between these two variables is positive, the decrease in the decision-by-resolution share generated by the hierarchic problem structure is mitigated by the simultaneous variation of the number of opportunities.
A further insight for a research using the GCM is understanding what determines the share of decisions by resolutions.All in all, our analysis shows that since the presence of a problem is the defining difference between a decision by resolution and a decision by oversight, it is problem hierarchy that leads to the greatest decline in decisions by resolution.An organization that wants to avoid a proliferation of problem-free decisions should take care to spread problems out widely and avoid having lots of problem-poor, solution-rich decision makers looking for work.
To learn a broader lesson from these results, our sensitivity analysis approach can suggest where further modeling efforts could be fruitful.In our example, given the importance of problem structure, one may think of building sub-models that simulate access of problems to choice opportunities at a finer granularity than the anarchy-hierarchy procedures in the current GCM implementation.

Discussion
We now discuss in depth a number of possible issues that we mentioned in the previous sections.

The case of computationally heavy models
Thanks to the fast execution time of the GCM, we were able to explore all input combinations, but our approach may help the analyst also when this is not possible.On the one hand, consider an analyst wishing to apply our design to a time consuming model.Given a budget of available simulator runs, she can reduce the number of model evaluations in several ways, for instance by considering groups of inputs rather than individual assumptions (Saltelli et al. 2008); as an alternative, she might apply a so-called fractional factorial design, which allows the exploration of the model input space at a lower number of locations than in a full factorial design.In using fractional factorial designs, the recommended procedure would be that the investigator pre-identifies potential interactions of interest (for instance, all second order interactions), and then chooses the design that allows the desired level of resolution.Here, the analyst has available a variety of choices based on orthogonal arrays (Morris et al. 2008) or supersaturated designs (Lin 2000).Moreover, (or alternatively) the analyst may fit the model response with an emulator, and focus on parametric sensitivity in a predetermined region of the model input space.If, according to statistical performance measures, the emulator fit is accurate, then the original (and time consuming) mathematical model can be replaced by the fast emulator, eliminating the computational burden restrictions (see Appendix 1 on works applying emulation in agent-based modeling).
On the other hand, the methodological part of our approach can provide guidance to analysts dealing with models for which computational burden does not allow one to apply the proposed design.Given high resource constraints, it becomes even more delicate to carefully examine the triplets element/goal/method to find the combination(s) that could maximize the insights on model behavior and aid managerial intuition.The analyst could apply our sensitivity analysis iteratively, in conjunction with model building.This would make it possible to avoid computationally heavy assumptions, if cheaper assumptions lead to similar results with respect to the goals of interest.

When is it sensitivity analysis and when is it a new model?
In Sect.3.1, we defined principles as those elements of an ABM that characterize it: changing principles would lead to a new model.As it could be difficult to tell principles from assumptions, especially when considering alternative procedures, this might lead the analyst to the border between exploring a new model versus performing a sensitivity analysis of the same model.If the use of different procedures corresponds to the use of different theories, then, indeed, the analyst might be interpreting results that are in between a sensitivity analysis of a model and the comparison between two alternative models.To be concrete, consider the choice in the GCM between anarchy and hierarchy for, say, participant structure.Both anarchic and hierarchic participant structures are consistent with the principles of the GCM listed in Sect.3.4.However, a procedure that computes the optimal assignment of participants in a way that the number of decisions by resolution is maximized would be in contrast with the first principle from Sect.3.4, and thus lead to a new model.Our investigation would suggest that it is not possible to draw a sharp line for all ABM elements and all ABMs.Telling whether an assumption is in conflict with a principle depends on domain knowledge, and is outside the scope of sensitivity analysis.

Selection of elements
Choosing which elements are subjected to sensitivity analysis and which are ignored depends on many aspects.First, the analyst could disambiguate between elements that are central to the research question, which we call elements of interest, and elements that are needed for internal consistency of the model, which we call incidental.Under this distinction, the analyst could focus the sensitivity analysis on elements of interest.For instance, in the GCM, anarchy vs. hierarchy is certainly an element of interest, while the assignment of ability/efficiency/difficulty (our example element in class E in Sect.3.4) is an incidental element, which we do not vary.
Other classifications of elements are possible.In terms of parameters, Smith and Rand (2018) suggest four possible criteria to select the ones appropriate for a robustness test.Their analysis suggests that model-altering parameters and parameters about whose value the research is uncertain should be part of a sensitivity analysis.Additionally, it can be possible to add controlling parameters, which have policyintervention potential, and environmental parameters.The researcher decides which subset of these parameter groups should be included in the analysis, depending on what the research question is.
Finally, some combinations of elements may yield identical results simply through the mechanics of the model.For example, in the GCM, for values of grid size and number of agents around the current assignment, it turns out that it is equivalent to increase grid size or decrease the number of agents of each type, as what determines the number of potential decisions is the "density" of agents.-Weactually performed additional experiments increasing the number of agents from 25 to 40 and keeping the density fixed, and results confirmed this assumption.-However,testing this effect for a larger number of agents might result in different insights: for instance the model might evidence the appearance of self-organizing behavior.More generally, understanding these equivalences requires a detailed knowledge of the simulator, which may not be possible for large and complicated ABMs.

Assigning values to input variables
To carry out a sensitivity analysis, the analyst is actually asked to make decisions on what element(s) to vary, and on what values to assign to the element(s) under scrutiny.This is a critical step, because the indications of the model as well as the results of the sensitivity analysis depend on the quantitative assumptions.However, this is a common problem for any sensitivity analysis method.We can distinguish alternative situations.In a modeling phase in which information collection on the inputs is partial, the analyst may be willing to assign wide or conservative (albeit subjective) ranges to the inputs, to gain preliminary insights about, among others, the correctness of the model behavior.In a second phase, the researcher might hone in on parameter ranges in which the model changes its behavior, e.g., when a previously found result switches its sign.For complex ABMs that attempt to replicate real world phenomena, relevant parameter ranges might be suggested by field experts.Varying procedures might further lead to the observation that the domain of procedures is actually infinite and that while an infinite domain for a parameter is tractable, such a domain is intractable for procedures.In this respect, our approach can even help the analyst to recognize that she is not in a position to make decisions concerning numerical assignments.Then, the quantitative part of our approach would not be applicable.It would remain an open research question whether such an impasse would be pointing the analyst towards using an alternative sensitivity analysis approach or towards additional modeling efforts or information collection, which are needed before a sensitivity analysis can be fully informative.Also, the assignment of values cannot be separated from the structure of the model: constraints may impose to bind together certain elements, limit the magnitude of relative changes or even their possibility of varying individually.In certain situations one might not be able to disentangle individual effects by varying each input separately.These structural constraints then impact the design that can be chosen for a meaningful sensitivity analysis.
Our approach also opens further research questions.First, while sensitivity analyses on parameters mostly lead to further data collection, sensitivity analyses on procedures and behavioral rules can guide further research on how agents interact or behave.In this respect, our method could increase the synergy between lab experiments and ABMs, highlighted in the recent work of Smith and Rand (2018).Second, the issue about the border of an ABM, i.e., whether changing an element leads to a new model, actually raises the question of what is the border of sensitivity analyses.

Conclusion
This work contributes to the use of ABMs by studying their sensitivity analysis.We have proposed a general conceptual structure that classifies the main elements of an ABM and a six-step approach that allows the researcher to introduce an elementmethod oriented analysis.We have studied a design that enables the variation of potentially all ABM elements simultaneously.Our approach borrows from the literature on the design of computer experiments and adds the calculation of a variety of sensitivity measures; in particular, the design allows the randomized evaluation of finite change sensitivity indices for determining individual and interaction effects.For direction of change, we have introduced a convenient graphical representation modifying the well known ICE plots to account for the stochastic nature of the ABM response.
To illustrate our method, we have carried out a thorough sensitivity analysis of the GCM, varying three procedures and seven parameters.Among other results, the analysis reveal that the most important element is non-parametric and interaction occurs between a parameter and a non-parametric assumption; this interaction would have been overlooked if the sensitivity had focused only on parameters.Limitations of the approach and direction of future research have been underlined in the discussion section.concerning agent-based modeling in the field of decision and management science during the period from 2000 to 2009.In the following decade (2010-2019), the same searching criteria reveal 1261 articles, with an increase of +295%.However, sensitivity analysis remains a challenging step that is often missing or performed in a sparse and unsystematic way.In this section, we review the literature on sensitivity analysis of ABMs according to three criteria: the goal (also called setting), the technique or method adopted, and the element under investigation.
Regarding goals or settings, the most frequently adopted setting in the literature is robustness.In a robustness analysis the modeler is interested in gaining an understanding of how stable a result is, given plausible variations of the model inputs/assumptions.In practice, this is done by checking whether parameter values different from those in a baseline scenario impair the main qualitative/theoretical conclusions of the paper.For instance, Smith and Rand (2018) study the dynamics of a network of workers that either narrow or widen their social network upon being fired, and test whether this choice has an effect on wealth inequality.They show that it does, and that this finding is robust to different unemployment rates and initial network conditions.Similarly, in the context of a network of innovation diffusion, Stummer et al. (2015) study the adoption dynamics of a second-generation biofuel in Austria.The authors perform a robustness analysis on the size of the population and on several network parameters, showing that the path of the biofuel market share is similar for all distinct parameter values (the only exception is the case of very sparsely connected networks).Other examples of this use of sensitivity analysis can be found in Lenox et al. (2006), Puranam and Swamy (2016), Jaspersen and Peter (2017), Keuschnigg and Ganser (2017), Clement and Puranam (2018), and Hernandez and Menon (2018).
While only a minority of authors consider other goals beyond robustness, another relatively popular goal of sensitivity analysis is input prioritization (Dancik et al. 2010;Jiang et al. 2016;Broeke et al. 2016).For example Delre et al. (2016) study the effects of shared consumption on the life cycle of new product launches and, using sensitivity analysis, they show that shared consumption has the largest impact on the output.A similar but opposite goal for sensitivity analysis is the identification of least important parameters, which can then be fixed to save computational cost and reduce output complexity (Ligmann-Zielinska 2018).Another setting that can be found in the literature is direction of change.For instance, Levinthal (1997) shows that a decrease in the smoothness of the fitness landscape negatively affects the ability of organizations to survive in changing environments.Other examples of studies that frame the sensitivity analysis around this setting are Lenox et al. (2006) and Sauvageau and Frayret (2015).An important but often overlooked sensitivity setting is interaction quantification.While this analysis is carried out in few works, an exception is represented by the study of Dancik et al. (2010), where the interaction effects between parameters are thoroughly explored and the most important interactions are reported.Other studies instead rely on the exploration of the behavior of the model by plotting a two-dimensional surface obtained varying two arbitrarily chosen inputs, while keeping all others fixed (e.g., Troitzsch 2012).Usually researchers are interested only in two-element interactions, since higher-level interactions are computationally expensive to evaluate.We note that authors adopts sensitivity analysis to answer one of the mentioned settings, without a-priori stating the goal of the analysis.
Regarding techniques, the most widely adopted method is the variation of oneparameter-at-a-time.This is not surprising, since it is the most common technique through which the robustness of a model is tested.It is used for example in works such as Levinthal (1997), Zhao and Ma (2016), Hassani-Mahmooei and Parris (2013) and Lenox et al. (2006).A further widely adopted method is scenario analysis.An illustrative example is the work of Leitner et al. (2017): the authors combine variations in 5 parameters to obtain 1296 scenarios on which they evaluate their model.The scenarios are then analyzed graphically, without a formal quantification of sensitivity measures.A similar exercise, but with a smaller number of scenarios, is performed in He et al. (2017) and Zhao and Ma (2016).Global sensitivity methods are applied in fewer studies.Segovia-Juarez et al. ( 2004) employ uncertainty and sensitivity analysis with Latin hypercube sampling and the partial rank correlation coefficient.In their study, Marino et al. (2008) extend the sensitivity analysis of the Segovia-Juarez et al. ( 2004) model by performing variance-based sensitivity analysis.For computationally expensive models, we find works that use metamodelling techniques before performing global sensitivity analysis.An example is represented by Dancik et al. (2010), who approximate their model with Gaussian Processes to perform a variance-based sensitivity analysis.Gaussian Process/Kriging emulation is also used in Parry et al. (2013) and Dosi et al. (2018), while Happe et al. (2006) adopt a regression-based metamodel.More recently, Fadikary et al. (2018) propose the use of quantile Gaussian processes as emulators to deal with the stochasticity of the ABM response.Support vector machines are used as a metamodelling tool in Fonoberova et al. (2013).
A third way to analyze the literature is to consider the elements of the ABM that are subjected to investigation.(See Sect. 3 in the main text for the conceptual map).In the majority of cases, the sensitivity analysis concerns only parameters, while all other elements of the model are fixed (e.g., Segovia-Juarez et al. 2004;Fonoberova et al. 2013;Leitner et al. 2017).The parameters that are varied can fall into several categories: some parameters are attributes of agents (e.g., the number N of decisions that organizations face in the NK model), other parameters are properties of the environment (e.g., the number of locations in a spatial model), yet other parameters quantify how many agents of each class exist (Restocchi et al. 2018), etc.Some parameters characterize a behavioral rule, thus impacting a non-parametric element (Riggs et al. 2008).For example, Laurie and Jaggi (2003) study the effect of a vision parameter in the moving rule of a modified version of the Schelling model.Finally, some studies go beyond parameters and change non-parametric elements using scenario analysis, where two or more scenarios are compared while keeping all other elements fixed.For instance, to optimize the distribution of initial alternative fuel vehicles refueling stations, Zhao and Ma (2016) compare four scenarios derived from different targeting strategies.We observe that investigations that involve both parameters and procedures are rarely done, possibly due to the lack of a systematic framework.
We finally note two works that review the use of sensitivity analysis in agentbased modeling before us, namely Lee et al. (2015) and Broeke et al. (2016).Lee et al. (2015) review alternative local and global sensitivity methods and explain the challenges/limitations related to their use.They highlight several technical aspects and the challenges of carrying out sensitivity analyses, given the stochastic, spatial and dynamic nature of the output of ABMs.Broeke et al. (2016) compare one-factor-at-a-time methods with global methods (regression-based and variance-based).Both works focus on methods and on sensitivity analysis of ABMs to variations in input parameters only.

Sensitivity methods: quantitative details
In this appendix we provide background information on Sect. 4. We mathematically define sensitivity analysis methods and give more details about visualization tools.Here we only review tools already existing in the literature; our contributions and some specific details of our analysis are described in Appendix 3. As a guide throughout this section, the reader can refer to Table 2.

Preliminaries
For the purposes of sensitivity analysis, we regard any model as a black box inputoutput mapping that relates inputs X to output Y.Given the scope of our work, we call this black box a "simulator".Here and in what follows, we use capital letters to indicate the generic name of variables, and small letters to indicate specific values that those variables can take.We write the simulator as Y = g( ) + ( ), where g ∶ X par ↦ ℝ , X par is the parametric input space (the subscript par stands for parameter), in general a subspace of ℝ d , d is the number of simulator inputs, the symbol denotes the inputs vector = {X 1 , X 2 , ..., X d } .The term ( ) is a stochastic error term which, if present, makes the response of the simulator stochastic.The simulator is called deterministic otherwise.In case the simulator is deterministic, the output or quantity of interest is Y = G = g( ) .If the simula- tor is stochastic, then we consider as quantity of interest Y = [G| ] , or some other function of the simulator output distribution, such as a given quantile.If the simulator inputs are uncertain, we denote by (X par , B(X par ), ℙ ) the corresponding probability space and by F the simulator input (cumulative) distribution.

Local methods
Let 0 and 1 denote two possible values of the simulator inputs.The corresponding values of the quantity of interest are y 0 = g( 0 ) and y 1 = g( 1 ) .We can regard 0 and 1 as two scenarios, one corresponding to all inputs at the base case and the other with all inputs at the sensitivity case.If the analyst shifts the inputs across these two scenarios, she obtains a variation in the model output response equal to y = y 1 − y 0 .In the literature, methods have been developed to quanti- tatively explain the change y , apportioning it to the main (individual) effects of each input, to the interactions of pairs, triplets of inputs, etc.. Formally, several works have shown that y can be expanded as (Li et al. 2001)   where In the above equation, the symbol x 1 i ∶ 0 −i denotes the point in X par obtained by setting x i at the sensitivity case, while keeping all remaining inputs at the base case.Similarly, the symbol x 1 i,j ∶ 0 −i,j denotes the point in X par obtained by setting x i and x j at the sensitivity case, while keeping the remaining inputs at the base case.
The first order terms in Eq. 1), denoted by 0→1 i , are called main effects of the model inputs Li et al. (2001).These effects are obtained by what is known as a "one-at-a-time" experimental design, in the sense that they are obtained by varying one input only.Main effects can be visualized through Tornado Diagrams (Eschenbach 1992).The normalized version of main effects, is the Newton quotient of g at 0 for the change = x 1 i − x 0 i .When considered as a function of x i , the numerator in Eq. (3) contains the univariate function , which is the one-way sensitivity function of g with respect to x i .The second order terms in Eq. (1), denoted by 0→1 i,j , quantify the residual interaction between x i and x j .These terms are interaction effects Kleijnen (2015).Dividing the second order effects by (x 1 i − x 0 i )(x 1 j − x 0 j ) , we obtain the second order Newton quotients: If inputs are denominated in different units, the second order Newton quotients are not directly comparable.In order to homogenize the scale on the horizontal axes of the graphs in Fig. 5, we reported the distribution of (1) .
The sign of i,j (or  ⋆ i,j ) delivers information about the sign of the pairwise interaction between X i and X j at 0 .A positive i,j implies a synergistic interaction, a nega- tive i,j implies an antagonistic interaction.A similar interpretation holds for higher interactions and for higher order terms in (1).
In the literature, one also defines the total finite change effect of x i (Li et al. 2001).This effect is the sum of all the terms in Eq. ( 1) that contain index i: The total effect 0→1 i is the fraction of y that can be apportioned to x i , due to its individual effect and to all its interactions with the remaining inputs.A visualization that makes it possible to see main, and interaction effects is the Generalized Tornado Diagram (Borgonovo and Smith 2011).In order to compute all terms in the decomposition in Eq. ( 1), one needs to evaluate the model at all vertices of the hypercube joining 0 and 1 .Thus, one is effectively considering a full factorial design between the two model inputs (see below).

Towards global methods: experimental designs
In a global sensitivity analysis the model is explored at N locations, x 0 , x 1 , … , x N .If parameters are discrete then a finite (albeit large) N may suffice in exploring the entire input space, while if parameters are continuous there are infinitely many points to inspect.For space limitations, we present intuition about some designs, referring the reader to Chapter 2 in Saltelli et al. (2008) and to Barton (2015) for a more extensive treatment.We start from discrete input spaces.A standard choice is to consider all possible combinations of the values that model inputs can takethis is known as a full factorial design (Montgomery 2000).Of course, running the simulator at all these input points can be extremely expensive, so computationally cheaper solutions have been proposed.For example, one may consider fractional factorial designs (Montgomery 2000).These make it necessary to run the simulator at only a subset of the points of the full factorial design, selected so as to be able to calculate interaction effects up to a certain order.In principle, the use of an experimental design does not require that inputs are assigned a probability distribution.In the case inputs are assigned a probability distribution, then points are sampled from this distribution.Several sampling techniques are available, among which crude Monte-Carlo, Quasi-Monte Carlo and Latin Hypercube Sampling (Le Guiban et al. 2018).As mentioned in the main text, some sampling techniques ensure that the input space is explored uniformly, without oversampling a portion of the input space and undersampling another part.These sampling techniques are called lowdiscrepancy sequences; a popular one is the Halton sequence (Saltelli et al. 2008). (5

Global methods: variance-based
Among global sensitivity methods, variance-based methods play an important role.Let us assume, for the moment, that inputs are independent.One then has the possibility to decompose the simulator output variance in an equation analogous to Eq.
(1).As proven in Efron and Stein (1981), we can write where In the above equation, 2 Y∶i is the portion of the variance attributed to X i alone, 2 Y∶i,j the portion of the variance attributed to the residual interaction between X i and X j , etc.The function F denotes the cumulative distribution with respect to its arguments, and the functions g i (x i ), g i,j (x i , x j ) are obtained by taking partial expectations of Y through the following equations: We observe that Y|X i = x i is the conditional regression function of Y with respect to X i .In practice, explicitly computing the integrals in Eq. ( 8) and the expected values in Eq. ( 9) is unfeasible, so practitioners using variance-based methods use techniques such as partitioning the input space Saltelli et al. (2008).We discuss in Appendix 3 how we compute variance-based measures for our design.
In the literature, one considers the normalized version called first order Sobol's index, giving the relative contribution of input i to total variance.One also defines the total sensitivity index as (Homma and Saltelli 1996): giving to the relative contribution to variance of input i, taking into account its interactions with other inputs.The total effect 2 i has been proposed as a sensitivity measure in the operations research field by Wagner (1995).Let u ⊆ {1, 2, … , n} be a subset of indices, say Y∶i,j = ∫ g i,j (x i,j ) 2 dF i,j (x i , x j ).… is a global measure of the interaction among the parameters in group u.
Specifically, when u = {i, j} , is a global measure of the strength of the pairwise interaction between x i and x j .We recall that when inputs are dependent, first order variance-based sensitivity measures S i remain well posed.However, the interpretation of higher and total order indexes becomes less clear (Kucherenko et al. 2012).

Global methods: distribution-based
A further class of sensitivity measures whose definition remains well posed under input dependence is the class of "distribution-based" or "moment-independent" sensitivity measures.These sensitivity measures are obtained considering the entire distribution of model output instead of a particular moment (e.g., variance).In general, a distributionbased sensitivity measure can be written as where d(⋅, ⋅) is a distance or divergence between distributions, F Y is the marginal distribution of the simulator output and F Y|X i is the conditional distribution of the simulator output given X i .Depending on the choice of d(⋅, ⋅) , one obtains alterna- tive ways of measuring importance.For instance, if d is set equal to the Kuiper distance between cumulative distributions functions, one obtains the global importance measure (Baucells and Borgonovo 2013) Alternatively, if Y is absolutely continuous, one can choose the L 1 -norm between densities f Y and f Y|X i , calculating the so called −importance measure Borgonovo et al. ( 2014): We observe that, in the case the model output distribution is unimodal, then Ku i and i coincide [see Borgonovo et al. (2014)].It is important to note that sensitivity measures 2 i , Ku i and i possess the nullity-implies-independent property.That is, a null value of any of these sensitivity measures occurs if and only if Y is independent of X i .First order variance-based sensitivity measures ( S i ) do not possess this property. (12) Finally, other global methods (which we cannot describe due to space limitations) include regression-based methods Storlie et al. (2009), Shapley values Owen (2014) and others Ghanem et al. (2016).

Global methods: direction of change
All methods listed so far provide indications concerning factor prioritization.However, global methods exist for direction of change and interaction quantification as well.For example, a standard method for direction of change is based on integrating g over the marginal distributions of all inputs but X i .This leads to partial depend- ence functions Friedman (2001), a family of well known indicators used in machine learning.The expression of a partial dependence function is This function is visualized through a partial dependence plot.We note that h i (x i ) in Eq. ( 16) is the expectation of the functions g(x i ; −i ) with respect to F −i ( −i ) .Con- sider a given point 0 −i ; the function z i (x i ) = g(x i ; 0 −i ) is the one-way sensitivity func- tion of X i at 0 −i [see, among others (Borgonovo (2017, Ch. 4)].Note that g(x i ; 0 −i ) is termed an individual conditional expectation (ICE) function in Goldstein et al. (2015).Thus, a partial dependence function is the average of one-way sensitivity functions (or of ICE functions).

Our design, sensitivity measures and S-ICE plots
This section presents the proposal for a design that allows the simultaneous variation of several (possibly all) elements of an ABM and the estimation of global sensitivity measures.

Input space
All local and global sensitivity methods listed in Appendix 2 are defined on numeric input spaces.To include non-parametric elements of an ABM into the sensitivity analysis toolbox, consider the categorical variable X p with support X p = {0, 1, ...x p } .In the case X p = x p , the simulator is run under the p th + 1 available procedure.For example, in the case of the GCM, X p could be problem structure.We would have X p = {0, 1} , corresponding to anarchy and hierarchy. 5We can then enlarge the input space to account for all combinations of parameter values and possible assigned procedures.More generally one can assign categorical variables to any non-parametric element of an agent-based simulator, considering the enlarged input space X ABM = X par × X non−par .The analyst can then assign any desired distribution on the (16) measurable space (X ABM , B(X ABM )) .This opens the door to using any fully fledged global sensitivity analysis method and design.

Definition of the lattice
We aim at a scheme that leads to a full factorial design.To do so, we propose to take as support for the parameters the set } , where x low i and x high i are the lowest and upper boundary values assigned to X i , and x mid i is an intermediate value.Then, X ABM becomes a lattice, whose points are all combinations of the ABM input elements.In case the simulator is expensive to run, one can use cheaper experimental designs, such as the ones listed in Appendix 2.

Local effects at multiple locations
Evaluating the simulator on this lattice, the analyst is effectively evaluating the model at all points necessary to obtain the full finite change decomposition in Eq. (1).Not only, but the lattice allows the analyst to obtain several replicates of first order effects and higher order finite change effects for each input.Evaluating many local effects is by no means novel in sensitivity analysis.For example, it is at the basis of the classical method of Morris (1991).By that method, the analyst computes main effects at a certain number of randomly sampled points in the lattice, and then obtains global measures as appropriate averages of the random main effects.However, the Morris method only makes it possible to address the factor prioritization Fig. 6 The points of the experimental design used in this work.Left graph: two inputs, x 1 and x 2 with two values, 0 and 1, and one level for each value.Right graph: two inputs, x 1 and x 2 ; x 1 with two values, 0 and 1, and one level ; x 2 with three values, 0, 1 and 2, and two levels.Blue arrows: individual effects; red arrows, interaction effects.(Color figure online) goal.Our design instead makes it easy to also address direction of change and interaction quantification.
As a reference, consider Fig. 6.The left graph displays a two-dimensional design, with two inputs and two values (one level) per input.We have four possible combinations, displayed as vertices of the rectangle.The arrows evidence shifts of the inputs from one point to another on the grid.Movements parallel to the vertical and horizontal axis are "one-at-a-time shifts"; diagonal movements are instead "interaction shifts".Evaluating the model at all points in the lattice allows us to obtain: a) the complete decomposition of the finite change A→C y = g(x 1 1 , x 1 2 ) − g(x 0 1 , x 0 2 ) , with the main effects A→C 1 . But we can also find the main effects B→D . These effects are associated with the arrows displayed in the graph (blue for individual and red for interaction effects).The evaluation of the model at all points of the lattice, by the above rearrangement, allows us to find two main effects and two finite change interaction indices.(These are all the distinct effects, however one could consider additional effects or alternative versions of these effects by flipping the sign-inverting the direction of the arrows in Fig. 6.)The right graph displays a two dimensional design in which x 1 has one levels while x 2 has two levels.In this case, the total number of points in the grid is six, and combination of these six model evaluations allows us to find six individual effects for each input and six second order interaction effects.
In general, letting N i denote the number of main effects ( 0↦1 i ) associated with x i , we have a lattice of N = ∏ d l=1 k l points, where k i is the number of levels of each varia- ble.The number of first order effects associated with input The number of second order effects for the pair x i , x j is N i,j = One can verify these formulas for the examples in Fig. 6.For the case in the left graph, both variables x 1 and x 2 have only one level ( k 1 = k 2 = 1 ).Thus, the number of first order effects for each variable is N 1 = N 2 = 1∕2 ⋅ 2 ⋅ 2 = 2 , giving a total of four first order effects.The total number of interactions is N 1,2 = 1∕2 ⋅ 2 ⋅ 2 = 2 .In the right graph, variable x 1 has again one level, while x 2 has two levels ( giving a total of nine first order effects.The total number of interactions is N 1,2 = 2∕2 ⋅ 2 ⋅ 3 = 6.
For the GCM analysis in Sect.6, in the notation of this section, we have k 1 = ⋯ = k 4 = 1 and k 5 = ⋯ = k 10 = 2 , so that N = 5832 , from which we compute N 1 = ⋯ = N 4 = 5832 and N 4 = ⋯ = N 10 = 11664 first order indices and a grand total of 664, 848 second order interaction effects.

Sensitivity measures from the design
Once these effects are available, we calculate sensitivity measures.So far, we have been defining global sensitivity measures in an abstract way.However, when we are dealing with a simulation model with no explicit functional form for g, it is impossible to evaluate the integrals in Eq. ( 8).Several methods have been proposed for this task, with corresponding experimental designs (Saltelli et al. 2008).Given the discrete nature of the input space implied by our experimental design (Appendix 3), we use a method that obtains global measures (such as the total index 2 i ), from local measures (such as finite change sensitivity indexes 0→1 i ), evaluated at many points in the input space.Consider N + 1 locations x 0 , x 1 ,..., x N and compute a series of indexes 0→1 i , 1→2 i , … , N−1→N i .It is possible to see that the following equation holds Gamboa et al. (2016): where i is a random variable denoting the population of first order effects (that is, i is a random variable whose realizations are 0→1 i , 1→2 i , ..., N−1→N   i   ).This means that calculating the variance of the available first order effects and dividing by two leads to a proxy (an estimator) of the total index of the corresponding simulator input.
Specifically, from the main effects it is then possible to compute: (a) N i Newton quotients 0→1 i for each input and (b) one total variance-based effect for each input from Eq. ( 17).From interaction effects, one can calculate second order Newton quotients, gaining insights on the sign and magnitude of interactions.
Moreover, the design yields a dataset of N × d input-output realizations.This input-output sample can be used to calculate the total order variance-based index using Eq.17, and one can also estimate the global sensitivity measures S i (first order variance-based), Ku and through a so-called "given data approach" Plischke et al. (2013).[Specifically, we use the betaKS3.mfunction available at https:// zenodo.org/ record/ 88533 2#.XgoB-kdKiUk.]

S-ICE plots
An ICE plot is made of two main constituents: first, the graph of a partial dependence function h i (x i ) in Eq. ( 16), second the graphs of a number of ICE functions.The simultaneous use of these two constituents was suggested by Goldstein et al. (2015) to enrich the information provided by partial dependence functions.In fact, the marginal expectation at the basis of h i (x i ) produces a smoothing that hides patterns real- ized for specific values of 0 −i .Instead, such patterns are potentially revealed if one plots also the g i (x i ; 0 −i ) functions, as Goldstein et al. (2015) observe.In our case, regarding the first constituent, we observe that marginal expectation coincides with conditional expectation in our design.Denoting the realizations of X i with (x 1 i , x 2 i , … , x K i ) , the partial dependence plot represents the set of pairs ( x k i , [Y|X i = x k i ] , k = 1, 2, ..., K ).We shall use big black circles on the vertical axis to represent [Y|X i = x k i ] .Regarding the graph of ICE functions, consider pair of points of the type (x 0 i , g(x 0 i ; k −i )) and (x 1 i , g(x 1 i ; k −i )) .Note that the two values g(x 0 i ; k −i ) and g(x 1 i ; k −i ) belong to the univariate function z i (x i ) , whose graph is represented by two points on the vertical axis (one at x 0 i and the other at x 1 i ) in the case x i is binary.To evidence pairs of corresponding points (x 0 i , g(x 0 i ; k −i ) and (x 1 i , g(x 1 i ; k −i )) , we join such pairs with a line.We use blue if the output change is positive, red if the output ( 17) change is negative and gray if it is non-statistically significant.In fact, we recall that g(x 0 i ; k −i ) = [Y| = (x 0 i ; k −i )] is the result of an expectation over simulation repli- cates.That is, the quantities reported in the graph are ĝ(x 0 i ; k −i ) = ̂ [Y| = (x 0 i ; k −i )] .That is, due to numerical noise, we may register ĝ(x 0 i ; k −i ) − ĝ(x 0 i ; k −i ) ≠ 0 , albeit small, even if in reality it is g(x 0 i ; k −i ) − g(x 0 i ; k −i ) = 0 .Then, visualizing the change with a blue or red line would be misleading, because, in fact, the change would not be statistically significant.We then use the Student t-test for differences in the two conditional means ̂ [Y|X i = x 0 i ] and ̂ [Y|X i = x 1 i ] to assess the statistical significance of the change.In the case the test provides a "not significant" report, a gray line is reported.
The prototypical representation of an ICE plots involves continuous functions with continuous inputs on an interval.In this case, the graphs of h i (x i ) and of each g(x 0 i ; k −i ) are continuous lines.Implementation-wise, the graph is obtained by discre- tizing the values of X i into a dense enough grid and then reporting the corresponding values of h i (x i ) and of each g(x 0 i ; k −i ) on the y/axis to render the visual impression of continuity.Nonetheless, the definitions of ICE functions as well as partial dependence plots hold unchanged when inputs are discrete.The difference in the visual representation is that on the vertical axis one registers a sequence of points (dots) rather than a continuous line.To illustrate, in the case of binary inputs with values x 0 i and x 1 i , a partial dependence plot reduces to two points (x 0 i , h i (x 0 i )) and (x 1 i , h i (x 1 i )) ; the plot of any individual ICE function is a pair of points of the type (x 0 i , g(x 0 i ; k −i )) and (x 1 i , g(x 1 i ; k −i )).We also note that ICE plots deliver information about whether g responds additively (or not) to variations in x i .In fact we can prove that if g is additive in X i , and X i is binary, then all slopes in ICE functions are parallel.Proposition If model input x i is binary and the output is additive in x i , then all lines in the direction plot of x i are parallel.
Proof If x i is binary, let x − i and x + i denote its two possible values, with x − i < x + i .Let a i (x − i ) , a i (x + i ) denote the values of the univariate function a i (x i ) at x − i and x + i , respectively.If g is additive in x i then it can be written as g( ) = a i (x i ) + a −i ( −i ).Consider next the one-way sensitivity function z i (x i ) = g(x i ∶ k −i ) .At any point (x + i , k −i ) we have: g(x + i ∶ k −i ) = g i (x + i ) + g −i ( k −i ) and at any point (x − i , k −i ) , we have: . Therefore, we find Note that the above equation implies the constancy of the Newton quotients, q.e.d.◻ Consequently, one or more changes in slope in the S-ICE plot of a binary variables implies that the input of interest is involved in interactions with the other inputs.

NetLogo implementation
In this section, we take a look at the NetLogo interface of the Garbage Can Model by Fioretti and Lomi6 (Fig. 7).This particular software implementation includes some extensions of the traditional GCM, such as the possibility to backpass difficult problems and to have a competence/incompetence system to assign parameters for ability, efficiency and difficulty.To remain as close as possible to the original theory of the GCM, we do not consider these extensions in our study.However, it would be straightforward to introduce these elements in the sensitivity analysis to gain more insights into their role within the model.
In our simulation, we adopt a 15 × 13 grid where 25 agents of each class move in a random fashion.The number of agents is proportional to the size of the grid so that there is a non-trivial amount of interactions between agents; note that, once the density of agents is fixed, the share of decisions of each kind is not affected by simultaneous changes in the number of agents and the size of the grid.Participants are represented on the grid with the man icon, solutions with a triangle, problems with a circle and choice opportunities with a square.Ability, efficiency and difficulty are shown at the bottom of the correspondent icon.A run of the simulation is composed of 1000 time steps, as this is enough to have stable values for the share of the different decision styles.Participants stay the same during the course of the simulation, while solved problems and adopted solutions leave the grid and are replaced with new ones with different difficulty and efficiency.Ability, efficiency and difficulty are assigned independently of the ranking of each agent within the organization.

Fig. 2
Fig. 2 Six steps for sensitivity analysis of agent-based models

Fig. 4
Fig. 4 S-ICE plots for the share of decisions by resolution as a function of the elements.Red: negative change.Blue: positive change.Grey: change is not statistically significant.(Color figure online)

Fig. 5
Fig. 5 The five strongest interactions for decisions by resolution

Table 2
Mapping from goals to methods, depending on the scale (local, global, and our design) Numbers refer to the equation in which each method is

Fig. 7
Fig. 7 NetLogo interface of the Garbage Can Model implementation

Table 1
Inputs for the GCM sensitivity analysis of their respective ranges-see Table1.Ranges are selected to avoid overlap of the minimum and the maximum values of attributes.