Spatial Causal Tensor Completion for Multiple Exposures and Outcomes: An Application to the Health Effects of PFAS Pollution

Spatial Causal T ensor Completion for Multiple Exp osures and Outcomes: An Application to the Health Effects of PF AS P ollution Xiao dan Zhou Departmen t of Statistics, North Carolina State Univ ersity Brian J Reic h Departmen t of Statistics, North Carolina State Univ ersity Sh u Y ang ∗ Departmen t of Statistics, North Carolina State Univ ersity Marc h 18, 2026 Abstract P er- and polyfluoroalkyl substances (PF AS) are t ypically encountered as mixtures of distinct c hemicals with distinct effects on multiple health outcomes. Estimating join t causal effects using spatially-dep enden t observed data is challenging. W e propose a spatial causal tensor completion framew ork that jointly mo dels multiple exp osures and outcomes within a lo w-rank tensor structure, while adjusting for observed confounders and laten t spatial confounders. This metho d combines a lo w-rank tensor represen tation to p o ol information across exp osures and outcomes with a sp ectral adjustmen t step that incorp orates graph-Laplacian eigen vectors to approximate unmeasured spatial confounders, implemented via a pro jected-gradient descen t algorithm. This framework enables causal inference in the presence of unmeasured spatial confounding and p erv asiv e missingness of p otential outcomes. W e establish theoretical guaran tees for the estimator and ev aluate its finite-sample p erformance through extensiv e sim ulations. In an application to national PF AS monitoring data, our approac h yields more conserv ativ e and credible causal relationships b et ween PF OA and PF OS exp osure and 13 c hronic disease outcomes compared with existing alternativ es. K eywor ds: Unmeasured Spatial Confounding; Mixture Mo deling; Poten tial Outcomes; Prop ensit y Score; In verse Probabilit y W eigh ting. ∗ The authors thank Dr. Ana Rapp old of the United States Environmen tal Protection Agency for her help accessing and interpreting the data. This w ork w as partially supp orted by National Institutes of Health gran ts R01ES031651 and R01ES036270. 1 1 In tro duction P er- and p olyfluoroalkyl substances (PF AS) are a class of synthetic c hemicals widely used since the 1950s due to their resistance to water, oil, and heat ( Gaines 2023 ). Because of their p ersistence in the environmen t, they are often referred to as ‘forev er chemicals,’ raising substan tial public health concerns. PF AS c hemicals rarely o ccur in isolation; drinking water con tamination typically inv olves mixtures of comp ounds such as PFO A and PF OS. Multiple PF AS c hemicals hav e b een asso ciated with a range of adv erse outcomes, including altered imm une and th yroid function, liv er disease, lipid and insulin dysregulation, kidney disease, adv erse reproductive and dev elopmen tal outcomes, and cancer ( F enton et al. 2021 ). Y et despite growing atten tion, evidence for many health outcomes remains inconsisten t and inconclusiv e ( NASEM et al. 2022 ), p oten tially due to the complexity of the problem. In this w ork, w e study the effects of PF AS exp osures on multiple health outcomes. This naturally gives rise to a multiple-exposure, m ultiple-outcome setting. The analysis is complicated b y t w o ma jor c hallenges: first, PF AS are typically encountered as mixtures of distinct c hemicals with different sources and health effects, making it difficult to isolate the contributions of individual comp onen ts; second, b oth exp osures and outcomes exhibit spatial patterns, introducing the p ossibility of unmeasured spatial confounding ( Reich et al. 2021 ). These complexities motiv ate the dev elopment of new statistical approac hes tailored to this problem. A large statistical literature has developed to quantify the effects of exp osure mixtures ( Joub ert et al. 2022 , Kang et al. 2023 ). Widely used approac hes include w eighted quantile sum regression ( Carrico et al. 2015 ), Bay esian k ernel machine regression ( Bobb et al. 2015 ), and quan tile-based g-computation ( Keil et al. 2020 ). How ever, they generally fo cus on a single outcome and do not accoun t for unmeasured spatial confounding. Univ ariate-outcome spatial metho ds ha v e been prop osed to adjust for spatially structured confounders ( Thaden 2 & Kneib 2018 , Papado georgou et al. 2019 , Osama et al. 2019 , Keller & Szpiro 2020 , Schnell & Papadogeorgou 2020 , Marques et al. 2022 , Dup ont et al. 2022 , Guan et al. 2023 , Gilb ert et al. 2021 , Wiec ha et al. 2025 ) but these approaches are limited to one exp osure and one outcome at a time, making them unsuitable for analyzing PF AS mixtures with m ultiple health endp oin ts. P arallel to these developmen ts, existing w ork has reframed causal effect estimation as a matrix or tensor completion problem, where coun terfactual outcomes are reco vered under lo w-rank structures ( P oulos et al. 2021 , Auerbac h et al. 2022 , Zhen & W ang 2024 , Mao et al. 2024 ). This line of researc h enables principled imputation of missing p oten tial outcomes and extends naturally to m ultiv ariate settings. Ho wev er, these metho ds typically treat eac h exp osure indep enden tly , o verlooking latent relationships among exp osures. More recen t approaches in tegrate exposures in to a single tensor represen tation, allo wing for joint mo deling of exp osure effects and improv ed efficiency , as demonstrated in Abadie et al. ( 2024 ), Agarwal et al. ( 2020 ), Gao et al. ( 2024 ), and Gao et al. ( 2025 ). T wo recent parallel w orks also address m ultiple exp osures and outcomes with spatial confounding, though with differen t strategies. Prim et al. ( 2025 ) propose a spectral confounder adjustmen t that uses a three-w ay tensor ov er exp osure, outcome, and spatial scale to model causal effect, assuming that confounding bias dissipates at more lo cal spatial scales. Our approac h is more general in that we use a tensor for the p oten tial outcomes, rather than using tensors as comp onen ts of a linear mo del. W u & F ranks ( 2025 ) dev elop a latent factor panel approach for spatiotemp oral data that lev erages rep eated temp oral observ ations to identify causal effects under a factor confounding assumption, where unmeasured confounders are captured b y lo w-dimensional latent factors that ev olv e ov er time, com bined with partial interference assumptions that limit the degree of spatial spillo ver. By contrast, w e prop ose a spatial causal tensor completion framework sp ecifically designed 3 for cross-sectional spatial data with m ultiple binary exp osures and multiple outcomes. Our metho d com bines a low-rank T uc ker tensor representation ( K olda & Bader 2009 ) to p o ol information across exp osures and outcomes with a sp ectral adjustmen t step that incorp orates graph-Laplacian eigen vectors to appro ximate unmeasured spatial confounders, implemented via a spatial pro jected-gradient descen t algorithm. Unlik e panel approaches that exploit temp oral v ariation, our framework relies on spatial smo othness assumptions to enable causal inference in the presence of spatially patterned unmeasured confounding and p erv asiv e missingness of potential outcomes. W e establish rigorous theoretical guaran tees, including F rob enius error b ounds for the estimated p otential outcomes tensor that decomp ose in to in terpretable comp onents, and provide uncertain t y quan tification of a verage treatment effects through doubly robust v ariance estimation. W e demonstrate the framework’s adv antages through extensiv e simulations and an application to an national PF AS monitoring data. Our analysis cov ers 5,495 public w ater systems across 47 states and 13 chronic disease outcomes, reflecting typical con temp orary drinking-w ater exp osure levels rather than extreme con tamination scenarios. W e find that most previously reported PF AS–health asso ciations are substantially atten uated after adjustment for spatial confounding. Ho wev er, PFOS exp osure remains significantly asso ciated with elev ated risk of hypertension, to oth loss, asthma, and ob esit y—asso ciations that p ersist after adjustmen t and thus app ear less likely to reflect geographic artifacts alone. These findings underscore the imp ortance of accounting for spatial confounding in en vironmen tal health studies and suggest that targeted concerns ab out PF OS w arrant con tinued attention. This work makes three primary con tributions. First, w e develop a unified causal inference framew ork that jointly mo dels m ultiple exp osures and outcomes through low-rank tensor completion, enabling efficient estimation of factorial treatmen t effects in environmen tal mixtures studies. Second, we in tegrate graph-Laplacian sp ectral adjustmen t into the tensor 4 completion pro cedure to account for unmeasured spatial confounding—a p erv asive challenge in environmen tal epidemiology . Third, w e establish theoretical guaran tees for the prop osed estimator and v alidate its finite-sample p erformance through sim ulations. Our application to national PF AS data pro vides a large-scale, spatially-adjusted assessmen t of drinking-water PF OA and PFOS exp osure on chronic disease outcomes, yielding more conserv ative and credible effect estimates than existing metho ds. The remainder of the pap er is organized as follows. Section 2 describ es the motiv ating dataset. Sections 3 and 4 presen t the prop osed metho dology and its theoretical prop erties, with additional empirical ev aluations rep orted in Supplemen t D . Section 5 applies the prop osed approac h to the motiv ating data. Section 6 concludes. 2 PF AS Data Description The EP A’s PF AS Analytics T ools ( US EP A 2024 ) provide national-lev el PF AS data in Public W ater Systems (PWS). UCMR5 rep orts m ultiple PF AS chemicals. Among these, PF OA and PFOS are the most frequen tly detected comp ounds, while the remaining chemicals w ere detected at only a small n umber of lo cations. Therefore, we fo cus on PFO A and PFOS in this study . The EP A establishes minim um rep ortable limits (MRL) as the lo w est measurable concen- tration of a con taminant that is achiev able b y the ma jorit y of participating lab oratories, indep enden tly of the con taminant health effects. F or this reason, most often only the presence or absence of the chemical ab ov e the MRL is used as the exp osure measure. In addition, the health-based Maxim um Con taminant Lev el Goals (MCLGs) defined b y EP A for PFO A and PF OS are both zero ppt (also expressed as ng/L). Therefore, the exp osure is taken to b e the binary indicator of detectable PF O A and/or PF OS, given four p ossible exp osure levels: Neither PFO A or PF AS, only PF O A exp osure, only PFOS exp osure, and 5 b oth PF OA and PFOS exp osures. W e use health outcome data provided b y the Centers for Disease Control and Prev en tion (CDC) PLA CES dataset ( CDC 2023 ), whic h contains mo del-based census tract-level es- timates of disease prev alence among the adult p opulation. The 13 health outcomes are: arthritis (AR THRITIS), hypertension (BPHIGH), cancer (non-skin) or melanoma (CAN- CER), asthma (CASTHMA), coronary heart disease (CHD), c hronic obstructiv e pulmonary disease (COPD), diagnosed diab etes (DIABETES), high cholesterol among adults who ha v e ev er b een screened (HIGHCHOL), chronic kidney disease (KIDNEY), depression (MHL TH), ob esit y (OBESITY), stroke (STROKE), and all to oth loss among adults aged ≥ 65 y ears (TEETH). Demographic co v ariates were obtained from the 2020 US Census Bureau via the tidycensus R pac kage ( W alker et al. 2021 ); sp ecific cov ariates are describ ed in Section 5 . Health outcomes and demographic co v ariates were link ed to PWS using area-weigh ted a veraging to aggregate tract-lev el v ariables to the PWS level ( Supplemen t C.1 ). This yielded 5,495 PWS samples across 47 states in the con tinental United States. The T able in Supplemen t C.2 summarizes PWS-lev el exp osure groups and outcomes. 3 Spatial Causal T ensor Mo del In this section, w e in tro duce our spatial causal tensor framew ork and estimation strategy . Section 3.1 in tro duces notations, Section 3.2 formalizes the p oten tial-outcomes setup, and Section 3.3 states the required causal assumptions. Section 3.4 then outlines a three-step estimator with a spatial pro jected gradient descen t (S-PGD) algorithm. 6 3.1 Notations and T ensor Op erations Lo wercase letters (e.g., x, y ) denote scalars; b old lo wercase letters (e.g., x , y ) denote v ectors; b old upp ercase letters (e.g., A , B ) denote matrices; calligraphic upp ercase letters (e.g., X , Y , Z ) denote tensors. A tensor X ∈ R N 1 × N 2 × N 3 is a 3-w ay array indexed b y ( i 1 , i 2 , i 3 ) , where i k ∈ [ N k ] , and [ N k ] = { 1 , ..., N k } . The F robenius norm of X is defined as ∥X ∥ 2 F = P i 1 ,i 2 ,i 3 X 2 i 1 i 2 i 3 . The sp ectral and n uclear norms of tensors are defined formally in Supplemen t A . The inner pro duct is ⟨X , Y ⟩ = P i,l,o X iℓo Y iℓo . The mo de- k unfolding of a tensor X ∈ R N 1 ×···× N d rearranges the tensor into a matrix b y aligning the mo de- k fib ers, a one-dimensional slice of a tensor obtained b y fixing all indices but one, as columns of the resulting matrix, X (1) ∈ R N 1 × ( N 2 N 3 ) , X (2) ∈ R N 2 × ( N 1 N 3 ) , X (3) ∈ R N 3 × ( N 1 N 2 ) . The mo de- k tensor-matrix pro duct, denoted X × k M , with M ∈ R r × N k , yields a tensor of dimension where mo de k is replaced b y r . F or example, X × 1 M ∈ R r × N 2 × N 3 . This op eration inv olv es unfolding X along mo de k , applying matrix multiplication, and refolding to the tensor form. F or a comprehensive introduction to tensor operations, see Kolda & Bader ( 2009 ). 3.2 T ensor F orm ulation of P oten tial Outcomes Our prop osed metho d is a no v el attempt to impute the missing p otential outcomes using lo w-rank tensor completion under the p oten tial outcomes framework. Supp ose there are N observ ational units, eac h exp osed to a K -dimensional binary exp osure vector A i ∈ { 0 , 1 } K , whic h defines L = 2 K p ossible exp osure com binations. W e index each exp osure combination b y its decimal equiv alen t, denoted ℓ i ∈ { 1 , . . . , L } . F or clarit y , let a ( ℓ ) ∈ { 0 , 1 } K represen t the binary exp osure vector corresp onding to exp osure level ℓ . F or example, when K = 2 , if A i = ( A i 1 , A i 2 ) , then ℓ i = A i 1 + 2 A i 2 . In the motiv ating example, N = 5,495 denotes the n umber of geographic units (PWS) and K = 2 represents tw o PF AS exp osure v ariables (PF O A and PF OS). 7 Under the p otential outcomes framew ork ( R ubin 1974 ), each unit i has a set of p oten tial outcomes {Y ilo : ℓ = 1 , . . . , L, o = 1 , . . . , O } , one for every p ossible exp osure com bination ℓ and outcome dimension o . In the PF AS analysis, O = 13 outcomes are considered, as listed in Section 2 . Only the outcome corresp onding to the realized exp osure is observed, while the others remain coun terfactual. W e organize these outcomes in to a tensor Y ∈ R N × L × O (left of Figure 1 ). Similarly , the treatment assignment is represen ted by the tensor A ∈ R N × L × O , with A iℓo = 1 if unit i receiv ed exp osure ℓ , and A iℓo = 0 otherwise. Note that A iℓo do es not v ary across the third outcome-mo de dimension, so it effectiv ely reduces to an N × L matrix. W e main tain the tensor representation for consistency with the outcome tensor and to simplify presen tation. Each unit receives one exp osure, i.e. P ℓ A iℓo = 1 . W e assume that the observed outcome tensor follows a lo w-rank T uck er decomp osition (Righ t of Figure 1 ): Y = Y ∗ + E = G × 1 U 1 × 2 U 2 × 3 U 3 + E , E iid ∼ sub-Gaussian (0 , σ 2 ) . (1) In ( 1 ), G ∈ R r 1 × r 2 × r 3 is the core tensor and U 1 ∈ R N × r 1 , U 2 ∈ R L × r 2 and U 3 ∈ R O × r 3 are factor matrices corresp onding to spatial, exp osure, and outcome v ariation, resp ectiv ely . T o capture b oth measured and unmeasured sources of spatial v ariation, w e mo del U 1 = Z η Z + S η S , where Z ∈ R N × r 1 are measured co v ariates with effect η Z ∈ R r 1 × r 1 , and S ∈ R N × r 1 are unmeasured spatial comp onen ts with effect η S ∈ R r 1 × r 1 . By unmeasured spatial confounder, w e refer to latent v ariables that exhibit spatial structure and simultaneously influence exp osures and outcomes. The measured cov ariates Z and unmeasured spatial confounding S join tly affect b oth the outcomes and the treatmen t assignment mec hanism, as describ ed b elo w. In this representation, r 1 , r 2 and r 3 are less than or equal to N , L , and O , resp ectiv ely . The factor matrices U can b e view ed as the principal comp onen ts in eac h mo de, and the en tries 8 Figure 1: Illustration of the spatial causal tensor decomp osition. Left: a tensor of coun terfactual outcomes, where observ ed en tries (y ellow) and missing en tries (gra y) are arranged b y spatial units (PWS), ex- p osures, and outcomes. Right: T uck er de- comp osition into spatial ( U 1 ), exp osure ( U 2 ), and outcome ( U 3 ) factors with a lo w-rank core G . The spatial factor U 1 is further expressed as a com bination of measured cov ariates Z and latent spatial comp onen ts S . The illustration omits the noise term E for clarit y . in G en tail the level of interaction b et w een differen t comp onen ts. The lo w-rank assumption on Y ∗ can b e justified in the follo wing wa y . First, the heterogeneity among lo cations ma y b e attributed to r 1 laten t spatial factors; e.g., so cio demographic trends. Second, the v ariation from the PF AS exp osure may b e summarized into r 2 k ey factors, such as c hemicals with shared sources or physical, c heminformatic prop erties. Finally , the m ultiple health resp onses ma y b e related due to r 3 hidden biological path wa ys. Our parameter of interest is the av erage treatmen t effect (A TE) defined for any exp osure ℓ (relativ e to the reference level L ) and outcome o by: θ ∗ ℓo = E  P N i =1 ( Y iℓo − Y iLo ) / N  . 3.3 Causal Assumptions T o formalize the exp osure mec hanism, the exp osure probability for unit i (with cov ariates Z i and unobserved spatial factor S i ) to receiv e exp osure ℓ is π ( ℓ | Z i , S i ) = P ( A i = a ( ℓ ) | Z i , S i ) . (2) W e adopt the following classical causal assumptions to justify our estimating strategy . Assumption 1 (SUTV A; Stable Unit T reatment V alue Assumption) (1) the p o- 9 tential outc omes for any unit do not vary with the tr e atment assigne d to other units; (2) ther e ar e no differ ent versions of e ach tr e atment level le ading to differ ent p otential outc omes. Assumption 1 rules out in terference b et ween geographic units and requires that each exp osure condition is well-defined. The no-interference comp onent requires that PF AS exp osure in one geographic unit do es not affect health outcomes in another. This is plausible in our setting b ecause exp osure is measured through drinking w ater systems that serve distinct p opulations. Poten tial violations could nonetheless arise through p opulation mobilit y or en vironmental contamination spreading b etw een neigh b oring regions. W e interpret our results as direct effects of lo cal drinking water exp osure, ac knowledging that spillo v er effects, if presen t, are not captured. The consistency comp onen t requires careful in terpretation: our binary exp osure indicator aggregates heterogeneous concentration lev els ab ov e the detection limit. Consequen tly , estimated effects represent a v erages ov er the observed distribution of exp osure in tensities rather than effects of a precisely defined exp osure level. Assumption 2 (Latent Ignorability) A i ⊥ {Y iℓo : ℓ, o } | ( Z i , S i ) . In other wor ds, Z i and S i ac c ount for al l c onfounders influencing tr e atments and outc omes. Assumption 2 extends the standard no-unmeasured-confounding condition by allo wing for unobserved spatial confounders, suc h as unmeasured industrial activit y or correlated en vironmental exp osures that v ary smo othly ov er space. W e mo del these through low- dimensional latent spatial comp onents approximated via sp ectral adjustmen t using the graph Laplacian eigenbasis. This approac h is effectiv e for confounders v arying at spatial scales captured by the selected eigen vectors, but may fail for confounders v arying at finer scales than our adjacency structure can resolve or for confounders that lack spatial structure altogether. The main threats to this assumption include: non-spatial confounders not captured by our measured cov ariates (such as b ottled water); and temp oral misalignment, 10 as our exp osure measurements (2023–2024) and health outcomes (2017–2018) span different p erio ds, though the long en vironmental p ersistence and biological half-liv es of PF AS partially mitigate this concern. Assumption 3 (Latent Positivit y) Ther e exists a c onstant p min ∈ (0 , 1 2 ) such that for every unit i and exp osur e level ℓ ∈ { 1 , . . . , L } , p min ≤ P  A i = a ( ℓ ) | Z i , S i  ≤ 1 − p min . That is, e ach unit has a nonzer o pr ob ability of r e c eiving any exp osur e c ombination, ruling out deterministic assignment by c ovariates or lo c ation. Assumption 3 requires that eac h exp osure combination has nonzero probability of o ccurring, ruling out deterministic assignment b y geograph y or cov ariates. Limited ov erlap is a kno wn c hallenge in en vironmen tal studies, particularly when con tamination is highly lo calized. W e empirically assess this assumption b y examining the distribution of estimated prop ensit y scores in our application ( T able 3 in Supplement C.2 ). The results indicate mo derate o v erlap: only limited units ha v e prop ensit y scores b elow 0.01, though a notable fraction of exp osed units (10 to 48%) ha ve scores b elow 0.05, with the PFO A-only and PFOS-only groups exhibiting the most limited supp ort. 3.4 A Three-Step Estimating Pro cedure W e pro ceed in three steps. Step 1 uses an unw eighted S-PGD fit with a Laplacian eigen basis to learn the laten t spatial comp onen t b S (0) while estimating the T uck er factors from observed en tries. Step 2 fits the exp osure mo del b π ( ℓ | Z i , b S (0) i ) and constructs inv erse-probability w eigh ts c W (0) that accoun t for exp osure-induced missingness. Step 3 re-runs S-PGD with the weigh ts to re-estimate the spatial comp onen t and complete the full p oten tial-outcome tensor, yielding b S and b Y , and th us the final A TE estimates. 11 3.4.1 Step 1: Sp ectral adjustmen t for spatial confounding Let Q b e the graph Laplacian built from the spatial adjacency matrix, with eigendecom- p osition Q = Φ Λ Φ ⊤ . Low-frequency eigen vectors in Φ capture broad spatial trends. W e assume the unmeasured spatial comp onen t S follo ws a Gaussian Marko v random field that can b e well-appro ximated b y the span of the first k lo w-frequency eigenv ectors, S ≈ Φ 1: k β , U 1 ≈ Z η Z + Φ 1: k β , w e denote [ Φ 1 , Φ 2 , ..., Φ k ] as Φ 1: k , with β ∈ R k × r 1 are the corresponding co efficien ts. W e obtain b S (0) b y minimizing the unw eighted reconstruction loss arg min G , U 2 , U 3 , η Z , β ( ||A ◦ ( Y − G × 1 ( Z η Z + Φ 1: k β ) × 2 U 2 × 3 U 3 ) || 2 F ) . using the S-PGD solv er describ ed in Section 3.4.4 . This step provides an initial estimate of the latent spatial comp onen t. 3.4.2 Step 2: Prop ensit y score and weigh t estimation By Assumptions ( 1 )–( 3 ), in verse-probabilit y w eigh ting creates a pseudo-p opulation in which exp osure assignment is indep endent of the p oten tial outcomes. Consequently , the exp osure- induced missingness in Y is ignorable after weigh ting, enabling un biased tensor completion under the T uck er model. Sp ecifically , we estimate the exp osure probabilities b π ( ℓ | Z i , b S (0) i ) using multinomial logistic regression, and further construct in v erse probability w eigh ts c W (0) , b π ( ℓ | Z i , b S (0) i ) = P ( A i = a ( ℓ ) | Z i , b S (0) i ) , c W (0) iℓo = 1 ( A i = a ( ℓ ) ) b π ( ℓ | Z i , b S (0) i ) , ∀ ℓ, where the indicator function 1 ( A i = a ( ℓ ) ) ensures that weigh ts are only applied to observed exp osure-outcome pairs. Since inv erse w eigh ting can b ecome unstable in settings with many 12 treatmen t categories, calibration w eigh ting may b e emplo yed as a more stable alternative ( Y ang 2018 , Gao et al. 2023 ). 3.4.3 Step 3: W eigh ted T ensor Completion In Step 1, the initial spatial comp onen t b S (0) is learned from an unweighte d tensor fit that do es not yet account for the exp osure mechanism. Using the w eights from Step 2, we re-estimate b oth the spatial comp onent and the tensor to correct for exp osure-induced selection, yielding refined b S and b Y . Sp ecifically , w e refit the tensor with a w eigh ted loss, arg min G , U 2 , U 3 , η Z , β ( ∥ q c W (0) ◦ A ◦ ( Y − G × 1 ( Z η Z + Φ 1: k β ) × 2 U 2 × 3 U 3 ) ∥ 2 F ) , using the S-PGD. This rew eighted tensor fit pro vides the completed p otential-outcome tensor b Y , which can directly feed into the av erage treatment effector estimators with standard errors describ ed in Section 4.3 . 3.4.4 Spatial Pro jected Gradient Descent Algorithm Classical PGD implicitly assumes a fixed pro jection space for factors in any mo de. In con trast, in the S-PGD algorithm–whic h w e use in Steps 1 and 3 abov e–up dates the pro jection space b y selecting Laplacian eigen v ectors in a step wise fashion. At eac h step, we expand the pro jection b y one candidate eigenv ector, using w arm-start from the previous solution, and retain the expansion only if it improv es a BIC-type criterion. This allo ws the algorithm to adaptiv ely capture spatial structure rather than committing to a pre-sp ecified set of eigenv ectors. In the ab o ve, the T uck er ranks r 1 , r 2 , r 3 are assumed to b e known. In practice, they are selected via cross-v alidation in Step 1 of the estimating pro ceduce. Complete algorithmic details, including BIC criterion, warm-start initialization scheme and stopping criteria, are pro vided in Supplemen t B . 13 4 Theoretical Prop erties 4.1 T ensor Structural Assumptions Before presenting our main theorem, we introduce additional structural assumptions on the laten t spatial v ariables and on the lo w-rank tensor mo del. These assumptions provide the tec hnical foundation needed to establish the asymptotic b ounds for our weigh ted spatial tensor completion estimator. Assumption 4 (Unmeasured Smo oth Spatial V ariables) L et S ∈ R N × r 1 denote r 1 smo oth gr aph signals over a c onne cte d, undir e cte d gr aph G = ( V , E ) with | V | = N . L et Q b e the normalize d gr aph L aplacian with eigende c omp osition Q = ΦΛΦ ⊤ , wher e Φ = [ ϕ 1 , . . . , ϕ N ] , wher e Φ has orthonormal c olumns, and 0 = λ 1 ≤ λ 2 ≤ · · · ≤ λ N . A ssume that e ach c olumn S ( ℓ ) of S admits the sp e ctr al exp ansion: S ( ℓ ) = P N j =1 c ( ℓ ) j ϕ j , with | c ( ℓ ) j | ≤ C · j − α/β , for some c onstants C > 0 and α > β > 0 , uniformly over ℓ = 1 , . . . , r 1 . Assumption 4 formalizes the idea that unmeasured spatial confounders v ary smo othly o ver geograph y . By requiring that their eigen-expansion co efficients deca y at a p olynomial rate, this assumption ensures that most of the spatial v ariation is captured b y a small num b er of the leading eigen vectors of the Laplacian. Assumption 5 (Incoherent T ensor P arameter) A ssume U 1 , U 2 and U 3 have ortho g- onal c olumns, and ther e exists some c onstant µ 0 and µ 1 such that Y ∗ ∈ C ( r 1 , r 2 , r 3 , µ 0 , µ 1 ) , wher e max  N r 1 ∥ U 1 ∥ 2 2 , ∞ , L r 2 ∥ U 2 ∥ 2 2 , ∞ , O r 3 ∥ U 3 ∥ 2 2 , ∞  ≤ µ 0 , max k ∈{ 1 , 2 , 3 } ∥M ( k ) ( G ) ∥ ≤ µ 1 v u u t N LO µ 3 / 2 0 ( r 1 r 2 r 3 ) 1 / 2 , wher e for matrix V , define ∥ V ∥ 2 2 , ∞ = max i ∥ V i : ∥ 2 and V i : is the i -th r ow ve ctor of V . Assumption 5 is common in the matrix/tensor completion literature ( Candes & Plan 2010 , 14 Candes & Rec ht 2009 ). Incoherent matrices represent a class where eac h en try con tains a comparable amoun t of information, making it feasible to complete the matrix from a small num b er of observ ations. In practice, this assumption states that information ab out outcomes is reasonably distributed across geographic units, exp osures, and diseases. 4.2 Asymptotic Bounds In this section, w e establish the asymptotic consistency of the estimated parameter tensor b Y obtained by the prop osed spatial tensor completion metho d. W e present a simplified version of our main result here; the complete theorem statement with full tec hnical conditions is pro vided in Supplemen t A.2 . Theorem 1 (F rob enius Error Bound - Simplified V ersion) Under A ssumptions ( 1 )– ( 5 ), with pr ob ability at le ast 1 − 4 ( N + L + O ) − 2 , our estimator satisfies: ∥ ∆ ∥ 2 F N O ≤ max ( A 1 , A 2 , A 3 , A 4 , A 5 ) , wher e A 1 = µ 2 1 p min log( N + L + O ) N , A 2 = C ( r 1 ∧ r 2 ) 1 / 2 r 3 / 4 3 ( r 1 r 2 ) 1 / 4 p max p 2 min σ √ r 1 k − τ 2 s L log( N + L + O ) O , A 3 = C ′ ( r 1 ∧ k ) r 2 r 3 max { r 1 ∧ k , r 2 , r 3 } + µ 2 1 ( r 1 ∧ r 2 ) r 3 / 2 3 µ 3 / 2 0 ( r 1 r 2 ) 1 / 2 p 2 min ! p 2 max p 4 min σ 2 log( N + L + O ) O , A 4 = C ′′ p max p 3 min σ 2 √ r 1 k − τ 1   1 + C 5 s log( N + L + O ) N O   , A 5 = C ′′′ ( r 1 ∧ k ) r 2 r 3 ( r 1 ∨ r 2 ∨ r 3 ) 2 max { r 1 ∧ k , r 2 , r 3 } 1 p 2 min ( N ∨ O ) log ( N + L + O ) N O , wher e a ∧ b , a ∨ b denote the minimum and maximum of a and b r esp e ctively. 15 W e pro vide detailed in terpretations of eac h error comp onent in Supplemen t A.2 . 4.3 A v erage T reatmen t Effect After tensor completion, we estimate the a verage treatment effect θ ∗ ℓo = E [ N − 1 P N i =1 ( Y iℓo − Y iLo )] comparing exp osure ℓ to reference exp osure L for outcome o . Given the imputed p oten tial outcomes, the outcome imputation (OI) estimator is giv en by: b θ (OI) ℓo = 1 N N X i =1 ( b Y iℓo − b Y iLo ) . By Cauc h y-Sch warz inequality and Theorem 1 , this estimator satisfies | b θ (OI) ℓo − θ ∗ ℓo | ≤ ∥ b Y − Y ∗ ∥ F / √ N . Ho wev er, the OI estimator do es not ac hiev e √ N -consistency because tensor completion conv erges at a slow er rate than √ N determined by the terms in Theorem 1 . T o ac hiev e √ N -consistency and enable v alid asymptotic inference, we construct the augmen ted in v erse probabilit y w eigh ting (AIPW) estimator ( Bang & Robins 2005 ): b θ (AIPW) ℓo = 1 N N X i =1 " b Y iℓo + 1 ( A i = ℓ ) b π ( ℓ | Z i , b S i )  Y ( A i ) io − b Y iℓo  ! − b Y iLo + 1 ( A i = L ) b π ( L | Z i , b S i )  Y ( A i ) io − b Y iLo  !# . The AIPW estimator enjoys rate double robustness: it achiev es √ N -consistency when the outcome model and the prop ensit y score mo del are consistently estimated, with the pro duct of their estimation errors b eing of smaller order than N − 1 / 2 , which p ermits the use of flexible mo deling approaches ( Chernozh uko v et al. 2018 ). In our setting where b oth mo dels share the spatial comp onen t, double robustness applies conditionally on adequate spatial appro ximation b S and dep ends on which term in Theorem 1 dominates. Detailed conditions for √ N -consistency are provided in Supplement A.3 . V ariance is estimated via the empirical influence function, and 100(1 − α )% confidence in terv als are constructed 16 as b θ (AIPW) ℓo ± z α/ 2 q b V ℓo / N , where b V ℓo = N − 1 P N i =1 c IF 2 i with details in Supplemen t A.3 . 4.4 Sim ulation Summary W e conducted a simulation study to ev aluate finite-sample p erformance; full details are pro vided in Supplement D. Using a 20 × 20 spatial grid with tw o binary exp osures and ten outcomes, we compared our prop osed spatial tensor metho d against alternativ es includ- ing non-spatial tensor completion, regression, and kriging-based metho d. A cross settings v arying in tensor complexity , strength of spatial confounding, outcome noise, and exp osure o verlap, our metho d achiev ed the low est or near-lo west mean squared error while main- taining appropriate co verage for confidence interv als. These results confirm that join tly mo deling exp osures and outcomes within a low-rank tensor framew ork, combined with sp ectral adjustment for spatial confounding, yields efficiency gains ov er outcome-b y-outcome approac hes. 5 PF AS Health Effect Analysis 5.1 Study Design W e apply our metho d to examine the effects of PF AS exp osures on adult c hronic disease outcomes, fo cusing on PF OA and PFOS. Outcomes include disease rates of 13 common c hronic conditions, eac h prepro cessed by shifting to ensure p ositivit y , applying a logarithmic transformation, and scaling the entire tensor to ha v e zero mean and unit v ariance. The same transformations are applied to all comp eting metho ds. Our study pro vides broader p opulation co v erage than most prior PF AS health in v estigations, encompassing 5,495 public w ater systems serving appro ximately 24.8 million p eople across 47 states. Imp ortan tly , the exp osure levels w e analyze are relatively lo w (mean concentrations 17 of 0.00373 ng/L for PF OA and 0.00498 ng/L for PF OS among exp osed systems), reflecting t ypical con temp orary U.S. drinking-water conditions rather than highly contaminated sites. This con trasts with earlier studies fo cused on heavily p olluted regions with substan tially higher exp osures. F or example, Biggeri et al. ( 2024 ) examined a p opulation exp osed to mark edly elev ated drinking-w ater PF AS concentrations in the V eneto region of Italy , while Barry et al. ( 2013 ) inv estigated the Mid-Ohio V alley cohort using directly measured serum PF AS concen trations, whic h pro vide more precise exp osure quantification than en vironmen tal w ater measurements. Our contribution lies in conducting a large-scale assessment of low- lev el, p opulation-represen tativ e drinking-water exp osure. Ho w ever, this approach necessarily relies on environmen tal monitoring data rather than individual biomark er measurements, whic h limits our ability to capture v ariation in individual absorption, metab olism, and cum ulativ e exposure from non-water sources. Figure 2 depicts the conceptual structure of the study . Our primary fo cus is PF AS exp osure through drinking water, as captured b y ‘UCMR5’, which is affected by both observ able so cio- economic and urbanization factors as well as unmeasured spatial confounders. These factors ma y directly influence health outcomes, in tro ducing potential confounding. Although our analysis fo cuses on w aterb orne PF AS, the diagram also marks other exp osure routes–such as o ccupational sources, PF AS in air, and additional p ollutan ts–that are unobserv ed and appro ximated through their assumed spatial structure. The diagram situates our target path wa y within this broader causal context and underscores the role of b oth measured and unmeasured determinan ts in shaping asso ciations b etw een PF AS exp osure and c hronic disease outcomes. Our final data tensor had dimensions (5 , 495 , 4 , 13) , corresp onding to N = 5 , 495 geographic units, L = 4 exp osure combinations (binary indicators for PF OA and PF OS), and O = 13 outcomes. Spatial distance was defined using regional cen troids, with adjacency constructed 18 Figure 2: Conceptual diagram of PF AS health effect analysis. via a k = 4 nearest-neigh b ors algorithm symmetrized by union of directed edges. The a v erage n um b er of neigh b ors p er region is therefore sligh tly greater than 4. In Supplemen t C.2 visualizes the PF AS distribution and disease distribution. In b oth the prop ensity score mo del for c hemical exp osure and the outcome tensor mo del, w e incorp orate cov ariates related to industry employmen t (p ercen tage of p opulation emplo yed in agriculture, construction, man ufacturing, wholesale, retail, transp ortation, information, finance, professional, educat ion, arts, other services, public administration) as well as a range of so cio-economic indicators, including: median home v alue, the p ercentage of Hispanic/Latino residen ts, the p ercentage of foreign-b orn residents, median income, the p ercen tage of individuals aged 19 to 64 with health insurance, and the percentage of residen ts holding a bac helor’s degree or higher. Population densit y is included as a pro xy for urbanization. The T uck er ranks are tuned via 5-fold cross-v alidation ov er the grids r 1 ∈ { 1 , . . . , 10 } , r 2 ∈ { 1 , . . . , 3 } , and r 3 ∈ { 1 , . . . , 5 } . T o impro ve the robustness of missingness imputation in the tensor, w e implement a 5-fold cross-fitting pro cedure. 19 5.2 Results 5.2.1 Comp eting Metho ds W e compare five estimation approaches, eac h com bining an outcome mo del with a prop en- sit y score mo del under the augmented inv erse probabilit y w eighting (AIPW) framework as describ ed in Section 4.3 . Spatial-T ensor : Our prop osed metho d using spatial tensor com- pletion for outcomes and eigen vector-adjusted multinomial logit for prop ensities. T ensor : T ensor completion without spatial adjustment for outcomes, with standard multinomial logit for prop ensities. Spatial-PS : Standard regression fit separately p er outcome, with eigen v ector-adjusted m ultinomial logit for prop ensities ( Davis et al. 2019 ). Standard Re- gression : Separate regressions p er outcome with standard m ultinomial logit for prop ensities. A ugmen ted Regression : Separate regressions p er outcome, with standard multinomial logit for propensities. 5.2.2 F ull F actorial Exp osure Effects Figure 3 presents estimated o dds ratios across the full factorial exp osure structure, comparing our prop osed spatial tensor metho d with four alternativ e approaches. Regression-based metho ds pro duced inflated and often implausible asso ciations, with o dds ratios exceeding 3.0 for sev eral outcomes. In con trast, tensor-based metho ds pro vided more conserv ativ e and in ternally consistent estimates, with odds ratios typically b elo w 1.05. The stark differences b et ween metho ds reflect distinct approaches to handling missing p oten tial outcomes and confounding. Regression-based estimate eac h outcome indep enden tly , neither b orro wing strength across the 13 disease endp oin ts nor imp osing structure on the coun terfactual outcome surface. T ensor-based approaches, by contrast, exploit lo w-rank structure to p o ol information across outcomes and exp osures sim ultaneously , and regularize estimates through the assumption that a small num b er of laten t factors explain most 20 outcome v ariation, leading to impro v ed efficiency . The spatial tensor metho d further incorp orates graph Laplacian eigenv ectors to approximate unmeasured spatial confounders, remo ving spurious geographic asso ciations while preserving efficiency (substantially narro w er confidence interv als) gains from joint mo deling. The non-spatial tensor iden tified significant adv erse asso ciations for PF OS with four outcomes and one for PFO A, while com bined PF O A+PFOS was asso ciated with elev ated odds for nearly all outcomes. Interestingly , PFO A app eared protectiv e for asthma, whereas PFOS had the opp osite effect. After adjusting for unmeasured spatial confounders via 10 selected eigen vect ors, the spatial tensor attenuated most asso ciations: only tw o adverse health effects for PF OS (tooth loss, ob esity) and one protectiv e effect for PF O A+PF OS remained significan t (COPD). Sensitivit y analyses in Supplemen t C.2 demonstrate that including additional eigen vec- tors progressiv ely attenuates effect estimates, indicating that unadjusted mo dels conflate geographic confounding with causal effects. Cross-v alidation confirms the spatial tensor ac hiev es sup erior out-of-sample prediction, v alidating improv ements in b oth confounding con trol and predictiv e accuracy . 5.2.3 Marginal Exp osure Effects T o summarize exp osure effects across the factorial structure, w e define marginal effects that a verage ov er the distribution of the other exp osure–quan tifying the exp ected impact of eac h PF AS chemical regardless of whether the other is present or absen t, and relev ant for regulatory p olicy . Figure 4 displays the resulting marginal o dds ratios with 95% confidence in terv als, comparing the spatial-tensor and non-spatial tensor mo dels. Detailed calculations are provided in Supplemen t C.2 . After adjusting for unmeasured spatial confounders via the spatial-tensor approach, PF O A show ed no significant asso ciations with an y outcome 21 except for a mo dest protective effect on asthma (OR = 0.995, 95% CI: 0.991-0.999). In con trast, PFOS retained significan t adv erse asso ciations with hypertension (OR = 1.026, 95% CI: 1.001-1.052), to oth loss (OR = 1.019, 95% CI: 1.005-1.033), asthma (OR = 1.007, 95% CI: 1.003-1.012), and ob esity (OR = 1.023, 95% CI: 1.001-1.046). 5.2.4 Laten t F actor Structure Figure 5 (top panel) visualizes the latent geographic structure captured b y the spatial tensor mo del through the factor matrix U 1 . W e applied K-means clustering to the factor loadings to identify eigh t (the selected rank of U 1 ) distinct spatial patterns. The dominan t pattern (dark green) is concen trated throughout m uch of the Eastern United States, and the second cluster (orange) app ears concentrated in parts of the South-Central region, while smaller clusters capture localized geographic features suc h as the P acific Coast, Upper Midw est. Figure 5 (b ottom panel) displa ys the latent disease structure enco ded in factor matrix U 3 , whic h captures shared patterns of outcome v ariation across the 13 chronic diseases. Comp onen t 1 (red) reflects a dominan t shared morbidit y factor, with nearly all diseases load- ing p ositiv ely , representing common risk path wa ys that affects m ultiple c hronic conditions sim ultaneously . Comp onen ts 2 and 3 capture more n uanced con trasts across disease groups. In terpreting these finer structural patterns in terms of sp ecific biological mechanisms remains an imp ortan t direction for future research, as they may reveal how chemical exp osures align with distinct path w ays of chronic disease etiology . These factors demonstrate the tensor mo del’s ability to capture and exploit outcome structure when estimating exp osure effects. 6 Discussion T o estimate the effects of PF AS chemicals on a range of health outcomes, w e prop ose a spatial causal tensor completion framework for multiple binary exp osures and outcomes, 22 PFOA PFOS PFOA+PFOS Spatial−T ensor T ensor 0.98 1.00 1.02 1.04 1.06 0.98 1.00 1.02 1.04 1.06 0.98 1.00 1.02 1.04 1.06 OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH PFOA PFOS PFOA+PFOS Spatial−PS Standard Regression Augmented Regression 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH Odds Ratio Figure 3: Estimated o dds ratios for 13 disease outcomes under PFO A, PFOS, and combined exp osures across differen t mo dels. Black indicates statistically significan t asso ciations; gra y indicates non-significant results. Note that the x-axis scales differ across panels, so effect sizes are not directly comparable b et ween the top and b ottom rows. 23 PFOA (marginal) PFOS (marginal) Spatial−T ensor T ensor 0.98 1.00 1.02 1.04 1.06 0.98 1.00 1.02 1.04 1.06 OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH OBESITY KIDNEY ASTHMA STROKE CANCER CHD COPD MHL TH DIABETES HIGHCHOL TEETH ARTHRITIS BPHIGH Figure 4: Estimated marginal o dds ratios for 13 disease outcomes under PFO A and PFOS exp osures across differen t mo dels. Black indicates statistically significan t asso ciations; gra y indicates non-significant results. in tegrating a lo w-rank tensor structure to p o ol information and sp ectral adjustment to appro ximate unmeasured spatial confounders, em b edded within a pro jected-gradient descen t algorithm. This design enables causal inference in the presence of spatial confounding and p erv asiv e missingness. Theoretically , we establish the F rob enius-error b ounds for the completed outcome tensor, showing estimation error dep ends on spatial appro ximation error, prop ensit y-score estimation error, tensor rank, and others. Extensiv e simulations corrob orate these findings, showing our estimator achiev es efficiency gains against comp eting approac hes. Applied to national PF AS monitoring data linking PFO A/PFOS exp osures with 13 chronic disease outcomes, our metho d yields substantially more conserv ative o dds ratio estimates compared to existing alternativ es. The substantial attenuation following spatial adjust- men t demonstrates the imp ortance of accoun ting for latent geographic confounding. The alternativ es identified m ultiple significan t adv erse asso ciations (b ottom panel of Figure 3 ). After incorp orating tensor structure and laten t spatial confounding, our mo del (top panel) yielded a far sparser set of significan t and attenuated effects. This systematic shift from 24 −0.75 −0.50 −0.25 0.00 0.25 0.50 TEETHLOST BPHIGH OBESITY ARTHRITIS HIGHCHOL DIABETES MHL TH COPD CHD STROKE ASTHMA KIDNEY CANCER Loadings Component 1 Component 2 Component 3 Figure 5: Laten t factors in the spatial tensor model. T op: Clustering of geographic units based on the lo cation factor c U 1 . Bottom: Loadings in outcome factor c U 3 . 25 near-univ ersal asso ciations to selectiv e findings suggests that conv entional approaches ma y substantially inflate effect estimates when spatially-structured confounders remain uncon trolled. Our findings hav e implications for interpreting the broader PF AS epidemiologic literature. Recen t meta-analyses hav e demonstrated fairly consistent adverse asso ciations b etw een PF OA/PF OS exp osure and hypertension ( Xiao et al. 2023 ), blo o d lipid levels ( Liu et al. 2023 ), and particular t yp es of cancer ( Bartell & Vieira 2021 ). Our spatially-adjusted results for PFOS and h yp ertension are consisten t with this meta-analytic evidence. Ho wev er, for man y other health outcomes, the epidemiologic evidence remains inconsistent, as exemplified b y studies of diab etes ( Gui et al. 2023 , Sun et al. 2018 ) and ob esit y ( F rangione et al. 2024 ). Our results suggest that inadequate con trol for geographic confounding ma y hav e inflated effect estimates. The heterogeneity across studies may thus reflect not only differences in exp osure lev els, population c haracteristics, but also different spatial confounding bias. Our analysis is sub ject to several limitations. First, disease outcomes are deriv ed from mo del- based prev alence estimates rather than clinical measurements ( CDC 2023 ), in tro ducing measuremen t error that likely atten uates observ ed asso ciations. Second, exp osure assessment relies on public water system monitoring and do es not capture individual-lev el exp osure v ariation or non-water exp osure routes. Third, the cross-sectional linkage of 2023-2024 w ater monitoring data with 2017-2018 health outcomes limits causal inference, though the long en vironmen tal p ersistence and biological half-liv es of PF AS partially mitigate temp oral misalignmen t. F ourth, while our metho d adjusts for smo oth spatial confounding through Laplacian eigenv ectors, it cannot address all forms of unmeasured confounding, particularly individual-lev el factors uncorrelated with geograph y . Finally , our factorial design treats exp osure as binary indicators rather than mo deling dose-resp onse relationships at v arying concen tration levels. 26 Sev eral promising directions exist for extending this framew ork. First, incorp orating con tinuous exp osures through functional tensor representations would enable mo deling of dose-resp onse relationships rather than binary exp osure indicators. Second, while the tensor framework naturally pro duces unit-level treatment effect estimates, developing v alid inference pro cedures for heterogeneous and individualized causal effects remains an imp ortant metho dological c hallenge. Third, extending the approach to longitudinal settings with time-v arying exp osures and outcomes would accommo date dynamic exposure patterns and enhance its applicabilit y in environmen tal health and related domains. SUPPLEMENT AR Y MA TERIAL Supplemen t A: F ull Theorem, proof and interpretation. Supplemen t B: Detailed algorithm. Supplement C: Additional results with PF AS data analysis. Supplement D: Sim ulation study . References Abadie, A., Agarwal, A., Dwiv edi, R. & Shah, A. (2024), ‘Doubly Robust Inference in Causal Latent F actor Mo dels’, arXiv pr eprint arXiv:2402.11652 . Agarw al, A., Shah, D. & Shen, D. (2020), ‘Syn thetic interv entions’, arXiv pr eprint arXiv:2006.07691 . A uerbac h, J., Slawski, M. & Zhang, S. (2022), ‘T ensor completion for causal inference with m ultiv ariate longitudinal data: A reev aluation of covid-19 mandates’, arXiv pr eprint arXiv:2203.04689 . Bang, H. & Robins, J. M. (2005), ‘Doubly robust estimation in missing data and causal inference mo dels’, Biometrics 61 (4), 962–973. 27 Barry , V., Winquist, A. & Steenland, K. (2013), ‘P erfluoro o ctanoic acid (pfoa) exposures and inciden t cancers among adults living near a chemical plant’, Envir onmental he alth p ersp e ctives 121 (11-12), 1313–1318. Bartell, S. M. & Vieira, V. M. (2021), ‘Critical review on pfoa, kidney cancer, and testicular cancer’, Journal of the A ir & W aste Management A sso ciation 71 (6), 663–679. Biggeri, A., Stoppa, G., F acciolo, L., Fin, G., Mancini, S., Manno, V., Minelli, G., Zamagni, F., Zamboni, M., Catelan, D. et al. (2024), ‘All-cause, cardiov ascular disease and cancer mortalit y in the p opulation of a large italian area con taminated by p erfluoroalkyl and p olyfluoroalkyl substances (1980–2018)’, Envir onmental He alth 23 (1), 42. Bobb, J. F., V aleri, L., Claus Henn, B., Christiani, D. C., W right, R. O., Mazumdar, M., Go dleski, J. J. & Coull, B. A. (2015), ‘Bay esian k ernel machine regression for estimating the health effects of multi-pollutant mixtures’, Biostatistics 16 (3), 493–508. Candes, E. J. & Plan, Y. (2010), ‘Matrix completion with noise’, Pr o c e e dings of the IEEE 98 (6), 925–936. Candes, E. & Rec h t, B. (2009), ‘Exact matrix completion via conv ex optimization’, Com- munic ations of the A CM 55 (6), 111–119. Carrico, C., Gennings, C., Wheeler, D. C. & F actor-Litv ak, P . (2015), ‘Characterization of weigh ted quan tile sum regression for highly correlated data in a risk analysis setting’, Journal of agricultur al, biolo gic al, and envir onmental statistics 20 , 100–120. CDC (2023), ‘Places: Lo cal data for b etter health, census tract data 2023 release’ . Chernozh uk ov, V., Chetv eriko v, D., Demirer, M., Duflo, E., Hansen, C., New ey , W. & Robins, J. (2018), ‘Double/debiased mac hine learning for treatment and structural parameters’ . Da vis, M. L., Neelon, B., Nietert, P . J., Hunt, K. J., Burgette, L. F., Lawson, A. B. & 28 Egede, L. E. (2019), ‘A ddressing geographic confounding through spatial prop ensity scores: A study of racial disparities in diab etes’, Statistic al Metho ds in Me dic al R ese ar ch 28 (3), 734–748. Dup on t, E., W o o d, S. N. & Augustin, N. H. (2022), ‘Spatial+: A no vel approach to spatial confounding’, Biometrics 78 (4), 1279–1290. F enton, S. E., Ducatman, A., Bo obis, A., DeWitt, J. C., Lau, C., Ng, C., Smith, J. S. & Rob erts, S. M. (2021), ‘Per- and p olyfluoroalkyl substance toxicit y and human health review: Curren t state of kno wledge and strategies for informing future researc h’, Envir on- mental toxic olo gy and chemistry 40 (3), 606–630. F rangione, B., Birk, S., Benzouak, T., Ro driguez-Villamizar, L. A., Karim, F., Dugandzic, R. & Villeneuv e, P . J. (2024), ‘Exp osure to p erfluoroalkyl and p olyfluoroalkyl substances and p ediatric ob esity: a systematic review and meta-analysis’, International Journal of Ob esity 48 (2), 131–146. Gaines, L. G. (2023), ‘Historical and curren t usage of p er- and p olyfluoroalkyl substances (pfas): A literature review’, A meric an Journal of Industrial Me dicine 66 (5), 353–378. Gao, C., Chen, H., Zhang, A. R. & Y ang, S. (2025), ‘Causal inference on sequen tial treatmen ts via tensor completion’, arXiv pr eprint arXiv:2511.15866 . Gao, C., Y ang, S. & Kim, J. K. (2023), ‘Soft calibration for selection bias problems under mixed-effects mmdels’, Biometrika 110 (4), 897–911. Gao, C., Zhang, Z. & Y ang, S. (2024), ‘Causal customer c hurn analysis with low-rank tensor blo c k hazard model’, arXiv pr eprint arXiv:2405.11377 . Gilb ert, B., Datta, A., Casey , J. A. & Ogburn, E. L. (2021), ‘A causal inference framew ork for spatial confounding’, arXiv pr eprint arXiv:2112.14946 . 29 Guan, Y., P age, G. L., Reic h, B. J., V entrucci, M. & Y ang, S. (2023), ‘Sp ectral adjustment for spatial confounding’, Biometrika 110 (3), 699–719. Gui, S.-Y., Qiao, J.-C., Xu, K.-X., Li, Z.-L., Chen, Y.-N., W u, K.-J., Jiang, Z.-X. & Hu, C.-Y. (2023), ‘Asso ciation b etw een p er- and p olyfluoroalkyl substances exp osure and risk of diab etes: A systematic review and meta-analysis’, Journal of Exp osur e Scienc e & Envir onmental Epidemiolo gy 33 (1), 40–55. Joub ert, B. R., Kioumourtzoglou, M.-A., Chamberlain, T., Chen, H. Y., Gennings, C., T uryk, M. E., Miranda, M. L., W ebster, T. F., Ensor, K. B., Dunson, D. B. et al. (2022), ‘Po wering researc h through innov ative metho ds for mixtures in epidemiology (prime) program: No vel and expanded statistical metho ds’, International Journal of Envir onmental R ese ar ch and Public He alth 19 (3), 1378. Kang, S., F ranks, A., Audirac, M., Braun, D. & An tonelli, J. (2023), ‘Partial identification and unmeasured confounding with m ultiple treatments and multiple outcomes’, arXiv pr eprint arXiv:2311.12252 . Keil, A. P ., Buckley , J. P ., O’Brien, K. M., F erguson, K. K., Zhao, S. & White, A. J. (2020), ‘A quan tile-based g-computation approach to addressing the effects of exp osure mixtures’, Envir onmental he alth p ersp e ctives 128 (4), 047004. Keller, J. P . & Szpiro, A. A. (2020), ‘Selecting a scale for spatial confounding adjustmen t’, Journal of the R oyal Statistic al So ciety Series A: Statistics in So ciety 183 (3), 1121–1143. K olda, T. G. & Bader, B. W. (2009), ‘T ensor decomp ositions and applications’, SIAM r eview 51 (3), 455–500. Liu, B., Zhu, L., W ang, M. & Sun, Q. (2023), ‘Asso ciations b et w een p er- and p olyflu- oroalkyl substances exp osures and blo o d lipid lev els among adults—a meta-analysis’, Envir onmental he alth p ersp e ctives 131 (5), 056001. 30 Mao, X., W ang, H., W ang, Z. & Y ang, S. (2024), ‘Mixed matrix completion in complex surv ey sampling under heterogeneous missingness’, Journal of Computational and Gr aphic al Statistics 33 (4), 1320–1328. Marques, I., Kneib, T. & Klein, N. (2022), ‘Mitigating spatial confounding by explicitly correlating gaussian random fields’, Envir onmetrics 33 (5), e2727. NASEM et al. (2022), Guidanc e on PF AS exp osur e, testing, and clinic al fol low-up . Osama, M., Zachariah, D. & Schön, T. B. (2019), Inferring heterogeneous causal effects in presence of spatial confounding, in ‘In ternational Conference on Machine Learning’, PMLR, pp. 4942–4950. P apadogeorgou, G., Choirat, C. & Zigler, C. M. (2019), ‘Adjusting for unmeasured spatial confounding with distance adjusted prop ensit y score matc hing’, Biostatistics 20 (2), 256– 272. P oulos, J., Albanese, A., Mercatan ti, A. & Li, F. (2021), ‘Retrosp ectiv e causal inference via matrix completion, with an ev aluation of the effect of europ ean in tegration on cross-b order emplo ymen t’, arXiv pr eprint arXiv:2106.00788 . Prim, S.-N., Guan, Y., Y ang, S., Rapp old, A. G., Hill, K. L., T sai, W.-L., Keeler, C. & Reic h, B. J. (2025), ‘A sp ectral confounder adjustmen t for spatial regression with multiple exp osures and outcomes’, arXiv pr eprint arXiv:2506.09325 . Reic h, B. J., Y ang, S., Guan, Y., Giffin, A. B., Miller, M. J. & Rapp old, A. (2021), ‘A review of spatial causal inference metho ds for environmen tal and epidemiological applications’, International Statistic al R eview 89 (3), 605–634. R ubin, D. B. (1974), ‘Estimating causal effects of treatments in randomized and nonran- domized studies. ’, Journal of e duc ational Psycholo gy 66 (5), 688. 31 Sc hnell, P . M. & P apadogeorgou, G. (2020), ‘Mitigating unobserv ed spatial confounding when estimating the effect of sup ermark et access on cardiov ascular disease deaths’, arXiv pr eprint arXiv:1907.12150 . Sun, Q., Zong, G., V alvi, D., Nielsen, F., Coull, B. & Grandjean, P . (2018), ‘Plasma concen trations of p erfluoroalkyl substances and risk of t yp e 2 diab etes: A prosp ectiv e in v estigation among us women’, Envir onmental he alth p ersp e ctives 126 (3), 037001. Thaden, H. & Kneib, T. (2018), ‘Structural equation mo dels for dealing with spatial confounding’, The A meric an Statistician 72 (3), 239–252. US EP A (2024), ‘Epa’s pfa analytics to ol’ . Accessed: 2024-08-12. URL: https://e cho.ep a.gov/tr ends/pfas-to ols W alker, K., Herman, M., Eb erw ein, K. & W alker, M. K. (2021), ‘P ackage ‘tidycensus”, MIT . Wiec ha, N., Hoppin, J. A. & Reich, B. J. (2025), ‘T w o-stage estimators for spatial con- founding with point-referenced data’, Biometrics 81 (3), ujaf093. W u, J. & F ranks, A. (2025), ‘A latent factor panel approach to spatiotemp oral causal inference’, arXiv pr eprint arXiv:2509.10974 . Xiao, F., An, Z., Lv, J., Sun, X., Sun, H., Liu, Y., Liu, X. & Guo, H. (2023), ‘Asso ciation b et w een p er- and p olyfluoroalkyl substances and risk of hypertension: A systematic review and meta-analysis’, F r ontiers in public he alth 11 , 1173101. Y ang, S. (2018), ‘Propensity score weigh ting for causal inference with clustered data’, Journal of Causal Infer enc e 6 (2), 20170027. Zhen, Y. & W ang, J. (2024), ‘Nonnegativ e tensor completion for dynamic counterfactual prediction on co vid-19 pandemic’, The A nnals of A pplie d Statistics 18 (1), 224–245. 32