Counterfactual reasoning in space and time: Integrating graphical causal models in computational movement analysis

Movement analysis is distinguished by an emphasis on understanding via observation and association. However, an important component of movement from the human and computer modeling perspective is the processes that bring about movement behavior in the first place. This article contextualizes the graphical causal modeling framework (for association, intervention, and counterfactual causal analysis) in GIScience, and more specifically within movement analysis studies. This is done by modeling the movement behavior of football players, applied to spatiotemporal data generated by an agent‐based simulation. The movement dataset is thoroughly analyzed to infer the statistical associations among its variables, to estimate the effect of an intervention on some of those variables, and to answer a few counterfactual questions from the observations. We conclude that causal graphs (i.e., directed acyclic graphs), if implemented correctly, can assist analysts in infering causal relations from movement data. This research suggests the integration of causal graphs and agent‐based paradigms as one solution for computational movement analysis.


| INTRODUC TI ON
Spatial positions are such fundamental elements of geospatial dynamic phenomena that modeling their movement has become a prominent activity in GIScience (Dodge et al., 2020;Goodchild & Glennon, 2008;Hornsby & Yuan, 2008;Long et al., 2018;Miller et al., 2019).Of specific interest for enhancing utility of such a space-time resource are movement-related queries (e.g., Would the pedestrian flow have been different, had the city council not intervened on that intersection?Why do the cows move to the corner of the paddock in spring: is the grass better or do they feel safer there?Did the team-formation cause the first goal in that football match?).These examples are distinguished by seeking an understanding beyond mere observations and associations, a level of understanding that entails searching for causal relations.This degree of knowledge often permits scientific explanations of phenomena, along with the capability to understand responses to system changes (Griesmaier, 2006;Lipton, 2010).This level of analysis is not advocated by the data-centric view, which ultimately leads to descriptive or predictive models (Han et al., 2011), in GIScience and Computational Movement Analysis (CMA).The introduction of such data-centric approaches (i.e., data mining and knowledge discovery) into GIScience has clearly led to significant progress in organizing unstructured and messy movement data, but not modeling the processes that bring about movement behaviors.However recent studies have acknowledged the obligation to embrace causal reasoning in CMA (Rahimi et al., 2021;Xu & Kwan, 2020).
Following the "causal revolution" in various sciences in the late-20th-century (Pearl & Mackenzie, 2018), and given the accessible resources (i.e., data and computational power), there is emphasis and impetus in GIScience to address causes and effects of movement behaviors (Bleisch et al., 2014;Rahimi, 2021).However, a rigorous theoretically discussed and practically tested methodology is required to encourage developments beyond pattern recognition and visualization.
A considerable number of scientific investigations of causality are founded on David Lewis's (1973Lewis's ( , 1979Lewis's ( , 1986) "counterfactual" theory of causation.Simply put, the idea is to explore what would have happened in the most similar "parallel world" if, contrary to fact, a certain aspect of the perceived reality had been different (Rahimi et al., 2021).A promising methodology comes from a new stream of causal analysis in geospatial science: statistical-based approaches in which analysts try to model the cause-effect relationships among parameters underlying a spatiotemporal mechanism (Christiansen et al., 2020;Lozano et al., 2009;Luo et al., 2013;Zhu et al., 2017).This methodology is mainly associated with the counterfactual causal analysis framework and graphical causal models (GCMs), largely articulated by Judea Pearl (2000Pearl ( , 2009)).
Visual representation of cause-and-effect relations in the form of graphs have become a common accepted practice in causal modeling (Glymour et al., 2019;Glymour & Greenland, 2008;Greenland et al., 1999;Morgan & Winship, 2015;Robins, 2011).GCMs represent the analyst's causal hypotheses about how the world works and how the observations should be collected and analyzed (Greenland, 2010).Among different types of graphs (e.g., directed graphs, directed cyclic graphs, etc.), Directed Acyclic Graphs (DAGs) are used most commonly (Suzuki et al., 2020;VanderWeele & Robins, 2009).DAGs consist of a set of vertices connected with paths or edges to form a graph that is called directed if each path is directed from one vertex to another, and acyclic if it has no directed cycles.It is noteworthy that DAGs do not necessarily follow a tree structure.
DAGs communicate the experts' qualitative information about associations among variables.However, due to their nonparametric nature, they make no claim about the functional form of such dependencies.There are various methodologies to bridge between the qualitative causal claims in DAGs and quantitative association in data (Robins & Richardson, 2011).Structural equation modeling (SEM) is a well-established language for causal and counterfactual analysis (Bollen, 1989;Bollen & Pearl, 2013;Bullock et al., 1994;Christ et al., 2014;Pearl, 1998Pearl, , 2012)).GCMs, DAGs, and SEM are popular in various scientific fields (Kunicki et al., 2023;Li et al., 2018), while they have rarely been considered and discussed in GIScience and CMA.This article aims to investigate the applicability of GCMs in GIScience and to adopt some of their fundamental concepts in movement studies.This is done by providing an applied movement-related example.A detailed description of GCMs and DAGs can be found in the literature (Elwert, 2013;Hitchcock, 2018;Pearl, 1985Pearl, , 1988Pearl, , 1995Pearl, , 2000Pearl, , 2009)).For a discussion on the theoretical basis of causal models and the implementations of DAGs in GIScience, see Rahimi et al. (2021).Section 2 introduces a causal model to represent the movement mechanism and its constituent variables in the context of a football game.Section 3 explains how a causal model must be analyzed and presents results of intervention and counterfactual queries on the movement dataset.Finally, Section 4 wraps up the results of the analysis through a detailed discussion and presents a few suggestions for future study, followed by the conclusion.

| A N UMERI C AL E X AMPLE: ANALY ZING C AUS E S OF MOVEMENT IN A S IMUL ATED G AME
In testing the applications of DAGs in CMA, an example that is spatiotemporally constrained, well-expressed with procedural knowledge, with abundant and dense spatial, temporal and attribute data, is needed for ease of modeling and validation.Here we focus on the movement behavior of football players.Football games have recently attracted scientific attention (Brillinger, 2007;Herold et al., 2019;Lucey et al., 2013;Rein & Memmert, 2016;Sarmento et al., 2014).Movement analysis of football players has also been studied in GIScience (Andrienko et al., 2016(Andrienko et al., , 2017(Andrienko et al., , 2019)).

| Observations
It is initially important to know the mechanism that has generated the data, to examine the credibility and plausibility of DAGs in CMA, allowing the results of the causal analysis to be validated.The generated dataset from a simulation with known parameters, reported by Rahimi et al. (2022), meets this criterion.In this data-generating agent-based football simulation, each game can be run for 27,000 time intervals (each interval is 100 ms), which is one half of a game.In each time interval, players' attributes, their coordinates (referenced relative to the center of the pitch), and their agent-defined objectives of movement are logged.The objectives of movement are in fact a list of seven general tasks (see the next sub section), one of which a player can choose to execute at each time interval.The data also includes the position and the possessor of the ball.The ball possessor is recorded for each individual player, in each time interval, as one of four ("nobody," "opponent," "me," or "my teammate") categories.
In CMA, the movement behavior of an object is described by a set of parameters.Table 1 illustrates a list of movement parameters together with a few other items of complementary information that are chosen from the generated observations to describe the behavior of players.These variables and their corresponding ranges are specific to one of the agents, Player 10 (P 10 ), that is selected for analysis.This individual plays a central role for Team A in the simulation, which can execute various offensive (e.g., dribbling or shooting the ball) and defensive (e.g., marking an opponent, tackling the ball) actions (Rahimi et al., 2022).

| Model
An autonomous movement decision involves selecting a goal, direction, and speed.According to the model configuration in Rahimi et al. (2022), simulated players choose direction and speed regarding their objectives (A-1 to A-7).To select their objectives, players first try to find out "who possesses the ball."If it is not them or one of their teammates, they find a location at which they can intercept (collide or co-locate with) the ball.This is a function of their abilities and the velocity of the ball.If this location is inside their allocated role-area, they attempt to get possession of the ball (A-3).Otherwise, they decide to closely move with (or mark) the nearest opponent (A-2).
When players themselves possess the ball, they choose either to carry (A-4) or shoot (A-6) the ball towards the opponent's Goal, or to pass it forward to one of their teammates (A-5).To do so, they consider their distance to other opponent players, the opponent's Goal, and the center of their role-boxes along with their current abilities.
But if one of their teammates possesses the ball, players try to open-up the space (A-7) by moving forward and towards one of the sidelines of the pitch.Lastly, to incorporate a degree of freedom into the game, players have been designed to move randomly (A-1) around the center of their role-area.Players execute this action when their distance to the ball is more than a certain length (80 pixels or 40 m).In this model configuration, at each time interval, players perceive their environment and consider their pace, agility, and shooting abilities together with their current level of energy, which is reduced over time as a function of stamina, to determine what to do and where to go. Figure 1 illustrates the above causal mechanism in a DAG G1, on the variable set V, described in Table 1.While this model is clearly a simplification of reality, it captures many of the key decisions that a football player would make in a real-world setting.

TA B L E 1
A list of all parameters and their corresponding range set for P 10 .

Parameter
Variable Range (R:) Objective of movement (category a ) V obj {A-1, A-2, A-3, A-4, A-5, A-6, A-7} Ball possessor (category a ) V b_possess {nobody, opponent, me, teammate} Distance to the center of the role-box (pixel b ) Distance to the opponent's goal (pixel) Distance to the ball (pixel) Distance between the collision location c with the ball and the center of the role-box (pixel) Azimuth e to the center of the role-box (angle) Azimuth to the opponent's goal (angle) Azimuth to the ball (angle) Azimuth to collision location with the ball (angle) Azimuth to Player 13 (angle) Azimuth to Player 17 (angle) a Objective of movement and ball possessor are both nominal variables.A value has been assigned to each category but there is no intrinsic ordering to the categories.
b Each pixel in the simulation represents 0.5 m 2 in real world.
c Where P 10 can co-locate with the ball, given its maximum speed and the velocity of the ball.
d Player 13 and 17 are on average the closest opponent players to P 10 , which based on the model configuration (Rahimi et al., 2022) have, among other players, the highest impact on its behavior.
e Azimuth from point A to point B is the clockwise angle between the north reference ray on A to Ray AB.
f The model assumes that players always head forward, so their heading and direction is the same angle.

| Specification
Although DAG G1 communicates our qualitative information about associations among variables, it makes no claim about the functional form of such dependencies.As discussed in the introduction, this article implements SEM as a mathematical tool to fill this gap.SEM consists of a set of equations representing a system of cause and effect among a set of variables (Hitchcock, 2018).They are classically based on the assumption of normal distributions of variables and continuous outcomes as in conventional linear models.However, working with DAGs often requires an expanded set of models (e.g., generalized models).Pearl (2009) describes that in causal interpretation of SEMs, to accommodate models involving discrete variables and effect heterogeneity, the path coefficients must be redefined as a general capacity to transfer changes among variables of interest.Therefore, a DAG can be represented with the nonparametric SEMs that make no claim about the distributions or functional form.Modelers may specify a parametric statistical model of their DAGs, often considering the type of variables and dependencies in the target dataset (Elwert, 2013).Structural Equations (1-4) represents a generalized linear model specification of DAG G1. (1) representing the causal mechanism underlying the movement decisionmaking process of P 10 in the simulation (see Table 1 for an explanation of each variable).Note that U energy and U obj , respectively, represents all unobserved variables that may affect V energy and V obj .According to the model, U energy only includes stamina, while U obj includes pace, agility, and shooting.
where, the path coefficients i , i , i stand for the mean effect of variable i on energy level (V energy ), objective of movement (V obj ), and speed (V speed ), respectively.While the coefficient obj,j denotes the mean effects of variable j , moderated by V obj , on direction (V heading ).A notable detail is the existence of two unobserved exogenous variables (U energy , U obj ) in DAG G1 that do not appear in the right-hand side of Equations ( 2) and ( 4).This is because although U energy and U obj vary across different players, they are assumed constant for an individual player during the entire game (see Rahimi et al., 2022).This implies that the distribution of V energy and V obj for an individual player (here P 10 ) can be determined by their direct observed causes, thus U energy and U obj can be removed from Equations (2-4).

| Effect heterogeneity
Before moving to the next section, it is important to discuss a hidden causal mechanism, known as effect heterogeneity, in DAG G1.In DAGs, a variable may partially or entirely mediate the effect of a variable on another: in Figure 1, V energy is on both paths directed from V time to V speed , while V obj is only on one path.In addition to the effect mediation, a variable may modify the effect of a variable on another: in Figure 1, V obj modifies the effect of all other variables (V a_role , V a_op_goal , V a_ball , V a_collision , V a_p13 , and V a_p17 ) pointing to V heading .In explaining this effect modification, although V a_role , V a_op_goal , V a_ball , V a_collision , V a_p13 , and V a_p17 dictate the P 10 's heading, their effect magnitude and direction can vary based on the value of V obj .In conventional SEM, the modification effect is represented by an arrow pointing from the modifier variable to another arrow that connects the cause to the effect variable.This effect heterogeneity is often illustrated differently in DAGs and must be identified from the structural equation that represents the causal graph.Here, the framework suggested by VanderWeele and Robins ( 2007) is followed to encode this direct effect modification, which can be read off Equation (4).Pearl and Mackenzie (2018) articulate the counterfactual causal inference process into a conceptual framework -the ladder of causation-within which they distinguish three association, intervention, and counterfactual levels of causal analysis (Rahimi et al., 2021).To implement this framework in CMA, as a numerical example, the causal effect of variable V energy on the movement behavior of P 10 is estimated through Sections 3.1-3.3.Also, after specifying a DAG (Equations 1-4) and prior to estimation of the causal effects, analysts must ensure that their DAG is identified (Christ et al., 2014).This is a critical step in counterfactual causal analysis.To maintain the flow of the article, the identification analysis of DAG G1 is presented in Appendix A.

| Associations
Perceiving ("seeing") the environment and finding associations among these observations, framed by a mental model, is the first step towards causal reasoning.The association level in statistical terms involves computing the probability distribution of variables of interest given their hypothetical cause(s) within passive observations.In this stage, we aim to compute the dependencies claimed in Equations (1-4) by fitting a generalized structural equation model (GSEM) on the observations using maximum likelihood estimation.Fitting GSEMs often involves various pre-and post-estimation steps, such as "exploratory and confirmatory factor analysis," obtaining the parameter estimates, or "investigating goodness of fit" (Deng et al., 2018;Skrondal, 2004).Such processes are generally important but omitted here for simplicity.
The emergent behavior of autonomous decision-makers in agent-based models, particularly in agent-based movement models, is always known as being much more intricate than the sum of their attributes and the characteristics of their environment.Given such complexities, the fitted model seems to fairly estimate the association between a few observed parameters and the movement behavior of P 10 .Figure 2 compares the observed and estimated distributions of four endogenous variables (V energy , V obj , V speed , and V heading ), in Equations (1-4).
Figure 2a shows that Equation (1) has almost perfectly modeled P V energy | V time .This level of accuracy is expected as this variable has only one predictor.Figure 2b also presents a high accuracy in estimating A small difference between the modeled and observed probabilities is detectable for V obj = A-1 (walk randomly).It is because of a restrictive rule in the F I G U R E 2 Observed and estimated distributions of endogenous variable V energy (a), V obj (b), V speed (c), and V heading (d) in the fitted GSEM, represented in Equations (1-4).Variables V energy , V obj , V speed , and V heading represent P 10 's energy level, objectives of movement, speed, and heading.
simulation model that dictates P 10 to execute A-1 when its distance to the ball (V d_ball ) is more than 80 pixels (40 m).
The value of this crisp threshold is hard to estimate, which seems to be responsible for this slight misestimation.
The probability distribution of V obj nevertheless has been correctly modeled in more than 95% of occurrences.
However, according to Figure 2c, the model has difficulty estimating P V speed | V energy , V obj .Here, the continuous variable V speed has two explanatory variables in Equation (3), one of which (V obj ) is a categorical variable.
Without detailing the discussions on regression analysis with categorical variables, the effect of V obj on V speed is estimated by comparing the average value of V speed over different strata of V obj .This often makes estimating the upper and lower values of the dependent variable difficult, which appears to be the case for V speed .It is worth mentioning that the average estimated value of V speed over different strata of V obj has almost remained the same as in the observations.Modeling P V heading | V obj , V a_role , V a_op_goal , V a_ball , V a_collision , V a_p13 , V a_p17 also proves to be as complicated as it looks in Equation ( 4), with the modification effect of the categorical variable V obj .Figure 2d compares the modeled and observed distribution of V heading .This figure points out several important issues.First, according to Rahimi et al. ( 2022), P 10 heads towards the opponent goal while possessing the ball.On its way, if P 10 encounters an opponent player and decides to dribble pass them, it changes its heading with the maximum possible turning angle.This angle appears to be towards north (0 and 360°) or south (180°) on most occasions.Thus, the observations contain two clusters of frequently occurring values, around 180 and 360° that seem to average out in the estimation.Second, P 10 executes A-7 (open-up the space), when one of its teammates possesses the ball and its distance to the ball is less than a certain length.In such situations, P 10 heads forward towards one of the corners of its allocated rolebox, which are located around 70 and 115 azimuth degrees relative to its position in most occasions.Therefore, two other less frequently occurring values have been observed around 70 and 115°.However, the fitted model instead has predicted higher probabilities around 90°.This has a simple logical explanation.When itself or one of its teammates is in possession of the ball, P 10 aims to move towards the opponent's goal that on average is in its 95° azimuth angle during the observation.This scenario happens very often (nearly 20% of times during the simulations), hence the model estimates a high coefficient for V a_op_goal on V heading that ultimately results in estimating a higher chance of heading towards the opponent goal.

| Intervention
Investigators of causation need to look beyond mere recognition of patterns and regularities to model the outcomes of interventions.In observational causal analysis, researchers are interested in obtaining the probability distribution of a variable when another variable is forced to take on a specific value to examine the causal relations between them.This section is particularly focused on measuring the average amount of change, causal effect, that an intervention on V energy imposes on the behavior of P 10 , using the "calculus of intervention."Formally speaking, it aims to compute: The above postintervention probabilities are obtained to estimate the average causal effect of V energy on each V obj , V speed , and V heading .This estimation is based on the most general nonparametric criterion, called "do-calculus" or "calculus of intervention" (Pearl, 2009;Shpitser & Pearl, 2008).This criterion is written as modification in the current model that permits replacing Equation (1) with a constant value.A common practice, for a model that is identified and permits causal interpretation (see Appendix A for a detailed identification analysis of DAG 1), suggests the following three simple steps (Mahoney et al., 2013): • Set V energy to the minimum in its range, or other relevant values (e.g., 1 standard deviation below the mean), do v energy = 21.5 ,while keeping all control variables as they have been observed, then estimate the distribution of V obj , V speed , and V heading .
• Set V energy to the maximum in its range, do v energy = 100 , and again estimate the distribution of V obj , V speed , and V heading .
• Compute the differences between the estimated distribution of V obj , V speed , and V heading over the two above states of V energy .
This procedure allows assessment of the potential impact of a change in V energy on the behavior of P 10 .Figure 3 is a graphical representation of do v energy = 100 .Figure 3a compares pre-and post-intervention distributions of V obj .The first noticeable effect is the decreasing likelihood of V obj = A-2 (watch an opponent) opposed to the increasing likelihood of V obj = A-3 (try to get the ball).Another effect relates to the A-5 (pass the ball).When P 10 is in possession of the ball, it may consider executing this action for one of the two following reasons.It is either about to cross the boundary of its dedicated role-box area, or at least one of the opponent players is confronting while P 10 's energy level is not enough for dribbling past the opponents.Again do v energy = 100 makes P 10 to decide more often to perform A-4 (carry the ball), when it is challenged, which decreases the likelihood of V obj = A-5.This effect may not be observable in Figure 3a.However, the effect of such intervention is presented in Table 2.This figure also shows that the probability of the other three actions have remained almost the same, meaning that V energy does not have a substantial direct effect on their occurrences.

Figure 3b
illustrates the distribution of V speed .This variable has two parents, V energy and V obj , when V energy is forced to take on the constant value 100, V speed only varies over the variation of V obj .According to Rahimi et al. (2022), P 10 sets its speed to 30, 40, 90, or 100% of its current maximum speed, 1 depending on its objective of movement (V obj ). Figure 3c compares the pre-and post-intervention distribution of V heading .The likelihood of P 10 facing towards the opponents' goal (V obj = 95) increases when its energy level is fixed to 100.This is aligned with the changes in the distribution of V obj , where the chance of V obj = A-4 (carrying the ball towards the opponent's goal) increases by do v energy = 100 .
Table 2 numerically presents the estimated value of V obj , V speed , and V heading prior to and after do v energy = 21.5 and do v energy = 100 .This table also presents both marginal (direct) and total effect of V energy on those variables.
The effect sizes (marginal and total) are the difference between the fitted and predicted mean probabilities or values of the response variable (e.g., V speed ), when the explanatory variable (here V energy ) is set to a fixed value and all other variables are set to their mean (StataCorp, 2019).First, recall that according to DAG G1, V energy is expected to affect V obj only through the direct causal path V energy → V obj .This is apparent through the same marginal and total effect magnitudes in Table 2, over each class of V obj .These results now have a causal interpretation.The effect magnitude value for V obj _A-3, for example, implies that 1 unit increase in V energy , on average, increases the chance of V obj = A-3 by 0.0007.These interventions however show that V energy has no direct causal effect on V obj _A-6 (shoot the ball) and V obj _A-7 (create space).According to the model specification (Rahimi et al., 2022), the former is a function of V b_possess and V d_op_goal , while the latter is only caused by V b_possess .
Table 2 shows a considerable difference between pre-and postintervention estimation of V speed .Thus, V energy imposes a direct and an indirect effect on V speed .The direct effect is on average equal to 0.0077 and the indirect effect is 0.0005.Each unit increase in V energy , therefore, on average increases the value of V speed by 0.0082.
Different from other two dependent variables, V energy only effects V heading through the indirect (mediated) path V energy → V obj → V heading .Thus, while the direct effect is equal to zero, the total effect is considerable and on average around −0.027.

| Counterfactual
Unlike in the intervention stage, where modelers are interested in average changes that a variable can impose on another, the counterfactual analysis considers an individual event or value of a variable.The former describes, for example, how much the likelihood of V obj = A-3 (try to get the ball) changes if V energy is set to 100.The latter examines whether A-3 would have happened, either probabilistically or deterministically, in a most similar possible world where only V energy was.An answer to such a speculative scenario is often turned into a counterfactual F I G U R E 3 Pre-and post-intervention probability distribution of variable V obj (a), V speed (b), and V heading (c).Variables V energy , V obj , V speed , and V heading represent P 10 's energy level, objectives of movement, speed, and direction.
statement whose antecedent is hypothetical.If X and Y are two variables in a structural causal model M (a model that is identified and permits causal interpretation) the counterfactual question "what would have been the value of Y, had X been x?" can be formalized as follows: where, Y x (u) is the solution for Y (in unit, individual, u), had X been replaced by the constant value x, in the modified structural model M x (Pearl, 2000).Now, let us examine if the behavior of P 10 , in some certain situations, is counterfactually dependent on its energy level.For this, the values of V obj , V speed , and V heading are computed at a few critical scenarios that has led to scoring a goal in the simulated game.Restating it formally, the objective is to estimate the following counterfactuals: where, MG1 (represented in Figure 4) is a modified sub-model of DAG G1 within which the path between V time and V energy has been removed (i.e., time does not have a depleting effect on P 10 's energy).In other words, Equation ( 1) is replaced by V energy = 100 in a scenario to investigate the causal relations between P 10 's energy level and its movement behavior at a specific point in space and time.The following is a short description of two assessed scenarios related to the counterfactual questions (summarized in Table 3): would P 10 have behaved differently, in the exact same situation, had (only) V energy been 100?Scenario 1: According to the data, P 10 shoots the ball (A-6) at time 4501, while V energy = 87.This attempt ended up with loss of possession of the ball and receiving a goal against shortly after, at time 4742.The model accurately estimated the same value for V obj at the association level.The observed and estimated value of P 10 's heading are also both the same (V heading = 52).The estimated value of P 10 's speed (V speed = 1.5), however, is slightly smaller than its observed value (V speed = 1.6) due to a negligible inaccuracy of the model described in the previous sub-section.The counterfactual column in Table 3 shows the solutions for V obj , V speed , and V heading in a hypothetical world, where P 10 has its full energy level.The results show that this player would have behaved the same, had V energy been 100, instead of 87.

Marginal (direct) Total
Post-intervention estimation Effect magnitude

Post-intervention estimation Effect magnitude
V energy = 21.5 V energy = 100 V energy = 21.5 V energy = 100 Assessing two different scenarios for both A-5 and A-6 helps to be more confident in concluding that V energy is counterfactually a cause of the movement behavior of P 10 Also, although V energy is a cause of A-4 and A-5, there is no evidence that P 10 's A-6 is counterfactually related to V energy .

| D ISCUSS I ON
The current abundance of observations has encouraged spatial scientists to go beyond pattern recognition applications, aiming to discover causal narratives that underlie spatiotemporal phenomena.This article has presented one approach to accommodate counterfactual causal analysis (i.e., graphical causal modeling framework) and adopt its concepts in a GIScience context, and more specifically within CMA.This popular approach in general involves a three-level process for inferring quantitative causal evidence from passive observations: • Computing association between the assumed causal (X) and effect (Y) events (fitting a SEM the data); • Estimating the average change that an intervention on X could impose on Y (setting X to a fixed value x, estimating the probability distribution of Y, and calculating the differences between pre-and postintervention probabilities); • Quantitatively representing what would have happened to Y, in a similar hypothetical world, had X been x (i.e., a counterfactual scenario: estimating the value or the probability of Y at a single point in space and time).
For practicality, this research worked on a generated dataset from an agent-based movement simulation reported by Rahimi et al. (2022).Along with its introductory purposes, analyzing this simple movement dataset served an important purpose in highlighting inherent problems in more realistic movement causal modeling.The following outlines a few pitfalls that movement researchers should consider: • The first, and perhaps the most important, issue is the causal assumptions (e.g., causal Markov condition and causal faithfulness condition), implied by DAGs.A causal hypothesis should be founded on a set of testable assumptions that analysts must be willing to defend on scientific grounds.It is these causal assumptions that permit a causal interpretation of estimated associations among observations.But some of these assumptions can be easily violated in spatiotemporal causal analysis due to the complexity involved.Causal Sufficiency, for example, could be a restrictive causal assumption in movement analysis, where there are so many temporal, spatial, socio-economic, and psychological factors driving movement behaviors.An implicit communication of such assumptions and making sure that they can hold against data is critical in counterfactual analysis.
• Another closely related difficulty that a movement researcher must address is the identification analysis to see if all (un)conditional (in)dependencies claimed in a model can hold against the data.Identification analysis of a DAG with all environmental, individual, and socio-economic factors involved in human movement decisions can be quite challenging and time-consuming.Additionally, space and time can modify or mediate the effect of almost all factors involved in a movement decision.Such an effect heterogeneity also substantially adds to the difficulty of identification analysis.
• Movement analysts also must decide how to interpret their DAGs and estimate the causal effects.Depending on the problem at hand, they must choose between a parametric or nonparametric, and linear or nonlinear specification that have their own pros and cons.This again seems problematic in CMA, as movement data may come from different sources with variables in different binary, ordinal, nominal, or continuous formats.
Although many useful models (e.g., generalized linear models) have been developed to work with such datasets, they usually come with their own difficulties (e.g., non-convergence, high-dimensionality). Researchers must be ready to try various models when they encounter such problems.
• A complexity issue arises around the number of possible causal models that can hypothetically explain a spatiotemporal phenomenon.In a real-world study, researchers may have several candidate models.These models must be filtered in several steps: (1) in the exploratory data analysis, (2) in identification analysis, and finally (3) in association analysis step (e.g., using the goodness of fit and other relevant criteria).Conducting all these steps offers a comprehensive methodology that makes use of both top-down (expert knowledge) and bottom-up (knowledge extracted from data) in selecting the best candidate model.
• Another concern may arise around the prediction capability of such models that sometimes may be more desirable than their explanatory credibility.It seems that even if DAGs were be capable of inferring causes of movement, predicting the future state of a moving object with such models may not be an easy task.This is however not a restrictive problem for DAGs, and GCMs in general, as they are primarily implemented where the focus is on understanding causality in observational settings.A number of other machine learning methods (e.g., random forest, deep learning) have been proven effective for prediction purposes in GIScience.Therefore, the main purpose of a study (i.e., prediction or causality) must be decided prior to engaging in a complex and time consuming graphical causal modeling practice.
Among various definitions of causality, this article narrowed its attention down on the counterfactual dependence account of "token" causation that seeks relations between two events at particular points in time (Herd & Miles, 2019).DAGs, if implemented correctly, appear to be quite useful (in a limited sense) for inferring causes of movement at specific points in time and space.Although such cause-and-effect relations can have their applications in some studies, it seems too simplistic to reduce a cause of movement to a specific event in space and time: spatiotemporal movements are often products of a chain of events within a complex system of interrelationships and spillover effects.GCMs of course are not to blame for this limited example.Such a simplistic case was presented as a bridge between movement and counterfactual causal inference practices.Future studies must be conducted to examine the application of GCMs for inferring more general (i.e., type causation) causal relations in movement studies (e.g., why does the central defender of Team A mostly remain close to their own goal area during the entire game?).
These views are not completely unknown in GIScience and movement studies.Agent-based modeling, for example, is a prevalent methodology for investigating more general causal relations that is now popular in spatial sciences (Haklay et al., 2001;Mansona et al., 2020;O'Sullivan, 2008;Schelhorn et al., 2005).A short-term solution to some of the above issues could be found in the integration of GCMs with the agent-based modeling paradigm into a data-driven agent architecture (Rahimi et al., 2021).Integrating sophisticated statistical methods and agent-based models is also not unfamiliar in the broader causal analysis context, see for example (Kocabas & Dragicevic, 2013;Maes et al., 2007).
Using validated DAGs to define the dynamic mechanism of a simulation is in fact a pre-development validity assessment, making the final verification and validation of the agent-based movement model more reliable.This is a process that is usually skipped in agent-based modeling practices due to the lack of appropriate statistical tests.
Implementing a graphical causal reasoning model as the inferential means for autonomous moving agents can help to examine the extent to which a causal claim can regenerate the observed movement behaviors.

| CON CLUS ION
This article has examined if the counterfactual causal inference process, offered by observation-based methodologies (i.e., graphical causal modeling), is compatible with computational movement analysis.This is a necessary effort to bridge between movement and causal analysis that had either been previously implicit or even ignored altogether.Following this objective, we applied directed acyclic graphs to the movement of a tightly constrained scenario in space, time, possible actions, and behaviors (a simulated football game).This preliminary attempt led to the conclusion that directed acyclic graphs, if implemented correctly, can assist analysts to infer causal relations from movement data.It is now important to put this methodology in practice on real movement observations.
P(A = a| do(B = b)), where P is the (postintervention) probability distribution of variable A if variable B equals value b.The notation do(B = b) indicates B was set to b by an intervention to analyze its causal effect on A. Intervening on variable V energy means a minimal