Transformer-Based Spatial-Temporal Counterfactual Outcomes Estimation
Authors: He Li
Haoang Chi
Mingyu Liu
Wanrong Huang
Liyang Xu
Wenjing Yang
Venue: ICML 2025 poster (2025)
Keywords: Counterfactual Outcomes Estimation
Potential Outcome Framework
Spatial-Temporal Data
Primary area: general_machine_learning->causality

Abstract: The real world naturally has dimensions of time and space. Therefore, estimating the counterfactual outcomes with spatial-temporal attributes is a crucial problem. However, previous methods are based on classical statistical models, which still have limitations in performance and generalization. This paper proposes a novel framework for estimating counterfactual outcomes with spatial-temporal attributes using the Transformer, exhibiting stronger estimation ability. Under mild assumptions, the proposed estimator within this framework is consistent and asymptotically normal. To validate the effectiveness of our approach, we conduct simulation experiments and real data experiments. Simulation experiments show that our estimator has a stronger estimation capability than baseline methods. Real data experiments provide a valuable conclusion to the causal effect of conflicts on forest loss in Colombia. The source code is available at this [URL](https://github.com/lihe-maxsize/DeppSTCI_Release_Version-master).

Transformer-Based Spatial-Temporal Counterfactual Outcomes Estimation
He Li Haoang Chi
Mingyu Liu
Wanrong Huang Liyang Xu' Wenjing Yang'
Abstract
know, time and space are closely entangled and inseparable (Einstein, 1915; 1922), and the physical world naturally possesses spatial-temporal attributes. Thus, it is crucial to perform counterfactual outcomes estimation with spatial- temporal attributes.
The real world naturally has dimensions of time and space. Therefore, estimating the counterfac- tual outcomes with spatial-temporal attributes is crucial problem. However, previous methods are based on classical statistical models, which still have limitations in performance and gener- alization. This paper proposes a novel frame- work for estimating counterfactual outcomes with spatial-tempora attributes using the Transformer, exhibiting stronger estimation ability. Under mild assumptions, the proposed estimator within this framework is consistent and asymptotically nor- mal. To validate the effectiveness of our approach, we conduct simulation experiments and real data experiments. Simulation experiments show that our estimator has a stronger estimation capability than baseline methods. Real data experiments pro- vide a valuable conclusion to the causal effect of conflicts forest loss in Colombia. The source code is available at this URL.
However, the classical causal inference frameworks (Pearl, 2010; Rubin, 1974) cannot be directly applied estimate counterfactual outcomes with spatial-temporal attributes. A few studies Christiansen et al. (2022); Papadogeorgou et al. (2022) have initially explored causal inference meth- ods for spatial-temporal data based on classical statistical models. Christiansen et al. (2022) propose to use the struc- tural causal models for apatial-temporal data. However, in their setting, current outcomes are mainly affected by cur- rent treatment. Thus, their estimands primarily focus on the causal effects simultaneously, which does not fully address the issue of temporal carryover effects. Papadogcorgou et al. (2022) extend the potential outcomes framework to the spatial-temporal setting. Although they propose the spatial- temporal potential outcomes framework, they don't 't explic- itly propose a method for the computation of propensity scores with spatial-temporaa attributes, which is essential for outcomes estimation. Besides, they use classical statis- tical methods, such as kernel methods, which could suffer from complicated data patterns. The kernel methods rely on correctly specified kernel functions and smoothing parame- ters to ensure optimal performance. Thus, previous works still have limitations in performance and generalization.
1. Introduction
Causal inference plays a vital role in various fields, such as epidemiology (Lawlor et al., 2008; Robins et al., 2000) and economics (Baum-Snow & Ferreira, 2015; Dague & Lahey, 2019). Understanding and utilizing causality helps us reveal the mechanisms of the physical world, predict the occurrence of events in the real world, and manipulate their outcomes. Counterfactual outcome prediction (Pros- peri et al., 2020; Wang et al., 2025) is a promising direction in causal inference. In practice, some causal problems have spatial-temporal attributes, such as the well-known "But- terfly Effect"": butterfly flapping its wings in Brazil may eventually lead to tornado in the United States. As we
In this work, we propose a novel framework for the esti- mation of counterfactual outcomes with spatial-temporal attributes using deep learning. An overview of the studied problem is shown in Figure 1, the studied ppatial-tempore data is the time series of point patterns. Our objective is to estimate the expected number of occurrences of the out- come events within a specific region under a counterfactual treatment assignment strategy. The counterfactual treatment assignment strategy refers to a strategy that does not exist in the observational data. In other words, we aim to investi- gate "what will happen to outcome events when we employ other treatment strategies?". Accordingly, we propose deep- learning- based estimators within a spatial-temporal causal inference framework motivated by Papadogeorgou et al. (2022). The proposed estimators account for the effects of past treatments on current outcomes, thus capturing the
Eqqual contribution 'College Computer Science and Tech- nology, National University of Defense Technology, Changsha, China *Intelligent Game and Decision Lab, Academy of Mili- tary Science, Beijing, China. Correspondence to: Wenjing Yang <wenjing. yang @ nudt.edu.en.>.
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Wransformer-Basee TpaporalTemporal Counterfactual Outcomes Estimation
Outcome occurred location
Treated location
Spatial points
current covariates
Now
Xı(s)
Time history treatment
history outcome
Zz-1(s)
Yg-1(8)
Past
Zt(s)
Y((s)
current treatment
current outcome
(a) Studied
data.
(b) Causal structure of the studied problem
Figure 1. Overview of the studied problem. The left figure demonstrates the studied temporal data. Each layer in the cube represents spatial point pattern in time step. The blue blocks represent the outcomes that occurred, while the orange blocks indicate the treated locations. The point patterns at multiple time steps from the past to the present constitute spatial-temporal data, which can be viewed high-dimensional tensors. The right figure illustrates the causal structure of the studied problem, where refers to spatial location and f represents time.
temporal causal effects. Besides, under mild assumptions, our estimators are based on inverse probability weighting and two innovations. First, based on convolutional neural networks (CNNs), we propose an efficient method to com- pute propensity scores when both treatment and covariates are high-dimensional, such point pattern series. Sec- ond, employ Transformer-based model to estimate the intensity functions of point processes that characterize spatial-temporal data. Moreover, our deep-learning -based estimators within this framework exhibit excellent statistical properties, such as consistency and asymptotic normality.
approach by both simulated and real experiments. The real data experiments study the cause-and-effect of con- flicts on forest loss in Colombia, revealing a valuable conclusion: longer and more intense conficts lead to an increase in forest loss.
2. Related Works
2.1. Temporal Causal Inference
Temporal causal inference aims to investigate the causal re- lationships among time series data. Several studies have fo- cused on estimating counterfactual outcomes or treatment ef- fects of time-varying treatments. Examples include Granger Causality (Granger, 1969), Marginal Structural Models (MSMs) (Robins et al.. 2000), and Recurrent Marginal Structural Networks (RMSNs) (Lim, 2018). However, time- varying confounders may lead to bias in the estimation. To address the problem of the time-varying confounders, several deep learning-based methods, such as Causal Trans- former, were proposed (Bica et al., 2020a;b; et al., 2020; Liu al., 2020; Melnychuk et al., 2022; There exist methods targeting other types of sequential data, such natural language (Chi et al., 2024) and temporal omics data (Zhang et al., 2024). While these methods pro- vide strong theoretical foundations for temporal causal in- ference, their models do not include spatial dimensions.
To validate the effectiveness of our approach, we conduct simulations and real data experiments. Simulation experi- ments demonstrate that deep-learning-bee -basee estimator exhibits lower estimation bias than the four baseline meth- ods. Real data experiments study the causal effect of con- flicts on forest loss in Colombia from the year 2002 to the year 2022. The results conclude that the longer duration and more intensity of conflicts will cause more forest loss in response.
al., 2021).
Overall, our contributions summarized as follows:
We study counterfactual outcomes estimation with the spatial-tempores attribute, a more general setting, and propose an effective deep-learning-baat solution.
We propose efficient CNN-based method address the calculation of propensity scores when both treat- ment and covariates are high- sigh--imensional, such as point patterns. In addition, we use the Transformer to model the intensity functions of point processes, used to characterize spatial-temporal data.
2.2. Temporal Counterfactual Prediction
Temporal counterfactual prediction refers to performing counterfactual outcome prediction under time-varying set- tings. Seedat et al. (2022) propose TE-CDE, a neural- controlled differential equation approach for counterfac- tual outcome estimation in irregularly sampled time-series
We empirically demonstrate the effectiveness our
2

---

Transformer-Baser Spatial-Temporal Counterfactual Outcomes Estimation
data. Wu al. (2024) introduce conditional gener- ative framework for counterfactual outcome estimation under time-varying treatments, addressing challenges in high-dimensional outcomes and distribution mismatch. El Bouchattaoui et al. (2024) propose an RNN-based ap- proach for counterfactual regression over time, focusing on long-term treatment effect estimation.
often characterized using point process models (Baddeley et al., 2008). The rectangle layer with colored blocks in Figure Ia illustrates an example of spatial point pattern. We employ spatial point patterns to describe the treatments and outcomes in time step.
Spatial Point Process and Intensity Function. The spa- tial point process is a statistical model that describes the spatial distribution of the locations of events in a given re- gion. Formally, let V-)) denote the counting measure, and D, c R is measurable subset of IR A spatial point process can be defined as collection of random variables (M()) c D,}, where (A) represents the number of events occurring in a specific spatial region A. An im- portant characteristic of the spatial point process is the ex- pected number of events occurring within a region A, i.e., [A)] It can be calculated by an intensity function, which represents the average density of events occurring within unit region. Let s denote location and ds denote a small region located at location with area v(ds), the definition of the intensity function X()) is:
2.3. Spatial Causal Inference
Spatial data is commonly found in the environment, such as air quality data a region and the incidence rate of a disease in a region (Cressie & Wikle, 2015). Spatial data often includes spatial correlation and heterogeneity, which challenge treatment effect estimation (Akbari et al., 2023; Ning et al., 2018). To address the problem of spatial cor- relation. Jarner et al. (2002) eliminate the effect of known covariates by matching the treatment group and the control group and estimating the latent spatial confounding under this design. Concerning the heterogeneity treatment, Causal Forest (Wager& & Athey, 2018) is non-parametric random forest-based method for heterogeneity treatment effects esti- mation. The above methods ignore the temporal dimension, while we focus on spatial-tempol data.
X(5) = lim E[W(ds)] o(de)-0 v(ds)
Then we have B[W(A)] (a A D.. Intensity functions are crucial to characterize spatial point processes, which can be utilized to generate spatial point patterns.
2.4. Opatial-Temporaa Causal Inference
Spatial-temporal data refers to the data consisting of spa- tial and temporal information (Cressie & Wikle, 2015). In the real world, spatial-tempone data, such as temperature variation over time in a region, is ubiquitous, and many of them are time series of remote-sensing images. However, there has been limited research on the causal inference of Pratial-temporal data until now. Christiansen al. (2022) extended the structure causal model to make it applicable to spatial- -entialleempoal data. However, in their settings, cur- rent outcomes are mainly affected by the current treatments, which does not fully solve the temporal carryover effect. Papadogeorgou et al. (2022) extended the potential outcome framework on spatial-temporal data and employed kernel functions to estimate the treatment effects. However, they didn't 't explicitly propose method to compute propensity scores with spatial-tempora attributes, which is essential for treatment effects estimation. In practice, causality is effective for better understanding multi-modal data, such as video (Li et al., 2023; 2022) and market data (Li et al., 2024).
Spatial Poisson Point Process. An important spatial point process is the spatial Poisson point process, which is also the primary focus of this paper. We continue to use X(s) to denote the intensity function of a Poisson point process. An important property of a Poisson point process is ((A) Poisson(AA). and > = ʃA ((s)ds. Poisson(.) denotes the Poisson distribution. According to the properties of the Poisson distribution, E[ = AA. The detailed proof in Appendix E.
4. Method
We first introduce the potential outcome framework for spatial-temporal data. Based on this framework, we describe the estimands of interest. Finally, we derive the estimators and implement them with deep learning models.
4.1. Potential Outcome Framework for Spatial-Temporal Data
3. Background
Before formalizing our problem setting and estimands, we introduce the potential outcome framework for spatial- temporal data (Papadogeorgou et al., 2022). Specifically, denote Zj(s) as a binary treatment variable, and 2z(s) = 1, Ze(s) = 0 represent the presence and absence of treatment at time t and location s, respectively. Let = (1,2...., denote the time index set, and Q denote spatial region.
Now we introduce the important concepts used in this frame- work and estimators.
Spatial Point Patterns. The point patterns describe the locations where events occur within given region and are
3

---

Transformer-Basce Spatial-Temporal Counterfactual Outcomes Estimation
Then Z, = (+)ls Q) is the collection of binary treat- ment variables of all the locations in the region Ω2 Ze can be regarded spatial point pattern generated by spatial point process, with Z, (s)= = indicating the presence of a point and (s) = O indicating its absence. Z = [Zelt E 7) is the collection of treatments in time index set y. By taking a subset of Z, we construct a historical treatment set up
{21, 22, FF, (z)),, to intervene the treatment. Similarly, let (n(.) = Fh, (.) x (., x... x Fiuz (.) denote a joint distribution with M independent ones, and 2 M 1.1] (2--M+1: 2 , 24) denote the last M elements of Zea. We can apply (.) to Z[t M + 1.1] to obtain 2g1 (Fn) == {Z1,..., 2,, Fh,, (2+-M++1), .1.,. ..... Fn (z))}. It is noted that the proposed treatment intervention method is an im- provement of the method in (Papadogeorgou et al., 2022). Specifically, Papadogeorgou et al. (2022) assume the in- tervention distributions within FFr are all the same (Fi), while we allow the distributions in FØr to be diverse and change with time (i.e., Fi.)). Recall that FHH denotes the intervened distributions of previous treatments. Therefore, the estimands in our setting are more general than those in (Papadogcorgou et al., 2022).
to time t, i.e., Ze:= {21, ,,....,Z,, . By assigning Zst as Zgt, we have the historical treatment realization žgt { 22,,,2,.... .,,,,........................... 21, where 21, 22,..., be are the assignments of 21, 2 Zz Then, we denote S.,- {s (+) = 1) as the set of locations that have received treatment in time f.
Similarly, we can formalize the outcome. Specifically, denote (s) as a binary outcome variable, and (s) I, (s) 0 represent the occur- rence and non-occurrence of the target event at time t and position 8 under treatment Zei. respectively. Let Y((Z.)) = (Yr(Zz))
Estimands.
Based on the treatment intervention method, we define the expected number of outcomes in time / and region w as the estimand of interest N2 (P'H). Under the intervention distribution (.),), we have
(2} denote the potential out- come in a region Ω given the historical treatment Zgt, by as- signing Zei as Zse. and we can distinguish (zs)) from the notion of observed outcome y Yz( (Z, -zz-= (z<z) represents a spatial point pattern generated by indicating the presence of point and (s) = 0 indicating its absence. Likewise, Yse and Y represent the potential outcome and observed outcome up to time I, respectively. Sym"(=..) = (s e (2x)) = 1} is the collection of
(FH) = San (Syp)(es((PP)) nwldFu(2) M 1,4.)
(1)
spatial point process, with (s)
Furthermore, taking the mean values of No (FFFH) over time t M, +,...,,TT T,
we obtain,
1
T
N=(FH) T T-M+14 Σ (FH..
(2)
the locations whose target outcomes occur in time t.
It is noted that the treatments of M periods follow the intervention distribution Fa we specified, which does not necessarily exist in the observable data. Thus, the proposed estimands are counterfactual. To better demon- strate the estimands, we provide a running example that walks through each term of the estimands in Appendix B.
Let X, denote the covariates at time t and Xar (X1, X.) denote the covariates up to time f, and Ie 13 22,.,, is its realization. Finally, by combin- ing the treatment, outcome, and covariates, we obtain the historical information. Let hige = {zge, Yal (gse} be the observed historical information up to time t, and Hst = (Zez, Yst, ) be the potential historical information up to time t. Clearly, heg c Hst. We summarize all notations in Appendix A.
4.3. Assumptions
To identify the above estimands, we introduce the necessary mild assumptions.
4.2. Counterfactual Outcomes
Assumption 1. (Unconfoundedness)) We assume that con-
dition on hgg. Ze is not dependent Hga 7, 1l H Hgelhst- This assumption is also known as ignorability (Rubin, 1978), which is commonly used in causal inference, indicating the absence of unmeasured covariates.
We aim to estimate the expected number of outcome- occurring locations in an area under counterfactual treat- ment assignment mechanism. First, we introduce the method for treatment intervention.
Assumption 2. (Poisson Assumption) We assume that Z., ((sz)) and Zilhgt-1 are generated from spatial Poisson point processes. This assumption reasonable since the spatial Poisson point processes widely used to character- ize the distribution of discrete events in region (Cressie & Wikle, 2015).
Treatment Intervention. We consider stochastic treat- ment intervention, making treatment follow a specific in- tervention distribution. Let F.) denote a spatial point pattern distribution drawn from a spatial Poisson point process with intensity function h. Then (x)) repre- sents we assign treatment realization zt following the distribution F.).). For historical treatment realization Zge (*1, 22, 2 we can apply Fx(-) to the last element of Eạt (the latest treatment), i.c., Zes (F).) =
Additional assumptions and their explanations can be found in Appendix C. We provide the proof of identifiability for the estimands in Appendix D.
4

---

Transformer-Bant Spatial-Tempores Counterfactual Outcomes Estimation
4.4. Estimators
4.5. Theoretical Properties
With the above assumptions, we derive estimators for the corresponding estimands. The proposed estimators are based on inverse probability weighting (IPW). As key component of IPW, we first introduce the propensity score within the framework used in this work.
Now we introduce the important theoretical properties of the propensity score and estimator.
Proposition 1. The propensity score is a balancing score. For any t E 7
2, ll lgg)
Propensity Score. The propensity score within this frame- work is the probability of treatment Z₁ given historical in- formation lig.11 denoted as (z+) = P(Z.= ze(/gge-).
Proposition 2. Dimensional reduction property of the propensity score: if Ze 1l Hgrhhs then Z. 11 Hz)les(z)..
Counterfactual Probability. We define the counterfactual probability of treatment Z, as (x.)), Px.)) fn (Z, 2). fr(-) denotes the probability density function of the intervention distribution (.) A detailed explanation of why (*)))) is "counterfactual" is Appendix H.
(The consistency and asymptotic normal- ity of the estimator) Let Var denote the Variance, (V, denote the number of the observed outcome Y, in region w. If all assumptions hold, and T co, we have that
Proposition 3.
VT(N -N((H)) 4 N(O,)),
Intensity Function. In section 4.1, we have defined (za))) (Y, (s))) 22). (z is discrete spatial point pattern generated by spatial Poisson point process. We denote the intensity function of this latent spatial Poisson point process as (ex. this intensity function is assumed to be continuous, serving as a spatial smoothing of (Est).
where = limT -0 T-M+1 El-M os and Up
V Var [ar[[]]---+) Ph (N,( N (Y.)
(,,,)
The proofs of the above propositions are in Appendix C.
4.6.
Based Realization
Inverse Probability Weighting Estimator. With the above components, the inverse probability weighting estimator 8)(Fn,)) is derived as follows:
This subsection introduces how we employ neural networks to realize the above estimators. Figure shows the full architecture of model.
))(H,,))= I (z)( j-t-M+1 ez(2,) (x))
4.6.1. CALCULATE THE PROPENSITY SCORE
Challenges. We first discuss the challenges that high- dimensional data brings and why the classical classification- based method cannot address this problem. In the classical setting, treatments take limited values (e.g., binary treat- ments). Thus, we can utilize covariates to train a classi- fier of treatments and employ this classifier to predict the propensity scores. However, in our setting, the treatment Z, is high-dimensional. Specifically, 7, (7{(s)ss e 2} con- tains all the binary treatment variables of locations in region Ω. If there are 100 locations s, the Z, can take 2100 values. Then we need to train a classifier with 2100 classes, which is unacceptable. To address the challenges, we propose the following method.
Y: FF, s) can be viewed as the intensity function of a spa- tial Poisson point process that generates (z (FrH According to the definition of intensity funetion in Section 3, the expected number of outcome-occur points in a re- gion w is obtained by integrating (F(FH, over region w. Therefore, the estimator of No (FF))) is shown as follows:
(FH) J. (FH,
(4)
Then we obtain the estimator of (FH), (FH)) 1 Ef,r (FH).
Dimension Reduction. Since it's hard to compute the propensity scores of high-dimensional treatments directly, we utilize a dimension reduction map to project treatments into a low-dimensional space. Let R(-) denote the dimen- sion reduction map. We define it as follows:
T-M+1
We briefly outline the derivation of our estimator. It is built upon the Inverse Probability Weighting (IPW) approach from Marginal Structural Models (MSMs) (Robins et al., 2000). In particular, the term (z.) (8) in Eq. (3) denotes the intensity function of the observed outcome at time and
R(Z,) (Z)(s)) 1, {
(5)
Ph. (=z)
(2,), serves as a weighting factor analogous to that used in MSMs. The method extends the IaWeb weighting mechanism from
location s, while the product -
R R(Z,) maps Z, to low-dimensional space, representing the count of treated locations in time t. We focus on the count of treated points when constructing estimators. There- fore, this dimension reduction map preserves the informa-
MSMs to a spatial-temporal framework.
5

---

Transformer-Bane Based Spatial-Tempon Counterfactual Outcomes Estimation
Phy (zz) eパ(2))
Estimator: (F(F..)) =
(2z))
Propensity Score: eạ(2))
Intensity Function: Zyp (Est)
(s)
Hss-1
CNN
= ENN
Farward
Decoćés (MLP
Aype (2gr) E
historical information
Hat-1 st-1. 1i Y st- x R2x3((-1)
Matsilooa Alcckion
dimension reduction
coordinates s. sy (2xc) (5..
Transformer Eincoder
Latent Process
R(2,) N
Poisson(21)
ej(2,) E (0,1)
Counterfactual
Probability: Pny (2z)
average
22 E EN
RZZ)EN = Polsson(Az)
(2)) (0.1)
samples of (z,,)
10.
R: real number
2, parameters of poisson distribution
Figure Full model architecture.
tion in the treatments while facilitating efficient computation of propensity scores.
Consequently, according to the properties of the Poisson distribution, F R (hgge) synonymous with the pa- rameter A that governs R (1)))----1- Detailed derivation can be found in Appendix E. Utilizing the Poisson distribu- tion, the propensity score be calculated by the following equation:
After the dimension reduction of treatment, we express the propensity score in the following form:
cx((2z)) P(R(Z,) R(2))
(6)
R((e) A
Then we specify the distribution ((R()))= = ((z))1
e((R(z)))
R(zz)!
to calculate the propensity score. According Assump- tion 2, he is generated by a spatial Poisson point process. With the properties of spatial Poisson point process in section 3, we have (( ~ Poisson(Ai). Thus, the task of calculating the propensity score is transformed into estimating the parameter A₁ of the Poisson distribution Poisson(A)
4.6.2. CALCULATE COUNTERFACTUAL PROBABILITY
Analyzing the probability distribution of high-dimensional variables directly challenging due to the "curse of dimen- sionality". which leads to data sparsity, increased complex- ity of relationships, and higher computational costs (Ver- leysen & François, 2005). Thus, utilize the dimension reduction map to project Ze into low-dimensional space. Consequently, the transformed counterfactual probability takes the form:
CNN-based Regression. We formulate the estimation of into regression problem. Specifically, employ the high- dimensional series hig-1 to regress scalar R( (())). Accord- ing to the definition of regression, such regression model outputs ERR( Z, (hgr-1) Since CNNs efficiently extract local features from high-dimensional tensors. They cap- ture spatial hierarchies and patterns, making them ideal for data like 3D structures. We employ CNNs the regression model and the MSE loss during training. Specifically, denote the output of CNNs as output(-), then MSE ((output(hgg-)) =² A detailed structure of the used CNNs is in Appendix K.
Ph (R(2)) =
R(=)).
(8)
With the properties of spatial Poisson point processes in section R(Z,) follows Poisson distribution given by: R ((Z)) Poisson (22). Then, the objective is transformed into estimating the parameter X2.
As shown in Figure 2, to achieve this goal, we draw sam- ples from the Poisson point process (2,)) and estimate E E(R(Z)) using the sample mean of R(Z,). With the prop- erty of Poisson distribution, E. ((R()))) equivalent X2,
Reegression models can be seen modeling the probability distribution of the dependent variable, with the predicted output corresponding to conditional expectation.
6

---

Seansformer-Basee Temporll Counterfactual Outcomes Estimation
the detailed derivation is in Appendix E. Leveraging the defi- nition of Poisson distribution, the counterfactual probability is computed by the following equation:
from time-series images acquired by the Landsat satellite (Hansen et al., 2013). We consider forest loss events in Colombia from 2002 to 2022 and treat them as outcomes.
A22 R(=.)
UCDP Georeferenced Event Dataset. The UCDP Georef- erenced Event Dataset is a comprehensive dataset document- ing global conflicts, providing details on parameters such as the temporal aspects and geographical coordinates of conflict events (Croicu & Sundberg, 2015). We exclusively examine conflict events transpiring within the territorial boundaries of Colombia and consider them as treatments. The chronological span of our analysis encompasses the years from 2002 to 2022.
Pa (R(2)) = R(z))!
A2.
(9)
4.6.3. TRANSFORMER- BRASSE INTENSITY FUNCTION ESTIMATION
As for the estimation of the intensity function, we draw inspi- ration from the maximum likelihood estimation (MLE). The the neural network to model (scx) (s), and train this network to maximize the likelihood function. We define net: : Ω - R as the output of this neural network. The training objective is to minimize the following function:
core idea is
A brief overview and detailed description of the real-world data are in Appendix L.
18 (yphe Σ In(net(s:))+ + S n (L((lpp) (10)
5.2. Baselines, Metrics, and Implementations
Baselines. The chosen baselines are state-of-the-art litera- ture on the counterfactual outcomes estimation (Lim, 2018; Robins al., 2000), and the heterogeneous treatment ef- feets estimation (Wager & Athey, 2018). These are: MSMs (Robins et al., 2000), RMSNs (Lim, 2018), Causal Forest (Wager & Athey, 2018). We also employ a linear regression- based method (LR) as an additional baseline for comparison. Since these baselines cannot handle the high-dimensional series directly, we transformed all our synthetic data into scalar series adapt the baselines to our setting, Appendix N details the baselines and adaptation.
where s; E (y.'(sg.) is the coordinate in (gc) q is the data distribution obtained from the Transformer encoder and Gaussian sampling, and is a prior standard Gaussian distribution. (L(9lpp) is the Kullback-Leible divergence between two distributions q and p. Note that (n(ne((s))) and Sa (s) are the components of the likelihood function of the spatial Poisson point process. A detailed derivation of the training objective is in Appendix F. An overall archi- tecture of the network is depicted in Figure 2.
Why choose the Transformer? During training, we input coordinates sequence in (zz)(=z))) into the network. Since the Transformer can capture long-term and high-order de- pendencies and meanwhile enjoy computational efficiency (Zuo et al., 2020), the ability capture such dependencies creates more powerful models than RNNs, which facilitates the estimation of intensity function (Zhou et al., 2022; Zuo et al., 2020). Thus, we employ the Transformer to realize this network. In Section 5.3, through experiments, we fur- ther demonstrate the superiority of the Transformer model.
Metrics. For synthetic experiments, we compute the true counterfactual outcomes as ground truth and employ the relative error rate (RER) between the estimation and the ground truth as a metric. Details of the ground truth are in Appendix M.2. The formula of RER is shown below:
RER |sttmated Value True Valuel
|rue Valuel
where ||| denotes the absolute value.
For real data experiments, since the true counterfactual out- comes of real data are unknown, we consider the consistency of our conclusions with existing literature on the effects of human conflicts on natural resources.
5. Experiments and Results
We evaluate our method on both synthetic and real datasets.
5.1. Datasets
Error Bars. For each experimental setup, we conduct 20 independent runs. Reported values are the means, and error bars indicate the standard deviation (+o) over the 20 runs.
5.1.1. SYNTHETIC DATA
We employ the intensity functions of spatial Poisson point processes to generate synthetic spatial-temporal data. Gen- eration details are in Appendix M.
Implementations. For CNNs in section 4.6.1, the details are epoch learning rate 0.001, batch size 64. For the Transformer-basedd neural network in Figure 2, the de- tails are as follows: The number of blocks in the Trans- former encoder is 8, the number of attention heads per layer is 16, the number of layers in the Multi-Layer Perceptron decoder (MIP) 8, epoch = 300, learning rate = 0.0001.
5.1.2. REAL WORLD DATA
Forest Change Data. The Global Forest Change Data is an annually updated global dataset capturing forest loss derived
7

---

Transformer-Bant Spatial-Temporal Counterfactual Outcomes Estimation
5.5. Ablation Studies
Table 1. Real data experiments results. The data in the table are the estimated numbers of forest loss events.
Replace Transformer with RNNS. In Section 4.6.3, we discuss the reason for choosing the Transformer the back- bone to estimate the intensity function. To further demon- strate the Transformer's superiority over RNNs, we replace the Transformer with RNNs and replicate the synthetic ex- periments. We present the results of T=64, M=1 in Figure 4, other results and details of RNNs are in Appendix J. Ac- cording to Figure 4, results show that the estimation ability of RNNS is inferior to that of the Transformer.
3
c= 5
c= 6
c=
M = 20.6 = 2.3 20.5 2.2 20.7 2 20.5 1.9 20.8 + 2.0
M 3 21.5 + 1.4 21.6 24 22. 1.9 23.0 = 1.3 23.3 1.9
M 22.4 a 2.3 22.9 a 1.8 23.6 1.7 24.2 = 1.2 24.7
M = 7 24.7 + 1.3 2 1.7 26.7 + 1.2 27.2 a 1.5 28.0 a 2.1
5.3. Synthetic Experiments
Relax the Poisson Assumption To validate the robustness of our methods over Assumption 2 (Poisson Assumption), we relax and replicate the synthetic experiments. Specifi- cally, we add the Gaussian kernel to the intensity functions of synthetic data, thus breaking the standard setting of the Poisson point process. We compare the results with the standard Poisson setting. We present the results of setting T=64, M=1 in Figure 5. The implementation details and other results are in Appendix J.2. According to Figure 5, relaxing the Poisson assumption does not lead to significant degradation in the estimation performance of our method, demonstrating its robustness to data distribution.
Treatment Intervention. Section 4.2, we introduced the crucial estimands (FFn, which need to determine the intervention distribution of treatments over the M pe- riods from t M to t. Let Xz, denote the intensity function of the spatial point process generating Zj. We re- place Xz, with the form: h2 x log (^z,) ). C is a constant, j E (t- (t- M + 1, M 2....... t}, and log (.) represents the natural logarithm. The constant c is introduced to control the magnitude of the intensity function. A larger e corresponds the larger values of the intensity fune- tion, indicating increased average density of treatments. Under our apatial-tempone setting, parameter M is cru-
cial for regulating the duration of the intervention. A larger M signifies an extended period of intervening in the treatments. In synthetic experiments, c (3, .... 7} and M e {1,3). Details of treatment intervention are in Appendix G.
5.6. Computation Efficiency
Now we introduce the computing platform and computation efficiency. All experiments were conducted on an NVIDIA RTX 4090 GPU and an Intel Core 14700KF processor. In synthetic experiments, the synthetic data dimension could be (100, (100,100,192) 190)) In addition, the time for one experimental setting was less than 10 minutes.
Results. We generate three synthetic datasets with time lengths of 32, 48, and 64 (i.e., T = 32, 48, 64.). We conduct experiments on each dataset, and in the main text, we present the experimental results of M = 1. Additional experiment results are Appendix For M = 1, we calculate different e {3,4, ... 7). gradually increasing the intensity of the treatment. The results are presented in Figure Figure 3 visually illustrates results. Across different time lengths, our method consistently outperforms the baselines.
6. Conclusion
In this work, we study the estimation of counterfactual out- comes for spatial-temporal data and develop estimators with deep learning. Our experiments synthetic datasets demon- strate the superior estimation capability of our estimator over four baselines, and experiments on real-world datasets pro- vide a valuable conclusion to the causal effect of conflicts on forest loss in Colombia. For future work, we look to consider more general types of spatial-temporal data, mov- ing beyond discrete point process data to encompass more general spatial-temporal stochastic processes.
5.4. Real Data Experiments
In real data experiments, we consider conflicts in Colombia as the treatment variable Z and forest loss in Colombia as the outcome variable Y. We set M E {1, 3,,,,5, 7}, and e
(3. (3.4...,,7 7), representing a gradual increase in the intensity and duration of conflicts.
Acknowledgements
Results. The results are illustrated in Table Table 1 il- lustrates that the estimated forest loss increases as M and increase. Thus, our results suggest that longer and more intense conflicts in Colombia may increase forest loss. Pre- vious studies De Jong et al. (2007): Eniang et al. (2007); Garzón & Valánszki (2020); Kanyamibwa (1998) conclude that conflicts may harm natural resources like forests, which consistent with our conclusion.
This work was partially supported by the National Natural Science Foundation China (No. 62372459).
8

---

Wransformer-Basel Spatial-Temporal Counterfactual Outcomes Estimation
Ours
Causal Forest
LR
MSMs
RMSNs
T- 64, M-1
T=48, M=1
T=32, M=1
(Reletive Error Pate)
0.6
0.6
0.6
0.2
0.2
0.0g
0.0
0.0
Values
Values of
Values of
Figure 3. Synthetic experiments results of M = 1. The horizontal axis represents the values of c, while the vertical axis represents the relative error rate (RER). Lower lines the graph correspond to methods with higher estimation accuracy.
T-64, M-1
References
RER (Relative Error Rate)
Transformer
Akbari, K., Winter, S., and Tomko, M. Spatial causality: systematic review on spatial causal inference. Geograph- ical Analysis, 55(1):56-89, 2023.
RNN
0.4
0.-
Baddeley, al. Analysing spatial point patterns in r. Technical report, Technical report, CSIRO, 2010. Version 4. Available www. csiro. ..... 2008
0.1
Baum Soou-Snow, N. and Ferreira, F. Causal inference in urban and regional economics. In Handbook of regional and urban economics, volume 5, pp. 3-68. Elsevier, 2015.
Values of
Figure Comparison results of Transformer and RNN.
Bica, I., Alaa, A., and Van Der Schaar, M. Time series deconfounder: Estimating treatment effects over time in the presence of hidden confounders. In International conference on machine learning, pp. 884-898 PMI.R, 2020a.
64,
RER (Relative Error Rate)
0.150
Poisson Gaussian
0.125
0.100
Bica, I., Alaa, A. M., Jordon, J., and van der Schaar, M. Estimating counterfactual treatment outcomes over time through adversarially balanced representations. arXiv preprint arXiv:2002.04008 2020h.
0.075
0.050
0.025
0.000
Chi, H., Li, H., Yang. W., Liu, F,, Lan, L., Ren, X., Liu, T., and Han, B. Unveiling causal reasoning in large lan- guage models: Reality or mirage? Advances in Neural Information Processing Systems, 37:96640-96670, 2024.
Values of
Figure 5. Results of relaxing the Poisson assumption.
Christiansen, R., Baumann, M., Kuemmerle, T., Mahecha, M. D., and Peters, J. Toward causal inference for spatio- temporal data: conflict and forest loss in colombia. Jour- nal of the American Statistical Association. 117(538): 591-601, 2022.
Impact Statement
In this work, propose framework to estimate the coun- terfactual outcomes with apatial-temporaa attributes, which has positive societal impacts of helping predict the coun- terfactual outcomes of spatial-temporal data. Besides, all data used this work are publicly available, and our newly released assets are the source code of our paper, which does not contain unsafe images or text.
Cressie, and Wikle, C. K. Statistics for spatio-temporal data. John Wiley Sons, 2015.
Croicu, M. and Sundberg, R. Ucdp georeferenced event dataset codebook version 4.0. Journal of Peace Research, 50(4):523-532, 2015.

---

Transformer-Basee Spatial-Temporal Counterfactual Outcomes Estimation
Dague, L. and Lahey, N. Causal inference methods: Lessons from applied microeconomics. Journal of Pub- lic Administration Research and Theory, 2933)- 2019.
Li, J., Niu, L.,, and Zhang. L.. From representation reason- ing: Towards both evidence and commonsense reason- ing for video question-answering. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 21273-21282, 2022.
De Jong, W., Donovan, D., and Abe, K.-i. Extreme con- flict and tropical forests, volume 5. Springer Science & Business Media, 2007.
Li, J., Niu, I.., and Zhang, L.. Knowledge proxy inter- vention for deconfounded video question answering. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 2782-2793, 2023.
Einstein, A. Die feldgleichungen der gravitation. Sirzungs- berichte der Königlich Preußischen Akademie der Wis- senschaften, pp. 844-847, 1915.
Li, J., Zhou, Z., Zhang. J., Cheng, D., and Jiang, C. Hfter- net: Hicrarchical fusion transformer for interbank credit rating and risk assessment. IEEE Transactions on Neural Networks and Learning Systems, pp. 1-15, 2024.
Einstein, A. Die grundlage der allgemeinen rela- tivitätstheorie, volume 49. JA Barth, 1922.
El Bouchattaoui, M., Tami, M., Lepetit, B., and Cournède, P.-H. Causal contrastive learning for counterfactual re- gression over time. Advances in Neural Information Pro- cessing Systems, 37:1333-1369, 2024.
Li, R., Shahn, Z., Li, Lu, M., Chakraborty, P., Sow, D., Ghalwash, M., and Lehman, L.-w. H. G-net: deep learning approach to g-computation for counterfactual outcome prediction under dynamic treatment regimes. arXiv preprinl arXiv:2003. 10551, 2020.
Eniang, E. A., Haile, A., and Yihdego, T. Impacts of land- mines on the environment and biodiversity. Envtl. Pol'y & L. 37:501, 2007.
Lim, B. Forecasting treatment responses over time using re- current marginal structural networks. Advances in neural information processing systems, 31, 2018.
Garzón, F. A. M. and Valánszki, I. Environmental armed conflict assessment using satellite imagery. Journal of Environmental Geography, 13(3-4):1-14, 2020.
Liu, R., Yin, C., and Zhang, P. Estimating individual treat- ment effects with time-varying confounders. In 2020 IEEE International Conference Data Mining (ICDM), pp. 182-391. IEEE, 2020.
Granger, C. W. Investigating causal relations by economet- ric models and cross-spectral methods. Econometrica: journal of the Econometric Society, pp. 424 438, 1969.
Melnychuk, V., Frauen, D., and Feuerriegel, S. Causal trans- former for estimating counterfactual outcomes. In Inter- national Conference Machine Learning, pp. 1 15293- 15329. PMLR, 2022.
Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A., Tyukavina, A., Thau, D., Stehman, S. V.. Goetz, S. J., Loveland, T. R., et al. High-resolution global maps of st-ceuuuy forest cover change. science, 342(6160):850-888 -53355058588-00883 2013.
Ning, B., Ghosal, S., and Thomas, J. Bayesian method for causal inference in spatially multivariate time series. arXiv preprint arXiv: /801.008000006888 2018.
Jarner, M. F., Diggle, P., and Chetwynd, A. Estimation of spatial variation in risk using matched case-control data. Biometrical Journal: Journal of Mathematical Methods in Biosciences, 44(8):936-945, 2002.
Papadogeorgou, G., Imai, K., Lyall. J., and Li, F. Causal in- data: estimating the effects
ference with of airstrikes on insurgent violence iraq. Journal of the Royal Statistical Society Series B: Statistical Methodol- ogy, 84(5):1969-1999) 2022.
Kanyamibwa, S. Impact of war on conservation: Rwan- dan euvironment and wildlife in agony. Biodiversity & Conservation. 7:1399-1406, 1998.
Pearl, J. An introduction causal inference. The interna- tional journal of biostatisties, 6(2), 2010.
Lavancier, F., Møller, J., and Rubak, E. Determinantal point process models and statistical inference. Journal of the Royal Statistical Society Series B: Statistical Methodol- ogy, 77(4):853-877, 2015. 2015.
Prosperi, M., Guo, Y., Sperrin, M., Koopman, J. S., Min, J. S.. He, X., Rich, S., Wang, M., Buchan, E., and Bian, J. Causal inference and counterfactual prediction in machine learning for actionable healthcare. Nature Machine Intelligence, 2(7):369-375. 2020.
Lawlor, D. A., Harbord, R. M., Sterne, J. A., Timpson, N., and Davey Smith, G. Mendelian randomization: using genes as instruments for making causal inferences in epidemiology. Statistics in medicine, 133))11333 2008.
Robins, J. M., Hernan, M. A., and Brumback, B. Marginal structural models and causal inference in epidemiology, 2000.
10

---

Transformer-Baser Spatial-Temporal Counterfactual Outcomes Estimation
Rubin, D. B. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of edu- cational Psychology, 66(5):688, 1974.
Rubin, D. B. Bayesian inference for causal effects: The role of randomization. The Annals of statistics, pp. 34-58, 1978.
Seedat, N., Imrie. F., Bellot, A., Qian, Z., and van der Schaar, M. Continuous-time modeling of counterfactual outcomes using neural controlled differential equations. arXiv preprini arXiv:2206.08. 08 2022.
Van der Vaart, A. W. Time series. VU University Amsterdam, lecture notes, 2010.
Verleysen, M. and François, D. The curse of dimensionality in data mining and time series prediction. In International work-conference on artificial neural networks, pp. 758- 770. Springer, 2005.
Vo, T. V., Wei, P., Bergsma, W., and Leong, T. Y. Causal modeling with stochastic confounders. In International Conference on Artificial Intelligence and Statistics, pp- 3025-3033, PMLR, 2021.
Wager, S. and Athey, S. Estimation and inference of hetero- geneous treatment effects using random forests. Journal of the American Statistical Association, 1136523):1228- 1242, 2018.
Wang, Q. C., Xiao, 7., Mao, Y., Qu, Y., Shen, J., Lv, Y., and Ji, X. Model predictive task sampling for efficient and robust adaptation, 2025. URL pettps://arxin org/abs/22 2501.511100 11039.
Wu, S., Zhou, W., Chen, M., and Zhu, S. Counterfac- tual generative models for time-varying treatments. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 2024.
Zhang, 7.,, Yang, W., Chen, 7., Chi, H., Wu, S., Zheng, W., Jin, R., Wang, B., Wang, Y., Huo, N., et al. A causal multiomics study discriminates the early immune features of adcevorrreed ebola vaceine recipients. The Innovation, 5(3), 2024.
Zhou, Z., Yang, X., Rossi, R., Zhao, H., and Yu, R. Neural point process for learning spatiotemporal event dynamics. In Learning for Dynamics and Control Conference, pp. 777-789. PMLR, 2022.
Zuo, S., Jiang, H., Li, Z., Zhao, T., and Zha, H. Transformer hawkes process. In International conference on machine learning, pp. 11692-11702. PMLR, 2020.

---

Supplementary Materials of Transformer-Based Spatial-Temporal Counterfactual Outcomes Estimation
A. Notations
A.1. Causal Model
Table 2. Notations for the causal model
Whole spatial region
One specific spatial region
7= (1, ..., T
Time index set
Spatial location Time
Number of time periods
Ze
Treatment in time Treatment up to time
Zst
70-M+1.1
Treatment between time t M and time t Treatment location in time
Sa
Y((Zs))
Potential outcome in time Observed outcome in time Outcome location in time t Covariates in time t
(=stt)
Sym"(zat)
Xt
2t
Covariates realization in time t
hg:
Observed historical information up to time Potential historical information up to time
Ha
A.2. Intervention
Table 3. Notation for treatment intervention
h
Intensity function of spatial point process
M
Duration treatment intervention Distribution spatial point patterns with intensity
Fn
Fi
Joint distribution of M independent Fn
Fi)(z)
Distribution of 2, is be assigned to Fh
2ct FH)
Distribution of 2[z-M+1.d in Zet is be assigned to Fи
A.3. Estimands
Table 4. Notation for interest estimands
(FH)
Expected number of outcomes that occur in time t and region w under intervention distribution Fø Average NY (Fn) over time M to time T
No(Fu)
12

---

Transformer-Bane Spatial-Temporal Counterfactual Outcomes Estimation
A.4. Estimators
Table 5. Notation for estimators
s) (Fn) (F(FH) (x(t)w))
Estimated intensity function of ye(2st (FR))
Estimator for N (F#)
Estimator for (FH
Intensity function spatial Poisson point process that gen- erates Y (Zse (F))
B. Running Example of the Estimands
Now we provide a running example for the estimands. For simplicity, consider the case of / = 8 and M = 3. Then the estimands
Ng (FH) Sgz (sss n (w)dE(((((s)))
represent the expected number of outcomes in t 8 and region w under distribution Fy. Next, we employ the Table 6 to elaborate on each term of the estimands.
Table 6. Interpretation of key terms the countcrfactual estimands.
Term
Interpretation
8= 33
Eyaluates outcomes at time 8, considering 3-times intervention persistence (times 6-8).
2[6.8]
Sequence of treatment variables over the time window [6,8].
(2(08))))
Joint probability distribution of 2[6.8] under counterfactual intervention strategy Fu. Observed outcome counts in region time under Fn-
(yg* nw 73
All possible values of 2[6.8]-
Expectation computation over all possible values of 2[6,8]-
C. Proofs
C.1. Assumptions
Assumption 1: Unconfoundedness. Papadogeorgou et al. (2022) proposed the unconfoundedness assumption for spatial- temporal data, in which they assume that given hgg. Ze is not dependent on any past or future potential outcomes and covariates. Their assumption quite strict, as the future is unlikely to affect the past, so we relax their assumption and assume that condition on hig., Ze dependent on Hz:
Ze 1l Healha
pk();) > 82. This assumption ensures that all treatment point patterns in any intervention distribution F₁ are also possible in the observed world. Papadogeorgou al. (2022) only use one treatment distribution FR. in their estimators. In contrast, used different treatment distributions in our estimators, so here we strictly assume unique constant 5, exists for any treatment distributions.
Assumption 2: Overlap. For any 2, E Zt, é and h, there exists unique constant 5, > 0, such that <<(*s)
Assumption 3.
Let denote the number of elements set, and Var denote the variance.
(a) There exists a constant 8y > such that (Syc (Ege)l < by for and Egt-
(b) Let Var[I] (x) 00.
(55) N (4)) there exists constant 0 such that T-M+I Σr P as
13

---

Transformer-Bsce Spatial-Tempont Counterfactual Outcomes Estimation
(c) (c) (Wy ASyphizad (s)ds NoYil < 3 and / is an infinitesimals.
(d) i 0l )
In Assumption 3, (a) means that the number of outcome events at any time has upper bound under any treatment. In our real-world scenario, we consider forest loss as the outcome event. Forest loss cannot be infinite, so (a) is reasonable. (b) assumes the convergence in probability of sequence. In (c), N₂ (Y,) is the actual number of outcome events in time and region w because our neural network smoothing can be viewed as an estimate the intensity function
point process.
We design so that the output values of the neural network are as large as possible where the outcome occurs, and as small as possible where the outcome does not occur, so (c) is reasonable. In (d), we assume that 3 tends to zero a faster rate than / i.e., BVT as T
C.2. Propositions
Proposition 1.
The propensity score a balancing score. For any / E 7
Ze (
Proposition 2.
Dimensional reduction property of the propensity score: if Z, II Hgllhgr then Z, a Hz)))
The consistency and asymptotic normality of the estimator. Let Var denote the Variance, No (Y)) denote the number of the observed outcome Y; in region W. If all assumptions hold, and T oo, have that
Proposition 3:
VT(N (FF) 4 N(O,)),
ΣT 2h where u limm- T-M+1
(25) and 0 Var[II] e,(o))
Nu Ho-)
C.3. Definition
Definition: Martingale Difference Series (Van der Vaart, 2010). Let (9,F,Pr) be probability space, filtration
F, (F,tt 0} is non-decreasing collection of a field on F, (e.g Fo c Ficc... c F; c... F). Let X, = (X,; 0} be a time series. A martingale difference series relative to a given filtration a time series X, such that, for any t:
(1) X₃ is Fr measurable.
(2) [[[ 0o.
(3) F =[
C.4. Theorem
Theorem: Central limit theorem for Martingale difference series (Van der Vaart, 2010). If X, is martingale difference series relative to the filtration Fe, and there exists constant > O, such that 4 ΣΣ E[ 2 v, and such for any O, then vnX, d NO,,)).
that Σバ EX 41xppe/
C.5. Proofs for Propositions
Proofs for Proposition We need show that IP 20le((2)), P(Z, 2))e((z))). Since es (2t) is a function of heg then PiZ. 2ele((1), hsa) (Z, 2g(hsr) c2(2)
P(Z, =x))(x))) E(P(Z, = =+)])))))
Based on the above, we prove the Proposition
14

---

Transformer-Basee Spatial-Temporal Counterfactual Outcomes Estimation
Proofs for Proposition 2. We need show that P(Z, 2, = (P =))e( Since ex(2)) is function of hgg and hige Hit then P(Z+ (e, (2 2x) and
P(Z. = 2, 21) > = P2 2n Hg. hgg))
IP( 2 kst)
(Unconfoundedness))
=((())
P( 2/ (=)))
Based the above, we prove the Proposition 2.
Proofs for Proposition 3. Let Err, = (F(FH) - No (FH) represent the estimation error at time t. We divide the estimation error Err, into two parts: the first part is Eic, which comes from the treatment allocation, and the second part is Ez., which comes from the spatial smoothing of our neural network. To be specific.
Erri= = II
pz)(2z) j-t-M-1 Cz(z;)
Soo ASyprcas (s)ds - (FH)) = Erc + E2
where
Phy (2g)
Eut= II
NYi -Y
j+t-M+1 cs(2,)
and
II (z)( Jo Xs Xsyppeen) (s)ds No (V.).
Ez1
M-t-M+1 eg(2;)
We will show that,
(i) VT( M-M+I ΣF FAt) 4. r),,
(ii) VT( T-M+I EL Ez) 0.
Proof of the asymptotic normality the first part of the estimation error. We use Definition C.3, the definition of the Martingale difference series, and Theorem C.4, the central limit theorem for the Martingale difference series to prove (i) ΣT M Eid) vO.
Lemma 1. Eit is martingale difference series with respect to the filtration F, = H M+1.
Proof for Lemma 1: We need to show that,
(1) En is HaL-M+1 measurable,
(2) F1e] 00,
(3) El [- E[E
It is easy to prove that (1) holds from the definitions of ELt and H From Assumption 2 and Assumption 3 (c) we have
(2)
Enl s
II
No. (0
co.
j-t-15+1 (2(2))
Therefore, E[[E (1 8. 0o, (2) is proved.
For the proof of (3), we will show that
E I (z)(2))
N. (.) N. (Fu),
j-t-M-1 ej(2,)
15

---

Transformer-Ban Baaeee Opatial-Temporal Counterfactual Outcomes Estimation
(,)
f I
Phy (2t) (Zs(((F))) P(Z--M-((Hs-M))
E II
No No(Y)) s+
j=t-M-1 c2(2))
j=t-M+1 eg(2g)
x
Ph, (2n) (x (Z5t
(2sc (FH))) M
J-f- Cg(2))
PP x 3F
f I (x
(Unconfoundedness))
j-t-M+1
= (Fr).
Then have
Phy No) EEN (FN)(FF))
[[En||asMM] = ET II
j-t-M+1 es(2,)
= N7 No
=0.
Based on the above, we prove the Lemma 1.
Lemma 2. T-M+I Σte -
0, for any 6 0.
Proof for Lemma 2. We will show
- as 7 - co. From Assumption 2 and Assumption 3 (6), we can obtain Eiel < dy 51 + 5-M). We show 61 + ) S M +
E' 8, M) VVV-MM -
(1 TT -M+++
M - dy 8, )|² sT,
let To - [1 (X when T To. we have öy (1 + M1 S M and Eid e M + 1. Therefore, have F1 >
- T Based the above, we prove the Lemma 2.
(=y(:)) (Y)) from Assumption 3(b), we have (,(2))
We have E[E V En [T
Σ
[[
2
T M +1
Combining Lemma 1, Lemma and the Central limit theorem for the Martingale difference series C.4, have
1
T Ev) d N(.v))
VT( T- M+
Σ
t=M
Based on the above, we prove the asymptotic normality of the first part of the estimation error.
Proof of the convergence in probability of the second part of the estimation error to zero. The second part of the estimation error represents the difference between the integral the outcome of the neural network smoothing and the actual number of outcomes. We will show
VT
E2) - O. f-
Σ
T
M+ +
16

---

Transformer-Basel Spatial-Temporal Counterfactual Outcomes Estimation
Ph (#s) (z,)
Let 04 [
then
T Σ t-M
Σ I=M
E2i = T M +1
arl (s (s)ds())ssss))) (s)ds-NN(Y)))
T M +
From Assumption 2, we obtain a, < 6, From Assumption 3(c), we obtain I As (s)ds ( Nw <3. Then we
have
T Σ E-M
Ezd
M=
1
ΣlS.
Σ
T
As
Asype (s) No
T M
T-M +1
EM
1
T
$ 8 M
Σ
ASypeies
T M
t-M
1
8 M
(T- M 1)3
T M
= 8 M B.
Because 3 can be arbitrarily small so 8, M 3 can also be arbitrarily small, combined with Assumption (d), we conclude that (T(T-M-1 Σt-M Ez.) O.
Combining the proof of the asymptotic normality of the first part of the estimation error and the proof of the convergence to - (FH)) 4 (0, v) as T + co. 0 of the second part of the estimation error, we prove (F))
D. Identifiability Proof
This section provides proof of identifiability for the estimands. Under the potential outcomes framework, identifiability refers to whether estimand can be represented using observable data. The formal definition is given below:
Definition D.1. (Identifiabilitty) parameter @ is said to be identifiable it can be expressed as function of the distribution of observed data under certain assumptions. The parameter 0 said to be nonparametrically identifiable if can be expressed as such function without any parametric assumptions the model.
In the above definition, the parameter 0 refers the estimation target, which could be causal effect counterfactual outcome. Identifiability is crucial in observational studies because an estimand must be identifiable in order to be estimable from observational data. If the estimation target depends on unobservable data, then it is non-identifiable and cannot be estimated. To derive the identifiability of the estimation target in this work, we introduce the following reasonable assumptions:
Assumption (Ignorability). We assume that, given the observed history hgg, the treatment variable Z. is independent the latent history Her:
Z. Hgl hige
(11)
Note that Papadogeorgou al. (2022) assume that Ze is independent of all past and future potential outcomes and covariates given hit- Their assumption is overly striet, the future cannot causally affect the past (due to the unidirectional nature of causality in time). Therefore, we relax their assumption by only requiring Z, to be independent of the latent history His.
Assumption (Poisson Assumption). We assume that Z₁, (Z))). and Z, ha-11 are generated by spatial Poisson point processes. This assumption reasonable because spatial Poisson point processes are widely used statistical tools to model the distribution of discrete events over spatial regions (Cressie Wikle, 2015).
Assumption (Consistency). At any time and in any region, if a treatment sequence øςr is assigned, then the observed outcome yo must equal the potential outcome under that treatment, i.e.,
y. Yi(Ze: 2ct).
(12)
This assumption ensures that the observed outcomes are consistent with the potential outcomes, enabling the estimation of counterfactual outcomes and causal effects.
17

---

Wransformer-Basee Spatial-Temporal Counterfactual Outcomes Estimation
Theorem D.2 (Identifiability of the Estimation Target). If Assumptions Ignorability, Poisson, and Consistency hold, then the
estimands No( (FF) is identifiable.
Proof. By the definition of N2 (FH), it suffices to show that (2s) is identifiabl..
E[y," (zsa (F#))] E[E[Y," (2tt
(2s( (FH)) har]]
(13)
(Fe hhets Zet ~ FH]]
(Fe))
(14)
=E[E FH)
hets Zet ~ FH]]
(15)
E[ legs Ass [H]],
(16)
Equation (13) follows from the tower property of conditional expectation, (14) is derived using the ignorability assumption (since y (zge (FH)) e Hst). and both (15) and (16) follow from the consistency assumption. The Poisson assumption ensures that the distributions of the potential outcomes and treatments are known, further reinforcing the identifiability of the estimation target.
The above establishes the identifiability of the estimands, which guarantees the solvability of the proposed estimation problem.
E. Derivation of Property of Poisson Point Process
Let X(s) denote the intensity function of a Poisson point process, A denote a spatial region, (AA) denote the number of events occurring within region A, and [[]( denote the expected number events occurring within region A.
According to the property of the Poisson point process, we have the following:
(A) Poisson(AA).
where Poisson(-) denotes the Poisson distribution and Xa its parameter.
Then
20 e ca As 8 Xi1 E[N((A] - ((((A))) i1) i=0
20
AsE
Σ F=
= dse = AA
i
(i-
According to the definition of the intensity function we have:
E[N(A)]= f A(s)ds XA.
F. Derivation of the Training Objective Function
The objective function Eq. (10) is mainly derived from the likelihood function of the Poisson point process and the objective function of the variational autoencoder (VAE). Specifically, we employ the idea maximum likelihood estimation (MLE), and the likelihood function of the Poisson point process is shown below:
Σ Loog(A (s:) (- (u)
MIE for the intensity function of the Poisson point process seeks the optimal intensity function X²(-) from the data that optimizes the above function. As for the (I((llp) in Eq. (10), it's a common component of VAE. Therefore, our method can seen as a combination of MLE and VAE. We use a neural network as the intensity function and seek the optimal intensity function by optimizing the objective function Eq. (10).
G. Explanation of the Treatment Intervention
In this section, we provide brief explanation of why the intervention treatment duration M has an impact on the results. In our problem setting, we assume treatments have delayed effects, i.e., temporal carryover. For example, an educational program results in poorer performances at first, as students adjust to new learning methods, but eventually, they may get long-term improvements. Therefore, under our spatial-temporal setting, the larger M may influence the estimation.
18

---

Transformer-Basee Spatial-Temporal Counterfactual Outcomes Estimation
H. Explanation of the term "Counterfactual Probability"
this section, we provide a brief explanation of why the counterfactual probability is "counterfactual.. Recall that pa (2t) the probability density function of intervention distribution (.) In this work, we aim to estimate the number of potential outcomes when the treatment follows an intervention distribution. We specify this intervention distribution, which does not necessarily exist in the observable data. For example, in the observable data, treatment may follow distribution Pa. and we want to find out what will happen to outcomes when treatment follows another distribution PB. Therefore, the counterfactual probability Ph (z²) is counterfactual and it's counter to the observable data or real-world.
I. Additional Experiments and Their Results
Ours
Causal Forest
LR
MSMs
RMSNs
T=64, M=3
T=48, M=3
T=32, M=3
RER (Relative Ereer Rate)
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0g
6
0.0g
Values of
Values of
Values of
Figure Experiments results of The horizontal axis represents the values of c, while the vertical axis represents relative error rate (RER). The lower lines in the graph correspond to methods with higher estimation accuracy. From left to right, the three columns respectively represent the experimental results with time lengths of 64, 48, 32 64, 48,32).
According to Figure 6, in the different experiment settings of M=3, the estimation errors of our method are relatively low in most cases, indicating the robustness of our method concerning time lengths and intervention settings.
J. Details of the Ablation Studies
J.1. Details of the RNNs
J.1.1. PARAMETERS OF THE RNNS
We employ the Gated Recurrent Unit (GRU) network. The input_size = 32, hidden_size = 32, num_layers = 3, dropout=0.1. batch first True.
J.1.2. COMPARISON RESULTS OF RNNS
64, M=3
T=48, M=3
0.7
RER (Re ative Error Rabe)
RER (Re attvettvv Error Rabe)
2.6
0.5
0.4
0.4
0.3
0.
0.1
5.0
Values of
Values
(a) Comparison results of T=64, M=3.
(b) Comparison results T=48
19

---

Transformer-Base Spatial-Tempores Counterfactual Outcomes Estimation
T-48, M=1
T-32, M-1
0.6
Transfomer
0.6
Transormer
RER (Relasive Erer Rate)
0.5
RNN
RER (Relative Error Rase)
RNN
0.4
0.2
0.1
0.3
0.0
Values of
Values of
(a) Comparison results of T=48, M=1.
(b) Comparison results of T=32, M=1.
T=32, M=3
RER (Relative Error Rate)
0.2
0.0
Values of
Figure Comparison results of T=32, M=3.
According to the above comparison results with RNNs, the Transformer backbone is superior to the RNNs in most settings.
J.2. Details of the Relaxation of Poisson Assumption
J.2.1. SETTINGS OF INTENSITY FUNCTIONS
To relax the Poisson assumption, we add the Gaussian kernels to the intensity functions (see Appendix M.1) that generate the
12
synthetic data. Specifically, replace the (+)( 222 (8) with Gaussian kernel of .iz(s) e -2v,
D7
02 (0-,1 (a) 20
with and (0) 2D (Dzg with
The is set constant of V2"
J.2.2. COMPARISON RESULTS
T=64, M=3
4=48, M=3
RER (Relative Error Rate)
Poisson
RER (Relative Error Rate)
Poisson
Gaussian
Gaussian
0.20
0.20
0.15
0.15
3.10
3.05
0.05
0.00
3.00
Values
Values c
(a) Relaxing results of T=64, M=3.
(b) Relaxing results T=48, M=3.
20

---

Transformer-Based Spatial-Temports Counterfactual Outcomes Estimation
T=48, M=1
T=32, M=1
0.34
Poisson Gaussian
RER (Relative Error Rate)
Poisson
0.22
Gaussian
0.12
0.30
RER (Relative
0.06
0.00
0.04
0.06
0.52
0.04
0.02
Values of
Values of
Relaxing results of T=48, M=1.
(b) Relaxing results of T=32, M=1.
RER (Relative Emro Rate)
0.25
Gauasian
0.15
0.05
Values of
Figure 12. Relaxing results of T=32, M-3.
K. The Design Details of Convolutional Neural Networks
3 MaxPool2d 64
3+3 Conv2d 64
3 Conv2d 64
3 Conv2d 128
3 Conv2d 128
3 Conv2d 256
3*3 Conv2d 256
3*3 Conv2d 256
3*3 Conv2d 256
3 Conv2d 256
input
output
PC
Figure 13. The design details of Convolutional Neural Networks.
L. Details of the Real Dataset
Conflicts in 2002
Forest Loss in 2002
Figure 14, An example of the real datasct,
21

---

Transformer-Bane Based Spatial-Temporal Counterfactual Outcomes Estimation
An example of the real dataset is shown in Figure 14, In Figure 14, the blue lines represent the roads in Colombia. The left part of the figure shows the conflicts in 2002 in Colombia, the red dots represent the conflict locations. The right part of the figure demonstrates the forest loss in 2002 in Colombia, the green dots denote the forest loss locations.
L.1. Forest Change Data.
We only consider the "Year of gross forest cover loss event (loss year)" section of the Global Forest Change dataset, which is raster data matrix (.tif file). The matrix elements are O, representing representing forest loss events occurring from 2001 to 2022 (Hansen et al., 2013). We selected a portion from the dataset that occurred within the territory of Colombia. We read the raster data matrix's 's elements from 2 to 22, calculating their latitude and longitude positions. As a result, we obtained the forest loss events that occurred Colombia from 2002 to 2022. Data is publicly available online from https: //g wel.g..aaaaaaaeee eerthengiiennee app/view www.vvwwwwwwlaaaa. global forest Forest-change change.
forest loss events, and values from to 22,
L.2. UCDP Georeferenced Event Dataset.
UCDP Georeferenced Event Dataset exists as tabular data (.xlsx file) (Croicu & Sundberg, 2015). We only consider the years in which conflicts occurred, the latitude and longitude coordinates of the conflict locations, and the country or region where the conflicts occurred, as contained in the tabular data. We directly filtered conflict events within Colombia from the tabular data, spanning from 2002 to 2022. Data publicly available online from https: /ucdp uu. se/ dwww/wwwlll oads/ index. htmlf i leredggobball
L.3. Road Data.
We select several major roads within Colombia from Google Maps. Specifically, for convenience, choose key points from the main roads of Colombia and connect these points to form polyline, representing the road. The coordinates of selected key points are shown below:
Roadl: (-77.66, 0.77), (-76.42, 3.09), (-76.06, 4.29).
Road2: (-76.06, 4.29), (-74.03, 4.66), (-73.04, 7.05), (-72.38, 7.75).
Road3: (-76.02, 4.31), (-75.54, 6.15), (-73.11, 7.09), (-74.22, 10.97), (-72.23, 11.33).
Road4: (-75.52, 6.15), (-74.69, 10.93).
Road5: (-75.52, 6.22), (-76.77, 8.42).
The first part of the coordinate is the longitude, and the second part is the latitude.
M. Synthetic Data Generating Process
We utilize reject sampling (Lavancier al., 2015) to sample spatial point patterns of treatments and outcomes based their intensity functions. Below, we introduce the setting for intensity functions.
M.1. Intensity Functions Setting
We consider four covariates X1 (s) and Xg(s) and X3(s) and Xa(s). Xi(s)= =3- +log(D2ss)) X2(s)= -3ds( Ds(s) is the distance from the location to the road on spatial area. D2(s) is the distance from s to the border of the spatial area. Ds(s) is the distance from s to the center of the spatial area. X x=x(s) (ft) and Xxx (1)+ are two intensity function,
the point patterns generated by them create (s) and X*().. Specifically, az at, ag and af. are bias constant. Let (s(s) represent the distance from to the closet points generated by (x) Da(s) represent the distance from to the closet points generated by X (s). Then Xz(s) Da(+), X Xa(s) Pa(s))
Based on all covariates, we determine the intensity function of the Poisson point process that generates treatment and outcome. Let (s) = (x) (x()). Xs(s), (s) denote the covariates vector, and Z;- (s) = e 2D2,, (+), (s) 2Dv, (+), Dz (8) ) the distance from s to the closet point in Sz,,, Dy (s) the distance from s the closet point in Sy The
22

---

Transformer-Bant Based Spatial-Temporal Counterfactual Outcomes Estimation
intensity function for the Poisson point process that generates Z, shown as follows:
Az, - x
(17)
In Eq. (17), Bo. 9x, 32, By, constant parameters. Let Z [t-3..t 2D p-s.n (6) Dz, 2..(() is the distance from
to the closet point in - Sz The intensity function for the Poisson point process that generates Y, is shown as follows:
Xy, (s) = 20 (())
(18)
In Eq. (18), 7o. 7x. 72. YY, are also constant parameters. the synthetic data generating process, we set 30 3x (1,1,1,1), 82 By 1, and fo 7x (1,1,1,1), y2 1, NY ag -0.2, af 2.3, a -0.2 and at 2.8.
M.2. Ground Truth Generation
For the computation of the true counterfactual outcomes, we first employ the method described in Appendix M.1 to calculate the intensity function of the spatial Poisson point process generating Y: (=st (F)). Subsequently, we utilize this intensity function to generate samples of Y, (Zsa(FH) Finally, the average of samples computed the true No (F#).
M.3. Synthetic Data Region
In synthetic data, the entire area rectangle with a length of and width of 1. For ease of computer processing, we divide this area into 100 squares. We draw some lines the rectangle, which consider roads. The synthetic data area is shown in Figure 15. In Figure 15, the red line represents the straight road, while the green dashed line represents the curved road. Based on the synthetic data region, we calculate the intensity function of the treatments and outcomes. While it is possible to generate an arbitrary number of intensity functions, due to computational limitations, we have generated 32, 48, and 64 intensity functions for the treatments and outcomes, corresponding to time lengths of 32, 48, and 64. Afterward, based on the intensity functions, we use rejection sampling to generate point patterns of treatments and outcomes. To be specific, twenty spatial point patterns are generated from each intensity function.
10
0.8
0.4
0.2
0.4
0.6
1.2
Figure 15. The synthetic spatial area. The red lines and green lines in figure represent synthetic roads.
N. Details of Baselines
N.1. Baselines Adaption
Now we introduce how we adapt baselines to our setting. Since baselines cannot directly handle high-dimensional data such as series of spatial point patterns, we transform all treatments (2)) and outcomes (Y₂) into the number of events contained in the treatments and outcomes (R(z)) and R(Y.)), R() is the dimension reduction map defined in Eq. (5) in the main text.
23

---

Transformer-Basee Spatial-Temporal Counterfactual Outcomes Estimation
Therefore, the data used for baselines are all scalar series. For baselines like MSMs, RMSNs, and Linear Regression (LR), fit scalar outcomes on treatments and covariates and make comparison with method. For Causal Forest, we first use linear regression model to fit scalar treatments outcomes, and then train the causal forest model to estimate the treatment effects. Finally, we combine the causal forest and regression model to build the counterfactual outcomes.
N.2. Marginal Structural Models (MSMs)
Marginal Structural Models (MSMs) (Robins et al., 2000) is the statistical tool used in observational studies to estimate causal effects, especially when dealing with time-varying exposures and confounders. MSMs employ the standard regression method as a base estimation model and adjust for these complexities using methods like inverse probability weighting (IPW) or g-estimation. In our configurations, we employ MSMs to outcomes on treatments and covariates and make a comparison with our method.
N.3. Recurrent Marginal Structural Networks (RMSNs)
Recurrent Marginal Structural Networks (RMSNs) (Lim, 2018) is a deep learning-basce method to forecast counterfactual outcomes. Different from the MSMs, RMSNs employ the RNN model to build eequence-t--sequence architectures for counterfactual outcome prediction. To be specific, (Lim, 2018) employs RNNS to construct propensity networks and prediction networks, and combines these modules form the inverse probability of treatment weighting (IPTW) estimation. In our configurations, we utilize RMSNs to fit outcomes on treatments and covariates.
N.4. Causal Forest
Causal Forest (Wager & Athey, 2018) is random forest-based method for heterogeneous treatment effects estimation, the treatment effects are estimated at the leaves of the random trees. We employ Python library called EconML https: / wwc.oma.....eeebb.o...ooo azurewebs i tes tee.ttts...eeeeeeeeeeee / to realize the Causal Forest. Although EconML. does not support the estimation of the counterfactual outcomes directly, for most estimators in EconML, we can combine a baseline predictive model with one estimator in EconML to construct the counterfactual outcomes estimation. Specifically, we use a regression model to fit treatments on outcomes, and then train the causal forest model to estimate the treatment effects. Finally, we combine the causal forest and regression model build the counterfactual outcomes.
N.5. Linear Regression (LR)
We develop a linear regression-basee method as an additional baseline. Since LR cannot directly handle matrix data such as spatial point patterns, we transform all treatments (2t) and outcomes (Y₂) into the number of events contained in the treatments and outcomes (R(z)) and R(Y.)), R() is the dimension reduction map defined in Eq. (5) in the main text. Therefore, the data used for LR are all scalar. Subsequently, we use R(21), R(22)..., ...., R(2M) regress R(YM). and the regression model obtained is denoted as E[R(YM) M [R(2M)) After obtaining the regression model, we replace the independent variables, R(21), R(22)..... R(2M) in the regression model E[R(Y,/ M) [R(ccM) ] with c * log( (og(R(z)) R(21)) )) C * log( ((s) )). and then input them into the model to obtain the predicted values No (Fn). Similarly, we can obtain estimates at time M 1. We use R(21), R(22).... R( ((2++)) to regress R(YM+1), and the regression model obtained is denoted as E[R( ([[((*+++))))) (((ss=)))(((=MM)))) After obtaining the regression model, we replace the independent variables, R(21), R(=2).... 1() in the regression model E [[R(Y, e* log R22++))), and then input them into the model to obtain the predicted values N (FF)). In this way, we repeat the process until we obtain the estimate at the final time T, NY (FF)). Then we can obtain estimates in the main text: No. (FH))= T-M-I ΣLM (FH). We use the torch.nn module in PyTorch to implement linear regression, with 10 epochs and a learning rate set to 0.00004.
e* log(
(R(Ya++))) 1) [(== with c * log( R( (()))), ), log( R( ((222)) )....
O. Limitations
In this section, we discuss the limitations of our work. We employ the unconfoundedness assumption to identify our estimands and prove the propositions. The unconfoundedness assumption excludes the influences of unmeasured covariates. In real-world scenarios, unmeasured covariates may exist thus this assumption may be violated. However, the unmeasured covariates are still an open problem (Pearl, 2010), and the unconfoundedness assumption is still widely used in causal inference. We should be open to the unconfoundedness assumption.
24


Reviews:
[{"review_id": "ixV1UwlKy3", "reviewer": "Reviewer_p6Zj", "summary": "The paper introduces an approach to estimate counterfactual outcomes in a spatial-temporal setting, where both treatment and outcome may be represented in a high-dimensional space. The proposed method adapts IPW to the spatial-temporal setting, leveraging propensity score. The approach is implemented in practice using a deep learning architecture, leveraging CNNs for computing propensity scores and transformers for encoding the intensity function. Finally, the method is compared empirically with traditional counterfactual outcomes estimation approaches and achieves significantly lower error in high-dimensional settings.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response to reviewer p6Zj \n\nWe greatly appreciate the reviewer's comments and suggestions. We would like to respond to each detailed point individually.\n\n> 1. The estimator estimates the expected number of outcomes in a specific region at a specific time. This estimator (Eq. 3) seems to be the main contribution of the paper, but it is unclear why this is the quantity of interest and how the estimator is derived.\n\nNow we explain why the expected number of outcomes in a specific region at a specific time is the quantity of interest. These estimands provide spatial-temporal decision-making evidence for policy makers. For example, during epidemic control, policy makers can evaluate the average number of infections in specific areas at specific times under different isolation strategies, thereby optimizing the allocation of isolation resources.\n\nNext, we briefly introduce how the estimator is derived. The estimator is derived based on the Inverse Probability Weighting (IPW) strategy in Marginal Structural Models (MSMs) [1]. Specifically, the $\\lambda_{Y_t^{ob}(z_{\\leq t})}(s)$ in Eq.3 is the intensity function of outcomes in time $t$ and location $s$, and $\\prod_{j=t-M+1}^{t} \\frac{p_{h_j}(z_j)}{e_j(z_j)}$ represents the weights similar to those in the MSMs. In summary, we adopt the weighting strategy of MSMs to the studied spatial-temporal setting. \n\n[1] Robins J M, Hernan M A, Brumback B. Marginal structural models and causal inference in epidemiology[J]. Epidemiology, 2000, 11(5): 550-560.\n\n---\n\n>2. The proofs seem to hold, but it seems generally unclear how to interpret these results and why they are important. Prop. 3 seems to be the most important as it seems to indicate that the estimator properly converges to the intended value, but there seems to be some subtleties in the result that could be elaborated.\n\n\nProp. 1 and Prop. 2 are two essential properties of the propensity score. Specifically, these propositions demonstrate that the propensity score can make the estimation more **accurate and efficient**. Therefore, by proving Prop. 1 and Prop. 2, we validate that the propensity score used for spatial-temporal data is \"reasonable\" and can assist in counterfactual outcome estimation.\n\n\nProp. 3 establishes the **theoretical reliability** of the estimator. It demonstrates that the expected value of the proposed estimator equals the estimands (i.e., unbiasedness) and that the variance of the estimation error converges to a finite value.\n\n---\n\n>3. I appreciate that the assumptions are clearly highlighted.\n\nWe are glad the highlighted assumptions are helpful.\n\n---\n\n>4. I think that the clarity of the paper is its greatest weakness. It would be much clearer if the paper walked through a running example, where each term of the estimand is clearly highlighted.\n\nWe appreciate the reviewer's suggestion on the running example. Following your advice, we provide an example that illustrates the estimands.\n\nFor simplicity, consider the case of $t=8$ and $M=3$. Then the estimands \n$$\nN^\\omega_8(F_H) =  \\int_{Z^3} |S_{Y_8^{ob}(z_{\\leq 8}(F_H))} \\cap \\omega| dF_H(z_{\\left[ \\text{6,8} \\right]}),\n$$\nrepresent the expected number of outcomes in $t=8$ and region $\\omega$ under distribution $F_H$. Next, we employ a table to elaborate on each term of the estimands.\n\n\n| Term | Interpretation |\n| --- | --- |\n| $t$=8, $M$=3 | Evaluates outcomes at time 8, considering 3-times intervention persistence (times 6-8). |\n| $z_{[6,8]}$ | Sequence of treatment variables over the time window [6,8]. |\n|$F_H(z_{\\left[ \\text{6,8} \\right]})$  | Joint probability distribution of $z_{[6,8]}$ under counterfactual intervention strategy $F_H$.|\n| &vert;$S_{Y_8^{ob}(z_{\\leq 8}(F_H))} \\cap \\omega$&vert; | Observed outcome counts in region $\\omega$ at time 8, under $F_H$.  |\n| $Z^3$ | All possible values of $z_{[6,8]}$.  |\n| $\\int_{Z^3}\\cdot dF_H(z_{\\left[ \\text{6,8} \\right]})$  | Expectation computation over all possible values of $z_{[6,8]}$. |\n\n---\n\n>5. Why is the dimensionality reduction in Eq. 5 possible? Are we not losing information?\n\nIn this work, we mainly estimate the outcome quantities. Therefore, we believe it is possible to reduce the studied variables to their quantity. We acknowledge that this strategy may result in some information loss; however, it is **useful for estimating the number of outcome events**. Moreover, this strategy **has its own theoretical significance** in terms of how to reduce the dimensionality for complicated problems.", "claims_and_evidence": "The paper makes three claims: (1) they propose an estimator for counterfactual outcomes in the spatial-temporal setting; (2) they use a deep-learning approach to perform this estimation, leveraging CNNs for propensity score calculation and transformers for modeling intensity functions; and (3) the approach is demonstrated to outperform baselines empirically. See below sections for the evaluation of these claims.", "methods_and_evaluation_criteria": "The data is separated into treatment $Z$, outcome $Y$, and covariates $X$. Importantly, each variable is collected across several time steps, and $Z$ and $Y$ are modeled as spatial point processes rather than traditional binary treatments and outcomes. An inverse probability weighting estimator is derived, composed of the propensity score, the counterfactual probability, and the intensity function of the outcome. The estimator estimates the expected number of outcomes in a specific region at a specific time. This estimator (Eq. 3) seems to be the main contribution of the paper, but it is unclear why this is the quantity of interest and how the estimator is derived.\n\nThe estimator is implemented in practice using a CNN to compute the propensity score and a transformer to compute the intensity function. The architecture seems to make sense for the task.", "theoretical_claims": "The paper proves three propositions in Sec. 4.5 about the estimator in Eq. 3. The proofs seem to hold, but it seems generally unclear how to interpret these results and why they are important. Prop. 3 seems to be the most important as it seems to indicate that the estimator properly converges to the intended value, but there seems to be some subtleties in the result that could be elaborated.", "experimental_designs_or_analyses": "The proposed architecture is tested in both synthetic and real data, compared to baselines which do not take the spatial-temporal nature of the data into account. The results seem to convincingly show that the proposed architecture works better in these data settings.", "relation_to_broader_scientific_literature": "The related works section adequately sets the background of the task.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "I appreciate that the assumptions are clearly highlighted.\n\nI think that the clarity of the paper is its greatest weakness. The paper uses a lot of convoluted notation and does not explain much of the math. It would be much clearer if the paper walked through a running example, where each term of the estimand is clearly highlighted.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "Why is the dimensionality reduction in Eq. 5 possible? Are we not losing information?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gk9L7BgdMZ", "reviewer": "Reviewer_pyND", "summary": "This paper introduces a Transformer-based framework for counterfactual outcome estimation in spatial-temporal data. It aims to improve causal inference in settings where treatments and outcomes are structured across both space and time. The authors propose a novel deep-learning-based estimator with a CNN-Based Propensity Score Estimation and Transformer-Based Intensity Function Estimation. The framework was supported by both theoretical insights and experiments on synthetic and real data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response to reviewer pyND\n\nWe are glad the reviewer found the theoretical framework rigorous. We would like to respond to each detailed point individually.\n\n>1. The theoretical evidence is adundant, whereas the empirical results lack sufficient ablation to support the claims such as the improvements from the CNN based propensity score estimation.\n\nWe really appreciate that the reviewer found the theoretical evidence abundant. Now we explain why removing the propensity score estimation is inappropriate in our setting.\n\nIn our work, the CNN-based propensity score estimation serves as a practical implementation of the propensity score in our estimator framework (Sec. 4.4). Therefore, if we remove it, we effectively eliminate the propensity score from the estimator framework. As a result, the estimator would be **theoretically incomplete** and unable to perform its intended function properly.\n  \n---\n\n>2. Methods do make sense, but the evaluation criteria is lacking.\n\nWe are glad the reviewer found the methods reasonable. Now we justify our choice of evaluation criteria. In the simulation experiments, we use the relative error rate as a direct measure of the estimator’s accuracy. In real-world data experiments, since the ground truth is unavailable, we assess our approach by comparing our findings with those of existing studies.\n\n---\n\n>3. Spatio-temporal causal inference is very important and broadly applicable.\n\nWe are glad the reviewer found the studied problem important. \n\n---\n\n>4. I would encourage including more references to the temporal counterfactual prediction literature as it is an active research area, including (but not limited to).\n\nWe really appreciate the reviewer's suggestions on the references. Following your advice, we **have cited and discussed** these excellent works in our paper. \n\n**Temporal Counterfactual Prediction**\nTemporal counterfactual prediction refers to performing counterfactual outcome prediction under time-varying settings. Seedat, N. et al. [1] propose TE-CDE, a neural-controlled differential equation approach for counterfactual outcome estimation in irregularly sampled time-series data. Wu, S. et al. [2] introduce a conditional generative framework for counterfactual outcome estimation under time-varying treatments, addressing challenges in high-dimensional outcomes and distribution mismatch. El Bouchattaoui, M. et al. [3] propose an RNN-based approach for counterfactual regression over time, focusing on long-term treatment effect estimation.\n\n[1] Seedat, N. et al. Continuous-time modeling of counterfactual outcomes using neural controlled differential equations. \n\n[2] Wu, S. et al. Counterfactual generative models for time-varying treatments.\n\n[3] El Bouchattaoui, M. et al. Causal contrastive learning for counterfactual regression over time. \n\n---\n\n>5. The details of the transformer architecture seems to be missing: what positional encoding are you using? what tokenization scheme is used?\n\nWe employ the Transformer to capture spatial information within a single time step. As a result, we do not apply fixed positional encoding.  \n\nFor the tokenization scheme, the Transformer processes discrete data points as input. Therefore, we treat each discrete coordinate as a token and map it to a high-dimensional embedding space.\n\n---\n\n>6. More ablation is needed to identify the contribution of each innovation, including the CNN-based propensity score estimator.\n\nWe appreciate the suggestion of the ablation of the CNN-based estimator. However, the ablation of the CNN-based estimator is inappropriate in our setting, for detailed reasons, please refer to response 1.\n\n---\n\n>7. It is unclear if the performance gain come from model backbone (i.e. transformer) or the proposed spatio-temporal causal inference framework. One way is to keep the model backbone simple for similar to other baselines, and adopt several ablation studies to demonstrate the gains brought by the two major innovations.\n\nCompared to previous baselines, our work primarily tackles a problem that these methods cannot handle. This is because the frameworks and implementations of these baselines are not applicable to spatial-temporal data. Therefore, we attribute the superior performance of our method to both its theoretical framework and model implementation.  \n\nWe appreciate the reviewer's suggestions on the ablation studies. However, it is not appropriate to maintain the same model backbone as other baselines. For example, we employ a CNN to process the high-dimensional tensor. If we replace the CNN with the simple logistic regression or linear regression used in MSMs, the propensity score in our work would not be estimable.", "claims_and_evidence": "The theoretical evidence is adundant, whereas the empirical results lack sufficient ablation to support the claims such as the improvements from the CNN based propensity score estimation.", "methods_and_evaluation_criteria": "Methods do make sense, but the evaluation criteria is lacking.", "theoretical_claims": "I didn't check all the details. The proposals and theorems look good.", "experimental_designs_or_analyses": "See below.", "relation_to_broader_scientific_literature": "Spatio-temporal causal inference is very important and broadly applicable.", "essential_references_not_discussed": "I would encourage including more references to the temporal counterfactual prediction literature as it is an active research area, including (but not limited to):\n\nSeedat, N., Imrie, F., Bellot, A., Qian, Z., & van der Schaar, M. (2022). Continuous-time modeling of counterfactual outcomes using neural controlled differential equations. arXiv preprint arXiv:2206.08311.\n\nWu, S., Zhou, W., Chen, M., & Zhu, S. (2024). Counterfactual generative models for time-varying treatments. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (pp. 3402-3413).\n\nEl Bouchattaoui, M., Tami, M., Lepetit, B., & Cournède, P. H. (2024). Causal contrastive learning for counterfactual regression over time. Advances in Neural Information Processing Systems, 37, 1333-1369.", "other_strengths_and_weaknesses": "Strengths\nOne thing that stands out is the consistent and rigorous theoretical framework of the work, which can potentially lay the foundation for future works on spatio-temporal causal inference.\n\nWeaknesses\nThe experiment section is lacking, especially:\n1. The details of the transformer architecture seems to be missing: what positional encoding are you using? what tokenization scheme is used? \n2. More ablation is needed to identify the contribution of each innovation, including the CNN-based propensity score estimator.\n3. Following my last point, it seems that when comparing to other baselines, it is unclear if the performance gain come from model backbone (i.e. transformer) or the proposed spatio-temporal causal inference framework. One way is to keep the model backbone simple for similar to other baselines, and adopt several ablation studies to demonstrate the gains brought by the two major innovations.", "other_comments_or_suggestions": "No", "questions_for_authors": "See above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "p5gwKJcMtT", "reviewer": "Reviewer_Uw2a", "summary": "The paper studied counterfactual outcomes with spatial-tempora attributes using transformers and proved consistency and aymptotic normality of the estimator. The authors also conducted synthetic experiments and studied forest loss in Colombia. This paper is generally well-written.\n\n*Edit after rebuttal: I changed score from 2 to 3.*", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response to reviewer Uw2a\n\nWe are glad the reviewer found our paper well-written. We would like to respond to each detailed point individually.\n\n>1. I failed to be convinced that there is significnat novelty on the theoretical front -- the proof in Appendix B2 is thorough, but Proposition 3 seems to be more like a routine proof. To be fair, the authors did not claim theoretical contributions, so this is fine by me. \n\nWe appreciate that the reviewer found our proofs thorough. Now, we elaborate on our key contributions. \nOne of our **key contributions** is leveraging deep learning to address the challenge of counterfactual outcome estimation for spatial-temporal data. Previous baselines are **unable** to achieve this. Specifically, our proposed method enables the estimation of propensity scores for spatial-temporal data. This **provides a feasible and efficient technical pathway** for causal inference in high-dimensional spatial-temporal settings. Besides, our theoretical analysis provides a guarantee for the proposed deep learning method.\n\n---\n\n>2. However, I think the lack of validation for the non-synthetic studies makes the story less convincing. I understand that there is non \"true answer\" for the non-synthetic studies, but the authors could perhaps perform sensitivity analysis or convergence analysis of their estimators.\n\nWe appreciate the reviewer's suggestions on sensitivity analysis and convergence analysis. Following your advice, we provide the theoretical sensitivity analysis in the **\"sensitivity.pdf\"** file in this anonymous [URL](https://anonymous.4open.science/r/DeppSTCI_Release_Version-master-3C91). \n\nWith regard to the convergence analysis, Proposition 3 demonstrates that as time progresses, our estimator converges to the estimands (ground truth).\n\n---\n\n>3. I think it is helpful that the authors discuss the computation cost but a comparison of computation cost with the baseline might shed more light on the comparison between the proposed method.\n\nWe appreciate the suggestion to compare the computational cost with the baselines. All baselines were also run on an NVIDIA RTX 4090 GPU with an Intel Core i7-14700KF processor. For detailed information, please refer to the table below.\n\nNote that these baselines were **not designed for spatial-temporal data** and cannot process such data. Therefore, we convert all spatial-temporal data into scalar values for processing by these baselines (Appendix K.1). As a result, although our method does not significantly outperform these baselines in computation time, it **addresses a problem that they cannot solve.**\n\n| Methods | Time required for one setting |\n| --- | --- |\n| Ours | about 10 minutes |\n| LR | about 1.13 minutes |\n| MSMs | about  1.45 minutes |\n| RMSNs | about 1.53 minutes |\n|Causal Forest| about 2.25 minutes |\n\n---\n\n>4. I think, despite this paper is generally well-written, the constribution seems to be a bit weak -- there is no significnat contributions on the theory front. However, the empirical studies and analyses are related limited compared with what I would expect from an application paper.\n\nWe appreciate that the reviewer found our work well-written. Our main contribution lies in **addressing the challenge** of counterfactual outcome estimation for spatial-temporal data, which previous baselines are **unable to solve.** We provide a **feasible and efficient technical pathway** for causal inference in high-dimensional spatial-temporal data. The theoretical section provides the **foundation** for the proposed method.\n\n---\n\n>5. In the figures/tables of empirical studies, perhaps adds error bars and std/confidence bands.\n\nWe appreciate the suggestion to include error bars and confidence intervals. However, in the additional experiments, we trained the neural networks 25 times with different training data, and the results were consistently superior. Therefore, we believe the additional experiments demonstrate the robustness of our method.\n\n---\n\n>6. Did you try to relax the Poisson assumption and see how well this method holds up in the synthetic setup?\n\nWe acknowledge that the Poisson assumption may be somewhat restrictive. In this work, we use the spatial Poisson point process as a fundamental hypothesis for modeling spatial-temporal data. Therefore, if the Poisson assumption is violated, our theoretical framework would need to be restructured, and its performance may be compromised. In statistics, modeling spatial data distributions without assuming a specific distribution remains an open problem, and the Poisson assumption is widely used to describe spatial events [1-2].\n\n[1] Cressie, N., & Wikle, C. K. (2011). Statistics for spatio-temporal data. John Wiley & Sons.\n\n[2] Cressie, N. (2015). Statistics for spatial data. John Wiley & Sons.", "claims_and_evidence": "The authors claimed that using transformers to estimate the intensity function under the spatial Poisson assumption outperformed baseline estimators.", "methods_and_evaluation_criteria": "The authors tested their method on a synthetic Poission point process as well as Colombia forest change data and UCDP georeferenced event dataset. The authors used relative error rate for the synthetic experiment and checked the consistency of their conlusions with prior work as well as human judgement for the non-synthetic sutdies.", "theoretical_claims": "This paper recalled known results on propensity scores, and proved a consistency result for the spacial estimator $N _ {\\omega} (Y_t)$. While I think recalling the known results is helpful for the structure, I failed to be convinced that there is significnat novelty on the theoretical front -- the proof in Appendix B2 is thorough, but Proposition 3 seems to be more like a routine proof. To be fair, the authors did not claim theoretical contributions, so this is fine by me.", "experimental_designs_or_analyses": "The authors are very clear on the baseline and evaluation metrics. However, I think the lack of validation for the non-synthetic studies makes the story less convincing. I understand that there is non \"true answer\" for the non-synthetic studies, but the authors could perhaps perform sensitivity analysis or convergence analysis of their estimators.\n\nI think it is helpful that the authors discuss the computation cost but a comparison of computation cost with the baseline might shed more light on the comparison between the proposed method.", "relation_to_broader_scientific_literature": "Satisfactory.", "essential_references_not_discussed": "Satisfactory.", "other_strengths_and_weaknesses": "I think, despite this paper is generally well-written, the constribution seems to be a bit weak -- there is no significnat contributions on the theory front, and this paper has a more application flavor. However, the empirical studies and analyses are related limited compared with what I would expect from an application paper.", "other_comments_or_suggestions": "In the figures/tables of empirical studies, perhaps adds error bars and std/confidence bands.", "questions_for_authors": "Did you try to relax the Poisson assumption and see how well this method holds up in the synthetic setup?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper addresses the challenge of estimating counterfactual outcomes from observational data in a spatial-temporal context, where treatment and outcome variables are both influenced by spatial and temporal factors. The proposed method operates under the assumption that there is no unmeasured confounding present in the observations. It utilizes propensity scores to adapt the well-known Inverse Probability Weighting (IPW) method to this spatial-temporal setting.

To implement this approach, the authors employ a deep learning architecture that combines Convolutional Neural Networks (CNNs) for computing propensity scores with transformers for encoding the intensity function. The effectiveness of this method is demonstrated through empirical comparisons with traditional counterfactual outcome estimation approaches, showing a reduction in error.

Overall, this paper expands on a classic IPW algorithm, applying it to high-dimensional, spatial-temporal data. It is particularly relevant for researchers and practitioners in reinforcement learning and causal inference. However, since this paper explores a novel domain, readers may find it challenging to visualize practical applications for causal effect estimation in spatial-temporal contexts. Therefore, it would be encouraged for the authors to include a running example that illustrates the proposed method throughout the paper.
Decision: Accept (poster)