---
title: "Estimator-focused diagnostics"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimator-focused diagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  out.width = "100%"
)
```

Most positivity diagnostics answer one question: does the assumption hold, and
where does support run thin? PoRT names the subgroups, effective data points
measure the support around each intervention, and the extrapolation check
locates units with no near neighbor in the opposite exposure group. The two
diagnostics in this vignette answer a different question. They take the
estimator you intend to use as given, and ask what a positivity violation does
to that estimator's answer. `check_eta_bias()` estimates the bias a positivity
violation inflicts on a chosen estimator under the data's own fitted mechanisms,
and the further bias a truncation rule would add on top of it.
`check_extrapolation()`, read through the lens of g-computation, shows which
units force the outcome model to predict outside the region where it was
trained.

## Setup

```{r setup, message = FALSE}
library(positively)
library(dplyr)
```

Both examples use `pos_violations`, the simulated dataset described in
`?pos_violations` and in the getting-started vignette. It carries two planted
problems: a structural violation, where no subject with `region == "b"` and
`x2 > 1` is ever exposed, and a practical near-violation, where exposure depends
steeply on `x1`, so the fitted propensity approaches zero in the lower tail of
`x1` and one in the upper tail while both exposure levels remain observed. Its
exposure is binary, so the sections on categorical and continuous exposures
simulate small examples of their own.

## ETA.Bias: the bias an estimator incurs

`check_eta_bias()` implements the ETA.Bias parametric-bootstrap diagnostic of
Petersen et al. (2012). The idea is to treat the data's own fitted mechanisms as
a known truth and measure how a candidate estimator behaves under it.

The procedure has three parts. First, it fits a propensity model for the
exposure and an outcome regression on the exposure and covariates. The
g-computation contrast from the outcome model,
averaged over the sample, becomes the target `truth`. This target depends on the
outcome model alone, so it is the same for every estimator and constant across
truncation levels. Second, it draws bootstrap datasets: each resamples the
covariate rows, draws a bootstrap exposure from the untruncated fitted
propensity, and draws a bootstrap outcome from the fitted outcome model. Third,
it refits both models on each bootstrap dataset and applies the candidate
estimator, bounding the fitted propensity into the truncation interval where the
estimator uses it. ETA.Bias is the mean of the bootstrap estimates minus the
target. Because truncation is a property of the estimator rather than of the
data-generating mechanism, the bootstrap exposure is always drawn from the
untruncated propensity, so the diagnostic isolates the bias the estimator plus
its truncation rule introduce.

`check_eta_bias()` needs an outcome, which `pos_violations` does not include, so
`check_positivity()` never runs it and you call it directly. We construct a
continuous outcome that depends on the exposure and the covariates. Its true
average treatment effect is one.

```{r outcome}
set.seed(2025)
eta_data <- pos_violations |>
  mutate(y = exposure + 2 * x1 + 0.5 * x2 + rnorm(n()))
```

We supply the covariates `x1` and `x2`. The steep dependence of the
exposure on `x1` is the source of the extreme propensities that inverse
probability weighting reacts to, and the outcome leans on `x1` more heavily than
on `x2`, so the violation has something to act on. We diagnose the
inverse-probability-weighting estimator, which divides by the fitted propensity
and is therefore the estimator most exposed to a positivity violation.

The bootstrap is random, so we set a seed in the chunk. We use `n_boot = 500`
draws, which takes about half a second with these covariates and the default
model formulas. The number of draws sets the Monte Carlo standard error attached
to each estimate, and that error has to be small relative to the bias you are
trying to read. The bias below is under a tenth of the target, so a few dozen
draws would leave it buried in Monte Carlo noise.

The headline measurement is the untruncated run. With no truncation the
estimator uses the fitted propensities as they stand, so the reported ETA.Bias
is the bias the positivity violation itself inflicts on inverse probability
weighting.

```{r eta-ipw}
set.seed(1)
eta_ipw <- check_eta_bias(
  eta_data,
  exposure,
  y,
  c(x1, x2),
  estimator = "ipw",
  n_boot = 500
)
eta_ipw
```

The message reports that the exposure was detected as binary. The printed output
reports the target truth, the estimator, and the ETA.Bias with its Monte Carlo
standard error. `tidy()` returns a table with one row per estimand term and
truncation level, carrying the columns `term` (the contrast the row estimates),
`truncation_lower` and `truncation_upper` (the bounds the estimator applied),
`bias`, `mc_se`, and `boot_mean` (the mean bootstrap estimate).

```{r eta-ipw-results}
tidy(eta_ipw)
```

The bias is positive and sits about three Monte Carlo standard errors above
zero, so it is a feature of the fitted mechanism rather than bootstrap noise.
The mechanism is the one the dataset was built to carry. Exposure depends
steeply on `x1`, so the fitted propensity approaches zero in the lower tail of
`x1` and one in the upper tail. Weighting by the inverse of those propensities
represents the two tails unevenly, and because the outcome depends on `x1`, the
weighted contrast tilts away from the target.

The `type = "bootstrap"` plot shows the distribution of the bootstrap estimates
with the target truth marked, one facet per truncation level. There is a single
level here, so there is a single facet. Its spread is the sampling variability
of the estimator under the fitted mechanism, and the offset between its center
and the dashed line is the ETA.Bias. The long left tail is the other half of the
picture: on the occasional draw a propensity lands close enough to zero or one
that one observation's weight dominates the estimate.

```{r eta-bootstrap, fig.alt = "A single-facet histogram of bootstrap estimates for the untruncated inverse-probability-weighting estimator, with a tall peak just to the right of a dashed line at the target truth and a long thin left tail reaching below minus two."}
autoplot(eta_ipw, type = "bootstrap")
```

The three estimators respond differently, and running them on the same settings
makes the contrast concrete. G-computation ignores the propensity altogether.
The augmented, doubly robust estimator uses it and then corrects with the
outcome model.

```{r eta-comparison, message = FALSE}
set.seed(1)
eta_gcomp <- check_eta_bias(
  eta_data,
  exposure,
  y,
  c(x1, x2),
  estimator = "gcomp",
  n_boot = 500
)

set.seed(1)
eta_aipw <- check_eta_bias(
  eta_data,
  exposure,
  y,
  c(x1, x2),
  estimator = "aipw",
  n_boot = 500
)

bind_rows(
  ipw = tidy(eta_ipw),
  gcomp = tidy(eta_gcomp),
  aipw = tidy(eta_aipw),
  .id = "estimator"
) |>
  select(estimator, bias, mc_se)
```

G-computation sits within one Monte Carlo standard error of zero, which is what
it should do: the target truth is itself a g-computation estimate, so the only
distance between the two is bootstrap noise. The augmented estimator also sits
within two standard errors of zero, its outcome-model correction absorbing most
of the weighting bias, and it carries more variance than g-computation does.
Inverse probability weighting is the estimator the violation reaches, which is
why we diagnose it here.

The second use of the diagnostic evaluates truncation as a response to the
violation. Bounding the fitted propensity away from zero and one before the
estimator divides by it is a common reaction to extreme weights, and it has a
cost of its own. A vector of lower bounds passed to `truncation_grid` becomes a
set of truncation levels `c(lower, 1 - lower)`, and every level reuses one
shared set of bootstrap draws and model refits, so differences across the sweep
reflect the truncation bound alone. A lower bound of zero applies no truncation,
so the sweep can carry the untruncated run as its first row.

```{r eta-sweep-run, message = FALSE}
set.seed(1)
sweep_result <- check_eta_bias(
  eta_data,
  exposure,
  y,
  c(x1, x2),
  estimator = "ipw",
  truncation_grid = c(0, 0.01, 0.05, 0.1),
  n_boot = 500
)
tidy(sweep_result)
```

The first row and the rows below it say different things. The first row is the
untruncated level, and its bias belongs to the positivity violation. Every row
below it adds bias that the truncation rule introduces, because the bound
replaces the extreme fitted propensities that carry the tail observations, and
those are the observations the weighting relies on to represent the tails of
`x1`. By the time the lower bound reaches 0.1, the bias the truncation rule has
introduced is an order of magnitude larger than the violation's own. The Monte
Carlo standard error moves the other way, shrinking as the bound tightens,
because bounding the propensity stabilizes the weights. That is the variance
side of the tradeoff: truncation buys stability and pays for it in bias.

The `type = "sweep"` plot shows ETA.Bias against the truncation lower bound with
a band of plus or minus two Monte Carlo standard errors. The dashed line at zero
marks the bias a positivity-robust estimator would achieve.

```{r eta-sweep, fig.alt = "ETA.Bias plotted against the truncation lower bound for the inverse-probability-weighting estimator, rising from just above zero at no truncation to about nine tenths at a lower bound of one tenth, with a two Monte Carlo standard error band that narrows as the bound tightens."}
autoplot(sweep_result, type = "sweep")
```

### A categorical exposure, term by term

A binary exposure has one contrast to report. A categorical exposure has one per
non-reference level, and a positivity violation seldom reaches them equally. The
example below plants three doses whose top level is nearly unreachable in the
lower tail of `x1`, while the middle level is well supported everywhere.

```{r dose-data}
set.seed(2025)
n <- 400
x1 <- rnorm(n)
x2 <- rnorm(n)

# The odds of the top dose climb steeply in x1, so subjects in the lower tail
# of x1 almost never receive it.
odds <- cbind(1, exp(0.5 * x1), exp(-1.5 + 3 * x1))
probs <- odds / rowSums(odds)
dose <- factor(
  apply(probs, 1, function(p) sample(c("none", "low", "high"), 1, prob = p)),
  levels = c("none", "low", "high")
)

dose_data <- tibble(
  dose = dose,
  y = as.numeric(dose == "low") +
    2 * as.numeric(dose == "high") +
    2 * x1 + 0.5 * x2 + rnorm(n),
  x1 = x1,
  x2 = x2
)
count(dose_data, dose)
```

Past two levels the treatment mechanism is a multinomial logit, fitted with the
nnet package and refit on every bootstrap draw. The chunks in this section are
gated on nnet being installed, so the vignette builds either way.

```{r dose-ipw, eval = requireNamespace("nnet", quietly = TRUE)}
set.seed(1)
dose_ipw <- check_eta_bias(
  dose_data,
  dose,
  y,
  c(x1, x2),
  estimator = "ipw",
  n_boot = 200
)
tidy(dose_ipw)
```

There is one row per non-reference level, and the two rows read differently.
The `low - none` contrast is within one Monte Carlo standard error of zero: both
levels are observed across the range of `x1`, so weighting has support to work
with. The `high - none` contrast is biased upward by about five Monte Carlo
standard errors, which is the planted violation showing up in the only estimand
it touches. Reporting a single number for the whole exposure would have averaged
the two together and understated the problem.

`truncation_upper` is missing on both rows because truncation past two levels
applies no upper bound. The paragraph on truncation below unpacks that.

G-computation, which never divides by the fitted probabilities, sits within two
Monte Carlo standard errors of zero on both terms.

```{r dose-gcomp, eval = requireNamespace("nnet", quietly = TRUE), message = FALSE}
set.seed(1)
dose_gcomp <- check_eta_bias(
  dose_data,
  dose,
  y,
  c(x1, x2),
  estimator = "gcomp",
  n_boot = 200
)

bind_rows(
  ipw = tidy(dose_ipw),
  gcomp = tidy(dose_gcomp),
  .id = "estimator"
) |>
  select(estimator, term, bias, mc_se)
```

The `type = "bootstrap"` plot gives each term its own row of facets, marked with
the truth that term is aimed at. The two distributions sit differently against
their dashed lines, which is the same reading the table gives.

```{r dose-bootstrap, eval = requireNamespace("nnet", quietly = TRUE), fig.alt = "Histograms of bootstrap estimates for the inverse-probability-weighting estimator, one row per contrast. The low versus none row is a tall peak centered on its dashed truth line, and the high versus none row sits to the right of its dashed line with a long left tail."}
autoplot(dose_ipw, type = "bootstrap")
```

Both terms are contrasts against a reference level, `none` here, which is the
first factor level unless `reference_level` names another. The reference is a
choice of parameterization rather than of estimand: the fitted mechanisms and
the bootstrap draws do not change with it, only which differences are reported.

Truncation changes shape along with the estimand. With two levels the fitted
propensity has one free coordinate, so a grid entry `g` bounds it into
`c(g, 1 - g)`. With more levels each row of fitted probabilities is a point on a
simplex, and a grid entry raises every coordinate below `g` up to `g` and
renormalizes the row to sum to one, matching `propensity::ps_trunc()`. Pinning
every level away from zero already bounds every weight denominator, so there is
no upper bound to apply, which is why `truncation_upper` is missing. A grid
entry has to sit below `1 / k` for an exposure with `k` levels, since `k`
coordinates pinned at `1 / k` already fill the simplex.

```{r dose-sweep, eval = requireNamespace("nnet", quietly = TRUE), message = FALSE}
set.seed(1)
dose_sweep <- check_eta_bias(
  dose_data,
  dose,
  y,
  c(x1, x2),
  estimator = "ipw",
  truncation_grid = c(0, 0.02, 0.05, 0.1),
  n_boot = 200
)
tidy(dose_sweep)
```

The rows walk the sweep once per term. The `high - none` bias more than
quadruples across the grid while the `low - none` bias stays under a tenth, so
the truncation rule charges almost all of its bias to the contrast that was
already the weakest. The Monte Carlo standard error falls on both terms, the
variance side of the same tradeoff the binary sweep showed.

```{r dose-sweep-plot, eval = requireNamespace("nnet", quietly = TRUE), fig.alt = "ETA.Bias against the truncation lower bound with one line per contrast. The high versus none line rises steeply from about three tenths to about one and a third, while the low versus none line stays close to zero across the whole grid."}
autoplot(dose_sweep, type = "sweep")
```

### A continuous exposure and its working model

A continuous exposure has no levels to contrast, so the estimand is the
coefficient set of a marginal structural model. `msm_formula` sets that working
model, and the default is linear in the exposure. The weights differ too:
inverse probability weighting for a continuous exposure divides by a density
rather than by a probability, weighting each observation by the stabilized
density ratio, the marginal density of the exposure over its conditional density
given the covariates.

That ratio is what a sharp treatment mechanism makes explosive. The example
plants one: `x1` explains nine tenths of the variance in the dose, so the
conditional density is far tighter than the marginal and the few observations
whose dose is unusual for their covariates carry enormous weight.

```{r infusion-data}
set.seed(2025)
n <- 400
x1 <- rnorm(n)
x2 <- rnorm(n)

infusion <- tibble(
  x1 = x1,
  x2 = x2,
  dose = 1.5 * x1 + 0.5 * rnorm(n)
) |>
  mutate(y = dose + 2 * x1 + 0.5 * x2 + rnorm(n))
```

```{r infusion-ipw}
set.seed(1)
infusion_ipw <- check_eta_bias(
  infusion,
  dose,
  y,
  c(x1, x2),
  estimator = "ipw",
  n_boot = 200
)
infusion_ipw
```

The term is named `dose`, the one non-intercept coefficient of the default
working model, and the truth is that coefficient under the observed fit: the
counterfactual mean outcome is evaluated over the deciles of the observed dose
and projected onto the working model by least squares. The bias is more than
three quarters of the truth, which is what a mechanism this sharp does to a
density-ratio weight. G-computation, reading the same fitted outcome model, sits
on zero.

```{r infusion-gcomp, message = FALSE}
set.seed(1)
infusion_gcomp <- check_eta_bias(
  infusion,
  dose,
  y,
  c(x1, x2),
  estimator = "gcomp",
  n_boot = 200
)
tidy(infusion_gcomp)
```

Because the working model defines the estimand rather than approximating one,
two analyses of these data that set different `msm_formula` values are aimed at
different targets and will report different biases. That is a property of the
estimand and not an inconsistency: the number reported is the bias of estimating
the working model the analysis actually uses, so it is worth reading beside the
formula that produced it.

Truncation is the third thing that changes with the type. There is no
probability to bound, so a grid entry `g` caps the stabilized weight at its
`1 - g` quantile, read once from the observed fit and held fixed across the
bootstrap draws. The augmented estimator is unavailable here, and
`truncation` and `reference_level` are errors rather than silent no-ops.

```{r infusion-sweep, message = FALSE}
set.seed(1)
infusion_sweep <- check_eta_bias(
  infusion,
  dose,
  y,
  c(x1, x2),
  estimator = "ipw",
  truncation_grid = c(0, 0.01, 0.05, 0.1),
  n_boot = 200
)
tidy(infusion_sweep)
```

`truncation_lower` is missing at every level, since there is no lower cap to
place on a density ratio, and `truncation_upper` holds the cap itself: infinite
where no truncation was applied, then the weight at each requested quantile. The
caps say how concentrated the ratio is, falling from about seven at the 0.99
quantile of the weights to under one at the 0.9 quantile. The bias says that
capping does not rescue the estimator. It grows slightly as the cap tightens
while the Monte Carlo standard error falls by a factor of four, so truncation
buys stability around an answer that is wrong by most of the target either way.

```{r infusion-sweep-plot, fig.alt = "ETA.Bias against the weight cap quantile for a continuous exposure, a single line rising gently from about eight tenths at no truncation to about one and a tenth at a cap quantile of one tenth, well above the dashed line at zero."}
autoplot(infusion_sweep, type = "sweep")
```

ETA.Bias captures the positivity, truncation, and sparsity component of bias and
excludes model-misspecification bias by construction. It reports a diagnostic
for a fitted mechanism and does not correct the estimate.

## Extrapolation as the g-computation view of positivity

G-computation predicts the outcome under each exposure value for every unit and
averages the contrast. It never divides by the propensity, so it does not
produce the extreme weights that a positivity violation creates for inverse
probability weighting. The violation does not disappear, though. It changes
form. Where one exposure group has no member near a unit's covariate profile,
g-computation must predict that unit's outcome under the missing exposure by
extrapolating the model beyond the data it was trained on. Under g-computation
the positivity problem becomes an extrapolation problem, and
`check_extrapolation()` measures it directly.

`check_extrapolation()` implements the extrapolation diagnostics of King and Zeng
(2006). For every observation it measures how well the opposite exposure group
covers that unit's position in covariate space, using Gower distances that
handle numeric and categorical covariates together. We pass all three
covariates, including the factor `region`.

```{r extrapolation}
extrapolation <- check_extrapolation(
  pos_violations,
  exposure,
  c(x1, x2, region)
)
extrapolation
```

The Gower distance between two rows is the mean over covariates of a
per-covariate distance: for a numeric covariate, the absolute difference divided
by the pooled sample range; for a factor, one when the categories differ and
zero when they agree. Distances lie in `[0, 1]`. The geometric variability is
one half the mean of the full distance matrix and sets the reference radius: an
observation counts an opposite-group unit as nearby when their Gower distance is
within one geometric variability, controlled by the `nearby` multiplier.

Three per-unit statistics are reported: `gower_min`, the distance to the nearest
opposite-group unit; `gower_mean`, the mean distance to the opposite group; and
`frac_nearby`, the fraction of the opposite group within the reference radius.

```{r extrapolation-results}
tidy(extrapolation)
```

The units with the lowest `frac_nearby` sit in the covariate tails and in the
structural subgroup, the same regions PoRT reports. The `type = "distribution"`
plot shows the distribution of `frac_nearby` within each exposure group. Most
units carry a third to a half of the opposite group nearby, and the thin left
tail is the handful with almost none, whose g-computation prediction under the
missing exposure rests on extrapolation.

```{r extrapolation-distribution, fig.alt = "Histograms of the fraction of nearby opposite-group units, faceted by exposure group, each rising from a thin left tail near zero to a peak around four tenths."}
autoplot(extrapolation, type = "distribution")
```

The convex-hull test is the sharper statement of the same idea. It asks whether
each observation lies inside the convex hull of the opposite group, solved as a
linear-programming feasibility problem on the numeric covariates. A unit outside
the opposite group's hull cannot be written as a weighted average of
opposite-group units, so any g-computation prediction for it under the missing
exposure is extrapolation by construction. The test runs automatically when
there are at most ten numeric covariates and the lpSolve package is
installed; lpSolve is a suggested dependency, so install it to enable the
hull test. `summary()` aggregates the results by exposure group, reporting the
fraction of each group inside the opposite group's hull.

```{r extrapolation-summary}
summary(extrapolation)
```

The `type = "hull"` plot shows how many units of each group fall inside the
opposite group's hull. The units outside it are the ones that require
extrapolation, and they include the structural subgroup that is never exposed.

```{r extrapolation-hull, fig.alt = "Stacked bars per exposure group showing the share of units inside versus outside the opposite group's convex hull."}
autoplot(extrapolation, type = "hull")
```

Read together, the two diagnostics describe the same positivity violation from
the estimator's side. ETA.Bias puts a number on what the violation costs the
inverse-probability-weighting estimator, and on what truncating the weights
would cost on top of that, and the extrapolation check shows the g-computation
estimator leaning on model predictions for the units in the tails and the
structural subgroup.

## Where to go next

The getting-started vignette, *Getting started with positively*, covers
`check_positivity()` end to end and the diagnostics that answer whether
positivity holds. The vignette *Density ratios and modified treatment policies*
covers `check_density_ratios()`, the diagnostic for stochastic and modified
treatment policies, with hand-built ratios and with a fitted longitudinal model.

## References

King G, Zeng L (2006). The Dangers of Extreme Counterfactuals. *Political
Analysis*, 14(2):131-159.
[doi:10.1093/pan/mpj004](https://doi.org/10.1093/pan/mpj004)

Petersen ML, Porter KE, Gruber S, Wang Y, van der Laan MJ (2012). Diagnosing and
responding to violations in the positivity assumption. *Statistical Methods in
Medical Research*, 21(1):31-54.
[doi:10.1177/0962280210386207](https://doi.org/10.1177/0962280210386207)