Estimating causal effects under spatial confounding and interference

The problem

In spatial observational studies two phenomena often occur together:

Spatial confounding (SC): an unmeasured spatial variable U(s) drives both the treatment Z and the outcome Y, so conditioning on the measured covariates X alone does not close the backdoor path.
Spatial interference (SI): a unit’s outcome depends on its neighbours’ treatment, summarised through a neighbourhood-exposure mapping E_i = Σ_j G_ij Z_j.

Methods that address only one of the two are biased when both are present. This package implements two estimators that handle them simultaneously.

library(spaci)

Simulate data

simulate_spatial_causal() draws data in which U(s) is an exponential Gaussian random field, treatment depends on U(s), and the outcome depends on both U(s) and neighbourhood exposure. The true ATT is 2.

sim <- simulate_spatial_causal(n = 250, delta_u = 2.0, tau_exp = 0.1, seed = 1)
str(sim, max.level = 1)
#> List of 7
#>  $ Y       : num [1:250] 4.86 7.78 4.33 5.85 7.15 ...
#>  $ Z       : int [1:250] 0 1 1 1 1 0 0 0 1 0 ...
#>  $ X       : num [1:250, 1:2] 0.1362 0.4072 -0.0697 -0.2477 0.6956 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ coords  : num [1:250, 1:2] 0.266 0.372 0.573 0.908 0.202 ...
#>  $ E       : num [1:250] 0.541 0.603 0.448 0.531 0.803 ...
#>  $ U       : num [1:250] -0.551 0.901 0.25 0.225 1.703 ...
#>  $ true_att: num 2

Estimate the ATT with every method

res <- spatial_ate(sim$Y, sim$Z, sim$X, sim$coords, tau = 0.1, seed = 1)
res
#>      Method      ATT        SE    Lower    Upper
#> 1  Naive PS 2.758583 0.1757337 2.414151 3.103015
#> 2      DAPS 2.611254 0.1795707 2.259302 2.963206
#> 3     iDAPS 2.289499 0.2152935 1.867531 2.711467
#> 4  recoverU 2.664288 0.1607680 2.349189 2.979388
#> 5 recoverU+ 2.251783 0.1285335 1.999862 2.503704

The naive propensity score ignores both SC and SI and is the most biased; iDAPS and recoverU+ adjust for both and sit closest to the true ATT of 2.

Visualising the comparison

plot_ate() draws a forest plot of the estimates and their confidence intervals, with a reference line at zero and (optionally) the true effect.

plot_ate(res, true_att = sim$true_att)

Forest plot of ATT estimates by method

iDAPS in detail

idaps() matches on the composite distance and reports the data-driven weights (π₁, π₂, π₃) on the propensity-score, spatial and interference components.

fit <- idaps(sim$Y, sim$Z, sim$X, sim$coords, tau = 0.1, seed = 1)
fit
#> Spatial causal effect estimate (iDAPS)
#> Estimand: average treatment effect on the treated (ATT)
#> 
#>   ATT   = 2.289 
#>   SE    = 0.2153 
#>   95% CI = [1.868, 2.711]
#> 
#>   Tuning weights: pi1=0.1, pi2=0.3, pi3=0.6 
#>   Matched treated units: 88 (dropped: 50 )
fit$weights
#> pi1 pi2 pi3 
#> 0.1 0.3 0.6

recoverU+ in detail

recoverUplus() recovers the spatial confounder from the residual Matérn field and augments the doubly robust estimator with it and the exposure term. The recovered confounder is returned for inspection.

fp <- recoverUplus(sim$Y, sim$Z, sim$X, sim$coords, tau = 0.1)
fp
#> Spatial causal effect estimate (recoverU+)
#> Estimand: average treatment effect on the treated (ATT)
#> 
#>   ATT   = 2.252 
#>   SE    = 0.1285 
#>   95% CI = [2, 2.504]
head(fp$extras$Uhat)
#> [1]  8.360562e-07  1.021347e-06 -5.595718e-07 -2.916725e-08  5.322667e-07
#> [6] -1.135186e-06

A small simulation study

Averaging over repeated data sets recovers the bias/MSE ordering reported in the paper (recoverU+ and iDAPS beat the naive comparator).

methods <- c("Naive PS", "DAPS", "iDAPS", "recoverU", "recoverU+")
nsim <- 200; true <- 2
store <- matrix(NA, nsim, length(methods), dimnames = list(NULL, methods))
for (s in seq_len(nsim)) {
  d <- simulate_spatial_causal(n = 250, delta_u = 2.0, tau_exp = 0.1)
  store[s, ] <- spatial_ate(d$Y, d$Z, d$X, d$coords, seed = s)$ATT
}
data.frame(Method = methods,
           Bias = round(colMeans(store - true, na.rm = TRUE), 3),
           MSE  = round(colMeans((store - true)^2, na.rm = TRUE), 3))