---
title: "Reproducing the published IPW workflow"
author: "Daijiro Kabata"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Reproducing the published IPW workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = requireNamespace("MatchIt", quietly = TRUE) &&
    requireNamespace("WeightIt", quietly = TRUE) &&
    requireNamespace("survey", quietly = TRUE) &&
    requireNamespace("ranger", quietly = TRUE) &&
    requireNamespace("rpart", quietly = TRUE) &&
    requireNamespace("xgboost", quietly = TRUE)
)
```

## Introduction

The worked example in Kabata, Stuart & Shintani (2024) uses the model-averaged propensity score for *inverse probability weighting* (IPW): the averaged score $\bar e(X)$ is converted into weights, and the treatment effect is estimated as a weighted regression of the outcome on treatment using `survey::svyglm()`, which supplies design-based (sandwich) standard errors. This vignette reproduces that workflow with `psAve`, for both the ATT (the paper's primary estimand) and the ATE (added in the paper's supplement), and shows the two equivalent routes to the weights: the ecosystem route through `WeightIt`, and the direct route through `weights()`.

We use the `lalonde` data from `MatchIt` throughout. See `vignette("psAve", package = "psAve")` for an introduction to `psave()` itself.

```{r, message = FALSE}
library(psAve)
data("lalonde", package = "MatchIt")
```

## ATT weighting

First, fit the model-averaged propensity score with the default settings — criterion `"prog"` (weighted ASMD of the model-averaged prognostic score) and `estimand = "ATT"`:

```{r}
set.seed(1234)
fit <- psave(treat ~ age + educ + race + married + nodegree + re74 + re75,
             data = lalonde, outcome = ~ re78)
fit
```

### The ecosystem route: `psave_weight()`

`psave_weight()` hands the averaged score to `WeightIt::weightit()` as a fixed, externally estimated propensity score (`ps = fit$ps`), reusing the formula and data stored in the `psave` object so no retyping (and no row-misalignment) can occur. The result is an ordinary `weightit` object, so the full `WeightIt`/`cobalt` toolkit applies:

```{r}
w <- psave_weight(fit)   # a genuine weightit object, estimand = "ATT"
w

cobalt::bal.tab(w, distance = data.frame(prog = fit$prog))
```

This is equivalent to the explicit call

```{r, eval = FALSE}
WeightIt::weightit(treat ~ age + educ + race + married + nodegree + re74 + re75,
                   data = lalonde, ps = fit$ps, estimand = "ATT")
```

### The direct route: `weights()`

The `psave` object already contains the IPW weights implied by the averaged score at the fitted estimand, and `weights(fit)` extracts them. For the ATT these are

$$
W_i = \begin{cases} 1 & A_i = 1 \\ \dfrac{\bar e(X_i)}{1 - \bar e(X_i)} & A_i = 0, \end{cases}
$$

the standard "weighting by the odds." They are identical (and unit-tested to be identical) to `WeightIt::get_w_from_ps(fit$ps, treat, estimand = "ATT")`.

### The paper's estimator: `survey::svyglm()`

The published example estimates the ATT by fitting a weighted generalized linear model of the outcome on treatment within a survey design, which provides robust (sandwich) standard errors:

```{r}
library(survey)

des.att <- svydesign(ids = ~1, weights = weights(fit), data = lalonde)
fit.att <- svyglm(re78 ~ treat, design = des.att)
summary(fit.att)
```

The coefficient on `treat` is the IPW-ATT estimate of the effect of the treatment on `re78` among the treated.

Two honest caveats, consistent with the paper. First, `svyglm()`'s sandwich standard errors treat the weights as *fixed*; they do not account for the estimation of the propensity score (including the selection of the mixing weights). This is standard practice in applied IPW analyses and is what the paper reports; if you want uncertainty that reflects the whole pipeline, bootstrap the entire procedure from `psave()` onward. Second, `svyglm()` always returns a `gaussian`-family fit for a continuous outcome; for binary or other outcome types, follow the effect-measure guidance in `WeightIt`'s and `MatchIt`'s estimation vignettes (which use the `marginaleffects` package) rather than reading off coefficients.

## ATE weighting

For the ATE, set `estimand = "ATE"` in `psave()`. This changes two things internally, both taken from the paper's supplement: the weights entering every balance criterion during the grid search, and the weights returned for estimation. The prognostic model is unchanged — it is always fit on untreated units only.

```{r}
set.seed(1234)
fit.ate <- psave(treat ~ age + educ + race + married + nodegree + re74 + re75,
                 data = lalonde, outcome = ~ re78, estimand = "ATE")
fit.ate
```

The ATE weights are the usual inverse-probability weights,

$$
W_i = \begin{cases} \dfrac{1}{\bar e(X_i)} & A_i = 1 \\ \dfrac{1}{1 - \bar e(X_i)} & A_i = 0, \end{cases}
$$

and the estimator is the same `svyglm()` call with the new weights:

```{r}
des.ate <- svydesign(ids = ~1, weights = weights(fit.ate), data = lalonde)
fit.ate.glm <- svyglm(re78 ~ treat, design = des.ate)
summary(fit.ate.glm)
```

Note that the ATE and ATT fits are *different `psave` objects*: because the criterion is computed with estimand-specific weights, the selected mixing weights $\lambda$ (and hence the averaged score itself) generally differ between estimands. Do not reuse an ATT-selected score for an ATE analysis.

The `weightit` route works the same way — `psave_weight(fit.ate)` passes `estimand = "ATE"` through automatically.

## Note on estimand-specific formulas

For reference, the conventions used by the selection criteria (details and derivations in `vignette("method-details", package = "psAve")`):

* **Weights.** ATT: $W = 1$ (treated), $\bar e/(1-\bar e)$ (untreated). ATE: $1/\bar e$ (treated), $1/(1-\bar e)$ (untreated). These enter both the balance criteria over the $\lambda$ grid and the returned `weights`.
* **Standardization.** The weighted absolute standardized mean difference (of covariates for `criterion = "smd"`, of the prognostic score for `criterion = "prog"`) is standardized by the *unweighted standard deviation in the treated group* — for **both** estimands. This follows the paper's supplement exactly; it differs from the pooled-SD convention some software uses for the ATE.
* **KS statistic.** `criterion = "ks"` uses the proper weighted empirical CDF in each arm, with the estimand-specific weights above.
* **Extreme weights.** Candidate propensity scores are clipped to `[0.01, 0.99]` (the `clip` argument) *before* averaging, which bounds all weights; the convex average never leaves the clipped range, so no re-clipping is applied afterward. Inspect `plot(fit.ate, type = "distribution")` and `summary(w)` for remaining extreme weights, which are a greater concern for the ATE than for the ATT.

Only the ATT and ATE are supported; the supplement's criterion formulas were validated for these two estimands, and `psave()` refuses others with an error rather than guessing.

## References

Kabata, D., Stuart, E. A., & Shintani, A. (2024). Prognostic score-based model averaging approach for propensity score estimation. *BMC Medical Research Methodology*, 24, 228. doi:10.1186/s12874-024-02350-y

Lumley, T. (2004). Analysis of complex survey samples. *Journal of Statistical Software*, 9(1), 1–19. doi:10.18637/jss.v009.i08