--- title: "Getting started with positively" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with positively} %\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%" ) ``` positively provides diagnostics for positivity violations and extrapolation in causal inference. Causal effect estimates rest on the positivity assumption: every unit must have a non-zero probability of receiving each level of the exposure within every stratum of the covariates. Where that probability is zero, or close to it, you face one of two problems: 1) for an exposure-based model, like Inverse Probability Weighting, you get a biased and unstable estimate or 2) for an outcome-based model, like a regression model, an estimate comes from the model extrapolating rather than from observed comparisons. Positivity violations come in two kinds, and they call for different responses. A structural violation is one where treatment is impossible for a subgroup, so no amount of additional data would fill the gap. A practical violation is one where treatment is rare in a region of covariate space, so both exposure levels remain observed overall but the propensity approaches zero or one in the tails. Structural violations usually mean the causal question has to be narrowed to the population where the exposure is possible. Practical violations are a finite-sample problem that shows up as unstable weights and inflated variance, and they can sometimes be managed by trimming, truncation, or a change of estimand. Distinguishing the two is the first task a positivity diagnostic faces, so the methods in positively locate the covariate regions responsible for a violation in addition to measuring its overall severity. ## Setup ```{r setup, message = FALSE} library(positively) library(dplyr) ``` positively includes a small simulated dataset, `pos_violations`, with two known problems planted in it so that a diagnostic can be checked against ground truth. It has a binary `exposure`, two numeric covariates `x1` and `x2`, and a two-level factor `region`. ```{r data} pos_violations ``` The two planted violations are documented in `?pos_violations`. The first is structural: no subject with `region == "b"` and `x2 > 1` is ever exposed, because the true propensity is zero there. The second is practical: exposure depends steeply on `x1`, so the fitted propensity drops below 0.01 in the lower tail of `x1` and rises above 0.99 in the upper tail, while both exposure levels stay observed across the sample as a whole. A good set of diagnostics should recover both, and distinguish the empty subgroup from the near-empty tails. ## Running check_positivity() `check_positivity()` detects the exposure type, or uses the one you declare in `exposure_type`, then checks positivity with the set of diagnostics appropriate to that type, composing the individual `check_*()` functions with their defaults. It collects each diagnostic's output into a container you can read and plot. You give it the data, the exposure column, and the covariates. The covariates are required, with no default, so that outcome columns are never swept in by accident. `vignette("continuous-exposures")` covers when declaring the type is the right call. ```{r check} check <- check_positivity(pos_violations, exposure, c(x1, x2, region)) ``` The message reports that the exposure was detected as binary. For a binary exposure the default set is three diagnostics: effective data points (`edp`), positivity regression trees (`port`), and an extrapolation check (`extrapolation`). Printing the container gives a report. ```{r print} check ``` The header states the exposure, its type, the sample size, and the covariates once for the whole run. Each diagnostic then gets a section, headed by the name you extract it with, reading what that diagnostic found. Diagnostics that found something come first, and the rest follow in the order they were requested. The `port` section reports how many covariate subgroups it examined and how many read as low support, and spells out the rule that produced the count rather than leaving you to assemble it from the parameters. The `extrapolation` section reports the geometric scale it judges nearness by, then how many units have an opposite-exposure unit within one of it and how many fall inside the opposite group's convex hull. The `edp` section reports the range of effective data points across the intervention values. The footer names what to do next. ## Reading the container `summary()` summarizes the container as a whole. It reports one row per statistic per diagnostic: the `value` the diagnostic computed and, for the statistics it pairs with a cut, the `threshold` behind that value. A threshold is stated in the units of the quantity it cuts, which are not always the units of the value beside it. The `port` count of low-support subgroups sits beside `beta`, a prevalence, because the cut was applied to each subgroup's prevalence in turn and the count is what came back. An `NA` means the overview has no one number to show, not that the diagnostic read nothing against anything. ```{r summary-container} summary(check) ``` `sniff_violations()` answers the other question: not what every diagnostic measured, but what any of them found. It returns one row per finding across the whole run, with the number behind each row in the `threshold` column where there is one, so a reader can see that a cut came from an argument they passed. A convex-hull test has no number behind it, so its row states `NA`. A diagnostic that found nothing contributes no rows, so an empty table is an answer rather than a failure. ```{r sniff} sniff_violations(check) ``` `scope` says what a row is about, and rows of different kinds are not comparable. A `subgroup` row names a covariate region you can act on, an `overall` row is a reading taken across the whole run, and a `unit` row is one row of a distribution whose aggregate is the reading. Unit rows are available on request, and they are not included by default because a diagnostic reporting per unit produces far more rows than one reporting per subgroup. ```{r sniff-unit} sniff_violations(check, scope = "unit") ``` To pull a single diagnostic out of the container, name it with `$` or `[[`. Its `tidy()` method returns that diagnostic's own tidy tibble, one row per unit or per subgroup depending on the method, and its `glance()` method returns the wide row of statistics the overview drew on. ```{r tidy-child} tidy(check$port) ``` `autoplot()` draws the whole container as a panel of each diagnostic's default view. The panel needs the patchwork package. Each view keeps the size the panel gives it, so a panel reads well for diagnostics whose default view is compact. PoRT reports a couple of hundred subgroups here, and its default chart is worth a figure of its own, so the panel below holds the other two. ```{r container-plot, eval = requireNamespace("patchwork", quietly = TRUE), message = FALSE, fig.width = 9, fig.height = 3.5, fig.alt = "Boxplots of effective data points at the two candidate exposure values, both spanning most of the range from zero to the largest observed value with points beyond their whiskers, beside histograms of the fraction of nearby opposite-group units faceted by exposure group, each rising from near zero to a peak around four tenths."} compact <- check_positivity( pos_violations, exposure, c(x1, x2, region), diagnostics = c("edp", "extrapolation") ) autoplot(compact) ``` Naming a diagnostic draws that one alone, and arguments go with the name, so `autoplot(check, "port", low_support_only = TRUE)` draws the PoRT chart restricted to its low-support subgroups. Naming a diagnostic composes nothing, so it works without patchwork. Each diagnostic can also be called directly. Calling `check_port()`, `check_edp()`, or `check_extrapolation()` directly returns the same diagnostic object that `check_positivity()` composed, and each object carries its own plot method through `autoplot()`. The sections below call the individual functions so that we can plot and interpret each diagnostic in turn. ## PoRT: locating the subgroups The positivity regression tree of Danelian et al. (2023) grows a shallow tree to find subgroups where one exposure level is rare or absent. It is the most direct diagnostic for the structural-versus-practical distinction, because it reports each problem as a human-readable rule over the covariates. ```{r port} port <- check_port(pos_violations, exposure, c(x1, x2, region)) port ``` The `low_support` column marks subgroups meeting the algorithm's published rule: a prevalence below `beta` or above `1 - beta` among subgroups at least `alpha` of the sample. Filtering to the low-support rows shows what PoRT found. ```{r port-low-support} tidy(port) |> filter(low_support) |> select(description, n, prevalence, low_support) ``` PoRT recovers both planted violations and describes them in the covariates. The rule `x2 >= ... & region = b` has an exposure prevalence of exactly zero: this is the structural subgroup, the cell that is empty by construction. The two `x1` rules, one in the lower tail with prevalence near 0.05 and one in the upper tail with prevalence near 0.96, are the practical near-violation, where exposure is almost but not entirely determined by `x1`. The `autoplot()` method draws the reported subgroups as a bar chart of prevalence, with dashed reference lines at `beta` and `1 - beta` and the low-support subgroups highlighted. The default draws every reported subgroup; `low_support_only = TRUE` restricts the chart to the low-support subgroups. ```{r port-plot, fig.alt = "Horizontal bar chart of exposure prevalence for the low-support PoRT subgroups, drawn near prevalence zero and one against dashed reference lines."} autoplot(port, low_support_only = TRUE) ``` ## Effective data points: how much support surrounds each intervention The effective-data-points diagnostic of Ring and Schomaker (2026) measures, for each unit moved to a candidate exposure value, how much of the observed sample lies nearby in covariate space. It sums a product kernel over the data, so the result runs from zero, meaning no observed support, up to the sample size. For a binary exposure the candidate values are the two levels, 0 and 1. ```{r edp} edp <- check_edp(pos_violations, exposure, c(x1, x2, region)) tidy(edp) ``` Effective data points are a relative measure. They are meaningful when compared across units or across candidate values for a fixed covariate set, not against a universal threshold. Units in a thinly supported region carry a low value at the exposure level they rarely receive, because few observed neighbors surround the intervened-on point. The `autoplot()` method reads effective data points by candidate value. The default view boxes them, with whiskers at the fifth and the ninety-fifth percentile following Ring and Schomaker (2026) and the values outside them drawn as points. `type = "histogram"` opens the same numbers up as a distribution faceted by candidate value; the left tail near zero is the set of units with little observed support under that intervention. ```{r edp-plot, fig.alt = "Histograms of effective data points faceted by candidate exposure value, each with a left tail of low-support observations near zero."} autoplot(edp, type = "histogram") ``` The planted regions make the comparison concrete. The following labels each unit by the planted problem it belongs to: the structural subgroup is `region == "b"` with `x2 > 1`, and the practical tails are the outer five percent of `x1`. The labels are keyed by `.id`, the row identifier carried by the per-unit results of `check_edp()` and `check_extrapolation()`, so they can be joined to both sets of results below. ```{r planted-regions} planted <- pos_violations |> mutate( .id = row_number(), planted = case_when( region == "b" & x2 > 1 ~ "structural", x1 < quantile(x1, 0.05) | x1 > quantile(x1, 0.95) ~ "x1 tail", .default = "elsewhere" ) ) |> select(.id, planted) planted |> count(planted) ``` Comparing effective data points at the exposed level across those labels shows the thin support directly. ```{r edp-planted} tidy(edp) |> filter(value == 1) |> left_join(planted, by = ".id") |> group_by(planted) |> summarize(edp = median(edp)) ``` The median unit elsewhere in the sample has about 73 effective data points under the exposed intervention. The planted regions carry a fraction of that: about 21 in the structural subgroup and about 11 in the `x1` tails, the same regions PoRT reported as rules. ## Extrapolation: distance and the convex hull The extrapolation check follows the logic of King and Zeng (2006): a unit whose covariate profile has no near neighbor in the opposite exposure group can only be compared by extrapolating the model beyond the observed data. It measures nearness with a Gower distance that handles numeric and categorical covariates together, and, when the covariate count is small enough, tests whether each unit falls inside the convex hull of the opposite group. ```{r extrapolation} extrapolation <- check_extrapolation(pos_violations, exposure, c(x1, x2, region)) tidy(extrapolation) ``` `summary()` aggregates the results by exposure group. It reports, for each group, the average distance to the nearest opposite-group unit, the share of units with a neighbor within one geometric variability, and the share inside the opposite group's hull. ```{r extrapolation-summary} summary(extrapolation) ``` The same planted labels trace the extrapolation results to the unit level. ```{r extrapolation-planted} tidy(extrapolation) |> left_join(planted, by = ".id") |> group_by(planted) |> summarize( frac_nearby = median(frac_nearby), in_hull = mean(in_hull) ) ``` The `x1` tails are the extreme case on both measures. The median tail unit has about a tenth of the opposite group nearby, against a third elsewhere, and about a fifth of tail units fall inside the opposite group's hull, against 86 percent elsewhere. The structural subgroup sits between the two: its units have fewer opposite-group neighbors than the rest of the sample, but the hull test runs on the numeric covariates alone, so a subgroup defined partly by `region` separates less sharply here than in PoRT, which searches the categorical covariate directly. 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. ```{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") ``` When the convex-hull test runs, the `type = "hull"` plot shows how many units of each group fall inside the opposite group's hull. Units outside the hull require extrapolation by construction. ```{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") ``` The three diagnostics tell one story about the planted problems, each at its own resolution. PoRT names the structural subgroup and the practical tails as explicit rules. Effective data points shrink to a fraction of the sample median in both regions under the exposed intervention. The extrapolation check places the `x1` tails far from any opposite-group neighbor and largely outside the hull, and it places the structural subgroup below the rest of the sample on nearby support, with the caveat that its hull test reads only the numeric covariates. ## Choosing diagnostics beyond the defaults The default set depends on the exposure type. A categorical exposure drops the extrapolation check; a continuous exposure adds hat values and highest density regions and categorizes the exposure for PoRT. You can also request an explicit subset with the `diagnostics` argument and tune any method through `args`, a named list of per-method option lists. ```{r subset} check_positivity( pos_violations, exposure, c(x1, x2, region), diagnostics = "port", args = list(port = list(alpha = 0.1)) ) ``` Two diagnostics are never run by `check_positivity()`, because each needs something it does not have. `check_eta_bias()` needs an outcome, and `check_density_ratios()` needs user-supplied density ratios. Call both directly. The other vignettes go deeper into each area: - *Positivity for continuous exposures* covers effective data points, hat values, and highest density region non-overlap, and when each one applies. - *Identifying Violations: PoRT and sPoRT* covers subgroup trees, tuning `alpha`, `beta`, and `gamma`, and the longitudinal and censoring variants. - *Estimator-focused diagnostics* covers ETA.Bias and g-computation extrapolation. - *Density ratios and modified treatment policies* covers `check_density_ratios()` with hand-built ratios and with a fitted longitudinal model. ## References Danelian G, Foucher Y, Léger M, Le Borgne F, Chatton A (2023). Identification of in-sample positivity violations using regression trees: the PoRT algorithm. King G, Zeng L (2006). The Dangers of Extreme Counterfactuals. *Political Analysis*, 14(2):131-159. Ring C, Schomaker M (2026). A Diagnostic to Find and Help Combat Stochastic Positivity Issues, with a Focus on Continuous Treatments.