---
title: "Introduction to hvtiRpropensity"
author: "John Ehrlinger"
date: today
format:
  html:
    toc: true
    toc-depth: 3
    number-sections: true
    embed-resources: true
    code-fold: false
vignette: >
  %\VignetteIndexEntry{Introduction to hvtiRpropensity}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 5
)
library(hvtiRpropensity)
```

## Overview

`hvtiRpropensity` provides propensity score methods for cardiac surgery
comparative-effectiveness research.  It ports the balancing logic from six SAS
templates into a tidy R API whose outputs feed
[hvtiPlotR](https://ehrlinger.github.io/hvtiPlotR/). The six templates ship
with the package, in `system.file("sas-templates", package = "hvtiRpropensity")`.

The package covers the full workflow from raw covariates to sensitivity analysis:

1. **Estimate a balancing score**: `ps_logistic()`, `ps_forest()`,
   `ps_ordinal()`, `ps_nominal()`, `bs_continuous()`, `bs_count()`.
2. **Balance treatment groups**: `ps_match()` (nearest-neighbour matching)
   or `ps_weight()` (IPTW).
3. **Diagnose balance and support**: the SMD tables, group counts and
   effective sample sizes in `$tables`, `ps_stddiff()` and
   `ps_stddiff_perm()`, and `ps_support()` for the common-support rules.
4. **Estimate the effect**: `ps_mw_var()` for a matching-weight difference in
   means, `ps_rmst()` for restricted mean survival.
5. **Check the assumptions**: `sa_overlap()` and `sa_trim_sweep()` for
   positivity and overlap, `sa_rosenbaum()` and `sa_evalue()` for unmeasured
   confounding.

All score, balance, sensitivity, fit, and validation objects share a common
`ps_data` S3 base class with four slots: `$data` carries patient-level outputs, `$meta`
records the declared contract and row accounting, `$tables` carries inference
and diagnostics, and `$models` retains every fitted model used to create the
result.

---

## Sample data

`sample_ps_data()` generates a synthetic cardiac-surgery dataset that mimics a
typical SAS export: two treatment arms (SAVR / TF-TAVR), several continuous and
binary covariates, a fitted propensity score (`prob_t`), and a `match` indicator
initialised to `0`.

```{r}
#| label: sample-data
dta <- sample_ps_data(n = 500, seed = 42)
head(dta)
```

```{r}
#| label: sample-data-summary
table(dta$tavr)
summary(dta$prob_t)
```

The two treatment groups differ on all covariates (younger patients with higher
ejection fraction were more likely to receive TF-TAVR), so the raw SMDs are
non-trivial and the propensity scores have good overlap.

```{r}
#| label: ps-hist
hist(
  dta$prob_t,
  breaks = 30,
  main   = "Propensity score distribution",
  xlab   = "Estimated probability of TF-TAVR",
  col    = "steelblue",
  border = "white"
)
```

---

## Outcome-model bundles

`fit_logistic()` fits binary, proportional-odds ordinal, and generalized-logit
nominal outcome models. The complete outcome levels and the event or reference
level are required. This makes the scientific contrast independent of factor
storage order.

```{r}
#| label: fit-logistic-complete
outcome_model <- fit_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data = dta,
  family = "binary",
  outcome_levels = c(0, 1),
  event_level = 1
)
outcome_model$tables$estimates
outcome_model$tables$fit_status
```

Every bundle records `bundle_version = 1L`, the formula, required predictors,
declared outcome levels, package versions, analysis counts, all model objects,
and four common inference tables:

| table | schema |
|---|---|
| `estimates` | `term`, `estimate`, `std.error`, `statistic`, `df`, `p.value`, `conf.low`, `conf.high`, `odds_ratio` (`rate_ratio` for the log-link count family), `pooled` |
| `covariance` | named coefficient covariance matrix |
| `by_imputation` | `imputation`, `term`, `estimate`, `std.error` |
| `fit_status` | `imputation`, `converged`, `n_input`, `n_analyzed`, `n_excluded` |

For stacked long-form imputations, the same model is fit independently within
each imputation. Predictions are averaged patient by patient, matching the
legacy `PROC SUMMARY mean=` workflow, while coefficient inference and
covariance use Rubin's rules. These answer different questions and are stored
separately. A missing patient, changed outcome or categorical predictor level,
aliased term, failed fit, or unusable covariance stops pooling.

```{r}
#| label: fit-logistic-mi
dta_mi <- do.call(rbind, lapply(1:3, function(imputation) {
  partition <- dta
  partition$age <- partition$age + (imputation - 2) / 10
  partition$imputation <- imputation
  partition
}))

outcome_model_mi <- fit_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data = dta_mi,
  family = "binary",
  id_col = "id",
  imputation_col = "imputation",
  outcome_levels = c(0, 1),
  event_level = 1
)
outcome_model_mi$tables$estimates
outcome_model_mi$tables$fit_status
```

`validate_logistic()` scores a validation cohort from a saved binary bundle
without refitting and without modifying the source object. It returns
equal-frequency calibration groups and a performance row containing observed
and expected event counts, O:E ratio, AUC, and Brier score. The first interface
is intentionally binary only.

```{r}
#| label: validate-logistic
validation <- validate_logistic(outcome_model, dta, groups = 10L)
validation$tables$performance
```

---

## Propensity score estimation

### Binary treatment: `ps_logistic()`

In practice a propensity score is estimated from covariates before matching or
weighting.  `ps_logistic()` fits `stats::glm(..., family = binomial())` and
appends several columns that mirror the original SAS template output:

| R column | SAS equivalent | Description |
|---|---|---|
| `prob_t` | `_p_` / `_propen_` | Propensity score (probability of treatment) |
| `logit_t` | `_logit_` | Log-odds: `log(p / (1 - p))` |
| `mt_wt` | `mt_wt` | Matching weight `min(p, 1-p) / (p·trt + (1-p)·(1-trt))` |
| `quintile` | `quintile` | Rank-based quintile (1–5) |
| `decile` | `decile` | Rank-based decile (1–10) |

```{r}
#| label: ps-logistic
dta_ps <- sample_ps_data(n = 500, seed = 42)
dta_ps$prob_t <- NULL
obj_ps <- ps_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data = dta_ps
)
print(obj_ps)
```

The scored dataset can be passed directly to `ps_match()` or `ps_weight()`.

```{r}
#| label: ps-logistic-cols
head(obj_ps$data[, c("id", "tavr", "prob_t", "logit_t", "mt_wt",
                      "quintile", "decile")])
```

```{r}
#| label: ps-logistic-to-match
m_from_ps <- ps_match(obj_ps$data,
                      score_col = obj_ps$meta$score_col,
                      seed      = 42)
print(m_from_ps)
```

### Multiply-imputed data

The SAS template fits the model separately on each imputed dataset
(`PROC LOGISTIC ... BY _IMPUTATION_`) and averages predicted probabilities
with `PROC SUMMARY mean=`.  Pass a stacked long-format MI data frame and set
`imputation_col` to replicate this.

```{r}
#| label: ps-logistic-mi
dta_base <- sample_ps_data(n = 300, seed = 10)
dta_base$prob_t <- NULL

# Simulate two imputed datasets stacked with an imputation index column
dta_mi <- rbind(cbind(dta_base, imp = 1L),
                cbind(dta_base, imp = 2L))

obj_mi <- ps_logistic(
  tavr ~ age + female + ef + diabetes + hypertension,
  data           = dta_mi,
  imputation_col = "imp",
  id_col         = "id"
)
print(obj_mi)

# Per-patient PS is averaged across imputations, matching PROC SUMMARY mean=
summary(obj_mi$data$prob_t)
```

### Random-forest score: `ps_forest()`

When you would rather not commit to a main-effects logistic model,
`ps_forest()` grows a classification forest with `randomForestSRC` and scores
each patient **out-of-bag**, using only the trees that did not see that
patient. It appends the same `prob_t`, `logit_t`, `mt_wt`, `quintile` and
`decile` columns as `ps_logistic()`, so everything below takes either.
`mt_wt` is the Li and Greene (2013) matching weight. Give a `seed`; the
forest is random.

```{r}
#| label: ps-forest
#| eval: !expr requireNamespace("randomForestSRC", quietly = TRUE)
op <- options(rf.cores = 2L)
obj_rf <- ps_forest(
  tavr ~ age + female + ef + diabetes + hypertension,
  data  = dta_ps,
  ntree = 200,
  seed  = 42
)
options(op)
summary(obj_rf$data$prob_t)
obj_rf$tables$smd
```

### Ordinal treatment: `ps_ordinal()`

For ordered treatment variables (e.g. NYHA functional class I / II / III–IV),
`ps_ordinal()` fits a proportional-odds cumulative logit model via
`MASS::polr()`, equivalent to the default `PROC LOGISTIC` behaviour for ordinal
responses in SAS.

The SAS template decomposes cumulative probabilities into marginal probabilities:
`p1 = col1; p2 = col2 - col1; p3 = 1 - col2`.  `ps_ordinal()` performs this
decomposition internally, appending one `prob_<level>` column per level.

```{r}
#| label: ps-ordinal
#| eval: !expr requireNamespace("MASS", quietly = TRUE)
dta_ord <- sample_ps_data_ordinal(n = 400, seed = 1)
obj_ord  <- ps_ordinal(
  nyha_grp ~ age + female + ef + diabetes,
  data = dta_ord
)
print(obj_ord)

# One marginal probability column per level
head(obj_ord$data[, c("id", "nyha_grp",
                       obj_ord$meta$score_cols)])

# Quintile distribution (ordered by P(highest level) as in SAS)
table(obj_ord$data$quintile)
```

### Nominal treatment: `ps_nominal()`

For unordered multi-level treatments (e.g. repair type: CE / COS / PER / DEV),
`ps_nominal()` fits a generalised logit model via `nnet::multinom()`, equivalent
to SAS `PROC LOGISTIC ... / link=glogit`.  The first factor level is the
reference category (`REF=first` in SAS); use `ref_level` to override.

```{r}
#| label: ps-nominal
#| eval: !expr requireNamespace("nnet", quietly = TRUE)
dta_nom <- sample_ps_data_nominal(n = 400, seed = 2)
obj_nom  <- ps_nominal(
  rtyp ~ age + female + ef + diabetes,
  data = dta_nom
)
print(obj_nom)

# One probability column per repair type (analogous to p_cos, p_per, etc.)
head(obj_nom$data[, c("id", "rtyp", obj_nom$meta$score_cols)])
```

---

## Balancing scores for continuous and count exposures

When the "treatment" is a continuous or count variable (e.g. units of RBC
transfused, nadir haematocrit), a propensity score is undefined.  The SAS
templates instead estimate a **balancing score**: the linear predictor from a
saturated model.  Patients with similar balancing scores have similar expected
exposure values, enabling within-stratum comparisons analogous to PS
stratification.

### Continuous exposure: `bs_continuous()`

Mirrors `tp.rm.continuous.balncing_score.sas` (PROC REG balancing score on
nadir HCT).  The fitted value from `stats::lm()` is used as the score.

```{r}
#| label: bs-continuous
dta <- sample_ps_data(n = 400, seed = 5)
obj_bs <- bs_continuous(
  ef ~ age + female + diabetes + hypertension,
  data   = dta,
  id_col = "id"
)
print(obj_bs)

# Balancing score and decile-based strata
head(obj_bs$data[, c("id", "ef", "bs", "quintile", "decile", "cluster")])

# Binary stratum indicators stra_1 ... stra_10 mirror the SAS %grp macro
table(obj_bs$data$cluster)
```

### Count exposure: `bs_count()`

Mirrors `tp.pm.count.balncing_score.sas` (PROC GENMOD dist=nb link=log).  The
**linear predictor on the log scale** (`xbeta`) is used as the balancing score,
not the fitted count. The log scale is continuous and unbounded, making it
more suitable for rank-based stratification.

```{r}
#| label: bs-count
#| eval: !expr requireNamespace("MASS", quietly = TRUE)
dta_cnt <- sample_ps_data_count(n = 400, seed = 6)

# Negative binomial to account for overdispersion (SAS dist=nb)
obj_cnt <- bs_count(
  rbc_tot ~ female + age + diabetes + hypertension,
  data = dta_cnt,
  dist = "negbin"
)
print(obj_cnt)

head(obj_cnt$data[, c("id", "rbc_tot", "bs", "cluster",
                       "stra_1", "stra_2")])
```

---

## Propensity score matching

### Build the matched object

`ps_match()` performs greedy 1:1 nearest-neighbour matching without
replacement.

Treated patients are matched in a shuffled order, so the pairs depend on
the random-number stream. Give `seed` on every call and the match is
reproducible whatever ran earlier in the session.

```{r}
#| label: match-build
m <- ps_match(dta, seed = 42)
print(m)
```

### Diagnostics

`summary()` prints the available diagnostic tables.

```{r}
#| label: match-summary
summary(m)
```

The `$tables` slot is a named list you can access directly:

```{r}
#| label: match-tables
m$tables$smd_before
m$tables$smd_after
```

### Extract the matched dataset

`$data` is the *full* dataset with a `match` column set to `1` for matched
pairs.  Filter to `match == 1` to obtain the matched subset.

```{r}
#| label: match-data
matched <- m$data[m$data$match == 1L, ]
nrow(matched)
table(matched$tavr)
```

### Caliper matching

A caliper restricts matches to pairs whose propensity scores differ by no more
than the supplied threshold, reducing the number of matched pairs in exchange
for tighter balance.

```{r}
#| label: match-caliper
m_cal <- ps_match(dta, caliper = 0.05, seed = 42)
print(m_cal)
```

### Matching on the logit with a caliper in SD units

A caliper of 0.05 on the probability scale means different things near 0.5
and near the tails. Austin (2011) recommends matching on the **logit** of the
propensity score with a caliper of 0.2 pooled standard deviations of that
logit. Pass the `logit_t` column that `ps_logistic()` appends, set
`scale = "logit"`, and give the caliper as `caliper_sd`. `ps_match()` turns it
into an absolute caliper, `caliper_sd * sqrt((var_t + var_c) / 2)`, and
records that value in `$meta$caliper`.

```{r}
#| label: match-logit
m_logit <- ps_match(obj_ps$data,
                    score_col  = "logit_t",
                    scale      = "logit",
                    caliper_sd = 0.2,
                    seed       = 42)
m_logit$meta[c("scale", "caliper_sd", "caliper", "n_matched", "n_unmatched")]
m_logit$tables$smd_after
```

The object records `scale = "logit"`, and the functions that read it
downstream (`sa_overlap()`, `ps_support()`, `ps_rmst()`) back-transform the
score with `plogis()` before they work on the probability scale.

---

## Inverse-probability-of-treatment weighting

### ATE weights (default)

`ps_weight()` computes IPTW weights and appends them to `$data`.  The default
estimand is the **average treatment effect** (ATE).

```{r}
#| label: weight-ate
w_ate <- ps_weight(dta, estimand = "ATE")
print(w_ate)
```

```{r}
#| label: weight-ate-summary
summary(w_ate)
```

The `iptw` column holds the stabilised weights.

```{r}
#| label: weight-ate-head
head(w_ate$data[, c("id", "tavr", "prob_t", "iptw")])
summary(w_ate$data$iptw)
```

### ATT and ATC weights

For the **average treatment effect on the treated** (ATT), treated patients
receive weight `1` and controls receive `ps / (1 - ps)`.

```{r}
#| label: weight-att
w_att <- ps_weight(dta, estimand = "ATT")
summary(w_att$data$iptw)
```

For the **average treatment effect on the controls** (ATC), the roles are
reversed.

```{r}
#| label: weight-atc
w_atc <- ps_weight(dta, estimand = "ATC")
summary(w_atc$data$iptw)
```

### Weight trimming (winsorisation)

Extreme weights can destabilise estimates.  Supplying `trim` winsorises weights
to the (`trim`, `1 - trim`) quantile range.

```{r}
#| label: weight-trim
w_trim <- ps_weight(dta, estimand = "ATE", trim = 0.05)
summary(w_trim$data$iptw)
```

---

## Balance statistics and effect estimates

The `$tables` of a match or a weight give a quick SMD per covariate. Three
functions port the CCF SAS macros for the fuller analysis, and tests compare
each with the macro's output.

### Standardized differences: `ps_stddiff()`

`ps_stddiff()` reproduces `%stddiff`. You say which formula each variable
gets (Gaussian, non-Gaussian ordinal, binary, categorical) instead of letting
the storage type decide. With `weight_col`, each group's weights are first
scaled to average 1, so stabilised and unstabilised weights give the same
answer. The result is a proportion of a standard deviation, not a percent.

```{r}
#| label: stddiff
sd_mw <- ps_stddiff(obj_ps$data,
                    gaussian   = c("age", "ef"),
                    binary     = c("female", "diabetes", "hypertension"),
                    weight_col = "mt_wt")
sd_mw$tables$stddiff
```

### A permutation reference: `ps_stddiff_perm()`

`ps_stddiff_perm()` reproduces `%stddiffci`. It shuffles the group labels and
reports where the observed difference falls among the shuffled ones. The
percentiles are not a confidence interval; they show what label shuffling
alone produces. A weighted run needs a `reweight` function that recomputes
the weights for each permutation, so this example is unweighted.

```{r}
#| label: stddiff-perm
perm <- ps_stddiff_perm(obj_ps$data,
                        gaussian = c("age", "ef"),
                        binary   = "female",
                        n_perm   = 200,
                        seed     = 1)
perm$tables$stddiff_perm
```

### Matching-weight difference in means: `ps_mw_var()`

`ps_mw_var()` reproduces `%mw_var`: the weighted difference in means between
groups, with a bootstrap standard deviation and percentiles. The weights are
held fixed in the bootstrap, as in the macro, so the SD ignores the
uncertainty in estimating the score. Here `ef` stands in for an outcome.

```{r}
#| label: mw-var
mw <- ps_mw_var(obj_ps$data,
                outcomes   = "ef",
                weight_col = "mt_wt",
                n_rep      = 200,
                seed       = 1)
mw$tables$mw_var[, c("outcome", "mean_0", "mean_1", "estimate", "sd",
                     "p2_5", "p97_5")]
```

### Restricted mean survival: `ps_rmst()`

For a time-to-event outcome, `ps_rmst()` integrates weighted Kaplan-Meier
curves to a horizon `tau` and reports the treated-minus-control difference
with a bootstrap interval. On a `ps_match()` object it analyses the matched
pairs only, unweighted, and resamples whole pairs. It refuses score weights
on a matched set, which would adjust for the same confounders twice. On a
scored or weighted object it runs the weighting schemes you name.

```{r}
#| label: rmst
#| eval: !expr requireNamespace("survival", quietly = TRUE)
# Synthetic follow-up in years, censored at 5
dta_s <- obj_ps$data
dta_s$t <- withr::with_seed(7, stats::rexp(nrow(dta_s), 0.2))
dta_s$e <- as.integer(dta_s$t < 5)
dta_s$t <- pmin(dta_s$t, 5)

m_s <- ps_match(dta_s,
                score_col  = "logit_t",
                scale      = "logit",
                caliper_sd = 0.2,
                covariates = c("age", "female", "ef", "diabetes", "hypertension"),
                seed       = 42)
rm_match <- ps_rmst(m_s, time_col = "t", event_col = "e", tau = 4,
                    n_boot = 100, seed = 1)
rm_match$meta$resampling
rm_match$tables$estimates[, c("estimator", "n", "diff_days",
                              "lo_days", "hi_days")]
```

---

## Sensitivity analysis

A propensity analysis rests on two assumptions the data cannot prove. The
first is positivity: every patient could have received either treatment, so
the two score distributions overlap. `sa_overlap()`, `sa_trim_sweep()` and
`ps_support()` check it. The second is no unmeasured confounding.
`sa_rosenbaum()` and `sa_evalue()` ask how strong an unmeasured confounder
would need to be to change the conclusion. Each `sa_*()` function returns a
`ps_data` object, with its results in `$tables` and `$meta`.

```{r}
#| label: sa-setup
dta  <- sample_ps_data(n = 500, seed = 42)
m    <- ps_match(dta, seed = 42)
w    <- ps_weight(dta, estimand = "ATE")
```

### Overlap and positivity: `sa_overlap()`

Before any analysis, check that the PS distributions overlap sufficiently.
`sa_overlap()` reports the common support region, the fraction of each group
outside it, and near-positivity-violation flags.

```{r}
#| label: sa-overlap
ov <- sa_overlap(m)
ov$meta$overlap_region
ov$tables$outside_overlap
ov$tables$positivity_flags
```

### Comparing support rules: `ps_support()`

There is more than one way to call a patient "off support". `ps_support()`
flags patients under several rules (outside a fixed trimming range, outside
the range both groups occupy, or any logical flag you supply) and measures
how well the rules agree, with Cohen's kappa and the Jaccard overlap of the
patients each rule drops.

```{r}
#| label: ps-support
sup <- ps_support(obj_ps, trim = c(0.1, 0.9))
sup$tables$counts
sup$tables$agreement
```

### Weight-trimming sensitivity sweep: `sa_trim_sweep()`

Extreme IPTW weights inflate variance and can dominate estimates.
`sa_trim_sweep()` sweeps a range of winsorisation thresholds and reports how
ESS and the weight distribution change.  Use this to choose a trimming level
that controls extreme weights without unacceptable loss of sample size.

```{r}
#| label: sa-trim-sweep
sweep <- sa_trim_sweep(w, trim_seq = seq(0, 0.10, by = 0.01))$tables$sweep
head(sweep)

# ESS at no trimming vs. 5% trimming
sweep[sweep$trim %in% c(0, 0.05), c("trim", "ess_control",
                                      "ess_treated", "max_weight")]
```

### E-values: `sa_evalue()`

The E-value (VanderWeele & Ding 2017) is the minimum risk-ratio that an
unmeasured confounder would need with *both* the treatment and the outcome to
fully explain away the observed association.  A large E-value indicates a more
robust finding. For an odds ratio or hazard ratio you must say whether the
outcome is rare (under about 15 percent by the end of follow-up) with `rare`;
the conversion to a risk ratio depends on it.

```{r}
#| label: sa-evalue-rr
# Suppose a matched analysis yields RR = 1.8 (95% CI 1.3 to 2.5)
ev <- sa_evalue(1.8, ci_lo = 1.3, ci_hi = 2.5, type = "RR")
ev$meta$evalue_estimate   # confounder must be this strong to explain the RR
ev$meta$evalue_ci         # must be this strong to explain even the CI lower bound
```

```{r}
#| label: sa-evalue-rd
# Risk difference of 0.08 with baseline outcome probability 0.15
ev_rd <- sa_evalue(0.08, ci_lo = 0.02, ci_hi = 0.14,
                   type = "RD", p0 = 0.15)
ev_rd$meta$evalue_estimate
```

```{r}
#| label: sa-evalue-or
# Odds ratio 1.6 (1.1 to 2.3) for a common outcome
ev_or <- sa_evalue(1.6, ci_lo = 1.1, ci_hi = 2.3, type = "OR", rare = FALSE)
ev_or$tables$evalue
```

### Rosenbaum sensitivity bounds: `sa_rosenbaum()`

For matched analyses, `sa_rosenbaum()` computes worst- and best-case p-value
bounds for the Wilcoxon signed-rank test under a hidden-bias parameter Γ
(Rosenbaum 2002).  Γ = 1 means no hidden bias; the **sensitivity value** is
the largest Γ at which the analysis remains significant.

```{r}
#| label: sa-rosenbaum
res <- sa_rosenbaum(m, outcome_col = "ef", gamma_max = 3, gamma_inc = 0.25)

# Summary: how sensitive is the result?
res$meta$sensitivity_value   # largest Gamma where upper p-value < 0.05
res$meta$n_pairs

# Full bounds table
head(res$tables$bounds)
```

---

## Downstream use with hvtiPlotR

`hvtiPlotR` draws the two standard figures. Both of its constructors return
a data object, and `plot()` turns that into a bare `ggplot2` figure you then
decorate. The chunks below run when `hvtiPlotR` is installed.

### Mirror histogram

`hv_mirror_hist()` reads the `$data` of a match directly. Its defaults are
the column names this package uses (`prob_t`, `tavr`, `match`). It plots on
the probability scale, so give it `prob_t` even for the logit-scale match
from earlier; `ps_logistic()` left both columns in `$data`. The treated group
plots above the axis, controls below, with matched patients shaded.

```{r}
#| label: mirror-hist-match
#| eval: !expr requireNamespace("hvtiPlotR", quietly = TRUE)
mh <- hvtiPlotR::hv_mirror_hist(m_logit$data,
                                score_col = "prob_t",
                                group_col = "tavr",
                                match_col = "match")
plot(mh)
```

For a weighted analysis, pass the weight column and the bars become weighted
counts.

```{r}
#| label: mirror-hist-weight
#| eval: !expr requireNamespace("hvtiPlotR", quietly = TRUE)
mh_wt <- hvtiPlotR::hv_mirror_hist(w_ate$data,
                                   score_col  = "prob_t",
                                   group_col  = "tavr",
                                   weight_col = "iptw")
plot(mh_wt)
```

### Covariate balance plot

`hv_balance()` does **not** read `$tables` as they come. It wants one long
data frame, one row per covariate and comparison, with columns `variable`,
`group` and `std_diff`, and the difference in **percent** (its reference
lines sit at plus and minus 10). The SMD tables here are proportions of a
standard deviation, one table per comparison, so stack them and multiply by
100. Here is the logit-scale match from earlier.

```{r}
#| label: balance-reshape
balance_long <- rbind(
  data.frame(variable = m_logit$tables$smd_before$variable,
             group    = "Before match",
             std_diff = 100 * m_logit$tables$smd_before$smd),
  data.frame(variable = m_logit$tables$smd_after$variable,
             group    = "After match",
             std_diff = 100 * m_logit$tables$smd_after$smd)
)
balance_long
```

For a `ps_weight()` object the two tables are `smd_unweighted` and
`smd_weighted`.

```{r}
#| label: balance-plot
#| eval: !expr requireNamespace("hvtiPlotR", quietly = TRUE)
cb <- hvtiPlotR::hv_balance(balance_long)
plot(cb)
```

---

## Session info

```{r}
#| label: session-info
sessionInfo()
```