Paper deep dive
Spatiotemporal Proximal Causal Inference under Hidden Confounding and Interference
Omar Faruque, Pavan Raj Ravi, Jianwu Wang
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Estimating causal effects from real-world spatiotemporal data is challenging due to hidden confounders and interference. Standard causal identification methods assume conditional exchangeability given observed covariates, which fails whenever hidden confounders affect both treatment and outcomes - a common setting in domains such as climate, environmental policy, epidemiology, and regional economics. In this paper, we propose a novel spatiotemporal proximal causal inference framework that extends proximal identification theory to spatiotemporal settings. The proposed method jointly captures local and neighborhood-level confounding information by introducing treatment- and outcome-inducing proxies, and we derive a spatiotemporal outcome confounding bridge function that identifies the potential outcome without requiring direct recovery of the hidden confounder. We establish the identifiability of this bridge function under proxy exclusion restrictions and a spatiotemporal completeness condition, and show that the resulting estimator recovers the outcome through a proximal generalization of the g-computation formula. To operationalize this identification result, we propose a neural architecture that learns proxies via transformer-based spatiotemporal encoders - coupled with a conditional mutual information critic to enforce exclusion restrictions and a moment-matching network to guarantee that the learned bridge function satisfies the underlying identifying equation. We further introduce a stabilized weighting scheme to address treatment support imbalance. Experiments on synthetic datasets demonstrate that our approach achieves comparable performance to baseline causal inference methods, while providing, to our knowledge, the first theoretically grounded outcomes for the hidden confounding in the presence of spatiotemporal interference through a proximal causal inference framework.
Tags
Links
- Source: https://arxiv.org/abs/2608.01352v1
- Canonical: https://arxiv.org/abs/2608.01352v1
Trouble viewing inline? Open PDF directly →
Full Text
83,770 characters extracted from source content.
Expand or collapse full text
Spatiotemporal Proximal Causal Inference under Hidden Confounding and Interference Omar Faruque, Pavan Raj Ravi, Jianwu Wang Abstract Estimating causal effects from real-world spatiotemporal data is challenging due to two compounding challenges: hidden confounders that evolve jointly over space and time, and interference, whereby a unit’s outcome is caused by its own and its neighbors’ treatment. Standard causal identification methods assume conditional exchangeability given observed covariates, which fails whenever unmeasured spatiotemporal confounders affect both treatment assignment and outcomes — a common setting in domains such as climate, environmental policy, epidemiology, and regional economics. In this paper, we propose a novel spatiotemporal proximal causal inference framework that extends proximal identification theory to spatiotemporal settings. The proposed method jointly captures local and neighborhood-level confounding information by introducing treatment- and outcome-inducing proxies, and we derive a spatiotemporal outcome confounding bridge function that identifies the potential outcome without requiring direct recovery of the hidden confounder. We establish the identifiability of this bridge function under a set of proxy exclusion restrictions and a spatiotemporal completeness condition, and show that the resulting estimator recovers the average treatment effect through a proximal generalization of the g-computation formula. To operationalize this identification result, we propose a neural architecture that learns proxies via transformer-based spatiotemporal encoders — coupled with a conditional mutual information critic to enforce exclusion restrictions and a moment-matching network to guarantee that the learned bridge function satisfies the underlying identifying equation rather than reducing to an unconstrained regression. We further introduce a stabilized weighting scheme to address treatment support imbalance under continuous, spatially correlated interventions. Experiments on synthetic datasets demonstrate that our approach achieves comparable performance to state-of-the-art spatiotemporal causal inference methods, while providing, to our knowledge, the first theoretically grounded outcomes for the hidden confounding in the presence of spatiotemporal interference through a proximal causal inference framework. 1 Introduction Estimating the causal effects of interventions from observational spatiotemporal data is critical across scientific and policy domains, ranging from evaluating environmental regulations on regional air quality to assessing the impact of public health interventions across neighboring communities. Unlike conventional predictive modeling, causal inference estimates how outcomes would change under hypothetical interventions, enabling reliable decision-making and policy evaluation. Most observational spatiotemporal systems evolve simultaneously over space and time, and treatments applied at one spatial location may influence outcomes not only at the treated unit but also at nearby units through interference, which violates the standard no-interference (SUTVA) assumption (Rubin 1980). Another fundamental challenge in these settings is unobserved (hidden) confounders, such as hidden environmental conditions, socioeconomic factors, or unmeasured system states, that simultaneously affect treatment assignment and outcomes. Standard identification strategies — matching, inverse propensity weighting, and the classical g-formula (Robins 1986)- rely on strong assumptions such as sequential ignorability or no hidden confounding. These assumptions are difficult to justify in spatiotemporal observational data, where unmeasured states frequently influence both treatment and outcomes simultaneously. Methods that explicitly address unmeasured spatial confounding — through distance-adjusted propensity matching (Papadogeorgou et al. 2019) or joint spatial-interference identification (Papadogeorgou and Samanta 2023) — have generally treated confounding and interference as separate problems. The deconfounder-style factor models address hidden confounding by estimating a substitute confounder from multiple causes (Wang and Blei 2019) and have been extended to longitudinal settings (Bica et al. 2020a). Khot et al. (2026) observed that spatial interference and hidden confounding are not independent nuisances: spatial interference provides exactly the multi-cause structure a deconfounder needs, and proposed the Spatial Deconfounder to recover a substitute spatial confounder via a conditional variational autoencoder. Oprescu et al. (2026) developed GST-UNet, a neural framework for spatiotemporal causal inference under time-varying confounding without any hidden confounder. However, the identifiability guarantees of factor-model-based deconfounders have been repeatedly questioned, showing that recovering a latent confounder from a factor model does not, in general, non-parametrically identify the interventional distribution P(Y(a))P(Y(a)) without additional untestable assumptions (Ogburn et al. 2019; D’Amour 2019). Proximal causal inference (PCI) provides a way to handle hidden confounding that does not require the confounder to be fully factorized or recovered. Following the negative-control literature (Lipsitch et al. 2010; Miao and Tchetgen Tchetgen 2017), Miao et al. (Miao et al. 2018) established identification of causal effects by exploiting a pair of proxies — a treatment-inducing proxy and an outcome-inducing proxy — each associated with the hidden confounder but satisfying exclusion restrictions and completeness assumptions. Tchetgen Tchetgen et al. (2024) extended this identification strategy to a full potential-outcomes framework and to longitudinal settings. Despite this progress, these methods neither account for graph-structured spatial dependence nor accommodate interference and spillover effects. In this paper, we propose a Spatiotemporal Proximal Causal Inference framework for estimating causal effects from observational spatiotemporal data in the presence of hidden confounding and interference. Our framework extends proximal causal inference theory for spatiotemporal systems that jointly accommodates interference and hidden confounding by spatiotemporal proxy exclusion restrictions and a completeness assumption. To operationalize the proposed theory, we develop an end-to-end neural architecture utilizing the spatiotemporal transformer that jointly learns the latent proxy variables and the outcome bridge function. We also apply proximal exclusion restrictions through conditional mutual information regularization and a bridge momentum objective to ensure identification of the potential outcome. Our contributions are summarized as follows: (i) We formalize the identification theory for spatiotemporal proximal causal inference and prove that a spatiotemporal outcome confounding bridge function exists under the treatment- and outcome-inducing proxy and exclusion restrictions. (i) We propose an end-to-end neural framework for learning proximal proxies and bridge functions, consisting of a transformer-based treatment-inducing proxy encoder with a denoising diffusion decoder, a spatial self-attention outcome-inducing proxy autoencoder, and a transformer-based bridge function network, trained jointly with a conditional-mutual-information critic enforcing the proxy exclusion restrictions. (i) We provide empirical validation of the proposed method. Through extensive experiments using synthetic benchmarks with known ground-truth causal structure, we demonstrate that the proposed approach yields comparable results to existing causal inference methods under hidden confounding and spatial interference. 2 Related Works We give a brief overview of the related literature here; see Section A for a detailed discussion. Our work sits at the intersection of three main branches: (i) proximal causal inference, (i) causal inference under interference, and (i) deconfounding methods with factor models. Proximal causal inference. This provides an identification framework by introducing treatment- and outcome-inducing proxy variables together with bridge functions that recover causal effects without directly observing the latent confounder. The foundational work of Miao et al. (2018) established nonparametric identification through proxy exclusion restriction under completeness assumptions. Tchetgen Tchetgen et al. (2024) extended this to a potential-outcomes framework and to longitudinal settings. Shi et al. (2020) develop multiply robust estimators for proximal causal inference with categorical unmeasured confounders; Cui et al. (2024) generalize this to a fully semiparametric framework with multiple robustness; and Shi et al. (2026) apply proximal identification to synthetic control methods for panel data. These methods are developed for i.i.d. or purely temporal data and assume no interference. Causal inference under interference. These studies relax the no-interference (SUTVA) assumption and construct a low-dimensional mapping of neighboring effects as interference, including partial-interference designs and spatial generalized-propensity-score (Hudgens and Halloran 2008; Aronow and Samii 2017; Tchetgen and VanderWeele 2012; Forastiere et al. 2021; Sobel 2006; Giffin et al. 2023). A separate line addresses unmeasured spatial confounding directly, through distance-adjusted propensity matching (Papadogeorgou et al. 2019), orthogonalization of spatial trends (Dupont et al. 2022) or joint Bayesian models of interference and latent spatial fields (Papadogeorgou and Samanta 2023). Recently, the GST-UNet (Oprescu et al. 2026) model combines a U-Net-based spatiotemporal encoder with iterative G-computation to handle interference with time-varying confounding, but requires that all confounders be observed rather than addressing hidden confounding. Deconfounding methods with factor models. A complementary line of research employs deep latent variable models to infer substitute confounders from multi-cause observed data (Wang and Blei 2019), extended to longitudinal data via recurrent factor models (Bica et al. 2020b). STCINet (Ali et al. 2024) uses a U-Net with double attention for spatial interference with an autoencoder-based factor model to reduce bias from hidden, time-varying confounding. Khot et al. (2025) observe that spatial interference supplies the multi-cause structure required for a deconfounder and propose the Spatial Deconfounder, estimating a substitute confounder via a conditional variational autoencoder with a spatial prior. Our work differs fundamentally from existing literature, as we extend proximal causal inference to the spatiotemporal setting with hidden confounding and interference, introducing neighborhood-based proxies and a spatiotemporal completeness condition under which the resulting bridge function identifies potential outcomes without directly reconstructing a substitute confounder through the factor model. 3 Problem Formulation We consider spatiotemporal random variables observed over a discrete gridded region Nx×NYN_x× N_Y for T time steps. Each spatial location is indexed by i∈[1,2,3,…,NX∗NY]i∈[1,2,3,...,N_X*N_Y]. For each spatial location i, let (i)∈[1,2,3,…,NX∗NY]N(i)∈[1,2,3,...,N_X*N_Y] denote its spatial neighborhood defined by adjacent grid cells or a distance threshold. At each time step t∈[1,…,T]t∈[1,...,T] we observe continuous treatment Ai,tA_i,t, outcome variable Yi,tY_i,t, covariates Xi,tX_i,t. Besides these variables, we assume there are unobserved confounders Ui,tU_i,t, affecting both the treatment and outcome variables, and unobserved confounders evolve over time: Ut=g(Ut−1,ξt)U_t=g(U_t-1, _t). Here each spatiotemporal variable constructs a 3D tensor of dimension T×NX×NYT× N_X× N_Y. The uppercase symbol with only a time index XtX_t represents the whole spatial region at time t, and Xt:t+τX_t:t+τ and X¯t X_t both denote its value over a time interval of length τ starting at t. In the spatiotemporal interference setting, every variable defined at a location i must also carry information about its neighborhood (i)N(i), so in this paper we let the location subscript i implicitly denote the pair (i,(i))(i,N(i)), unless stated otherwise. Specifically, Xi,tX_i,t should be interpreted as the (Xi,t,X(i),t)(X_i,t,X_N(i),t). For notational convenience, we therefore omit the explicit dependence on (i)N(i) throughout the remainder of the paper. Figure 1: DAG representation of treatment, outcome, and confounders; (a) DAG with hidden confounder U and (b) proximal causal structure with treatment-inducing proxy Z and outcome-inducing proxy W. Building on Rubin’s potential outcomes framework, our main goal is to learn potential outcomes for location-specific treatment assignment considering previous historical observations. More specifically, for applied intervention ai,ta_i,t on a predetermined spatial region i⊆(NX,NY)i (N_X,N_Y) at time step t, the goal is to estimate the potential outcomes Yt+τ[ai,t]Y_t+τ[a_i,t] at a future time step t+τt+τ for τ≥1τ≥ 1. E[Yt+τ[ai,t]|At,Xt,Ut]E[Y_t+τ[a_i,t]|A_t,X_t,U_t] (1) 4 Spatiotemporal Causal Identification We follow the potential outcome framework [Rubin, 1974], and for observational systems, identification of potential outcomes relies on standard causal inference assumptions. Besides these, we have to introduce some necessary assumptions for spatiotemporal and time-varying requirements. Assumption 1: Consistency and positivity. We assume that if the treatment AtA_t satisfies At=atA_t=a_t, then the potential outcome for the treatment is the same as the observed outcome Yt=Yt[at]Y_t=Y_t[a_t]. The probability of treatment assignment P(Ai,t=ai,t|Xt)>0P(A_i,t=a_i,t|X_t)>0 for all treatment values, and the treatment assignment is not confounded by future outcomes. Assumption 1 means that at each time step the probability of a treatment will be assigned is higher than zero. Existing methods also assume the standard conditional exchangeability assumption, where it is assumed that the treated and untreated units share the same distribution based on the sufficiently measured covariate set. In a spatiotemporal context, let Ai,tA_i,t be the treatment at location i and time t, and Yi′,t′(a)Y_i ,t (a) be the potential outcome at location i′i and time t′t (t≤t′)(t≤ t ) under treatment a. Conditional exchangeability is expressed as: Yi′,t′(a)⟂Ai,t|X1:t,A1:t−1,Y1:t−1Y_i ,t (a) A_i,t|X_1:t,A_1:t-1,Y_1:t-1 This assumption implies that by conditioning on observed covariates, we block all backdoor paths that cause confounding. This condition holds for systems without any hidden confounder, and this is not practical for real-world systems. Here we take advantage of the proximal causal inference method proposed by Tchetgen et al.(Tchetgen Tchetgen et al. 2024) to approximate the hidden confounder and remove the confounding bias by adjusting the conditioning set. In the spatiotemporal process, the hidden confounder is no longer a scalar nuisance, it typically evolves over time and space. Assumption 2: Markov dynamics of hidden confounding. We assume the hidden confounder is autoregressive Ui,t=f(Ui,t−1,ξt)U_i,t=f(U_i,t-1, _t), and the noise term ξ is independent of space and time. Assumption 3: Proximal causal structure. We assume that there is a treatment-inducing proxy Z and an outcome-inducing proxy W such that the conditional independence given below holds: Y⊥Z|A,U,XandW⊥(A,Z)|U,XY Z|A,U,X\;\;\;and\;\;\;W (A,Z)|U,X Following proximal causal inference, we assume the treatment-inducing proxy Z is the part of the hidden confounder U that does not directly affect the outcome (Figure 1(b)). According to Tchetgen et al. (Tchetgen Tchetgen et al. 2024), Z may or may not be a direct cause of A. In our case, we consider that Z is related to treatment A through its association with the hidden confounder U. Similarly, the outcome-inducing proxy W is part of the hidden confounder that is related to outcome Y, but not directly related to treatment A and treatment-inducing proxy Z. Z=fz(U,X)andZ⊥A|X,UZ=f_z(U,X)\;\;\;and\;\;\;Z\; \;A\;|\;X,U (2) W=fw(U,X)andY⊥W|X,U,AW=f_w(U,X)\;\;\;and\;\;\;Y\; \;W\;|\;X,U,A (3) Equation 2 states that Z is associated with U, and it does not itself drive treatment assignment beyond what is mediated by U. So, when the covariates X and hidden confounder U are fixed, Z has no causal effect on treatment. This is the spatial analogue of a classic negative control exposure (Miao et al. 2024; Miao and Tchetgen Tchetgen 2017). Likewise, W is related to U, but once X, A, and U are fixed, this carries no further information about the potential outcome. Equations 2 and 3 are also known as proxy exclusion restrictions (Cui et al. 2024; Ringlein et al. 2025). Assumption 4: Spatiotemporal completeness. For any square-integrable function v(⋅)v(·) E[v(Ui,t)|Ai,t,Zi,t,Xi,t,H¯t−1]=0⇒v(Ui,t)=0E[v(U_i,t)|A_i,t,Z_i,t,X_i,t, H_t-1]=0 v(U_i,t)=0 This is an extension of the completeness condition defined in (Tchetgen Tchetgen et al. 2024; Miao et al. 2024) in the spatiotemporal context; it requires that Z be rich enough to distinguish variation in the hidden-confounder vector U. This assumption guarantees that proximal causal learning can potentially approximate proxies Z and W to account for unmeasured confounding as long as the probability density field of U is no larger than that of either proxy (Shi et al. 2020; Miao et al. 2024). Intuitively, completeness requires that the proxy retains sufficient variation to distinguish every possible realization of the latent confounding process. Bridge Function In spatiotemporal settings, each unit’s potential outcome is simultaneously affected by local treatment, neighboring treatments through spatial interference and spillover, observed covariates, and hidden confounders. Following the proximal causal structure, we learn treatment-inducing proxy Z and outcome-inducing proxy W (assumption 3). As the spatial neighborhood resembles the multi-cause criteria mandatory for the deconfounder (Bica et al. 2020b; Hatt and Feuerriegel 2024a; Khot et al. 2026), we can use any factor model to learn the Z and W proxies, where Z and W carry both the local and neighborhood components. The target of proximal causal inference is to provide an outcome bridge function to enable identification of potential outcomes without explicitly recovering the unobserved process U. Definition 1: Spatiotemporal Bridge Function. A measurable function h(⋅)h(·) is called a spatiotemporal outcome confounding bridge function if it satisfies E[Yt|At,Zt,Xt,H¯t−1]=E[ht(At,Wt,Xt,H¯t−1)|At,Zt,Xt,H¯t−1]. splitE[Y_t|A_t,Z_t,X_t, H_t-1]=E[h_t(A_t,W_t,X_t, H_t-1)\\ |A_t,Z_t,X_t, H_t-1]. split (4) Intuitively, the bridge function does not estimate the hidden confounder itself; rather, it reproduces its effect on the conditional outcome distribution through the proxy W and Z. Unlike standard regression functions, the bridge function is therefore an identification object rather than a prediction model. Two important features of this h(⋅)h(·) function are specific to the spatiotemporal settings. First, the outcome-inducing proxy W must resemble the own-unit and neighborhood confounder U. Second, the arguments of h(⋅)h(·) include the neighborhood exposure to simultaneously encode the dependence of Y on both direct treatment and spillover exposure. An important condition is that the proxies Z and W share no direct association beyond what is mediated by the true confounding state. By utilizing the above bridge function, we conclude our identification results as follows: Theorem 1: Potential Outcome. Under the assumption of proximal causal structure, there exists a measurable function ht(At,Wt,Xt,H¯t−1)h_t(A_t,W_t,X_t, H_t-1) such that E[Yt|At,Zt,Xt,H¯t−1]=E[ht(At,Wt,Xt,H¯t−1)|At,Zt,Xt,H¯t−1] splitE[Y_t|A_t,Z_t,X_t, H_t-1]=E[h_t(A_t,W_t,X_t, H_t-1)\\ |A_t,Z_t,X_t, H_t-1] split Based on this function, the potential outcome for treatment a is defined as: ψ(a)=E[ht(At=at,Wt,Xt,H¯t−1)|At=at,Zt,Xt,H¯t−1]ψ(a)=E[h_t(A_t=a_t,W_t,X_t, H_t-1)|A_t=a_t,Z_t,X_t, H_t-1] Proof. The proof is provided in Appendix B.∎ Based on Theorem 1, for two different treatment values a and a′a^ at t, the average treatment effect (ATE) is: ATE(a,a′)=E[Yi,t(a)−Yi,t(a′)].ATE(a,a^ )=E[Y_i,t(a)-Y_i,t(a^ )]. (5) The expectation of the potential outcome can be identified by the bridge function through Theorem 1. So each expectation can be expressed as E[Yi,t(a)]=E[ht(at,Wt,Xt,H¯t−1)]E[Y_i,t(a)]=E[h_t(a_t,W_t,X_t, H_t-1)] E[Yi,t(a′)]=E[ht(at′,Wt,Xt,H¯t−1)]E[Y_i,t(a^ )]=E[h_t(a^ _t,W_t,X_t, H_t-1)] Now the ATE is computed by subtracting the second expectation term from the first. Since the distribution of covariates or proxies does not depend on the treatment argument, using linearity of the expectation, we get the following form of equation 5 for ATE computation. ATE(a,a′)=E[ht(at,Wt,Xt,H¯t−1)−ht(at′,Wt,Xt,H¯t−1)]ATE(a,a^ )=E[h_t(a_t,W_t,X_t, H_t-1)-h_t(a^ _t,W_t,X_t, H_t-1)] Connection to G-Computation The classical g-computation formula can be extended for a spatiotemporal setting from the above derivation. In the standard g-formula, the counterfactual mean is obtained by integrating the conditional outcome model over the covariate distribution, E[Y(a)]=∫E[Y|A=a,X=x]⋅P(x)E[Y(a)]= [Y|A=a,X=x]· dP(x) for treatment a and covariates X. In proximal causal inference, the outcome model is replaced with the identified outcome bridge function. Consequently, the proximal g-computation formula for the spatiotemporal settings with hidden confounding is E[Y(a)]=∫h(a,w,x)Fw,x(w,x).E[Y(a)]= h(a,w,x)dF_w,x(w,x). This expression remains identifiable even when the confounder is completely unobserved, provided that the proxy assumptions and completeness condition hold. 5 Implementation of the Proposed Framework The theoretical analysis above establishes the identifiability of the potential outcome in the form of the bridge function for spatiotemporal settings with hidden confounding and interference. In this section, we provide a neural implementation to learn the potential outcome following the procedure. The overall architecture of the proposed framework is illustrated in Figure 2. Figure 2: Architecture of the proposed proximal causal inference method. Model Architecture The proposed architecture instantiates the proximal identification results with three learnable modules: (i) Treatment-inducing Proxy Learning Module (fZ)(f_Z): a Spatiotemporal Transformer encoder followed by a conditional denoising diffusion decoder, which learns a representation of the treatment-inducing proxy Z, while satisfying the treatment proxy exclusion restriction. (i) Outcome-inducing Proxy Learning Module (fW)(f_W): a spatial self-attention autoencoder that learns a representation of the outcome-inducing proxy W, while remaining independent of treatment. (i) Outcome Bridge Module (hθ)(h_θ): a Transformer encoder that regresses Yi,tY_i,t from the learned proxy representations, observed covariates, and treatment, as established in Theorem 1. Treatment-inducing Proxy Learning Module. The treatment-inducing proxy Z should capture the part of the hidden confounder related to the treatment and independent of the potential outcome given the treatment and covariates. Here we adopt a transformer-based encoder and denoising diffusion module to learn the latent representation of Z Z. The spatiotemporal transformer encoder models both the spatial information, temporal dependencies, and non-linear interactions among covariates and treatment history, ei,tZ=Transformer(X¯,A¯)e^Z_i,t=Transformer( X, A). Self-attention is computed for each spatial location over the neighborhood for a time window [t−τ,t][t-τ,t]. Instead of reconstructing Z with a deterministic decoder, here we utilize a denoising diffusion module dθ(Zi,t|ei,tZ)d_θ(Z_i,t|e^Z_i,t) to capture its distribution better given the latent features. During the forward diffusion process, Gaussian noise is gradually injected into the latent encoder features z(k)=α¯kz(0)+1−α¯kϵ,ϵ∼(0,I),k=1,…,k,z^(k)= α_kz^(0)+ 1- α_kε,ε (0,I),k=1,...,k, with the standard noise schedule αk\ _k\. The backward diffusion process learns to recover the latent proxy by estimating the injected noise. The denoising process is trained with the standard simplified denoising score-matching objective: ℒdiff=Ek∼1,K,ϵ∼(0,I)[‖ϵ−ϵθ(z,k,e)‖22].L_diff=E_k \1,K\,ε (0,I) [ \|ε- _θ(z,k,e) \|^2_2 ]. Diffusion learning regularizes the proxy manifold and encourages the latent representation to approximate the distribution of the hidden confounder rather than overfitting the observed treatment. Outcome-inducing Proxy Learning Module. The outcome-inducing proxy W summarizes the hidden confounding information related to the outcome while remaining independent of treatment after conditioning on covariates and the outcome. Since the outcome exhibits strong spatial dependence due to interference, we construct the outcome proxy using a spatial self-attention autoencoder. The spatial self-attention encoder takes the history of observed covariates and the outcome (X¯,Y¯)( X, Y) for a time window [t−τ,t][t-τ,t]. The neighborhood information for each unit is embedded with the convolution kernel to aggregate temporal information, then multi-head spatial self-attention is applied to aggregate neighboring information, ei,tW=SpatialAttnEncW(Embed(X^,Y^))∈RdWe^W_i,t=SpatialAttnEnc_W(Embed( X, Y))∈ R^d_W. The attended representation ei,tWe^W_i,t is passed through an encoder-decoder architecture, and the bottleneck representation constitutes the learned outcome-inducing proxy W W. W^=EncW(ei,tW),(X^,Y^)=DecW(EncW(ei,tW)) W=Enc^W(e^W_i,t),( X, Y)=Dec^W(Enc^W(e^W_i,t)) Conditional Mutual Information Critic. The proximal causal structure assumption (Assumption 3) and proxy exclusion restrictions (equations 2 and 3) require the learned proxies to satisfy the conditional independence of treatment W⟂A∣X,UW A X,U and outcome Y⟂Z∣A,X,UY Z A,X,U. We estimate conditional mutual information (CMI) using a neural mutual-information critic Tϕ(⋅)T_φ(·), and approximate the operational implication at the representation level by penalizing the CMI. The variational critic for treatment-inducing proxy TϕZ(Y|Z^,X,A)=I^(Y;Z^|A,X)T^Z_φ(Y| Z,X,A)= I(Y; Z|A,X) and for the outcome-inducing proxy TϕW(A|W^,X)=I^(A;W^|X)T^W_φ(A| W,X)= I(A; W|X). The goal is to enforce these mutual information close to zero, I(Y;Z^|A,X)≈0I(Y; Z|A,X)≈ 0 and I(A;W^|X)≈0I(A; W|X)≈ 0, to satisfy the proximal exclusion restrictions. We adopt the Contrastive Log-ratio Upper Bound (CLUB) method for minimizing CMI (Cheng et al. 2020), because it provides an upper bound ICLUB(Y;Z^|A,X)≥I(Y;Z^|A,X)I_CLUB(Y; Z|A,X)≥ I(Y; Z|A,X), so that driving the estimate to zero provably drives the true CMI toward zero. For outcome-inducing proxy W W, let TϕW(a|w^,x)T^W_φ(a| w,x) be a variational approximation of the actual conditional density p(a|w^,x)p(a| w,x). The conditional CLUB bound is ICLUBW(A;W^|X)=Ep(a,w^,x)[logTϕ(a|w^,x)]−Ep(w^,x)p(a|x)[logTϕ(a|w^,x)], splitI^W_CLUB(A; W|X)=E_p(a, w,x) [log\,T_φ(a| w,x) ]-\\ E_p( w,x)p(a|x) [log\,T_φ(a| w,x) ] split, where the first term is estimated directly from matched samples (X,W^,X)(X, W,X) and the second term is approximated by ramdomly shuffling treatment A′A^ while keeping (W^,X)( W,X) fixed. The variational critic Tϕ(⋅)T_φ(·) is trained to maximize the log-likelihood of matched pairs, argmax1/N∑logTϕ(a|w^,xarg\,max\1/NΣ log\,T_φ(a| w,x\. So the CMI loss from the treatment-inducing proxy module and outcome-inducing proxy module are: ℒCMIW=−ICLUBW(A;W^|X)andℒCMIZ=−ICLUBW(Y;Z^|X,A).L^W_CMI=-I^W_CLUB(A; W|X)\,and\,L^Z_CMI=-I^W_CLUB(Y; Z|X,A). Outcome Bridge Module. Given the learned proxies (W^,Z^)( W, Z), the outcome bridge module estimates the confounding bridge function hθ(W^i,t,Ai,t,Xi,t,H¯t)h_θ( W_i,t,A_i,t,X_i,t, H_t). The bridge module consists of stacked transformer encoder layers that jointly model temporal dependency, spatial interactions, treatment history, and proxy representations. The transformer encoder computes TransformerEncoder(W^,A,X,H¯)TransformerEncoder( W,A,X, H) followed by a feed forward block to estimate the bridge function hθh_θ. hθ(W^i,t,Ai,t,Xi,t,H¯t)=FeedForward(TransformerEnc(W^i,t,Ai,t,Xi,t,H¯t)) splith_θ( W_i,t,A_i,t,X_i,t, H_t)=FeedForward(\\ TransformerEnc( W_i,t,A_i,t,X_i,t, H_t)) split (6) Self-attention is applied over the own/neighborhood token set, allowing hθh_θ to flexibly weight own-unit versus neighborhood-proxy contributions when predicting the outcome. For continuous intervention, treatment assignments are generally imbalanced across space and time. Certain treatment configurations occur frequently, whereas others are observed only sparsely, and the conditional density P(a∣X)P(a X) can be arbitrarily small or effectively zero. Consequently, directly minimizing the bridge prediction loss tends to bias the learned bridge function hθh_θ toward high-density treatment regions. We incorporate a stabilized weighting module into the optimization to overcome this issue. The idea is to compensate for the imbalance of treatment exposure by assigning larger weights for low-density treatment regions while preventing the bridge from being dominated by frequently observed treatments. This module estimates the density of the treatment given covariates p^(A|X) p(A|X), which is implemented using a transformer encoder and trained by maximum likelihood of the observed (A,X)(A,X) tuples: −logp^(A|X)-log p(A|X). The stabilized weights for each observation are computed as: ω^i,t=p^(Ai,t)/p^(Ai,t|Xi,t) ω_i,t= p(A_i,t)/ p(A_i,t|X_i,t). Here p^(Ai,t) p(A_i,t) is the marginal density of the treatment. Now, instead of using Y−hθY-h_θ, we use the weighted prediction objective ω^(Y−hθ) ω(Y-h_θ) to train the outcome bridge function. The bridge function of equation 6 is trained to match the observed outcome Y directly on (W,A,X)(W,A,X), rather than conditioning on treatment-inducing proxy ZZ. Without conditioning on ZZ, hθh_θ is just a simple regression E[Y∣W,A,X]E[Y W,A,X], which is not guaranteed to solve the proximal causal inference for hidden confounding effects. To enforce the conditional restriction on Z for hθh_θ, we include the bridge moment network (Mψ)(M_ψ) containing a conditional expectation of the bridge function residual, E[(Y−hθ(W,A,X))∣Z,A,X]E [(Y-h_θ(W,A,X)) Z,A,X ]. The bridge moment network contains transformer encoder layers following a feed-forward block to capture both spatial and temporal information. Assume the bridge function residual is R=Y−hθ(W,A,X)R=Y-h_θ(W,A,X). The bridge momentum network takes (Z,A,X)(Z,A,X) as input and estimates the residual, R^=E(R∣Z,A,X) R=E(R Z,A,X). Then, the treatment-inducing proxy Z is enforced into the outcome bridge estimation through the following bridge momentum loss, and the bridge momentum network tries to maximize this loss: ℒmoment=−ω^‖(Y−hθ(W,A,X))−Mψ(Z,A,X)‖2.L_moment=- ω||(Y-h_θ(W,A,X))-M_ψ(Z,A,X)||^2. This empirical momentum network pins hθh_θ to be the bridge function rather than an arbitrary regression on (W,A,X)(W,A,X) and explicitly enforces the identifying condition. Optimization The overall optimization objective of the proposed model comprises components from all three modules. To train the treatment-inducing proxy module, we enforce reconstruction loss of the input features (X,A)(X,A) with Mean Squared Error (MSE). Besides this, the denoising diffusion objective and KL divergence regularization is applied to the latent representation to capture the distribution of Z Z. Finally, the CMI constraint is integrated to make the learned proxy Z Z independent of the outcome to satisfy the proxy exclusion restriction. Hence, the overall treatment-inducing module loss becomes ℒfZ=||(A,X)−(A,X)^||2+ℒdiff+DKL(q(Z^)||P(Z^))+ℒCMIZL_f_Z=||(A,X)- (A,X)||^2+L_diff+D_KL(q( Z)||P( Z))+L^Z_CMI The outcome-inducing proxy module reconstructs the input features (X,Y)(X,Y) using the spatial autoencoder and learns the latent representation of W W. This module is optimized using the MSE loss of the reconstructed features and the CMI constraint of A and W W. ℒfW=‖(A,Y)−(A,Y)^‖2+ℒCMIWL_f_W=||(A,Y)- (A,Y)||^2+L^W_CMI The outcome bridge module is optimized by minimizing the empirical bridge equation with the observed outcome and the bridge momentum loss derived earlier. ℒhθ=‖Y−hθ(W,A,X)‖2+ℒmomentL_h_θ=||Y-h_θ(W,A,X)||^2+L_moment 6 Experiments We evaluate the proposed framework using synthetic spatiotemporal datasets with known ground truth that incorporate key challenges: hidden confounder, interference, spillover, confounding, and temporal carryover. Using the synthetic data generation process, we compare our proposed framework against a set of representative baselines. As ground truth potential outcomes for treatment intervention are never observable in real-world data, we follow standard practice to evaluate on synthetic data with a known process. Additional details on data generation, model hyperparameters, and validation on real-world data can be found in supplementary section D. Replication code is available at https://anonymous.4open.science/r/Spatiotemporal-Proximal-Causal-Inference-2BF6. Synthetic Data We construct a synthetic spatiotemporal environment on a 28×2828× 28 (NX×NY)(N_X× N_Y) grid and generate two synthetic datasets over T=5000T=5000 time steps to jointly exercise hidden confounding and interference. We implement the following data generation process: Ut=Ut−1+dt∇2Ut−1+dtαU¯,t−1+ϵt,U_t=U_t-1+dt∇^2U_t-1+dtα U_N,t-1+ _t, Xt=Xt−1+dt∇2Xt−1+dtβUt−1++ϵt,X_t=X_t-1+dt∇^2X_t-1+dtβ U_t-1++ _t, At=At−1+dt∇2At−1+dtγ1Ut−1+dtγ2Xt−1+ϵt,A_t=A_t-1+dt∇^2A_t-1+dt _1U_t-1+dt\; _2X_t-1+ _t, Yt=Yt−1+dt∇2Yt−1+dtω1Ut−1+dtω2Xt−1+dtω3At−1+dtω4A¯,t−1+ϵt, splitY_t=Y_t-1+dt∇^2Y_t-1+dt _1U_t-1+dt _2X_t-1+dt _3A_t-1\\ +dt _4 A_N,t-1+ _t, split where dtdt is the rate of change, ∇2∇^2 is the Laplacian operator for diffusion process, ϵε is independent noise, and A¯,t−1 A_N,t-1 is the lagged mean neighborhood (8 adjacent grid locations) treatment as interference. Here, U acts as the hidden confounder influencing all other variables and is never included in the dataset, X is the observed confounder, A is the treatment, and Y is the outcome variable. To generate counterfactual data, first, we generate treatment intervention (Acf)(A_cf) by nudging a specific region of the treatment variable A[nX1:nX2,nY1:nY2]A[n_X1:n_X2,n_Y1:n_Y2] and then following the diffusion and interference process. Based on this intervened treatment AcfA_cf, we generate counterfactual outcome YcfY_cf. Synthetic dataset 1 contains one hidden variable and one confounder, whereas synthetic dataset 2 contains two hidden variables and two confounders. Comparative Analysis We compare the proposed framework against five representative deep learning and spatiotemporal baselines. A convolutional-LSTM-based U-Net architecture and a spatiotemporal Transformer encoder-decoder architecture serve as fully naive baselines capturing spatial and temporal information, ignoring causal identification. STCINET(Ali et al. 2024), an interference-aware spatiotemporal model for direct and indirect effect estimation, and GST-UNet (Oprescu et al. 2026), which addresses time-varying confounding for spatiotemporal causal inference. However, these methods do not address hidden confounders. The Spatial Deconfounder (Khot et al. 2026) method generates a latent substitute for hidden confounders in the spatial setting using a C-VAE-based deconfounder, the closest comparator since it targets the hidden confounding and interference. All baselines and our method were trained on the factual data and evaluated on the counterfactual data. All models are implemented using PyTorch and trained using the Adam optimizer. Following standard practice in the literature we report the standardized RMSE of the estimated counterfactual outcome Y^cf Y_cf for the whole grid and the intervened region here, and additional analyses are reported in Section D. Table 1: Performance comparison of baseline methods on synthetic datasets. Method Synthetic-1 Synthetic-2 RMSEWRMSE_W RMSEIRMSE_I RMSEWRMSE_W RMSEIRMSE_I UNet 0.034 0.065 0.038 0.096 Tansformer 0.033 0.065 0.061 0.160 STCINET 1.008 1.466 0.359 0.910 GST-UNet 0.034 0.061 0.072 0.164 Spatial Decon. 0.341 0.304 0.099 0.027 Proposed 0.021 0.015 0.069 0.066 Table 1 reports the counterfactual prediction accuracy for the whole grid RMSEWRMSE_W and the region where we applied intervention RMSEIRMSE_I. The proposed method outperforms all baselines on dataset 1 and achieves better performance then basline causal models for dataset 2. The UNet and Transformer models applied attention-based spatiotemporal designs and provided comparable performance. GSTU-Net incorporates richer spatiotemporal representations with the attention gating mechanism and outperforms STCINet and spatial deconfounder models, but cannot eliminate bias introduced by unobserved confounding. Spatial Deconfounder achieves competitive performance on the 2nd dataset by learning substitute latent confounders through a conditional variational autoencoder. By contrast, the proposed method combines proximal identification theory with CMI restriction to estimate treatment- and outcome-inducing proxies, and the bridge network explicitly satisfies the proximal bridge moment equation. This theoretical grounding enables more accurate estimation of potential outcomes and treatment effects. Ablation Study To quantify the contribution of each proposed component, we perform a comprehensive ablation study. We consider the following variants of the proposed framework: (i) without diffusion in proxy learning, (i) without stabilized weight module, and (i) without bridge moment optimization. The details of these ablation versions are provided in Section E. The complete end-to-end framework consistently outperforms every ablation. From Table 2, we can see that removing the stabilized weight module substantially degrades the performance of outcome prediction due to the treatment imbalance. Then removing the bridge momentum objective also causes performance degradation, demonstrating that directly optimizing the outcome prediction loss increases bias. Table 2: Ablation analysis between proposed framework and its different variants. Method Synthetic-1 Synthetic-2 RMSEWRMSE_W RMSEIRMSE_I RMSEWRMSE_W RMSEIRMSE_I W/O Diffusion 0.022 0.016 0.098 0.091 W/O Stabilization 0.066 0.062 0.255 0.255 W/O Momentum 0.040 0.039 0.124 0.123 Proposed 0.021 0.015 0.069 0.066 7 Conclusion We proposed a Spatiotemporal Proximal Causal Inference framework for estimating causal effects from spatiotemporal data in the presence of hidden confounding and interference. We extended proximal causal inference from temporal settings to spatiotemporal systems by introducing treatment- and outcome-inducing proxies with explicit neighborhood mapping. Based on proxy exclusion and completeness assumptions, we proved that the resulting outcome confounding bridge function is identifiable. We operationalized the theoretical framework with a neural architecture combining transformer/diffusion proxy learners, CMI critics for proxy exclusion restrictions, and an adversarial moment-matching network enforcing the bridge equation, together with stabilized weighting to handle imbalanced treatment. For synthetic benchmarks with known ground truth, our proposed framework generates comparable results. Some limitations to mention are: the completeness condition becomes more erratic as neighborhood size grows, and counterfactual values far outside the observed treatment distribution may remain unreliable. We believe this work provides a framework for integrating proximal causal inference with deep spatiotemporal representation learning and offers a promising direction for future work. References S. Ali, O. Faruque, and J. Wang (2024) Estimating direct and indirect causal effects of spatiotemporal interventions in presence of spatial interference. In European Conference on Machine Learning and Knowledge Discovery in Databases (ECML-PKDD), Vol. 14943, p. 213–230. Cited by: Appendix A, Table 3, Appendix D, §2, §6. P. M. Aronow and C. Samii (2017) Estimating average causal effects under general interference, with application to a social network experiment. The Annals of Applied Statistics. Cited by: Appendix A, §2. A. Bennett, N. Kallus, and T. Schnabel (2019) Deep generalized method of moments for instrumental variable analysis. Advances in neural information processing systems 32. Cited by: Appendix A. I. Bica, A. Alaa, and M. Van Der Schaar (2020a) Time series deconfounder: estimating treatment effects over time in the presence of hidden confounders. In International conference on machine learning, p. 884–895. Cited by: Appendix A, Table 3, §1. I. Bica, A. Alaa, and M. Van Der Schaar (2020b) Time series deconfounder: estimating treatment effects over time in the presence of hidden confounders. In International conference on machine learning, p. 884–895. Cited by: §2, §4. P. Cheng, W. Hao, S. Dai, J. Liu, Z. Gan, and L. Carin (2020) Club: a contrastive log-ratio upper bound of mutual information. In International conference on machine learning, p. 1779–1788. Cited by: §5. Y. Cui, H. Pu, X. Shi, W. Miao, and E. Tchetgen Tchetgen (2024) Semiparametric proximal causal inference. Journal of the American Statistical Association 119 (546), p. 1348–1359. Cited by: Appendix A, Table 3, §2, §4. A. D’Amour (2019) On multi-cause causal inference with unobserved confounding: counterexamples, impossibility, and alternatives. arXiv preprint arXiv:1902.10286. Cited by: Appendix A, §1. N. Dikkala, G. Lewis, L. Mackey, and V. Syrgkanis (2020) Minimax estimation of conditional moment models. Advances in Neural Information Processing Systems 33, p. 12248–12262. Cited by: Appendix A. E. Dupont, S. N. Wood, and N. H. Augustin (2022) Spatial+: a novel approach to spatial confounding. Biometrics 78 (4), p. 1279–1290. Cited by: Appendix A, Table 3, §2. L. Forastiere, E. M. Airoldi, and F. Mealli (2021) Identification and estimation of treatment and interference effects in observational studies on networks. Journal of the American Statistical Association 116 (534), p. 901–918. Cited by: Appendix A, §2. A. Giffin, B. Reich, S. Yang, and A. Rappold (2023) Generalized propensity score approach to causal inference with spatial interference. Biometrics 79 (3), p. 2220–2231. Cited by: Appendix A, Table 3, §2. T. Hatt and S. Feuerriegel (2024a) Sequential deconfounding for causal inference with unobserved confounders. In Causal Learning and Reasoning, p. 934–956. Cited by: §4. T. Hatt and S. Feuerriegel (2024b) Sequential deconfounding for causal inference with unobserved confounders. In Causal Learning and Reasoning, p. 934–956. Cited by: Appendix A. M. G. Hudgens and M. E. Halloran (2008) Toward causal inference with interference. Journal of the american statistical association 103 (482), p. 832–842. Cited by: Appendix A, §2. A. Khot, M. Oprescu, M. Schröder, A. Kagawa, and X. Luo (2025) Spatial deconfounder: interference-aware deconfounding for spatial causal inference. arXiv preprint arXiv:2510.08762. Cited by: Appendix A, Table 3, §2. A. Khot, M. Oprescu, M. Schröder, A. Kagawa, and X. Luo (2026) Spatial deconfounder: interference-aware deconfounding for spatial causal inference. In International Conference on Machine Learning, Cited by: §1, §4, §6. M. Lipsitch, E. T. Tchetgen, and T. Cohen (2010) Negative controls: a tool for detecting confounding and bias in observational studies. Epidemiology 21 (3), p. 383–388. Cited by: Appendix A, Table 3, §1. L. Liu and M. G. Hudgens (2014) Large sample randomization inference of causal effects in the presence of interference. Journal of the american statistical association 109 (505), p. 288–301. Cited by: Appendix A. A. Mastouri, Y. Zhu, L. Gultchin, A. Korba, R. Silva, M. Kusner, A. Gretton, and K. Muandet (2021) Proximal causal learning with kernels: two-stage estimation and moment restriction. In International conference on machine learning, p. 7512–7523. Cited by: Appendix A. W. Miao, Z. Geng, and E. J. Tchetgen Tchetgen (2018) Identifying causal effects with proxy variables of an unmeasured confounder. Biometrika 105 (4), p. 987–993. Cited by: Appendix A, Table 3, §1, §2. W. Miao, X. Shi, Y. Li, and E. J. Tchetgen Tchetgen (2024) A confounding bridge approach for double negative control inference on causal effects. Statistical Theory and Related Fields 8 (4), p. 262–273. Cited by: §4, §4. W. Miao and E. Tchetgen Tchetgen (2017) Invited commentary: bias attenuation and identification of causal effects with multiple negative controls. American journal of epidemiology 185 (10), p. 950–953. Cited by: Appendix A, Table 3, §1, §4. E. L. Ogburn, I. Shpitser, and E. J. T. Tchetgen (2019) Comment on “blessings of multiple causes”. Journal of the American Statistical Association 114 (528), p. 1611–1615. Cited by: Appendix A, §1. M. Oprescu, D. Park, X. Luo, S. Yoo, and N. Kallus (2026) GST-unet: a neural framework for spatiotemporal causal inference with time-varying confounding. Advances in Neural Information Processing Systems 38, p. 17295–17322. Cited by: Appendix A, Table 3, Appendix D, §1, §2, §6. G. Papadogeorgou, C. Choirat, and C. M. Zigler (2019) Adjusting for unmeasured spatial confounding with distance adjusted propensity score matching. Biostatistics 20 (2), p. 256–272. Cited by: Appendix A, Table 3, §1, §2. G. Papadogeorgou and S. Samanta (2023) Spatial causal inference in the presence of unmeasured confounding and interference. arXiv preprint arXiv:2303.08218. Cited by: Appendix A, Table 3, §1, §2. G. V. Ringlein, T. Q. Nguyen, P. P. Zandi, E. A. Stuart, and H. Parikh (2025) Demystifying proximal causal inference. arXiv preprint arXiv:2512.24413. Cited by: §4. J. Robins (1986) A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect. Mathematical modelling 7 (9-12), p. 1393–1512. Cited by: §1. D. B. Rubin (1980) Randomization analysis of experimental data: the fisher randomization test comment. Journal of the American statistical association 75 (371), p. 591–593. Cited by: §1. X. Shi, K. Q. Li, M. Yu, W. Miao, A. K. Kuchibhotla, M. Hu, and E. T. Tchetgen (2026) Theory for identification and inference with synthetic controls: a proximal causal inference framework. Journal of the American Statistical Association 0 (0), p. 1–13. Cited by: Appendix A, §2. X. Shi, W. Miao, J. C. Nelson, and E. J. Tchetgen Tchetgen (2020) Multiply robust causal inference with double-negative control adjustment for categorical unmeasured confounding. Journal of the Royal Statistical Society Series B: Statistical Methodology 82 (2), p. 521–540. Cited by: Appendix A, Table 3, §2, §4. R. Singh (2020) Kernel methods for unobserved confounding: negative controls, proxies, and instruments. arXiv preprint arXiv:2012.10315. Cited by: Appendix A. M. E. Sobel (2006) What do randomized studies of housing mobility demonstrate? causal inference in the face of interference. Journal of the American Statistical Association 101 (476), p. 1398–1407. Cited by: Appendix A, Table 3, §2. E. J. T. Tchetgen and T. J. VanderWeele (2012) On causal inference in the presence of interference. Statistical methods in medical research 21 (1), p. 55–75. Cited by: Appendix A, Table 3, §2. E. J. Tchetgen Tchetgen, A. Ying, Y. Cui, X. Shi, and W. Miao (2024) An introduction to proximal causal inference. Statistical Science 39 (3), p. 375–390. Cited by: Appendix A, Table 3, §1, §2, §4, §4, §4. Y. Wang and D. M. Blei (2019) The blessings of multiple causes: rejoinder. Journal of the American Statistical Association 114 (528), p. 1616–1619. Cited by: Appendix A, Appendix A, Table 3, §1, §2. S. M. Woodward, M. Tec, and F. Dominici (2024) An instrumental variables framework to unite spatial confounding methods. arXiv preprint arXiv:2411.10381. Cited by: Appendix A. Appendix Appendix A Extended Literature Review In this section, we provide a detailed review of the literature surrounding our proposed method. Table 3 provides a comprehensive summary of different literature related to the proposed framework. Proximal Causal Inference In epidemiological research, negative controls serve as an essential diagnostic. Grounded in subject-matter knowledge, a negative control uses a variable that is known to have no causal relationship with either the exposure or the outcome under study, yet related to the same potential sources of confounding (Lipsitch et al. 2010). Consequently, detecting a non-zero association between the negative control and the variable of interest provides direct evidence of residual bias or unmeasured confounding in the observational design. Miao and Tchetgen Tchetgen et al. (2017) extend this into the distinction between a negative control exposure and a negative control outcome, the direct conceptual ancestor of the treatment- and outcome-inducing proxies used in the proximal causal inference. The paper (Miao et al. 2018) introduced the theoretical foundation of proximal causal inference, formalizing this into full nonparametric point identification. Instead of requiring ignorability, the authors proposed identifying causal effects given treatment- and outcome-inducing proxies satisfying exclusion restrictions and completeness assumptions. They established nonparametric identification by introducing an outcome confounding bridge function that solves a conditional moment equation. Tchetgen Tchetgen et al. (2024) reformulated proximal causal inference entirely within the Rubin potential outcomes framework and extended proximal identification to longitudinal, time-varying treatment settings, formalizing the bridge function as the solution to a Fredholm integral equation of the first kind and connecting proximal identification to the theory of ill-posed inverse problems. It clarified the assumptions required for identification and established the bridge-function perspective that has become central to modern proximal causal inference. In the line of estimation literature, Shi et al. (2020) developed multiply robust estimators for proximal causal inference for categorical unmeasured confounders that remain consistent if a subset of nuisance models is correctly specified. The paper also derived semiparametric efficiency theory and influence-function-based inference. Cui et al. (2024) generalized proximal estimation to a semiparametric framework, deriving the efficient influence function for the proximal causal effect. With this generalized semiparametric framework, flexible machine learning estimators can be used to estimate nuisance functions while preserving multiple robustness and asymptotic efficiency. On the nonparametric-estimation side, Singh (2020) and Mastouri et al. (2021) develop kernel-based two-stage regression estimators for the bridge function, by formulating this estimation as a nonparametric instrumental-variable-style problem in a reproducing kernel Hilbert space. On the other hand, Bennett et al. (2019) and Dikkala et al. (2020) develop adversarial/minimax estimators for the underlying conditional moment restrictions, training a bridge network against an adversarially trained test-function network, the direct methodological ancestor of the adversarial momentum matching network used to train our own bridge function. Some literature extended this proximal framework to different application domains of non-i.i.d. structure. Egami and Tchetgen Tchetgen (2023) apply proximal identification to peer effects on a single observed network. They used double negative controls to identify causal effects under unmeasured network confounding. This paper considers interference between units, though restricted to a static network without spatial or temporal confounder dynamics. Shi et al. (2026) connected proximal causal inference with synthetic control methods by incorporating bridge functions into panel-data causal estimation. The authors utilized untreated donor units as negative-control proxies for the latent factor driving both outcomes and treatment timing, generalizing classical linear-factor-model synthetic control to a nonparametric setting. Table 3: Summary of related works. Paper Research Branch Is Causal? Data Type Hidden Confounder Interference Lipsitch et al.(2010) Negative Control No epidemiological Yes No Miao et al.(2017) Proximal No Cross-sectional Yes No Miao et al. (2018) Proximal Yes I.I.D. Yes No Tchetgen Tchetgen et al. Proximal Yes I.I.D., Yes No (2024) Temporal Shi et al. (2020) Proximal Yes Cross-sectional Yes No Cui et al. (2024) Proximal Yes Cross-sectional Yes No Sobel (2006) Interference Yes Clustered No Yes Tchetgen Tchetgen et al. Interference Yes General No Yes (2012) Giffin et al. (2023) Interference Yes Spatial No Yes Dupont et al. (2022) Spatial confounding Yes Spatial Yes No Papadogeorgou et al. Spatial confounding Yes Spatial Yes No (2019) Papadogeorgou et al. Spatial confounding, Yes Spatial Yes Yes (2023) interference Oprescu et al. (2026) Interference, Yes Spatiotemporal No Yes time-varying confounding Wang et al. (2019) Deconfounder Yes Cross-sectional Yes No Bica et al. (2020a) Deconfounder Yes Longitudinal Yes No Ali et al. (2024) Deconfounder Yes Spatiotemporal No Yes Khot et al. (2025) Deconfounder Yes Spatial Yes Yes Proposed Framework Proximal CI Yes Spatiotemporal Yes Yes Causal Inference under Interference Classical causal inference methods assume the Stable Unit Treatment Value Assumption (SUTVA), which excludes interference between observational units. This branch of causal inference literature relaxes the no-interference (SUTVA) assumption. Developed design-based and randomization-based estimators or hypothesis tests under a known interference structure by summarizing the neighboring unit’s interference to a unit through a low-dimensional exposure mapping and mediation formulation (Hudgens and Halloran 2008; Sobel 2006; Tchetgen and VanderWeele 2012; Aronow and Samii 2017; Forastiere et al. 2021). Partial-interference designs assume interference is confined to disjoint clusters of units (Sobel 2006; Liu and Hudgens 2014). Giffin et al. 2023 extend generalized-propensity-score methods to continuous, distance-decaying spatial interference. All of these methods identify direct and spillover effects under design-based or observed-covariate exchangeability, but assume all relevant confounders are measured. Another group of literature addresses unmeasured spatial confounding directly, without modeling interference jointly. Spatial+ (Dupont et al. 2022) is a two-stage regression approach that reduces spatial confounding by removing spatial dependence from covariates. It works by replacing original covariates with residuals obtained after regressing away their spatial trends, thereby preventing collinearity with spatial random effects without altering the main outcome model. Papadogeorgou et al. (2019) develop distance-adjusted propensity score matching (DAPSM) to proxy smooth unmeasured confounders via geographic proximity, and instrumental-variable designs exploit exogenous spatial variation such as wind direction or policy boundaries (Woodward et al. 2024). Papadogeorgou et al. (2023) is, to our knowledge, the only method addressing spatial confounding and interference jointly prior to the deconfounder-based methods discussed below, via a Bayesian model of the latent spatial field combined with an interference structure. In the deep-learning-based design paradigm, Oprescu et al. (2026) propose GST-UNet, combining a U-Net spatiotemporal encoder with iterative regression-based G-computation to jointly handle interference, spatial confounding, and temporal carryover together with time-varying confounding. However, following the classical g-formula it builds on, this requires that all confounders be measured rather than hidden. Deconfounding via Factor Models Another major research direction proposes to estimate substitute confounders through latent-factor model fit to multiple simultaneously observed causes. Wang et al. (2019) introduce the deconfounder method, showing that when a unit receives multiple treatments, their joint distribution can estimate a substitute for a shared unmeasured confounder, restoring point identification of the average treatment effect for each cause. Bica et al. (2020a) extend this deconfounder method to longitudinal panel data via a recurrent factor model and propose the time series deconfounder. Hatt et al. (2024b) extend it further to general sequential and network settings. Two very recent works adapt this idea to the spatial/spatiotemporal setting specifically. Ali et al. (2024) propose STCINet, pairing a U-Net with double attention gating for spatial interference with an autoencoder-based factor model, to reduce bias from hidden, time-varying confounding, and introduce direct (DATE) and indirect (IATE) effect estimands. Khot et al. (2025) observe more explicitly that spatial interference itself provides the multi-cause structure required for the deconfounder, as a unit’s own and neighboring treatments are jointly shaped by the same latent spatial field. This paper proposed the Spatial Deconfounder, which reconstructs a substitute confounder via a conditional variational autoencoder with a spatial (Gaussian Markov random field) prior, proving identification of direct and spillover effects under a latent-field sufficiency assumption. These deconfounder-based methods estimate substitute confounders through representation learning but generally lack nonparametric identification guarantees. D’Amour (2019) and Ogburn et al. (2019), commenting directly on (Wang and Blei 2019), show that a factor model fit to the observed causes does not, in general, nonparametrically identify the interventional distribution without additional untestable assumptions. Since the same observed variables cannot simultaneously serve as measurements of the confounder and as the treatments whose effects are being estimated, identifiability in practice requires restrictive auxiliary conditions such as no single-cause confounders, or (in the spatial case) smoothness priors and outcome-model separability. The recent spatial deconfounder papers acknowledge explicitly this as a limitation rather than a resolution. Our proximal approach avoids this specific tension by using proxy variables that are structurally distinct from the treatment being intervened on, at the cost of requiring dedicated negative-control variables rather than reusing the treatment vector itself. Appendix B Outcome Bridge Function Proof of Theorem 1: Potential Outcome Proof. The treatment-inducing proxy Z contains some parts on unobserved confounder U (assumption 3). Now, according to the towers property of iterative expectation, conditioning on (At,Zt,Xt,H¯t−1)(A_t,Z_t,X_t, H_t-1), the expectation of outcome YtY_t is E[Yt|At,Zt,Xt,H¯t−1]=E[E[Yt|At,Ut,Xt,H¯t−1)]|At,Zt,Xt,H¯t−1] splitE[Y_t|A_t,Z_t,X_t, H_t-1]=E[E[Y_t|A_t,U_t,X_t, H_t-1)]\\ |A_t,Z_t,X_t, H_t-1] split Assumptions 3 and 4 state that the outcome-inducing proxy W reflects the characteristics of U and is independent of A. From the proximal causal inference, we found that there exists a function hth_t such that E[Yt|At,Ut,Xt,H¯t−1)]=E[ht(At,Wt,Xt,H¯t−1)|Ut,Xt,H¯t−1] splitE[Y_t|A_t,U_t,X_t, H_t-1)]=E[h_t(A_t,W_t,X_t, H_t-1)\\ |U_t,X_t, H_t-1] split Now replacing this function hth_t into previous conditional expectation gives, E[Yt|At,Zt,Xt,H¯t−1]=E[E[ht(At,Wt,Xt,H¯t−1)E[Y_t|A_t,Z_t,X_t, H_t-1]=E[E[h_t(A_t,W_t,X_t, H_t-1) |Ut,Xt,H¯t−1]|At,Zt,Xt,H¯t−1]|U_t,X_t, H_t-1]|A_t,Z_t,X_t, H_t-1] We can express this conditional expectation in the following integral form E[E[ht(At,Wt,Xt,H¯t−1)|Ut,Xt,H¯t−1]|At,Zt,Xt,H¯t−1]=E[E[h_t(A_t,W_t,X_t, H_t-1)|U_t,X_t, H_t-1]|A_t,Z_t,X_t, H_t-1]= ∫ht(At,w,Xt,H¯t−1)P(w|u,Xt,H¯t−1) h_t(A_t,w,X_t, H_t-1)dP(w|u,X_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)dP(u|A_t,Z_t,X_t, H_t-1) According to Fubini’s theorem of we can swap integrals in this form. ∫ht(At,w,Xt,H¯t−1)P(w|u,Xt,H¯t−1) h_t(A_t,w,X_t, H_t-1)dP(w|u,X_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)dP(u|A_t,Z_t,X_t, H_t-1) =∫ht(At,w,Xt,H¯t−1)∫P(w|u,Xt,H¯t−1)= h_t(A_t,w,X_t, H_t-1) dP(w|u,X_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)dP(u|A_t,Z_t,X_t, H_t-1) From assumption 3, we know W is independent of X and Z, so we can include these variables in the second integral without any effect, ∫ht(At,w,Xt,H¯t−1)∫P(w|u,Xt,H¯t−1) h_t(A_t,w,X_t, H_t-1) dP(w|u,X_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)=dP(u|A_t,Z_t,X_t, H_t-1)= ∫ht(At,w,Xt,H¯t−1)∫P(w|u,Xt,Zt,At,H¯t−1) h_t(A_t,w,X_t, H_t-1) dP(w|u,X_t,Z_t,A_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)dP(u|A_t,Z_t,X_t, H_t-1) Now applying the conditional law of total probability, we get ∫ht(At,w,Xt,H¯t−1)∫P(w|u,Xt,Zt,At,H¯t−1) h_t(A_t,w,X_t, H_t-1) dP(w|u,X_t,Z_t,A_t, H_t-1) dP(u|At,Zt,Xt,H¯t−1)=dP(u|A_t,Z_t,X_t, H_t-1)= ∫ht(At,w,Xt,H¯t−1)P(w|At,Zt,Xt,H¯t−1) h_t(A_t,w,X_t, H_t-1)dP(w|A_t,Z_t,X_t, H_t-1) We can write the expectation form of this integral as: E[Yt|At,Zt,Xt,H¯t−1]=E[ht(At,Wt,Xt,H¯t−1)E[Y_t|A_t,Z_t,X_t, H_t-1]=E[h_t(A_t,W_t,X_t, H_t-1) |At,Zt,Xt,H¯t−1]|A_t,Z_t,X_t, H_t-1] This eliminates the hidden confounder U from the conditional expectation entirely. By the completeness assumption, the function hth_t has a unique solution; hence, for treatment a we can say the potential outcome is ψ(a)=E[Yt(a)]=E[ht(At=at,Wt,Xt,H¯t−1)|At=at,Zt,Xt,H¯t−1] splitψ(a)=E[Y_t(a)]=E[h_t(A_t=a_t,W_t,X_t, H_t-1)\\ |A_t=a_t,Z_t,X_t, H_t-1] split ∎ Appendix C Consistency of Diffusion-Based Treatment Proxy Learning The identification results of the proposed method assume the existence of a valid treatment-induced proxy Z that satisfies the proximal assumptions. These results are independent of the particular estimation procedure used to obtain the proxy. In this section, we discuss the theoretical connection between the proposed diffusion-based proxy learner and the proximal identification framework by showing that the learned proxy asymptotically recovers the latent treatment-induced representation required by the identification theory. Let, Z⋆Z denote the ideal treatment-inducing proxy satisfying the proxy exclusion restriction and the spatiotemporal completeness assumption. The proposed diffusion network produces an estimate Z^=fZ(X,A) Z=f_Z(X,A), where (fZ)(f_Z) consists of a spatiotemporal Transformer encoder followed by a conditional denoising diffusion model. The Transformer encoder first extracts long-range temporal dependencies and spatial interactions among neighboring units, while the diffusion model learns the conditional latent distribution through iterative denoising. Unlike deterministic models, the diffusion process estimates an entire probability distribution over latent confounding representations, enabling the learned proxy to capture both the variability induced by hidden confounding and interference. Proposition C.1: Consistency of Diffusion-Based Treatment Proxy Learning Suppose • The Transformer encoder is sufficiently expressive to approximate the conditional representation of the observed treatment and covariates; • The diffusion model consistently estimates the conditional score function associated with the latent treatment distribution; • The optimization converges to the global minimum of the diffusion objective; then the learned proxy Z^=fZ(X,A) Z=f_Z(X,A) converges in probability to the latent treatment-inducing proxy with infinite training samples and steps, Z^→Z⋆. Z PZ . The previous proposition establishes that the diffusion network consistently estimates the latent treatment-inducing proxy Z Z. Now we need to validate whether replacing the ideal proxy by its learned estimate affects the identification results established earlier. The following corollary shows that asymptotic identification of the bridge function is preserved. Corollary C.1: Identification with Learned Diffusion Proxy If the conditions of the diffusion consistency proposition and the assumptions of the spatiotemporal proximal identification theorem hold and the proximal proxy assumptions are satisfied by the latent proxy Z Z, then replacing the ideal proxy Z⋆Z with the learned diffusion proxy Z Z preserves asymptotic identification of the outcome bridge function. Specifically, E[hθ(W,A,X)∣Z^,X,A]→E[hθ(W,A,X)∣Z⋆,X,A]E[h_θ(W,A,X) Z,X,A] PE[h_θ(W,A,X) Z ,X,A] Proof: From the previous proposition of consistency of diffusion-based treatment proxy learning, we get Z^→Z⋆ Z []PZ The outcome bridge function is identified under the spatiotemporal completeness assumption and is continuous with respect to its proxy argument. Therefore, by the Continuous Mapping Theorem, we have h(W,A,X;Z^)→h(W,A,X;Z⋆).h(W,A,X; Z) []Ph(W,A,X;Z ). Since the identified potential outcomes E[h(W,A,X)|Z,A,X]E[h(W,A,X)|Z,A,X] are continuous functions of the bridge function, Slutsky’s theorem implies that E[hθ(W,A,X)∣Z^,X,A]→E[hθ(W,A,X)∣Z⋆,X,A]E[h_θ(W,A,X) Z,X,A] PE[h_θ(W,A,X) Z ,X,A] Hence, replacing the ideal treatment-induced proxy with the learned diffusion proxy preserves asymptotic identification of all causal estimands established in the preceding sections. Appendix D Extended Experimental Details Synthetic Data We construct a synthetic spatiotemporal environment on a 28×2828× 28 (NX×NY)(N_X× N_Y) grid and generate two synthetic datasets over T=5000T=5000 time steps to jointly exercise hidden confounding and interference. For the synthetic dataset 1, we implement the following data generation process: Ut=Ut−1+dt∇2Ut−1+dtαU¯,t−1+ϵt,U_t=U_t-1+dt∇^2U_t-1+dtα U_N,t-1+ _t, Xt=Xt−1+dt∇2Xt−1+dtβUt−1++ϵt,X_t=X_t-1+dt∇^2X_t-1+dtβ U_t-1++ _t, At=At−1+dt∇2At−1+dtγ1Ut−1+dtγ2Xt−1+ϵt,A_t=A_t-1+dt∇^2A_t-1+dt _1U_t-1+dt\; _2X_t-1+ _t, Yt=Yt−1+dt∇2Yt−1+dtω1Ut−1+dtω2Xt−1+dtω3At−1+dtω4A¯,t−1+ϵt, splitY_t=Y_t-1+dt∇^2Y_t-1+dt _1U_t-1+dt _2X_t-1+dt _3A_t-1\\ +dt _4 A_N,t-1+ _t, split where dtdt is the rate of change, ∇2∇^2 is the Laplacian operator used to simulate the diffusion process, ϵε is independent noise, and A¯,t−1 A_N,t-1 is the lagged mean neighborhood treatment as interference with 8 adjacent grid locations considered as the neighborhood. Synthetic dataset 1 contains 4 simulated variables. Here, U acts as the hidden confounder influencing all other variables and is never included in the dataset for training and testing. X is the observed confounder, A is the treatment, and Y is the outcome variable. To generate counterfactual data (Acf)(A_cf), first, we apply an intervention to the treatment variable by increasing values a specific region A[nX1:nX2,nY1:nY2]=A[10:20,10:20]A[n_X1:n_X2,n_Y1:n_Y2]=A[10:20,10:20] of the treatment variable with 1.51.5 times and then following the diffusion and interference process to calculate the whole spatiotemporal field. Based on this intervened treatment AcfA_cf, we generate counterfactual outcome YcfY_cf following the data generation process. The synthetic dataset contains a total of eight simulated variables: two hidden variables, two confounders, two covariates, a treatment, and an outcome. The data generation process of this dataset is as follows: U1t=U1t−1+dt∇2U1t−1+dtα1U1¯,t−1+ϵt,U1_t=U1_t-1+dt∇^2U1_t-1+dt _1 U1_N,t-1+ _t, U2t=U2t−1+dt∇2U2t−1+dtα2U2¯,t−1+ϵt,U2_t=U2_t-1+dt∇^2U2_t-1+dt _2 U2_N,t-1+ _t, X1t=X1t−1+dt∇2X1t−1+dtα3U1t−1+dtα4U2t−1+ϵt,X1_t=X1_t-1+dt∇^2X1_t-1+dt _3U1_t-1+dt _4U2_t-1+ _t, X2t=X2t−1+dt∇2X2t−1+dtα5U1t−1+dtα6U2t−1+ϵt,X2_t=X2_t-1+dt∇^2X2_t-1+dt _5U1_t-1+dt _6U2_t-1+ _t, C1t=C1t−1+dt∇2C1t−1+dtα7U1t−1+dtα8U2t−1+ϵt,C1_t=C1_t-1+dt∇^2C1_t-1+dt _7U1_t-1+dt _8U2_t-1+ _t, C2t=C2t−1+dt∇2C2t−1+dtα9U1t−1+dtα10U2t−1+ϵt,C2_t=C2_t-1+dt∇^2C2_t-1+dt _9U1_t-1+dt _10U2_t-1+ _t, At=At−1+dt∇2At−1+dtα11U1t−1+dtα12U2t−1+dtβ1X1t−1+dtβ2X2t−1+dtγ1C1t−1+ϵt, splitA_t=A_t-1+dt∇^2A_t-1+dt _11U1_t-1+dt _12U2_t-1+\\ dt\; _1X1_t-1+dt\; _2X2_t-1+dt\; _1C1_t-1+ _t, split Yt=Yt−1+dt∇2Yt−1+dtα13U1t−1+dtα14U2t−1+dtβ3X1t−1+dtβ4X2t−1+dtγ2C2t−1+dtω1At−1+dtω2A¯,t−1+ϵt, splitY_t=Y_t-1+dt∇^2Y_t-1+dt _13U1_t-1+dt _14U2_t-1\\ +dt _3X1_t-1+dt _4X2_t-1+dt\; _2C2_t-1\\ +dt _1A_t-1+dt _2 A_N,t-1+ _t, split The symbols convey the same meaning as synthetic dataset 1. In the dataset, U1U1 and U2U2 are used as hidden variables, and X1X1 and X1X1 are common confounders for both treatment A and outcome Y. C1C1 is the covariate of A, and C2C2 is the covariate of Y, both are influenced by hidden confounders. The counterfactual treatment and outcome are simulated using the same procedure as the previous dataset; the treatment is increased by 1.51.5 times in the spatial region A[10:20,10:20]A[10:20,10:20], and then we generate the counterfactual AcfA_cf and YcfY_cf following the simulation process. Baseline Models We evaluate the proposed framework using synthetic spatiotemporal datasets with known ground truth that incorporate key challenges: hidden confounder, interference, spillover, confounding, and temporal carryover. Using the synthetic data generation process, we compare our proposed framework against a set of representative baselines from three categories:(i) purely predictive spatiotemporal architectures, (i) interference-aware causal estimators, and (i) hidden-confounding-aware causal estimators. The implementation code of these baselines is available at https://anonymous.4open.science/r/Spatiotemporal-Proximal-Causal-Inference-2BF6. From the first category, we used UNet: a convolutional-LSTM-based U-Net architecture and Transformer: a spatiotemporal Transformer encoder-decoder architecture. We considered GST-Unet (Oprescu et al. 2026), and STCINET (Ali et al. 2024) models from the second category and the Spatial Deconfounder method from the third category. The UNet baseline model is built using a ConvLSTM layer and an encoder-decoder structure using 2D convolution with attention gating, following the baseline used in GST-Unet literature. The ConvLSTM layer is used to aggregate temporal and spatial features from the input data, and then the encoder-decoder is used to learn latent features. The baseline Transformer model is designed using a spatiotemporal transformer architecture with a positional embedding to capture spatial details and temporal patterns from the input data. These models are trained using the AdamW optimization algorithm. The STCINET model is used from the official code repository published by the authors (https://github.com/iharp-institute/causality-for-arctic-amplification/tree/main/stcinet). This model contains an encoder-decoder architecture designed using ConvLSTM2D and Conv2D neural network layers, with attention gating to learn subtle features from the input spatial grid. Similarly, the GST-Unet (https://github.com/moprescu/GSTUNet/tree/main) and the Spatial Deconfounder (https://github.com/moprescu/Spatial-Deconfounder/tree/main) baseline models are also acquired from the official code repository shared by the authors. The GST-Unet model implements a U-Net encoder-decoder architecture with ConvLSTM2D and Conv2D, similar to the STCINET model. Additionally, GST-Unet contains G-computation heads implemented as a feed-forward network to predict outputs for future time steps. The Spatial Deconfounder model is implemented as a two-stage architecture. It utilizes a Conditional Variational Autoencoder (CVAE) to estimate a substitute for hidden confounders using an encoder-decoder module design using 2D Convolutional layers. Using the estimated substitute and observed variables, the outcome module predicts the target output. The outcome module is designed with a U-Net structure leveraging hierarchical Conv2D downsampling and upsampling operations. Table 1 reports prediction performance across all methods on the synthetic datasets. From the comparative analysis, we see that some baseline methods performed better than the proposed method for dataset 2, though these methods are not designed to handle the hidden confounders. To assess whether this quantitative gap corresponds to a qualitative difference in the spatial structure of the predicted counterfactuals — rather than only their aggregate accuracy — Figure 3 visualizes the predicted outcome fields for three methods at time step 1000. The visualization reveals a further distinction beyond aggregate accuracy: while U-Net and the Spatial Deconfounder achieve comparatively low pointwise error, their predicted counterfactual fields fail to preserve the spatial structure of the ground-truth outcome. This suggests these methods are fitting the marginal distribution of outcome values well without correctly capturing where treatment and spillover effects occur spatially, which pointwise error metrics alone do not penalize. In contrast, the proposed method’s predicted field closely tracks the ground truth’s spatial structure rather than treating each location’s outcome as an independent regression target. Figure 3: Visualization of ground truth counterfactual outcome (Ycf)(Y_cf) and predicted counterfactual outcomes by the proposed model and different baseline models at time step 1000. Hyperparameters To find the best hyperparameters for baseline methods, we started using the parameters suggested by the authors and gradually tuned those values to obtain better evaluation results. The results reported in the comparative analysis of the main article are obtained with tuned hyperparameters. The parameters used to generate evaluation results are given here. • Unet: batch_sizebatch\_size = 16, epochsepochs = 40; learning_ratelearning\_rate = 5e-4, weight_decayweight\_decay=1e-5, optimizeroptimizer= AdamW, and scheduler_patiencescheduler\_patience=5. • Transformer: batch_sizebatch\_size = 16, epochsepochs = 40; learning_ratelearning\_rate = 5e-4, weight_decayweight\_decay=1e-5, optimizeroptimizer= AdamW, and scheduler_patiencescheduler\_patience=5. • STCINET: batch_sizebatch\_size = 64, epochsepochs = 60; learning_ratelearning\_rate = 1e-3, weight_decayweight\_decay=exp(-0.1), optimizeroptimizer= Adam, scheduler_patiencescheduler\_patience=5, and loss_weightsloss\_weights = [0.25, 0.75]. • GST-Unet: batch_sizebatch\_size = 64, epochsepochs = 50; learning_ratelearning\_rate = 5e-4, learning_ratedecaylearning\_rate_decay=exp(-0.1), optimizeroptimizer= Adam, and scheduler_patiencescheduler\_patience=8. • Spatial Deconfounder: epochs_CVAEepochs\_CVAE = 50, epochs_HEADepochs\_HEAD = 30, radiusradius = 1, batch_sizebatch\_size = 16, optimizeroptimizer= AdamW, learning_ratelearning\_rate = 1e-5, and weight_decayweight\_decay=1e-3, • Proposed Method: batch_sizebatch\_size = 32, epochsepochs = [200, 400]; learning_ratelearning\_rate = 1e-4, optimizeroptimizer= Adam, beta_klbeta\_kl = 0.2, lambda_difflambda\_diff = 0.2, lambda_cmi_Zlambda\_cmi\_Z = 0.5, lambda_out_Wlambda\_out\_W = 1.0, lambda_cmi_Wlambda\_cmi\_W = 0.5, lambda_bridge_mselambda\_bridge\_mse = 1.0, and lambda_momentlambda\_moment = [0.5, 0.8]. Appendix E Ablation Study To assess the contribution of each component of the operationalization of the proposed proximal framework, we perform a comprehensive ablation study. We consider the following variants of the proposed framework: (i) without diffusion in proxy learning, (i) without stabilized weights module, and (i) without bridge momentum network. Figure 4 demonstrates the components of the proposed framework removed for different ablation variations. For the first variation of the proposed framework without the diffusion block in the treatment-inducing proxy learning, we just removed the diffusion module from the encoder and connected the output of the spatiotemporal transformer to the input of the decoder. Other building blocks of the proposed framework are unchanged and trained with the same settings as the original end-to-end framework. Next, we want to assess the contribution of the stabilized weights module introduced in the proposed framework to reduce the bias of highly frequent treatment values on the potential outcome prediction. To study this, we disconnect the stabilized weights block from the outcome bridge module, so each treatment value is weighted equally despite their density distribution. Finally, in the last ablation variation, we turned off the bridge momentum network from the outcome bridge module to check its contribution to identifying the outcome bridge function considering only outcome-inducing proxy W and without applying a condition on the treatment-inducing proxy Z. Table 2 shows that the complete end-to-end framework consistently outperforms every ablation variant. Removing the stabilized weight module substantially degrades outcome prediction, reflecting the treatment imbalance induced by the data-generating process: because treatment assignment includes spatial and temporal autocorrelation (Section D), the majority of grid locations and time steps take treatment values within a narrow, densely populated range, while a comparatively small fraction of the spatiotemporal field is exposed to sparse, extreme treatment values. Without reweighting, training gradients are dominated by this high-density region, and the bridge function h^θ h_θ underfits the sparse region, precisely the region an intervention analysis is often most interested in analyzing. Removing the bridge momentum network likewise degrades performance, demonstrating that directly optimizing the outcome prediction loss increases bias. Moreover, consistent with our identification argument that direct regression of Y based on (W,A,X)(W,A,X) alone, without conditioning on the treatment-inducing proxy Z, recovers E[Y|W,A,X]E[Y|W,A,X] rather than the true confounding bridge function (Section B). This confirms empirically that the momentum network is not merely a regularizer but is necessary for the learned h^θ h_θto satisfy the identifying bridge equation, matching the theoretical distinction we draw between the two objectives. Figure 5 visualizes the spatial structure of the predicted counterfactual outcomes of these ablation variants. Which demonstrate the spatial structure of the predicted output by the proposed framework close to the ground truth outcome. Figure 4: Visualization of components held out for different ablation variations of the proposed framework. Figure 5: Comparison of ground truth counterfactual outcome (Ycf)(Y_cf) and predicted counterfactual outcomes by different ablation variations of the proposed framework at time step 1000.