---
title: "Application: does SCR/SNCR technology reduce ambient ozone?"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Application: does SCR/SNCR technology reduce ambient ozone?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
data_path <- system.file("extdata", "analysis_dat.xlsx", package = "spaci")
have_data <- nzchar(data_path) && requireNamespace("readxl", quietly = TRUE)
have_geoR <- requireNamespace("geoR", quietly = TRUE)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = have_data)
```

## Background

Selective catalytic and non-catalytic reduction (SCR/SNCR) technologies are
installed at power-generating facilities to cut nitrogen-oxide emissions, a
precursor of ground-level ozone. The causal question is whether installing
SCR/SNCR **reduces ambient ozone**. This is a textbook setting for the two
problems this package targets:

* **Spatial confounding** — meteorology, terrain and socio-economic factors that
  affect both where the technology is installed and local ozone are largely
  unmeasured.
* **Spatial interference** — NO<sub>2</sub> and ozone are transported by wind and
  chemistry, so a facility's treatment affects ozone at *nearby* facilities.

The data (473 facilities, one binary treatment `SnCR`, 18 covariates and
longitude/latitude) were assembled by Papadogeorgou (2016) and previously
analysed by Papadogeorgou et al. (2019) and Pokal et al. (2023). This vignette
reproduces the analysis in Section 2.4 of the report with `spaci`.

> **Data.** A copy of the facility data ships with the package under
> `inst/extdata/analysis_dat.xlsx` and is loaded below via `system.file()`.
> It derives from Papadogeorgou (2016); please cite that source when using it.

## Loading and preparing the data

```{r prep}
library(spaci)
library(readxl)

dat <- as.data.frame(read_excel(data_path))

outcome     <- "mean4maxOzone"
treatment   <- "SnCR"
NO2         <- "totNOxemissions"        # a mediator, excluded from adjustment
coord_names <- c("Fac.Longitude", "Fac.Latitude")
covariates  <- setdiff(names(dat), c(outcome, treatment, NO2, coord_names))

dat <- dat[complete.cases(
  dat[, c(outcome, treatment, coord_names, covariates)]), ]

Y      <- dat[[outcome]]
Z      <- dat[[treatment]]
coords <- as.matrix(dat[, coord_names])
X      <- as.matrix(dat[, covariates])

c(n = nrow(dat), treated = sum(Z == 1), controls = sum(Z == 0))
```

### A note on units

The `mean4maxOzone` column is recorded in **parts per million**, whereas
Table 2.2 of the report is expressed in **parts per billion**. Rescale by 1000
to reproduce the reported numbers:

```{r units}
range(Y)
Y <- Y * 1000          # ppm -> ppb, to match the report's scale
```

## Estimating the effect with every method

```{r fit}
## use geoR (as in the original analysis) when available, else the built-in MLE
engine <- if (have_geoR) "geoR" else "mle"

res <- spatial_ate(Y, Z, X, coords,
                   tau = 0.2, caliper = 0.25,
                   matern_method = engine, seed = 1)
res
```

Alongside the values reported in Table 2.2:

| Method | Report ATT (95% CI) |
|---|---|
| Naive PS  | 1.98 (0.12, 3.84)  |
| DAPS      | 0.54 (−0.76, 1.83) |
| iDAPS     | −0.58 (−1.79, 0.63)|
| recoverU  | −0.15 (−0.84, 0.55)|
| recoverU+ | −0.20 (−1.08, 0.69)|

The **doubly robust** estimators (`recoverU`, `recoverU+`) reproduce the report
essentially exactly, because they do not depend on the random matching order.
`Naive PS` is also close.

```{r forest, fig.alt = "Forest plot of the estimated effect of SCR/SNCR on ozone", fig.width = 7, fig.height = 4}
plot_ate(res, main = "Effect of SCR/SNCR on ozone (ATT)")
```

## Why DAPS and iDAPS need a seed

The matching estimators pair each treated unit with its nearest available
control, processing treated units in **random order**; the point estimate
therefore depends on the RNG state. The original analysis script did not fix a
seed, so its Table 2.2 entries for DAPS and iDAPS are one particular draw.
Averaging over seeds shows the reported values sit inside the sampling
distribution (this loop is illustrative and not run at build time):

```{r seeds, eval = FALSE}
S <- 40
daps_att  <- vapply(1:S, function(s)
  daps(Y, Z, X, coords, caliper = 0.25, seed = s)$att, numeric(1))
idaps_att <- vapply(1:S, function(s)
  idaps(Y, Z, X, coords, tau = 0.2, caliper = 0.25, seed = s)$att, numeric(1))

c(DAPS_mean = mean(daps_att),  DAPS_range  = range(daps_att))
#> DAPS_mean 0.56   range [-0.40, 2.20]   (report: 0.54)
c(iDAPS_mean = mean(idaps_att), iDAPS_range = range(idaps_att))
#> iDAPS_mean 0.06  range [-0.58, 0.59]   (report: -0.58)
```

For a reproducible headline number, fix a seed (as above) or report the
seed-averaged estimate.

## Conclusion

Every method except the naive propensity score returns a confidence interval
that contains zero, and the sign flips from **positive under Naive PS** to
**negative once spatial confounding and interference are adjusted for**. As in
the report, there is **no evidence that SCR/SNCR installation reduces ambient
ozone** after accounting for both phenomena — and, importantly, ignoring them
would have led to the opposite (positive) conclusion.

## References

* Papadogeorgou, G. (2016). Data for *Adjusting for unmeasured spatial
  confounding with distance adjusted propensity score matching*.
* Papadogeorgou, G., Choirat, C. & Zigler, C. (2019). Adjusting for unmeasured
  spatial confounding with distance adjusted propensity score matching.
  *Biostatistics*.
* Pokal, S. et al. (2023). recoverU: a doubly robust estimator using partially
  recovered spatial confounders.
* Ogunsola, I., Johnson, O. & House, T. Unified methods for causal effect
  estimation: mitigating spatial confounding and interference concomitantly.