Package 'propensity'

Description

propensity provides tools for propensity score analysis in causal inference. Calculate propensity score weights for a variety of causal estimands, handle extreme propensity scores through trimming, truncation, and calibration, and estimate causal effects with inverse probability weighting. The package supports binary, categorical, and continuous exposures.

Weight functions

Calculate propensity score weights for different causal estimands:

• wt_ate(): Average treatment effect (ATE) weights

• wt_att(): Average treatment effect on the treated (ATT) weights

• wt_atu(): Average treatment effect on the untreated (ATU) weights (wt_atc() is an alias)

• wt_atm(): Average treatment effect for the evenly matchable (ATM) weights

• wt_ato(): Average treatment effect for the overlap population (ATO) weights

• wt_entropy(): Entropy weights for the entropy-tilted population

• wt_cens(): Censoring weights

Propensity score modifications

Handle extreme propensity scores before calculating weights:

• ps_trim(): Trim observations with extreme propensity scores

• ps_trunc(): Truncate (winsorize) extreme propensity scores

• ps_calibrate(): Calibrate propensity scores to improve balance

• ps_refit(): Re-estimate the propensity score model after trimming

Estimation

• ipw(): Inverse probability weighted estimator with variance estimation that accounts for propensity score estimation uncertainty

PSW class

The psw() class represents propensity score weights with metadata about the estimand and modifications applied:

• psw(), as_psw(), is_psw(): Create and test propensity score weights

• estimand(): Query the causal estimand

• is_stabilized(): Check if weights are stabilized

Author(s)

Maintainer: Malcolm Barrett [email protected] (ORCID) [copyright holder]

Authors:

• Malcolm Barrett [email protected] (ORCID) [copyright holder]

See Also

• vignette("propensity") for a getting started guide

• The package website for full documentation

Description

augment() works per observation rather than per estimate: the data an ipw() result was produced from, carried through in full, with the propensity score, the weights, the fitted values, and the residuals attached as dot-prefixed columns. Nothing is dropped and no column of the source frame is changed.

The added columns are read from the two models the result holds, so they describe the fit rather than the frame they are attached to. data says which frame to carry them on, which is how broom uses the argument: it names the data the fit was produced from, not new data to predict on. Supplying the modeling data rather than leaving the default is how covariates the outcome formula left out arrive beside the fit's own columns.

The weights are reported once, as .weights, which is a deliberate departure from broom. A broom augment() method carries the model frame's (weights) column through and adds no weights column of its own; here the weights are the central per-observation quantity of the analysis, and the psw() vector reports them with the estimand they target attached, which the model frame's plain numeric copy of the same values does not. The default frame therefore leaves (weights) out. A frame supplied through data is the caller's own and is carried as it arrives, so a weight column they hold stays where they put it, whatever it is named.

A column is an answer per observation only while the fit it is read from and the frame it is carried on describe the same observations, so three things are refused rather than reported. A result whose propensity score model and outcome model were fit to different rows is refused, because the propensity score of one observation would arrive beside the fitted value of another. Two models of the same data most often part this way over missing values, each dropping the rows a variable of its own is missing on. What is compared is the number of observations each model produced an answer for and, when the outcome model's frame names the exposure, the exposure values of the two model frames position by position, reading a factor and the numbers its labels spell as one encoding of one exposure. Exposure values that disagree prove the two models hold different observations; exposure values that agree prove nothing, since two different sets of rows are free to carry the same sequence of values, and two encodings that cannot be read onto each other, such as a factor of "a" and "b" against a recoding of it as 0 and 1, prove nothing either way. The check is one-sided in that direction on purpose: it refuses only what it can prove, so a fit whose models do line up is never refused over the labels the rows of either frame carry or the encoding its exposure is written in. A frame that already holds a column this call would add is refused as well, because the fit's column has nowhere to go: added beside the frame's, the two would be told apart by nothing, and written over it, the frame's data would be gone. Which names are in the way depends on the fit and the frame, .resid being among them only for a frame it would be added to and the propensity columns being the ones the propensity score model produced. Last, a result whose outcome model was fit without weights is refused under the class ipw() refuses the same fact with: .weights has nothing to report for such a fit.

Usage

## S3 method for class 'ipw'
augment(x, data = NULL, ...)

Arguments

Value

A tibble with one row per observation, holding every column of the source frame in its own order followed by:

.propensity

The propensity score, as the propensity score model predicts it: the probability of exposure for a binary exposure, and the conditional mean of the exposure for a continuous one. This is the propensity score the weight functions take, so feeding it back with the exposure returns the weights the outcome model was fit with.

⁠.propensity_<level>⁠

For a categorical exposure, which has a probability for every level rather than a single number, one column per level in place of .propensity, named for the level and ordered as the propensity score model orders its levels.

.weights

The weights the outcome model was fit with, as the psw() vector they were supplied as, so the estimand they record travels with them.

.fitted

The outcome model's fitted values, on the response scale.

.resid

The observed outcome less .fitted. Present only when the source frame holds the outcome the model was fit on, which the default frame always does. A factor or logical outcome is differenced on the 0/1 scale its fitted values are on.

See Also

ipw() for the estimator, tidy() for its estimates, and glance() for a one-row summary of the fit.

Examples

set.seed(123)
n <- 200
x1 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.5 * x1))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + 0.3 * x1))
dat <- data.frame(x1, z, y)

ps_mod <- glm(z ~ x1, data = dat, family = binomial())
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = quasibinomial(), weights = wts)
result <- ipw(ps_mod, outcome_mod)

augment(result)

# The modeling data carries through the covariates the outcome formula left
# out
augment(result, data = dat)

Description

Weights for a continuous exposure are a ratio of densities. Both densities are evaluated on the standardized residual $z_i = (A_i - \mu_i) / \sigma$, where $A_i$ is the exposure, $\mu_i$ the fitted conditional mean, and $\sigma$ the residual spread, so a specification describes a density on a standardized scale rather than on the scale of the exposure itself.

These constructors name the family and record the parameters that identify it. Pass one to the .density argument of wt_ate() or wt_cens(), which also accept the strings "normal", "laplace", and "kernel" for the families that need no parameters, and a bare function of one argument, which is wrapped with dens_fn().

• dens_normal() is the standard normal density, the default.

• dens_laplace() is the standard Laplace density, $\exp(-|z|) / 2$, which puts more mass in the tails than the normal. Its scale is estimated under the Laplace itself, where it is the mean absolute residual, or by the root mean square of the residuals, which is the spread every family without an estimator of its own is read at.

• dens_t() is Student's t density with df degrees of freedom, heavier tailed still, and heavier the smaller df is. Its scale is estimated under the t itself, by maximum likelihood, or by the root mean square of the residuals.

• dens_kernel() is a kernel density estimate of the standardized residuals, fit with stats::density() and interpolated to each observation. It assumes no family at all, at the cost of a density that is not a smooth function of the model's parameters: weights built from it have no closed-form standard error, and ipw() reports none for them at all. Bootstrap the whole fit by hand to put an interval around such an estimate.

• dens_fn() is a density you write yourself.

Usage

dens_normal()

dens_laplace(sigma_method = c("mle", "rms"))

dens_t(df, sigma_method = c("mle", "rms"))

dens_kernel(bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512)

dens_fn(f)

Arguments

Value

An object of class propensity_density: a list with the elements family, the name of the family; params, the parameters that identify it, which is an empty list for a family that takes none; and fn, the function that evaluates the density, which is NULL for dens_kernel(). dens_t() and dens_laplace() carry a fourth element, sigma_method, the name of the estimator their scale is read with.

The spread of the conditional density

Both densities are evaluated on a residual standardized by a spread, and for the conditional density that spread is estimated from the residuals of the propensity score model. A family with no scale estimator of its own takes the root mean square of those residuals, which is sigma_method = "rms" and is the maximum likelihood estimator of the spread of a normal density.

It is not the maximum likelihood estimator of the scale of a heavier-tailed family. The scale of such a family is smaller than its standard deviation, and the root mean square is pulled outward by the large residuals a heavy tail produces, which is what the family was chosen to accommodate. sigma_method = "mle" estimates the scale under the family itself, and both dens_laplace() and dens_t() take it by default, since the spread the density is divided by is the scale parameter of the family reading it. For dens_laplace() it is the mean absolute residual. For dens_t() it is the root of $\sum_i \left[(\nu + 1) r_i^2 / (\nu \sigma^2 + r_i^2)\right] = n$, where $r_i$ is the residual and $\nu$ is df. Each residual enters that sum through a bounded term, so a residual far out in the tail moves the estimate by less than it moves the root mean square. Prefer it when the residuals are heavy tailed, which is the case these families are for; the two estimators answer the same question, and answer it alike, as df grows and the t approaches the normal.

dens_normal() offers no such choice, because the root mean square is already its maximum likelihood estimator. Neither do dens_kernel() and dens_fn(): a kernel is fit to the standardized residuals and divided by the same spread, so the spread cancels up to the grid the estimate is interpolated on, and a density you write yourself names no family an estimator could be derived from.

The choice describes both densities of the ratio. The marginal density that stabilizes the weights is the exposure's own, read at the exposure's mean and at the spread the same estimator gives for it, so the two halves of the ratio are densities of the same width. A scale estimated by maximum likelihood is recorded by density_meta() as sigma = "mle", and ipw() estimates it alongside the propensity score model's coefficients, solving the equation it is the root of as part of its stacked system so that the standard errors account for it. Supplying a .sigma says the spread is a number of your own rather than one estimated from the residuals, so the two cannot be given together.

Examples

dens_normal()

dens_t(df = 4)

dens_t(df = 4, sigma_method = "mle")

dens_laplace(sigma_method = "rms")

dens_kernel(adjust = 1.5)

dens_fn(function(z) stats::dt(z, df = 4))

Description

Weights carry two records describing the exposure they were built for. exposure_type() returns the type of exposure the weight function was given, and density_meta() returns the record left by weights whose value is a ratio of densities.

Both describe the exposure rather than the units, so neither holds a length of its own: both survive subsetting, arithmetic, and anything else that changes the number of weights. Combining two sets of weights that record either of them differently drops the record they disagree on, with a warning of class propensity_metadata_conflict_warning. A result that no longer records an exposure type is not described as continuous, so it drops any density record along with it.

Usage

exposure_type(wt)

density_meta(wt)

## S3 method for class 'propensity_density_meta'
format(x, ...)

## S3 method for class 'propensity_density_meta'
print(x, ...)

Arguments

Details

exposure_type() returns "binary", "categorical", or "continuous". It returns NULL for weights that were built without an exposure to describe, such as those written with psw() directly.

density_meta() returns a record of what the ratio was built from, and NULL for weights that are not a ratio of densities, which is every set of weights for a binary or categorical exposure:

• density, the specification of the conditional density family, as built by dens_normal() and its relatives.

• numerator, what stabilized the weights: "marginal" for the marginal density of the exposure, "integrated" for the conditional density marginalized over the units, "model" for the conditional density a fitted model supplied to stabilize estimates, "score" for a stabilization_score the caller supplied, and "none" for weights that were not stabilized.

• numerator_model, that fitted model itself, and NULL under every other numerator. The model is kept whole rather than the numerator it evaluates to, since ipw() rebuilds that numerator at every value of its parameter vector, which takes the model's design and its coefficients. Two sets of weights stabilized on models that read different terms, or that were fit to different values of the same terms, record different numerators and are combined the way any other disagreeing records are.

• sigma, where the residual spread of the conditional density came from: "pooled" for the pooled residual root mean square, "mle" for a scale estimated under the family that reads it, which dens_t() and dens_laplace() take with sigma_method = "mle", and "supplied" for a .sigma the caller gave. Weights built from a dose model trimmed with ps_trim() record where the trim's spread came from.

• sigma_value, the single spread the weights were read at when it was not estimated by the weight function itself: one the caller supplied, or the spread a dose model trimmed with ps_trim() holds in its record, which is recorded whatever its source. It is NULL for a spread estimated from the residuals, by either estimator, and for one supplied per observation. A spread that is one number is a constant the weights can be rebuilt from, which is what ipw() needs of it; a spread that changes with the observation is not, so the record holds where it came from and nothing more.

Weights that carry a density record print it under their values, as the lines format() renders: the density family, the numerator, and the spread, followed by the formula of the model supplied to stabilize when a model estimated the numerator. Weights that carry none, which is every set of weights for a binary or categorical exposure, print exactly as they always have.

Value

• exposure_type(): a single string, or NULL.

• density_meta(): a list of class propensity_density_meta with the elements density, numerator, sigma, sigma_value, and numerator_model, or NULL.

• print() and format(): the record, invisibly, and a character vector of one line per choice the record describes.

See Also

psw() for the weight vector class and the rest of what it records, and wt_ate() for the functions that write these records.

Examples

set.seed(1)
ps <- runif(20, 0.2, 0.8)
trt <- rbinom(20, 1, ps)

# A binary exposure is weighted by propensity scores rather than by a
# density, so there is no density to record.
exposure_type(wt_ate(ps, trt))
density_meta(wt_ate(ps, trt))

dose <- rnorm(20)
mu <- 0.3 * ps
w <- wt_ate(mu, dose, exposure_type = "continuous")

exposure_type(w)
density_meta(w)

Description

glance() describes an ipw() result rather than its estimates: one row naming the estimand and counting the observations and the residual degrees of freedom of the system the standard errors came from. A fit reporting several effect measures, or several contrasts of a categorical exposure, still returns exactly one row.

Under M-estimation that system is the stacked estimating equations, which hold the propensity score model, the outcome model, and the effect measures at once. Its residual degrees of freedom are the observations it was solved on less the parameters it solves for, so a multinomial propensity score model leaves fewer of them than a binary one does on the same data. Linearization and the robust diagnostic stack nothing and record no parameter count, so the observations are the outcome model's and there is no count to subtract from them.

The columns and their types are the same on every route ipw() takes, so the rows of several results stack into one table.

Usage

## S3 method for class 'ipw'
glance(x, ...)

Arguments

Value

A one-row tibble with the columns:

estimand

The estimand the weights target, such as "ate".

nobs

The number of observations the outcome model was fit on, which is also what the stacked estimating equations are solved on under M-estimation. Reported by stats::nobs().

df.residual

The residual degrees of freedom of the stacked estimating equations, nobs less the number of parameters the system solves for. NA under linearization, which records no parameter count. Reported by stats::df.residual().

See Also

ipw() for the estimator, tidy() for its estimates, and augment() for its per-observation columns.

Examples

set.seed(123)
n <- 200
x1 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.5 * x1))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + 0.3 * x1))
dat <- data.frame(x1, z, y)

ps_mod <- glm(z ~ x1, data = dat, family = binomial())
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = quasibinomial(), weights = wts)
result <- ipw(ps_mod, outcome_mod)

glance(result)

Description

glance() describes a pool_ipw() result rather than its estimates: one row naming the estimand, counting the observations and the results that were pooled, and reporting the complete-data degrees of freedom the small-sample adjustment used. A pooled result reporting several effect measures, or several contrasts of a categorical exposure, still returns exactly one row.

Usage

## S3 method for class 'ipw_pooled'
glance(x, ...)

Arguments

Value

A one-row tibble with the columns:

estimand

The causal estimand every pooled result targeted.

nobs

The smallest number of observations any pooled result was estimated from. The imputed datasets are the same size, so this is that size unless a result was fitted on a subset of one. Reported by stats::nobs().

m

The number of results pooled.

dfcom

The complete-data degrees of freedom the Barnard-Rubin adjustment used.

See Also

pool_ipw() for the pooling, tidy() for the pooled estimates, and the Multiple imputation section of ipw() for the workflow the two belong to.

Examples

set.seed(2024)
n <- 150
x1 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.3 * x1))
y <- rbinom(n, 1, plogis(-0.4 + 0.9 * z + 0.5 * x1))
dat <- data.frame(x1, z, y)
dat$x1[rbinom(n, 1, 0.2) == 1] <- NA

imp <- mice::mice(dat, m = 2, print = FALSE, seed = 1)
fits <- with(imp, {
  ps <- glm(z ~ x1, family = binomial())
  w <- wt_ate(ps)
  om <- glm(y ~ z, family = quasibinomial(), weights = w)
  ipw(ps, om)
})

glance(pool_ipw(fits))

Description

The joint_wt_models method estimates the effects of a joint intervention on two treatments from the pair of fitted treatment models joint_wt_models() records and a weighted outcome model that reads both treatments. Standard errors are computed by M-estimation. Neither se_method = "linearization" nor se_method = "robust" is available here: neither path solves a stacked system, and the cell means and every contrast built from them are parameters of one.

The only supported estimand is "ate", which is what the product weights wt_joint() builds target.

Either treatment may be categorical rather than binary, fit with nnet::multinom(), in which case its block is the multinomial score that fit solves and the surface reports the whole crossing its levels make. Such a component may be stabilized like any other, its numerator estimated in a block of its own; see Joint exposures.

The second treatment may be a dose, whatever the first one is, in which case the surface is the marginal structural model's own coefficients rather than the cells of a crossing, with one level contrast per non-reference level of the first treatment and one dose slope at its reference level. The dose model is read through the same registry the single-treatment route reads, so an stats::lm(), a gaussian stats::glm() at an identity or a log link, and a MASS::rlm() fit with one of the psi functions MASS supplies are stacked, and the dose component's weights may be any density ratio wt_ate() builds. A dose model or a density the stacked system cannot differentiate is refused rather than resampled, since this route has no resampling method; see Joint exposures.

Either component may be stabilized, including on a numerator model the caller fit and passed to wt_ate()'s stabilize. Each component's numerator is estimated in a block of its own, so the standard errors account for it having been fitted; see Joint exposures for what each block holds and for what a numerator may condition on. A numerator the caller computed and passed as wt_ate()'s stabilization_score is read back off the product and held fixed, contributing no block, as it is for a single treatment.

Every score this route stacks is unweighted, so a treatment model fit with prior case weights is refused by the name it was recorded under in joint_wt_models().

This is the second of the two routes to a joint intervention. The other declares the crossing with causalgenerics::joint_exposure() and weights it with one multinomial model over the cells; see Joint exposures. Both report the same surface, so the choice between them is a modeling choice rather than a reporting one. Prefer this route when the two treatments call for different adjustment sets, or when the dependence of the second treatment on the first is what you want to model directly: it weights through the sequential factorization ⁠f(A | L) f(E | A, L)⁠, so each treatment gets its own model and its own covariates. Prefer the declared crossing when one model over the cells is the natural specification.

The factorization is validated when the container is built rather than here: joint_wt_models() refuses a second model that does not condition on the first treatment, and a pair that condition on each other. By the time a container reaches ipw() those questions are settled.

The multinom method estimates causal effects for a categorical exposure from a fitted nnet::multinom() propensity score model and a weighted outcome model. Standard errors are computed by M-estimation. Neither se_method = "linearization" nor se_method = "robust" is available for a categorical exposure, and both are refused with an error of class propensity_method_error: the stacked system has a sandwich for every fit this exposure type accepts, and the diagnostic sandwich describes two fitted cells that a categorical result does not report.

For a K-level exposure, the counterfactual mean at each level is reported first, under the effect label "mean", and then the effects for each non-reference level against the reference (first) factor level. The contrast column names the level each mean row belongs to and the pair of levels each effect row compares, as it does for a binary exposure.

Two fits are refused with an error of class propensity_model_family_error. A nnet::multinom() fit to two levels reports the single probability of a binary exposure rather than one probability per level; fit such an exposure with stats::glm() and family = binomial() and pass it to the glm method instead. A fit to a matrix response records no levels at all, since nnet::multinom() reads a matrix as counts; refit it with the exposure factor on the left-hand side.

The lm method estimates the causal dose-response effect for a continuous exposure from a fitted stats::lm() (or gaussian-family stats::glm()) propensity score model of the exposure and a weighted marginal structural outcome model. A MASS::rlm() fit with one of the psi functions MASS supplies and an mgcv::gam() of a gaussian family reach this method by inheritance and are stacked at their own scores, the additive fit with its smoothing parameters held at the values it reports. A gaussian model is read through its link, which may be identity or log; the other gaussian links, and every class the stacked system has no score for, are refused, as are weights built with a "kernel" density, which have no closed-form standard error. See Continuous propensity score models in ipw() for what each refusal says, what to do about it, and what the additive route conditions on. The only supported estimand is "ate".

Standard errors are computed by M-estimation; neither the linearization method nor the robust diagnostic is available for continuous exposures. The package offers no resampling method, so a fit the stacked system cannot write has no standard error here, and the refusal points at a bootstrap of the whole pipeline written by hand. See Continuous propensity score models in ipw().

Every term of the marginal structural model that reads the exposure must read the exposure and nothing else, and the reported effects are the coefficients those terms contribute, labeled "log(or)" for a logit-link outcome and "log(rr)" for a log-link one. A model with one exposure coefficient keeps the eight-column contract with no contrast column and labels its row "slope" at an identity link; a dose-response curve such as y ~ A + I(A^2) reports one row per coefficient under "coef", gaining a contrast column that names them.

A curve is also the one shape of fit that has a single reading. An exposure entering the outcome model through more than one design column, whether written as y ~ A + I(A^2) or as a basis such as poly(A, 2), splines::ns(A, 3), or splines::bs(A, 3), has no coefficient that is the effect of the exposure, because the dose response has a different slope at every dose. ipw() records the conditional reading for such a fit and reports it once as a message. Asking ipw() itself for effects = "marginal" is refused with an error of class propensity_ipw_effects_error; the result declares the conditional reading as the only one it supports, so causalgenerics::as_marginal(), stats::coef(), stats::vcov(), stats::confint(), generics::tidy(), and as.data.frame() refuse it from the result class instead, with an error of class causalgenerics_unsupported_reading_marginal. Marginalizing the curve over the observed doses is a separate estimand that this package does not compute; use the marginaleffects package on the conditional result: avg_slopes() for slopes, avg_comparisons() for contrasts, and avg_predictions() for causal dose-response functions. A basis term reading a covariate rather than the exposure is a covariate term however many columns it expands to, so it contributes no row and leaves the reading marginal.

A model frame records the term rather than the variables inside it, so the frame of a fit whose exposure enters through poly() or a spline holds no exposure column and .data must supply one. The propensity score design is then rebuilt through the model's terms object, whose predvars attribute records the basis the model was fit with, so the rebuilt columns are the fitted ones.

ipw() is a bring-your-own-model (BYOM) inverse probability weighted estimator for causal inference. You supply a fitted propensity score model and a fitted weighted outcome model; ipw() computes causal effect estimates with standard errors that correctly account for the two-step estimation process.

ipw() supports binary, categorical, and continuous exposures. For a binary or categorical exposure, a binary outcome returns the risk difference, log risk ratio, and log odds ratio, and a continuous outcome returns the difference in means. A continuous exposure is supplied through an stats::lm() or gaussian stats::glm() propensity score model with a weighted marginal structural outcome model whose link is identity, logit, or log, and ipw() reports the single exposure coefficient of that model. A MASS::rlm() fit with one of the psi functions MASS supplies and an mgcv::gam() fit of a gaussian family are read as the equations their own coefficients solve; every other subclass errors.

Usage

## S3 method for class 'joint_wt_models'
ipw(
  wt_mod,
  outcome_mod,
  ...,
  .data = NULL,
  .by = NULL,
  estimand = NULL,
  ps_link = NULL,
  conf_level = 0.95,
  se_method = c("mestimation", "linearization", "robust"),
  effects = c("marginal", "conditional")
)

## S3 method for class 'multinom'
ipw(
  wt_mod,
  outcome_mod,
  ...,
  .data = NULL,
  .by = NULL,
  estimand = NULL,
  ps_link = NULL,
  conf_level = 0.95,
  se_method = c("mestimation", "linearization", "robust"),
  .focal_level = NULL,
  effects = c("marginal", "conditional")
)

## S3 method for class 'lm'
ipw(
  wt_mod,
  outcome_mod,
  ...,
  .data = NULL,
  .by = NULL,
  estimand = NULL,
  ps_link = NULL,
  conf_level = 0.95,
  se_method = c("mestimation", "linearization", "robust"),
  effects = c("marginal", "conditional")
)

## S3 method for class 'glm'
ipw(
  wt_mod,
  outcome_mod,
  ...,
  .data = NULL,
  .by = NULL,
  estimand = NULL,
  ps_link = NULL,
  conf_level = 0.95,
  se_method = c("mestimation", "linearization", "robust"),
  effects = c("marginal", "conditional")
)

Arguments

Value

Methods of ipw() return an S3 object of class ipw. That result class is shared across packages and its components, its print() method, and its as.data.frame() method are documented at causalgenerics::new_ipw(). propensity's methods fill every component; three of them take propensity-specific values:

estimand

One of "ate", "att", "atu", "atm", "ato", or "entropy".

se_method

One of "mestimation", "linearization", or "robust".

fit

The fitted M-estimator object when se_method is "mestimation". Calling stats::vcov() or generics::tidy() on this object exposes the full stacked system of estimating equations, including the propensity score, outcome, and estimand parameters. It is NULL under "linearization" and "robust", neither of which solves a system.

The result answers the standard model accessors, which causalgenerics registers for the shared class: stats::coef() and stats::confint() for the reported effects, stats::vcov() for their covariance, stats::nobs() and stats::df.residual() for the counts describing the fit, and stats::weights() for the psw() vector the outcome model was fit with. Coefficients are named for the effect measure and the contrast together, the contrast being the level a counterfactual mean belongs to or the pair of levels an effect measure compares, and for the subgroup as well where .by reports one row per subgroup. Which surface stats::coef(), stats::vcov(), and stats::confint() report follows the presentation mode the result records, described under effects above.

Under se_method = "mestimation" the stored wt_mod and outcome_mod are the models supplied, carrying their block of the joint sandwich. Calling stats::vcov() on either reports a covariance that accounts for the whole system having been estimated from the same data, rather than the one the model's own fitting routine reports. The models supplied are left as they were fit. Linearization solves no such system, so its component models are stored unchanged.

Workflow

ipw() is designed around a three-step workflow:

1. Fit a propensity score model (e.g., logistic regression of exposure on confounders).

2. Calculate propensity score weights for your estimand (e.g., wt_ate()) and fit a weighted outcome model.

3. Pass both models to ipw() to obtain causal effect estimates with correct standard errors.

You are responsible for specifying and fitting both models. ipw() handles the variance estimation.

Estimands

The available estimands depend on the exposure type:

A continuous exposure supports only "ate"; any other request errors. For a categorical exposure, "att" and "atu" require a focal level, supplied through .focal_level or taken from the weights (see the multinom method).

A binary exposure takes no focal level: ipw() always treats the second sorted exposure level as the exposed group, which is the second level of a factor, TRUE for a logical, and the larger of two numeric values. Because "att" and "atu" weights are mirror images of each other, weights built with wt_att() or wt_atu() on the other level carry values that are correct for the reversed roles, and ipw() rejects them rather than correcting them as though the roles matched. To target the other level, relevel the exposure so that it sorts second and refit both models.

The estimand a weighted outcome model was fit for is not necessarily the one a tool that averages its per-row predictions reports, so see Standardizing model predictions in ps_tilt() before taking such an average of a model fit for anything other than "ate".

Effect measures

A binary or categorical exposure reports the counterfactual mean at each exposure level before the contrasts built from those means. Those rows carry the effect label "mean" and lead the table, one per level in level order with the reference level first. Each is the outcome model's prediction with every unit set to that level, averaged over the sample under the estimand's tilt, so a binomial outcome reports the marginal risk under each level and a continuous one reports the marginal mean. Reading a risk difference beside the two risks it is a difference of is what those rows are for: a difference of 0.05 is a different finding at risks of 0.02 and 0.07 than at 0.50 and 0.55.

For a binary exposure, the reported measures depend on the outcome model. A binary outcome (stats::glm() with family = binomial()) returns three measures:

• rd: Risk difference (marginal risk in exposed minus unexposed)

• log(rr): Log risk ratio

• log(or): Log odds ratio

A continuous outcome (stats::lm() or stats::glm() with family = gaussian()) returns only the difference in means (diff).

For a categorical exposure, the same measures are reported for each non-reference level against the reference (first) factor level, so a K-level exposure produces one block of measures per non-reference level.

The estimates table of either exposure carries a contrast column naming what each row is about: the level a mean row belongs to, and the pair of levels a contrast row compares, written "1 vs 0" or "b vs a" in the levels the models were fit on. A binary exposure names its one pair the same way a categorical one names each of its several, so the two tables are read alike.

For a continuous exposure, ipw() reports the exposure coefficients of the weighted marginal structural outcome model. A model with one exposure coefficient reports one row and no contrast column, and a dose-response curve such as y ~ A + I(A^2) reports one row per coefficient, with the contrast column naming the coefficient each row belongs to.

The effect label names the scale: log(or) for a logit link and log(rr) for a log link, since a coefficient there is a log odds ratio or a log risk ratio per unit of whatever column it multiplies. At an identity link the word depends on what the row claims. A model with one exposure coefficient reports slope, because that coefficient is the slope of the dose response everywhere. A model with several reports coef, because a curve has a different slope at every dose and no one of its coefficients is that slope; the row says which coefficient it is and leaves the dose response to be built downstream from coef() and vcov().

Every term of the marginal structural model that reads the exposure must read the exposure and nothing else. The requirement is on what a term reads rather than on how it is written, so A + I(A^2) and A + sin(A) are admitted and A * x1 is not: a term reading a covariate as well contributes a coefficient that depends on that covariate, and no row could name the effect it stands for. Read the whole coefficient vector from the returned fit object for a model this surface cannot report.

A basis matrix reads the exposure alone, so poly(A, 2), splines::ns(A, 3), and splines::bs(A, 3) are admitted on the same footing. One such term contributes one coefficient per basis column and each of them is a row, reported under coef at an identity link like any other multi-coefficient surface. The contrast column names the coefficient rather than the term, a distinction a basis makes and A + I(A^2) does not: the single term poly(A, 2) reports the contrasts ⁠poly(A, 2)1⁠ and ⁠poly(A, 2)2⁠. A basis reading a covariate is a covariate term and contributes no row, and one reading both, such as poly(A, 2):x1, is refused with every other mixed term.

A marginal structural model built on a basis of the exposure requires .data. A model frame records the term rather than the variables inside it, so the frame of such a fit carries the basis matrix and no column holding the exposure, which is where ipw() reads the exposure from when .data is absent. Omitting it errors, and the message names the exposure and directs to .data. Supplying it does not refit the basis: the outcome design is read from the fit, and a design rebuilt from .data is rebuilt through the fitted terms object, whose predvars attribute records the basis that term was fit with.

Use as.data.frame() with exponentiate = TRUE to obtain risk ratios and odds ratios on their natural scale.

Effect modification

.by names a modifier and asks for the effect within each of its levels alongside the effect for the sample as a whole. The reported rows come in four blocks, in this order:

• the overall rows, unchanged from the fit without .by, the counterfactual means and then the contrasts, under the group label "overall";

• the counterfactual mean at each exposure level within each level of the modifier, under a "var = value" label such as "sex = female";

• one block of contrasts per level of the modifier, under the same labels;

• one block of contrasts per non-reference level against the reference level, which is the modifier's first level, under a "var = value vs var = value" label.

The last block reports no means of its own. Its rows compare the effects in two strata, and a difference of two effects belongs to neither stratum's population, so there is no level whose mean it would be.

A stratum's effect is g-computation restricted to that stratum: the outcome model predicts each unit's outcome at each exposure level, and the predictions are averaged over the units of the stratum, weighted by the estimand's tilt. What a stratum row estimates is therefore the effect in that stratum's own covariate distribution, and the last block is a difference of two such effects. Those two populations differ whenever the covariates are distributed differently across the strata, so a difference between strata is a difference between two conditional effects and not, on its own, evidence that the modifier changes the effect. It is the estimand's population that the tilt fixes, and the modifier's stratum that the indicator fixes, and the row reports the effect where those two meet: an att fit reports the effect among the exposed units of a stratum, not among all of its units. For a categorical exposure the tilt an att fit standardizes over is the focal level's own propensity score, so a stratum row there reports the effect among the units of that stratum who look like the ones assigned the focal level.

A categorical exposure already reports one contrast per non-reference level, and .by crosses those contrasts with the strata: every block holds each contrast on each reported measure, so a K-level exposure with S strata reports the whole-sample block and then S + (S - 1) blocks of (K - 1) contrasts, beside S blocks of K means. A row is named by three things, the measure, the contrast, and the subgroup, and the estimates frame carries contrast and group in that order after effect. Within a contrast block the ordering is the one an ungrouped fit already uses, contrast-major and measure-minor, so a contrast's measures sit together; the blocks themselves run group-major.

The stratum rows and the stratum contrasts report the risk difference and the log risk ratio for a binary outcome, and the difference in means for a continuous one. They report no odds ratio. An odds ratio is noncollapsible, so a difference of two of them moves with the outcome distribution in each stratum whether or not the effect does, and reporting it as effect modification would invite exactly the reading it does not support. The overall rows keep all three measures.

Three requirements apply. The modifier must be a factor or a character vector with no missing values: a missing value names no subgroup, and a continuous variable names no fixed set of them, so cut it into groups first. A level no unit carries is dropped rather than refused, since the fits drop it too and an empty stratum has no mean to standardize. Every one of the modifier's remaining levels must hold every exposure level, since a stratum in which nobody took one of them identifies no contrast against it and the outcome model would extrapolate one from the strata that do; the remedy is a coarser modifier rather than a refit. Finally, .by requires se_method = "mestimation" and a binary or categorical exposure: the stratum means and their contrasts are parameters of the stacked estimating equations, which is what puts them in the same sandwich as everything else, and a continuous exposure reports its marginal structural model's own coefficient rather than a set of standardized means to modify. Each of those errors.

An outcome model with no term reading the exposure and the modifier together warns rather than errors. The stratum effects are g-computation on the model as specified, so such a model forces the same effect on every stratum up to the covariate distribution each one has, and the modification it reports is a property of the model rather than of the data. The fit still returns.

Stabilizing the weights on the modifier is worth doing here. A stabilized weight divides by the propensity score and multiplies by a numerator, and the numerator may condition on any variable the outcome model the effects are computed from already reads, without changing what is estimated. The modifier is such a variable by construction, so a numerator conditioning on it leaves the estimator consistent for the same stratum effects and changes only how variable the weights are. That is where it pays: the part of the exposure's variation the modifier explains cancels between the numerator and the propensity score, and the weights tighten by more than the default numerator tightens them.

The converse is the constraint, and it is the reason this is a statement about the model rather than a free choice. A numerator conditioning on a variable the outcome model does not read changes the estimand rather than the precision, because the pseudo-population it builds is one in which that variable still predicts the exposure, and confounding by it remains (Robins et al. 2000; Cole and Hernán 2008; Hernán and Robins 2020, Chapter 12). Condition the numerator on the modifier and on nothing the model does not hold. The same rule governs an exposure-by-modifier interaction reported without .by, since what matters is what the reported model reads and not which argument asked for the rows.

Build the numerator with wt_ate()'s stabilization_score, evaluated at the exposure each unit actually took; see Stabilization in wt_ate() for the recipe. The sandwich treats a supplied score as a known constant, so the standard errors do not account for the numerator model having been fitted. The default stabilizer is estimated in the stacked system instead, which is one parameter wider as a result. There is a third way: pass the fitted numerator model to stabilize itself, a stats::glm() of a binary exposure or an nnet::multinom() of a categorical one, and its own equations join the stack, which is what a supplied score has none of. That accounting has something to say only where the reported model is not saturated in the variables the numerator reads; a numerator on a modifier the reported model interacts with fully divides out of every cell, so the estimates and the standard errors are the ones the unstabilized weights give.

A stabilized fit reporting effects within the strata of a modifier, and carrying the default numerator, warns with class propensity_ipw_by_stabilizer_warning and returns: the default numerator conditions on nothing, which is consistent for the same stratum effects and leaves precision on the table. A numerator model the caller supplied is read the same way, since a fitted model carries the terms it read: one that never reads the modifier warns with the same class, naming the terms it was built on, and one that reads the modifier, on its own or through a transformation of it, is silent. A stabilization_score says nothing either way, being a vector of numbers that does not carry what it conditions on.

Every reported row is a parameter of one stacked system solved on one sample, so the covariance couples them: stats::vcov() reports nonzero entries between the overall rows, the stratum rows, and the stratum contrasts, and the standard error of a stratum contrast accounts for the two strata having been estimated together. Fitting each stratum on its own subset would report the same stratum effects with no covariance between them, leaving the difference between two strata untestable.

Joint exposures

causalgenerics::joint_exposure() crosses two discrete treatments into one categorical exposure and records the crossing on the vector. The result is a factor, so ipw() fits it exactly as it fits any exposure over the same cells. What the declaration changes is what gets reported: instead of each cell against the reference cell, the result is written in the two treatments.

• One row per cell, under the effect label "mean", carrying the counterfactual risk (or mean) had everyone been set to that cell. The cell is named in the contrast column and the group is "overall".

• The simple effects: each treatment's effect within a fixed level of the other, with the contrast naming the treatment and the level it is set to ("a: 1 vs 0") and the group naming the level the other treatment is held at ("e = 1"). These include the comparisons that are not against the reference cell, which the vs-reference reporting cannot express at all.

• The interaction: the difference between two of the first treatment's simple effects, keyed by the two levels of the second being compared ("e = 1 vs e = 0"). On the risk difference scale this is the additive interaction, and on the log risk ratio scale the multiplicative one, which is the log of the ratio of the two risk ratios.

A 2-by-2 crossing therefore reports fourteen rows for a binary outcome and nine for a continuous one. The interaction is reported once, under the first treatment's framing, rather than twice: interaction is symmetric in the two treatments, so the difference between the first treatment's two simple effects is the difference between the second treatment's two, and reporting both would be one quantity under two labels.

No odds ratio is reported anywhere on this surface. It is noncollapsible, so neither a simple effect reported beside one nor a difference of two of them says what it appears to.

Every row is a parameter of the same stacked system, so the covariance couples them: the interaction covaries with each simple effect it is built from, and with each cell mean underneath those. The cell means are the marginal-mean block the categorical path already estimates, and the simple effects and interaction rows are contrasts over that block, so an interaction equals the corresponding double difference of the cell means exactly rather than to within a solver residual.

Two restrictions apply. A declared crossing is reported for the "ate" estimand alone: every cell mean on this surface standardizes to one population, and a focal estimand standardizes each of them to a population the simple effects and the interaction are not defined over, so weights built for anything else error. And .by cannot be combined with a declared crossing: effect modification of a joint intervention is a three-way question, the interaction between two treatments within the levels of a third variable, which this surface neither answers nor sensibly projects. Either restriction is escaped the same way, by dropping the declaration with factor(x) and reporting each cell against the reference cell.

Two treatment models instead of a declared crossing

A joint intervention can also be weighted through the sequential factorization it really has, $f(A | L) f(E | A, L)$: one treatment model per treatment, recorded with joint_wt_models(), and the product of their weights, built with wt_joint(). Pass the container as wt_mod and an outcome model that reads both treatments as separate columns.

The reported surface is the one above, row for row and label for label: the same cell means, simple effects, and interaction, under the same effect, contrast, and group values. The two routes estimate the same estimand through different models of the same propensity score, so they agree as estimators rather than to the last bit, and they coincide when the two parameterizations are saturated in the same covariates.

Either treatment may have more than two levels, fit with nnet::multinom() and weighted with categorical wt_ate() weights. Such a component is stacked under the multinomial score its own fit solves, which reads no order into the pair, so a categorical treatment sits in either position. The surface grows with the crossing rather than changing shape: a 2-by-3 crossing reports twenty-four rows for a binary outcome and fifteen for a continuous one, being the six cell means, each treatment's simple effects within each level of the other, and the first treatment's interaction against each non-reference level of the second. A categorical component needs no stabilization, and a caller who asks for it anyway is answered: the marginal proportions of its levels, a stabilization_score, and a fitted nnet::multinom() numerator each get the block described under Stabilizing numerators, one per component. A categorical treatment may also be paired with a dose, whose vocabulary names one level contrast per non-reference level; see the dose section below.

A nnet::multinom() fit to two levels is a logistic regression solved by a different optimizer, and it is stacked as one. joint_wt_models() records such a component as "categorical", after the model's class, but the block written for it is the binomial score, the weight factor is the binary one wt_ate() already builds from such a fit, and the pair reports the same fourteen-row two-by-two surface a pair of binomial stats::glm() fits reports.

Prefer the two-model route when the two treatments call for different adjustment sets, or when the dependence of the second treatment on the first is what you want to model directly, since each treatment then gets its own model and its own covariates. Prefer the declared crossing when one model over the cells is the natural specification. The factorization itself is validated when the container is built rather than here: joint_wt_models() refuses a second model that does not condition on the first treatment, and a pair that condition on each other.

Both treatment models are stacked alongside the outcome model, so the weights entering the outcome score are rebuilt from both propensity score blocks on every evaluation and the reported standard errors account for having estimated both. The same two restrictions apply: the "ate" estimand only, and no .by. Standard errors come from M-estimation alone.

A joint intervention with a dose

The second treatment may be a dose rather than a second discrete treatment, recorded with any of the classes the single-treatment route reads a conditional density from, and weighted with a stabilized density ratio, which wt_joint() requires of a continuous component. A dose has no cells, so there is no crossing to report counterfactual risks over: the surface is the marginal structural model's own coefficients, which for an identity-link model are exactly the weighted fit's coefficients.

The first treatment may be binary or categorical. What the position settles is the dose: it is the second factor of the factorization, the treatment whose model conditions on the other, and nothing here carries the density of a first factor that is one, so a dose in the first slot and a pair of doses are both refused by type with an error of class propensity_ipw_exposure_error.

Which of two surfaces a fit reports is decided by the marginal structural model alone. A model written in bare treatment terms and fit with an intercept reports the vocabulary surface below, whose rows name the treatment each one varies and where it is evaluated. Every other treatment-reading model reports the coefficient surface, whose rows are named after the coefficients they report.

The vocabulary surface

For y ~ a * e with a binary and e a dose, three coefficients carry causal content and each is a row:

• the binary treatment's effect at a dose of zero, with the contrast naming the treatment and its levels ("a: 1 vs 0") and the group naming where it is evaluated ("e = 0");

• the dose's slope at the binary treatment's reference level ("e: per unit", "a = 0");

• the interaction, the change in the binary treatment's effect per unit of the dose, keyed by the one-unit step of the other treatment ("a: 1 vs 0", "e + 1 vs e").

The effect column names the scale, keeping the vocabulary each treatment type already uses: a level contrast is a diff and a per-unit change in the dose is a slope under an identity link, and both are log(or) under a logit link and log(rr) under a log link. An additive model, y ~ a + e, forces one effect for the binary treatment at every dose and one slope at either of its levels, so it reports two rows under the group "overall" and no interaction.

A categorical first treatment reports those same three readings with the level contrasts run out over the levels it has: one contrast row per non-reference level at a dose of zero, one slope row at the reference level, and one interaction row per non-reference level. For y ~ z * e with z taking "lo", "mid", and "hi", that is five rows, in the outcome design's column order:

• "z: mid vs lo" at "e = 0" and "z: hi vs lo" at "e = 0";

• "e: per unit" at "z = lo";

• "z: mid vs lo" and "z: hi vs lo", each at "e + 1 vs e".

A binary treatment is the one non-reference level of that rule rather than a case of its own, so the three rows above are this surface at two levels. An additive model reports the level contrasts and the slope under the group "overall" and no interaction, K rows for K levels. Every row is one coefficient of the weighted fit, so the surface grows by a row per level rather than by a table of contrasts, and per-arm dose slopes, which are not coefficients of this model, are built downstream from coef() and vcov(). No guard here reads how thin an arm is: a level the dose model rarely sees has its conditional density read far from any data, and the product weight carries that as a heavy tail. Check the weights' effective sample size and their spread within each level, as for any joint weight, before reading the rows.

Those three readings hold of a model that is linear in each treatment and carries an intercept, and of no other, which is why such a model reports this surface and no other model does. A model fit without one, written - 1 or + 0, reports the coefficient surface: dropping the intercept expands a factor treatment to an indicator for every level rather than for the non-reference levels alone, and where the columns survive it, as a 0/1 numeric treatment's do, the forced zero moves what the coefficients mean instead, so no row of this vocabulary would be true of the fit either way.

The columns are read as well as the formula: a factor treatment under a coding other than treatment contrasts leaves the terms bare while rescaling or recentering the column, so its coefficients are no longer the effects these rows name, and such a fit is refused rather than reported in this vocabulary. Refit the outcome model with the treatment as an unordered factor under treatment contrasts. A treatment with two levels also has a numeric coding whose bare term contributes that column, 0 for the reference level and 1 for the other; a treatment with more than two levels has none, so there it is the unordered factor alone. The refusal belongs to this surface rather than to the fit. A model the vocabulary has no reading for reports the coefficient surface instead, and a model written without an intercept is one such, so an ordered factor there contributes an a.L:e row named after the column it multiplies, claiming nothing about a coding.

An interaction of the treatment with the dose written where the dose has no term of its own, as y ~ a + a:e, is refused here too, as is a crossing carrying no treatment term either, y ~ a:e. R codes a factor in such an interaction with one indicator per level rather than with the contrasts it carries, so the design gains a reference-level column no row names and each column after it holds the level before it. Write the crossing as y ~ a * e, which gives the dose a term of its own and reports on this vocabulary.

The dose column is read the same way, against the dose the treatment models were fit to. A dose transformed in the outcome model's own formula, as I(e / 10) or scale(e), leaves that column alone and is no longer a bare term, so the fit reports the coefficient surface. A dose rescaled or recentered in the data after the treatment models were fit is refused on this surface instead: the models have to agree on what the variable is, and the weights the outcome model carries were built for the dose those models hold. Working on a rescaled dose is supported either by rescaling it before every model is fit or by writing the transformation into the outcome formula. Neither refusal reaches the coefficient surface, which reports whatever columns the model holds and names each row after its own column: a bare-term model on a rescaled dose fit without an intercept, y ~ a * e - 1, reports the coefficients and standard errors that reparameterization leaves.

The coefficient surface

A marginal structural model carrying a transformed or a basis treatment term, such as y ~ a * sin(e), y ~ a + e + I(e^2) + a:e, or y ~ a * splines::ns(e, 3), reports one row per treatment-reading coefficient. The contrast column carries the coefficient name exactly as the fit writes it, and there is no group column at all: the surface makes no claim about where a row is evaluated, because for a curve there is no one place. Build an interpretable dose response downstream from coef() and vcov(), whose covariance between the rows is what that takes. Marginalizing the curve over the observed doses is a separate estimand that this package does not compute, here as on the single-treatment route: use the marginaleffects package on the result, avg_slopes() for slopes, avg_comparisons() for contrasts, and avg_predictions() for causal dose-response functions.

Unlike the single-treatment route, this one announces nothing and refuses nothing for such a fit. There the surface a basis exposure leaves is the only reading the result has, so the reading is declared; here both readings are computed whatever the marginal structural model is, and effects is answered either way.

The effect label follows the outcome link, coef at an identity link and the same log(or) and log(rr) as everywhere else. coef is a word of its own rather than diff or slope, both of which carry an evaluated-at claim on the vocabulary surface.

Covariates and mixed terms

A covariate entering on its own is admitted on either surface and contributes no row, however many columns it expands to: it adjusts the marginal structural model and the surface has nothing to say about its coefficient. A term reading a treatment and a covariate together is refused on both, since its coefficient is a change in an effect per unit of that covariate and no row could name the effect it stands for. Nothing here is standardized over the covariates, which is what the declared crossing does instead.

Unlike the single-treatment route, this one needs no .data for a basis marginal structural model. Without one, each treatment and each treatment design is read off the model that owns it, so the outcome model frame is never asked for a column it does not hold.

A .data the caller does supply is the frame every design here is rebuilt from, on the terms the single-treatment routes rebuild theirs: both treatment columns, both treatment designs, each component's numerator design, and the counterfactual designs are built over the rows every model read and under the coding each fit recorded. A column one of those fits reads that is missing from .data errors with class propensity_columns_exist_error, and one supplied as another type errors with class propensity_ipw_data_error, naming the fits that recorded it, a treatment model by the treatment it fits. A set of fits made under na.action = na.exclude is reconciled with the frame they were given by restricting .data to the rows they analyzed.

The dose model contributes its coefficients and the conditional variance its density ratio divides by, all as parameters of the same stacked system. Each component's stabilizing numerator is estimated in a block of its own; see Stabilizing numerators, one per component.

Stabilizing numerators, one per component

Either component may be stabilized, and each component's numerator is estimated in its own block of the stacked system, so the standard errors account for that numerator having been fitted rather than reading it as a constant. What the block holds follows what the component's weights record. A binary component stabilized by the default numerator contributes the one marginal proportion that numerator is, a categorical one the $k - 1$ free proportions its levels leave, and a dose stabilized by the default numerator the exposure's two marginal moments. A component stabilized on a fitted model passed to wt_ate()'s stabilize contributes that model instead: one parameter per coefficient for a binary component, one run of them per non-reference level for a categorical one, and for a dose one per coefficient plus one for the spread its density is read at. A dose stabilized with numerator = "integrated" is built from the dose block and the data alone and contributes nothing, and an unstabilized discrete component contributes nothing either.

Those parameters are named ⁠stab_<component>_<parameter>⁠, where ⁠<component>⁠ is the name the component was recorded under in joint_wt_models(): a joint system carries two components, and a name saying only which role a parameter plays would not say whose. A categorical component composes that prefix with the multinomial layout its own block takes, ⁠<level>:<term>⁠ for a fitted numerator and the level alone for a marginal one, so a treatment z with levels lo, mid, and hi reports stab_z_mid and stab_z_hi for the marginal proportions and stab_z_mid:(Intercept), stab_z_mid:v, stab_z_hi:(Intercept), and stab_z_hi:v for a numerator of v.

A numerator model is read through the guards its own single-treatment route reads it through. A binary component's is read as an unpenalized binomial fit at a logit, probit, or cloglog link; a categorical component's as an nnet::multinom() fit that declares the treatment's levels in the order the component's own model declares them; and a dose's goes through the registry under Continuous propensity score models. A multinomial numerator fit under another level order, a model of a response other than the treatment it stabilizes, or one fit to a different set of observations errors with class propensity_ipw_numerator_error; one fit with case weights errors with class propensity_ipw_ps_weights_error. Every numerator design here is one of the designs a .data the caller supplies rebuilds, on the terms the single-treatment routes rebuild theirs, so a numerator that keeps no model frame and whose fitting call can no longer be evaluated errors with class propensity_ipw_data_error without one, and is rebuilt from .data with one.

Both components' numerators arrive in the same argument, one stabilize per component, so every refusal about one names the component it was built for beside that argument, stabilize for z. wt_mod here is the container the two treatment models arrived in rather than a model anyone could refit, so a refusal reading a numerator against its treatment model names that component too.

What a numerator may condition on is a statement about the reported model rather than a free choice, and it is the same statement here as for a single treatment. A numerator may condition on any variable the marginal structural model already reads, such as a modifier V: $P(Z | V)$ for a discrete treatment and $f(A | V)$ for a dose leave the estimator consistent for the same effects and only tighten the weights. A numerator conditioning on a variable the reported model does not read changes what is estimated rather than the precision, because the pseudo-population it builds is one in which that variable still predicts the treatment, and the estimates are then read with V or not at all. The accounting a stacked numerator adds has something to say only where the reported model is not saturated in the variables the numerator reads: a numerator that divides out of every cell of the reported model moves neither the estimates nor their standard errors.

The dose model and the ratio its weights are

The dose model goes through the registry described under Continuous propensity score models, so the classes and links this route stacks are the ones the single-treatment route stacks: an stats::lm(), a gaussian stats::glm() fit with an identity or a log link, a MASS::rlm() fit with one of the psi functions MASS supplies, and an mgcv::gam() of a gaussian family read at its own penalized score. A gaussian dose model fit with another link errors with class propensity_ipw_link_error, and a robust or an additive fit whose score the system cannot write raises the same conditions it raises for a single dose. The limitations of the additive route are the ones listed under Continuous propensity score models, and they apply here unchanged.

Every score this route stacks is the unweighted one, so a component model fit with prior case weights is not at the root of the block written for it and errors with class propensity_ipw_ps_weights_error, naming the component rather than the container the two models arrived in. Refit that component without weights.

The dose component's weights may be any ratio wt_ate() builds: a heavier-tailed .density, numerator = "integrated", or a single .sigma the caller fixed. wt_joint() records each component's ratio on the product, and the stacked system rebuilds the dose factor from that record, so what is differentiated is the weight the outcome model was fit with. An integrated numerator is built from the dose block and the data alone, so it adds no stabilization parameters; it is the conditional density averaged over the units the dose model was fit to and read at points spanning their dose, so it is rebuilt over those rows the way a single dose's is, and a dose model that keeps no model frame and whose fitting call can no longer be evaluated errors with class propensity_ipw_data_error under it, .data or no .data. A fixed scalar .sigma is a known constant, so the dose block is its coefficients alone and none of that number's uncertainty is carried. A spread supplied for each observation has no counterpart in the system and errors with class propensity_ipw_sigma_error.

A component built with a stabilization_score is rebuilt from the score itself. wt_joint() records each component's score on the product, and the sandwich treats it as a known constant here exactly as it does for a single treatment, so it adds no stabilization parameters and none of its uncertainty is carried. A product built by a version of this package that recorded no score, one assembled by hand, or one subset after it was built, which drops a score held per observation, says a component was stabilized without saying what by: the system stands that component's own marginal numerator in, and the weight-consistency check reports the difference, naming the components whose numerator was stood in for. Rebuild such weights with wt_joint() to have the score read back. A record left whole over fewer rows is refused instead: a model frame that drops incomplete rows re-attaches the record as it was, leaving a component's score at the length the product was built at, and a score that describes other observations than the ones being weighted errors with class propensity_ipw_stabilization_score_error.

Dose weights built with a "kernel" density are refused with class propensity_ipw_se_method_unavailable_error. The bandwidth is not a function of the parameters that the sandwich could differentiate, and the package has no second method to fall back on, so the refusal points at building the dose weights from a density the stacked system can write, or at bootstrapping the whole fit by hand.

Variance estimation

By default (se_method = "mestimation"), standard errors are computed by M-estimation. The propensity score model, the weighted outcome model, the marginal means, and the effect contrasts are stacked as a single system of estimating equations, and the empirical sandwich variance of that joint system is used. The stacked marginal means are standardized to the estimand's tilted target population, so a non-ate estimand with a covariate-adjusted outcome model reports the contrast for that target population rather than the full-sample average. Stacking the models accounts for the uncertainty introduced by estimating the propensity scores, avoiding the underestimated standard errors that arise from treating estimated weights as fixed. See Stefanski and Boos (2002) for the M-estimation framework.

M-estimation standard errors are available for all exposure types: binary (from a stats::glm() propensity score model), categorical (from a nnet::multinom() model), and continuous (from an stats::lm(), a gaussian stats::glm(), a MASS::rlm(), or an mgcv::gam() model; see Continuous propensity score models below). The atm weight pmin(e, 1 - e) is not differentiable at a propensity score of 0.5; deli's central finite difference straddles the kink there and averages the one-sided slopes. Its effect on the variance is negligible unless many observations sit at exactly 0.5.

A binary exposure's propensity score model is read as an unpenalized logistic fit, whatever the standard error method: the stacked system solves that score, and both the linearization and the robust diagnostic differentiate it. An mgcv::gam() fit of a binomial family is not at the root of that score but of a penalized one, so it is refused with class propensity_ipw_se_method_unavailable_error rather than described by a standard error for a model that was not fit. What is refused is the standard error and not the weights: wt_ate() reads the same fit and builds them from its fitted probabilities. Write the shape each smooth term carries out as columns of a stats::glm() formula, or bootstrap the whole fit by hand. An additive dose model is stacked at its penalized score instead, as described under Continuous propensity score models.

M-estimation standard errors for a categorical exposure allocate memory roughly linearly in the number of observations, on the order of 70 to 90 kilobytes per observation for a single fit, so a sample of 10,000 observations allocates on the order of hundreds of megabytes. The cost comes from the stacked estimating-equation machinery rather than any single term, and for very large samples you can expect long, garbage-collection-heavy fits. This is expected behavior.

Stabilizing numerators for a discrete exposure

A numerator estimated by a model the caller passed to wt_ate()'s stabilize is estimated here as well, whichever type the exposure is. A binary exposure's numerator model contributes one parameter per coefficient, named ⁠stab_<term>⁠, and a categorical exposure's contributes one per coefficient of its nnet::multinom() fit, named ⁠stab_<level>:<term>⁠, so that the level a coefficient belongs to is read off the name. The multinomial numerator's own score is stacked alongside the propensity score model's, written as the same multinomial estimating equation, so the standard errors account for the numerator having been fitted rather than reading it as a constant. The default stabilizer contributes what it estimates instead: the single marginal proportion stab_pi for a binary exposure, and the K - 1 free marginal proportions ⁠stab_<level>⁠ for a categorical one. A stabilization_score contributes nothing at all, being a number the caller fixed. What a dose's numerator contributes is described under Continuous propensity score models.

Reading an estimated numerator as known would be conservative for these weights rather than wrong, since stabilization is offered for the ate estimand alone (Lunceford and Davidian 2004; Kostouraki et al. 2024), and Shu et al. (2020) stack the numerator and the denominator of stabilized weights for that reason. Stacking the multinomial score therefore tightens the inference rather than rescuing it. Neither result isolates the numerator's own contribution, and Kostouraki et al.'s tutorial is written for a binary exposure, so what they give is the direction of the correction rather than a bound on its size.

A numerator model's design is one of the designs ipw() rebuilds when .data is supplied: it is rebuilt over the rows every model read and under the coding the fit recorded, so a column only the numerator reads is asked of .data and is held to the type the fit recorded for it. A .data missing such a column errors with class propensity_columns_exist_error, and one supplying it as another type with class propensity_ipw_data_error. Without .data the design is read off the fit itself, and a fit that keeps no model frame rebuilds one by re-evaluating its fitting call. A stats::glm() fit with model = FALSE is such a fit, so one made inside a function whose frame is gone errors with class propensity_ipw_data_error, and the remedy is to supply .data.

A multinomial numerator model is held to the propensity score model's levels in the propensity score model's own order. Everything its block contributes is positional: the coefficients are stacked level-major, and the softmax the numerator is rebuilt from reads the first level as the reference the others are contrasted with, so a fit made under another order is a different parameterization rather than a permutation of this one. Such a fit is refused with an error of class propensity_ipw_numerator_error, which is also what a model of a response other than the exposure, or one fit to a different number of observations, errors with. One fit with case weights errors with class propensity_ipw_ps_weights_error, for the reason a weighted propensity score model does. nnet::multinom() keeps no model frame, so the design its coefficients were fit over is recovered by re-evaluating the fitting call; a fit whose call can no longer be evaluated, such as one made inside a function whose frame is gone, errors with class propensity_ipw_data_error, and the remedy is to supply .data, which the numerator's design is rebuilt from.

A categorical exposure reports M-estimation standard errors alone, so a multinomial numerator model meets neither of the other two methods. The linearization's refusal of a fitted numerator model, below, and the robust diagnostic's reading of every weight as a known constant are both matters for a binary exposure.

Continuous propensity score models

A continuous exposure's weights are a ratio of densities centered on the conditional mean the propensity score model fits, so the stacked system carries that model's own score. Four classes supply one: an stats::lm(), a gaussian stats::glm() fit with an identity or a log link, a MASS::rlm() fit with one of the psi functions MASS supplies, and a gaussian mgcv::gam() fit with one of the same two links. A model fit with a family that is neither gaussian nor binomial is neither a dose model nor a binary propensity score model, and errors with class propensity_model_family_error naming the family it has. A gaussian model fit with any other link is refused by name, because the coefficients its iteration stops at are not a tight enough root to seed the stacked solve from; refit it with one of the two links or as an stats::lm(). Any other class errors, since solving it here would report a propensity score model that was never fit. Each class is read as the class it is rather than by what it inherits from, so a subclass of a supported fit errors as well.

A MASS::rlm() is stacked as the psi score its coefficients are the root of, clipped where the fit itself clipped. Each of the three psi functions MASS supplies has a loss of its own in the stacked system: MASS::psi.huber(), MASS::psi.bisquare(), and MASS::psi.hampel(). rlm() clips the residual divided by its scale estimate at the psi's own constants, and the stacked score clips the raw residual, so each constant is multiplied by the scale the fit settled on. Those constants are read off the psi the fit carries, which is where rlm() records a caller's choice of them, so a retuned psi is stacked at the values it was retuned to. That scale enters the location score as a known constant. It solves an equation of its own that the stacked system does not write, so none of its sampling variability is propagated, and the standard errors are those of a psi score read at a scale treated as fixed. The spread of the conditional density is a separate quantity, and comes from the estimator its family asks for, there as for every other class.

MASS::psi.bisquare() and MASS::psi.hampel() redescend, and that is worth knowing about the standard errors they get. A redescending psi's estimating equation has more than one root, so there is no single set of coefficients it identifies and the solve reports whichever root it is seeded at. The seed here is the fit's own coefficients, so what comes back is the root the user has, and the sandwich is the covariance of that root read locally. It carries no claim about the other roots, and it does not describe the sampling behavior of a procedure that might have converged to one of them instead. Use a bootstrap of the whole fit if that behavior is what you need.

Weights built with .sigma = fit$s, the robust scale, are read as the fixed spread they are, but that combination is usually a poor choice: the robust scale resists the very residuals the density is meant to spread over, so on contaminated data the conditional density is far narrower than the residuals are, and the ratio can collapse onto an effective sample size of about one. Leave .sigma unset unless you have a reason to hold the spread at a particular value.

Three robust fits are refused. One fit with a psi function of the caller's own has no loss here to write its score as, and one fit with method = "MM" reaches its coefficients through a high-breakdown start that decides which root of a redescending psi the fit finishes at and supplies the scale it clips at, neither of which is an equation this system writes. A sandwich read locally at those coefficients would therefore not describe how an MM fit behaves across samples. Both error with class propensity_ipw_robust_psi_error. One whose iteration stopped short of its own tolerance is the root of nothing, and errors with class propensity_ipw_convergence_error. The first points to a refit with one of the three psi functions or to a bootstrap of the whole fit written by hand, the second to the default method = "M", and the last to a larger maxit or a looser acc, which is what a fit needs to reach its own tolerance.

An mgcv::gam() is stacked as the penalized least squares it is, with its smoothing parameters held at the values the fit reports. Its design is the smooth basis it was built on rather than the columns its formula names, and its coefficients are the root of the least-squares score of that basis less the penalty those smoothing parameters define. An mgcv::bam() is the same fit computed for larger data and is read the same way. The additive part of an mgcv::gamm() fit is not: it carries the gam class alone, without the glm and lm classes ipw() reads a model through, and errors with class propensity_method_error. Refit it with mgcv::gam() or mgcv::bam(). A penalty attached to a parametric term, such as one from paraPen, belongs to no smooth term and is not one this route can place, so such a fit errors with class propensity_ipw_se_method_unavailable_error. A fit made under a smoothing floor from min.sp errors with the same class, because mgcv adds that floor to the penalty and leaves it out of the smoothing parameters the fitted object reports, so nothing in the fit records the penalty it was made under; refit without min.sp, which mgcv::bam() ignores in any case.

Four limitations come with the additive route, and each is a property of the method rather than of a particular fit.

1. The variance conditions on the smoothing parameters rather than estimating them, so none of the uncertainty of having chosen them is carried. A simulation put that at roughly three percent of the standard error at n = 500, roughly seven percent at n = 200, and less at larger samples. A fit handed its smoothing parameters and a fit that chose them report the same standard error here, which is the same limitation stated as an equality.

2. What the propensity score block reports is the frequentist variance of a penalized fit rather than the Bayesian covariance mgcv reports by default, so vcov() on the stacked fit reads like vcov(fit, freq = TRUE) and will not match vcov(fit).

3. The envelope is the one every other continuous fit has: a gaussian family, an identity or a log link, and no prior weights. A fit outside it is refused where every other model of the wrong family, link, or weighting is.

4. Coverage degrades as the weights grow heavier tails, and the interval is optimistic once the effective sample size falls much below 60 percent of the observations. That holds whatever the dose model is: the same simulation found it for a correctly specified stats::glm() on the same data.

Every density wt_ate() builds a ratio in is stacked as it was built: the normal, Laplace, and Student t families, and a density you wrote yourself, under either the marginal or the integrated numerator. ipw() reads the family, the numerator, and the spread off the record the weights carry and rebuilds the same ratio at each value of the parameter vector, so the weights it differentiates are the weights the outcome model was fit with. A marginal numerator contributes the exposure's two marginal moments to the parameter vector; an integrated one is built from the propensity score block and the data alone and contributes none. That numerator is the conditional density averaged over the units the propensity score model was fit to, read at points spanning their exposure, so it is rebuilt over those rows even where .data leaves the rest of the system standing on fewer of them. Rebuilding it needs the design those rows enter through, which .data cannot stand in for: a propensity score model that keeps no model frame, made inside a function whose frame is gone, errors with class propensity_ipw_data_error under numerator = "integrated", and the remedy is to refit it where its data is available. Weights carrying no such record, including any written with psw() by hand, are read as the normal ratio the package built before the record existed.

A numerator estimated by a model the caller passed to wt_ate()'s stabilize contributes that model: one parameter per coefficient, named with a stab_ prefix so the terms it shares with the propensity score model stay apart from it, and one for the spread its density is read at, sigma2_n. That parameter holds a squared scale rather than a variance, which are the same number only for a family whose scale is its standard deviation. Its score, and the equation that spread is the root of, which is the residual moment or the family's own score depending on what the density asks for, are stacked alongside the propensity score model's, so the standard errors account for the numerator having been fitted, which is what separates this from a stabilization_score: a score the caller computed is carried as a known constant and contributes no parameter at all. The numerator model is read through the registry above, so a class, family, or link whose score this system cannot write is refused there in the terms that registry refuses a propensity score model in, naming stabilize. A model of a response other than the exposure, or one fit to a different set of observations, errors with class propensity_ipw_numerator_error. One fit with case weights errors with class propensity_ipw_ps_weights_error, for the reason a weighted propensity score model does: the block written for it is its unweighted score, and its coefficients are not at the root of that.

A dose numerator model's design is one of the designs ipw() rebuilds when .data is supplied, on the terms a discrete exposure's is: it is rebuilt over the rows every model read and under the coding the fit recorded, so a column only the numerator reads is asked of .data and is held to the type the fit recorded for it. A .data missing such a column errors with class propensity_columns_exist_error, and one supplying it as another type with class propensity_ipw_data_error. Without .data the design is read off the fit itself, and a fit that keeps no model frame rebuilds one by re-evaluating its fitting call. A stats::lm() or stats::glm() fit with model = FALSE is such a fit, so one made inside a function whose frame is gone errors with class propensity_ipw_data_error, and the remedy is to supply .data.

One choice is refused for the standard error rather than for the model. A "kernel" density chooses its bandwidth from the residuals of the propensity score model, which is not a function of the parameters that the sandwich could differentiate, so weights built from it have no closed-form standard error and raise an error of class propensity_ipw_se_method_unavailable_error. The package offers no resampling method of its own, so the remedy the refusal names is a bootstrap of the whole pipeline written by hand: resample the rows, refit the propensity score model, rebuild the weights with wt_ate(), and refit the outcome model on each resample, then read the spread of the coefficient across the resamples. Doing it that way rather than inside ipw() keeps what each replicate re-estimates, such as a kernel bandwidth, under the writer's control. It is also how to carry a smoothing choice that the additive route conditions on: a bootstrap that refits the mgcv::gam() each time re-chooses the smoothing parameters along with everything else.

Standard errors by linearization

Setting se_method = "linearization" instead uses the influence-function linearization of Kostouraki et al. (2024). It is available only for a binary exposure and only for the ate, att, ato, and atm estimands; a categorical or continuous exposure, or the atu or entropy estimand, requires se_method = "mestimation". The linearization also supports only an outcome model whose formula contains the exposure alone. Its point estimates come from g-computation on the fitted outcome model, while its influence functions are those of the Hajek weighted-mean estimator; the two agree only for an exposure-only outcome model. A covariate-adjusted outcome model (any term beyond the exposure, including covariates, interactions, or transformed terms) errors on the linearization path and requires se_method = "mestimation", which stacks adjusted outcome models correctly. The conditional reading is restricted the same way: linearization stacks no system, so the outcome model it stores carries no block of one, and the standard errors that reading reports are refused with a classed error rather than replaced by the covariance the outcome model computed for itself, which treats the estimated weights as fixed.

Weights stabilized on a fitted numerator model are refused here as well, with an error of class propensity_ipw_se_method_unavailable_error. The influence functions read the stabilizer as a known constant, which is exact for the default stabilizer and for a stabilization_score the caller fixed and is not exact for a model fit to the same data, so the method is refused rather than a standard error reported for a numerator nobody fit. se_method = "mestimation" estimates the numerator model alongside everything else. se_method = "robust" reads every weight as a known constant by design, that being what the diagnostic is, so it is not refused.

The linearization outcome model must also carry an intercept, which is a stricter requirement than the M-estimation path imposes. Under a numeric coding such as y ~ z - 1 the linear predictor under no exposure is fixed at the link's zero point rather than estimated, so the g-computation marginal means no longer match the Hajek means the influence functions describe. The rejection is broader than that mechanism: a saturated factor coding such as y ~ 0 + zf does estimate both cell means, and the M-estimation path accepts it, but the linearization influence functions are derived for the intercept parameterization, so every no-intercept outcome model errors here. See Model requirements for the baseline contract both methods impose.

Standard errors as a diagnostic

Setting se_method = "robust" reports the sandwich the weighted outcome model computes for itself, with the weights entering as known constants. It is the linearization route with the correction for having estimated the propensity score dropped, so the point estimates are the linearization ones and only the standard errors move. What it reports for the two counterfactual means and the contrasts built from them is the delta-method reading of sandwich::vcovHC(outcome_mod, type = "HC0"), exactly: the influence variance is read as a plain mean of squares rather than as a sample variance, so what comes back is the HC0 sandwich itself rather than a finite-sample rescaling of it.

This is offered as a diagnostic and not as an inference. It ignores that the weights were estimated, and the correction it drops is generally what a propensity score model contributes to the variance, so the standard errors it reports generally understate it, and understate it most where the weights are least stable. Use "mestimation" or "linearization" to report a standard error. Reach for the diagnostic to see what the propensity correction is worth on a given fit, by reading it beside one of them.

A result fit this way carries the class "ipw_diagnostic_se" ahead of "ipw", and printing it writes one line naming the method, so the numbers are not read as accounting for the propensity score model. as.data.frame() and tidy() carry the same mark as an ipw_se_diagnostic attribute on the table. Every requirement of the linearization route applies here: a binary exposure, the ate, att, ato, and atm estimands, an outcome model of the exposure alone fit with an intercept, no .by, and no conditional reading. The categorical, continuous, and joint treatment routes refuse it as they refuse linearization.

The printed mark does not survive pool_ipw(): a pooled object is a summary of results rather than a result, and what says its standard errors are the diagnostic ones is the method it records, se_method = "robust".

Multiple imputation

With missing data, fit the whole analysis once per imputed dataset and pool the results. Everything the analysis needs is rebuilt inside a single expression, so each imputation gets its own propensity model, its own weights, and its own outcome model:

imp <- mice::mice(dat, m = 20, print = FALSE)
fits <- with(imp, {
  ps <- glm(z ~ x1 + x2, family = binomial())
  w <- wt_ate(ps)
  om <- glm(y ~ z, family = quasibinomial(), weights = w)
  ipw(ps, om)
})
pool_ipw(fits)

pool_ipw() is the recommended verb. It reads the results themselves rather than a tidied table of them, so the effect labels survive, a categorical result keeps the contrast each row reports, and the complete-data degrees of freedom fall back to the outcome models whenever a result records none of its own. That fallback is written for the condition rather than for one route: a fit under se_method = "linearization" or se_method = "robust" records none, and so does a result another package built on the same class from estimating equations that report no residual degrees of freedom. It also takes the smallest degrees of freedom across the pooled results, which is the conservative choice when they differ. tidy() and glance() report what it returns.

mice::pool() also works, for every exposure type, and is the right choice when a result has to travel through the same pipeline as other analyses. Four things are worth knowing about that route. It reaches the result through tidy(), which pools a categorical result grouped by term and contrast together, a grouping mice has supported since 3.15.0. mice::pool() groups correctly but its summary() prints only term, so a categorical pooled table shows each effect measure repeated with no contrast label beside it, and the labels have to be read off pooled$pooled instead. A result that records no complete-data degrees of freedom of its own, a fit under se_method = "linearization" or se_method = "robust" among them, leaves the pooled degrees of freedom missing unless dfcom is passed explicitly, as in mice::pool(fits, dfcom = df.residual(fits$analyses[[1]]$outcome_mod)). Every result of that kind has to be told, where pool_ipw() reads the outcome models for all of them unasked. That remedy is what the package's requirement of mice 3.18.0 or later is for: mice 3.17.0 introduced a regression in the dfcom argument of pool(), and 3.18.0 is the version that repairs it. And the exponentiate argument of summary() on a pooled mipo is not the one these methods take: a mipo records no scale for the rows it holds, so it exponentiates every one of them, returning a risk difference as its exponential and still labeling the row rd. tidy() and the pooled result's own frame exponentiate the rows reported on the log scale alone and relabel those.

To pool the conditional reading, record it inside the same expression with as_conditional(), and the pooled result reports the outcome models' coefficients rather than the causal contrasts.

Recording the reading beforehand is one route to it rather than the only one. pool_ipw() pools both readings of the results it is given, storing the one they record and keeping the other beside it, so as_marginal() and as_conditional() move a pooled result between them after the pooling and report a table the pooling already built rather than pooling a second time. tidy(), stats::coef(), stats::vcov(), stats::confint(), and the pooled result's own frame each take an effects argument naming a reading for the one call, which leaves the pooled result as it is. mice::pool() offers neither: it reaches a result through tidy() and returns a mipo holding the one table it asked for, so on that route the reading still has to be recorded before the pooling.

A reading the pooling could not compute is refused rather than reported under the other one's name. The conditional reading needs the covariance the joint estimation of the weights and the outcome implies, and a fit under se_method = "linearization" or se_method = "robust" stacks no such system and records no such block, so a set of those fits pools the marginal reading alone. The pooled result records why the other reading is missing, and asking for it, whether by flipping the result or by naming it for one call, errors with that recorded explanation under the classes causalgenerics_pool_missing_surface_conditional and causalgenerics_pool_missing_surface.

The analysis belongs inside with() rather than outside it because the propensity score model has to be estimated once per imputation. Weights built from propensity scores averaged across imputations are identical in every imputation, which leaves no between-imputation variance for the weights to contribute and no per-imputation estimation for the corrected standard errors to account for; both components of the variance are lost. Fitting within each imputation keeps the uncertainty of having estimated the weights inside each result, and the pooling adds the uncertainty the imputation itself contributed. Leyrat et al. (2019) compare the approaches and recommend this one.

Model requirements

ipw() cannot yet account for propensity scores that were trimmed (ps_trim()), truncated (ps_trunc()), or calibrated (ps_calibrate()) before weighting. The stacked estimating equations rebuild the weights from the propensity score model on every evaluation, so a weight that is no longer a deterministic function of that model breaks the sandwich variance. Supplying weights built from a modified score errors on either standard error method; refit the weights from the unmodified propensity score model. Density-trimmed weights for a continuous exposure (ps_trim() with method = "density" or method = "resid") carry the trim record and are refused as trimmed weights.

Weights bounded with wt_trunc() are refused as well, with their own error, for a different reason. A bound on the weights is not a smooth function of the propensity score model's parameters, so the stacked sandwich would report an interval for the truncated estimand without saying so. The four methods also differ in kind. A fixed bound ("wt", or "adaptive", which depends only on the sample size) is known before the weights are seen, while a percentile bound ("pctl", or "count", which bounds at an observed weight) is estimated from the weights it modifies and would need an estimating equation of its own. Only the fixed case is a candidate for support, and relaxing the refusal for it first requires a simulation that measures the coverage of the M-estimation interval with the bound inside the stack. Until then, rebuild the weights without wt_trunc() and report the M-estimation interval, or fit the outcome model with the truncated weights and report a fixed-weight sandwich, stating that it conditions on the bound.

An outcome model fit without weights also errors on either method. The outcome model must not carry an offset term on either method, since neither the stacked outcome score nor the linearization influence functions thread an offset; supplying one errors.

For a binary or categorical exposure, the outcome model formula must contain the exposure. The counterfactual designs are built by setting the exposure to each level in turn, so a model fit on covariates alone gives one identical design per level; it errors on either standard error method.

For those two exposure types, the outcome model must also be able to represent a baseline at every level, for example through an intercept or through a saturated factor coding of the exposure. A numeric no-intercept coding such as y ~ z - 1 errors on both standard error methods, for different reasons. On the M-estimation path the counterfactual design at the zero-coded level is identically zero, so the marginal mean there is fixed by the outcome link (0.5 under a logit or probit link, 0 for a linear model) rather than estimated from the data. On the linearization and robust paths the intercept is required outright, because without it the g-computation means stop matching the Hajek means the influence functions describe. A saturated factor coding such as y ~ 0 + zf is a reparameterization whose designs are the level indicators, so the M-estimation path accepts it and reproduces the with-intercept fit; the linearization and robust paths still require the intercept. A no-intercept model that keeps a covariate, such as y ~ z + x1 - 1, still estimates the marginal mean from that covariate and runs on the M-estimation path. A continuous exposure has no counterfactual designs, since ipw() reports the marginal structural model's own exposure coefficient, and is unaffected.

Three further requirements apply to the weights and the propensity score model. The outcome model weights must match the values implied by the propensity score model; a mismatch errors, on both standard error methods.

For a continuous exposure that requirement bears on the spread of the conditional density. ipw() stacks a single conditional spread alongside the propensity score coefficients, by the moment the pooled spread is the root of or, for a density built with sigma_method = "mle", by the score of the family it names, so the weights it rebuilds are the ones wt_ate() produces when no .sigma is given. A single .sigma is taken instead as a known constant: the weights are rebuilt at the number that was supplied, and the stacked system carries none of that number's uncertainty, which is what fixing a spread says. Weights built with an observation-level .sigma, such as influence(model)$sigma, are a different function of the data with no counterpart in the stacked system, and are refused before anything is solved, with an error of class propensity_ipw_sigma_error. Rebuild the weights with the pooled default or with one number to use ipw().

It bears on a stabilization_score in the same way. A score written per observation is one value per unit. Subsetting or reordering the weights with [ carries it along, but a slice that cannot place it, such as vctrs::vec_slice() or dplyr::arrange(), drops it, and a model frame that drops incomplete rows leaves it at the length the weights were built at. Either way the record still names a score as the numerator and the score in hand no longer describes the observations being weighted, so ipw() refuses before anything is solved, with an error of class propensity_ipw_stabilization_score_error. Rebuild the weights on the rows being analyzed, or stabilize on a single stabilization_score, which scales every weight and survives any restriction.

Weights stabilized on a fitted numerator model carry a further requirement, on outcome_mod rather than on the weights. A numerator fit on a covariate divides every weight by a quantity that varies with it, so the pseudo-population the weights build is one in which that covariate still predicts the exposure. Confounding by it remains there, and what the weighted fit identifies is the effect within its levels rather than the marginal one, which only a marginal structural model carrying it reports (Robins et al. 2000; Cole and Hernán 2008). ipw() compares the two: the variables on the right-hand side of the numerator model's formula against the variables outcome_mod's formula reads anywhere, and any the second does not read are reported once for the call, in a warning of class propensity_ipw_stabilizer_coverage_warning naming each of them. The comparison is of variables rather than of terms, so a model reading factor(v), poly(v, 2), or splines::ns(v, 3) reads v; a marginal structural model conditional on a transformation of a covariate is conditional on the covariate. The joint route asks the same question of each component's numerator and reports the answers together. The estimates are the ones the fit gives either way, so this is a warning and not a refusal: add the variables to outcome_mod, or refit the model given to stabilize without them. The default stabilizer, an intercept-only numerator model, and a stabilization_score all condition on nothing this can name and say nothing.

The propensity score model must not separate the exposure. A model whose covariates predict the exposure without error has no finite maximum likelihood estimate, and the propensity scores its coefficients imply reach exactly zero or one, leaving the corresponding weights undefined. For a binary exposure, whose propensity score model is a binomial glm(), both standard error methods reject such a fit at the same threshold: the fitted linear predictors, put through the link's inverse, saturate for at least one observation. Nothing short of saturation is rejected. A fitted model's predict() cannot show this on its own, because the inverse link glm uses is bounded away from zero and one, so a separated fit otherwise yields finite weights and a standard error that looks ordinary. For a categorical exposure, which runs on the M-estimation path alone, the threshold is narrower: only a probability of exactly zero at the level a unit was actually assigned is rejected, since that alone leaves the unit's own weight undefined, and a softmax column for a level the unit was not assigned may reach zero without the fit being refused. A continuous exposure has no saturating inverse link and is not checked.

The weight functions apply a stricter rule to the propensity scores they are handed: a categorical matrix holding an exact 0 or 1 anywhere is refused, whichever level the cell belongs to. A separated nnet::multinom() fit is therefore stopped where its weights are built rather than here. ps_trunc() repairs such a matrix, its categorical matrix method reading the closed interval so that a cell at an endpoint is bounded rather than refused, while ps_trim() holds the open interval as the weight functions do. Weights built from a trimmed or truncated score are their own refusal here, as are weights bounded with wt_trunc(), since ipw() cannot yet account for a modified propensity score or a truncated weight, so a separated fit is repaired for the weight functions rather than for this one. See Propensity scores at 0 and 1 in wt_ate() for the remedy.

The propensity score model must also be fit without case weights, since the stacked propensity score equations are unweighted and a weighted fit would not sit at the score root; that requirement applies to se_method = "mestimation" alone, because neither the linearization path nor the robust one restacks the propensity score model. Linearization still corrects for the uncertainty of estimating the propensity scores; it does so through the influence functions rather than through a stacked score. The robust diagnostic drops that correction, which is what makes it a diagnostic.

References

Stefanski LA, Boos DD. The calculus of M-estimation. The American Statistician. 2002;56(1):29–38. doi:10.1198/000313002753631330

Kostouraki A, Hajage D, Rachet B, et al. On variance estimation of the inverse probability-of-treatment weighting estimator: A tutorial for different types of propensity score weights. Statistics in Medicine. 2024;43(13):2672–2694. doi:10.1002/sim.10078

Lunceford JK, Davidian M. Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study. Statistics in Medicine. 2004;23(19):2937–2960. doi:10.1002/sim.1903

Shu D, Young JG, Toh S, Wang R. Variance estimation in inverse probability weighted Cox models. Biometrics. 2020. doi:10.1111/biom.13332

Robins JM, Hernán MA, Brumback B. Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;11(5):550–560. doi:10.1097/00001648-200009000-00011

Cole SR, Hernán MA. Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology. 2008;168(6):656–664. doi:10.1093/aje/kwn164

Hernán MA, Robins JM. Causal Inference: What If. Boca Raton: Chapman & Hall/CRC; 2020.

Leyrat C, Seaman SR, White IR, et al. Propensity score analysis with partially observed covariates: How should multiple imputation be used? Statistical Methods in Medical Research. 2019;28(1):3–19. doi:10.1177/0962280217713032

See Also

wt_ate(), wt_att(), wt_atu(), wt_atm(), wt_ato(), wt_entropy() for calculating propensity score weights.

ps_trim(), ps_trunc() for handling extreme propensity scores before weighting.

Examples

# Simulate data with a confounder, binary exposure, and binary outcome
set.seed(123)
n <- 200
x1 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.5 * x1))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + 0.3 * x1))
dat <- data.frame(x1, z, y)

# Step 1: Fit a propensity score model
ps_mod <- glm(z ~ x1, data = dat, family = binomial())

# Step 2: Calculate ATE weights and fit a weighted outcome model
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = quasibinomial(), weights = wts)

# Step 3: Estimate causal effects with correct standard errors
result <- ipw(ps_mod, outcome_mod)
result

# Exponentiate log-RR and log-OR to get RR and OR
as.data.frame(result, exponentiate = TRUE)

# Continuous outcome example
y_cont <- 2 + 0.8 * z + 0.3 * x1 + rnorm(n)
dat$y_cont <- y_cont
outcome_cont <- lm(y_cont ~ z, data = dat, weights = wts)
ipw(ps_mod, outcome_cont)

# Effect modification: the effect within each level of a factor, and the
# difference between them. The outcome model carries the exposure-by-modifier
# term that lets the effect differ across the strata.
dat$grp <- factor(rep(c("a", "b"), length.out = n))
ps_grp <- glm(z ~ x1 + grp, data = dat, family = binomial())
wts_grp <- wt_ate(ps_grp)
outcome_grp <- glm(
  y ~ z * grp + x1,
  data = dat,
  family = quasibinomial(),
  weights = wts_grp
)
ipw(ps_grp, outcome_grp, .by = grp)

# The same fit with the weights stabilized on the modifier. The numerator is
# the probability of the exposure each unit took given `grp`, which is the
# shape the default stabilizer has and conditions only on a variable the
# outcome model reads, so it tightens the weights without moving the estimand.
num_grp <- glm(z ~ grp, data = dat, family = binomial())
p_grp <- fitted(num_grp)
wts_stab <- wt_ate(
  ps_grp,
  stabilize = TRUE,
  stabilization_score = ifelse(dat$z == 1, p_grp, 1 - p_grp)
)
sd(wts_grp)
sd(wts_stab)
outcome_stab <- glm(
  y ~ z * grp + x1,
  data = dat,
  family = quasibinomial(),
  weights = wts_stab
)
ipw(ps_grp, outcome_stab, .by = grp)

# Continuous exposure: an lm propensity model of the dose on covariates,
# weights, which are stabilized by default for a continuous exposure, and a
# weighted marginal structural outcome model
a <- 0.5 + 0.8 * x1 + rnorm(n)
y_dose <- 1 + 0.6 * a + 0.3 * x1 + rnorm(n)
dat$a <- a
dat$y_dose <- y_dose
ps_cont <- lm(a ~ x1, data = dat)
wts_cont <- wt_ate(ps_cont)
msm <- lm(y_dose ~ a, data = dat, weights = wts_cont)
ipw(ps_cont, msm)

# A dose-response curve written as a basis matrix reports one row per basis
# coefficient. The model frame of such a fit holds the basis rather than the
# dose, so `.data` supplies the exposure.
msm_curve <- lm(y_dose ~ poly(a, 2), data = dat, weights = wts_cont)
ipw(ps_cont, msm_curve, .data = dat)


# Categorical exposure: a multinomial propensity model and per-level contrasts
set.seed(2)
n <- 300
x1 <- rnorm(n)
a <- factor(sample(c("a", "b", "c"), n, replace = TRUE))
y <- rbinom(n, 1, plogis(-0.5 + 0.4 * (a == "b") + 0.8 * (a == "c") + 0.3 * x1))
dat_cat <- data.frame(x1, a, y)

ps_cat <- nnet::multinom(a ~ x1, data = dat_cat, trace = FALSE)
ps_mat <- predict(ps_cat, type = "probs")
wts_cat <- wt_ate(ps_mat, dat_cat$a, exposure_type = "categorical")
outcome_cat <- glm(y ~ a, data = dat_cat, family = quasibinomial(), weights = wts_cat)
ipw(ps_cat, outcome_cat)

Description

is_ps_calibrated() tests whether x was calibrated, rather than whether it is still a ps_calib object: it is TRUE for a calibrated propensity score, for weights derived from one, and for scores that were calibrated before ps_trim() or ps_trunc() bounded them, each of which records the calibration in its own metadata.

Usage

is_ps_calibrated(x)

Arguments

Value

A single TRUE or FALSE.

See Also

ps_calibrate() to calibrate propensity scores.

Examples

ps <- runif(100)
exposure <- rbinom(100, 1, ps)

is_ps_calibrated(ps)

calibrated <- ps_calibrate(ps, exposure, smooth = FALSE)
is_ps_calibrated(calibrated)

Description

is_ps_trimmed() returns TRUE if x is a ps_trim object or a psw object created from trimmed propensity scores, and FALSE otherwise. This tests whether the object carries trimming information, not which individual units were trimmed; see is_unit_trimmed() for that.

Usage

is_ps_trimmed(x)

Arguments

Value

A logical scalar (TRUE or FALSE).

See Also

ps_trim() for trimming propensity scores, is_unit_trimmed() to identify which units were trimmed, ps_trim_meta() to retrieve full trimming metadata.

Examples

ps <- c(0.05, 0.3, 0.6, 0.95)
trimmed <- ps_trim(ps, method = "ps", lower = 0.1, upper = 0.9)

is_ps_trimmed(trimmed)
is_ps_trimmed(ps)

Description

is_ps_truncated() returns TRUE if x is a ps_trunc object or a psw object derived from truncated propensity scores. Use is_unit_truncated() to find out which observations were modified.

Usage

is_ps_truncated(x)

Arguments

Value

A single TRUE or FALSE.

See Also

ps_trunc(), is_unit_truncated(), ps_trunc_meta()

Examples

ps <- c(0.02, 0.3, 0.5, 0.7, 0.98)
is_ps_truncated(ps)

ps_t <- ps_trunc(ps, method = "ps", lower = 0.05, upper = 0.95)
is_ps_truncated(ps_t)

Description

is_refit() tests whether x is a ps_trim object whose propensity model has been refit on the retained (non-trimmed) observations via ps_refit().

Usage

is_refit(x)

Arguments

Details

On a psw vector built from trimmed propensity scores, the answer comes from the trimming record carried on the weights. is_refit() reads a single flag out of that record rather than a position, so unlike is_unit_trimmed() it answers from any record present, whatever length the weights have since taken on. It raises an error of class propensity_missing_meta_error only when weights marked as trimmed carry no record at all, which is what a subset of such weights leaves behind. Weights that were never trimmed have nothing missing and return FALSE.

Value

A single TRUE or FALSE.

See Also

ps_refit() to refit a propensity model after trimming, ps_trim() to trim propensity scores.

Examples

set.seed(2)
n <- 30
x <- rnorm(n)
z <- rbinom(n, 1, plogis(0.4 * x))
fit <- glm(z ~ x, family = binomial)
ps <- predict(fit, type = "response")

trimmed <- ps_trim(ps, lower = 0.2, upper = 0.8)
is_refit(trimmed)

refit <- ps_refit(trimmed, fit)
is_refit(refit)

Description

is_unit_trimmed() returns a logical vector indicating which observations were removed by trimming. This is a per-unit query, as opposed to is_ps_trimmed(), which tests whether the object has been trimmed at all.

The answer comes from the trimming record, which is written for a fixed set of observations and can both be lost and outlive them. On a ps_trim, vctrs::vec_slice() and c() drop it, and subassignment that grows the vector carries it across a length change; see ps_trim() for the whole contract. On a psw vector built from trimmed propensity scores, a subset drops it, while model.frame() re-attaches it to the shortened weights column of an outcome model fit on these weights.

is_unit_trimmed() therefore checks that the record covers the object it is given, and raises an error of class propensity_missing_meta_error when it does not, or when an object marked as trimmed carries no record at all, rather than name trimmed units at stale positions. Query the ps_trim object the record was written for instead.

That check counts observations, which a reordering does not change, so a ps_trim or psw reordered through vctrs rather than through [, by vctrs::vec_slice(x, 5:1) or dplyr::arrange(), drops the positions instead, and is_unit_trimmed() refuses the result. See ps_trim() and psw for the whole contract.

Usage

is_unit_trimmed(x)

Arguments

Value

A logical vector the same length as x, where TRUE marks a trimmed unit.

See Also

ps_trim() for trimming propensity scores, is_ps_trimmed() to test whether an object has been trimmed, ps_trim_meta() to retrieve full trimming metadata.

Examples

ps <- c(0.05, 0.3, 0.6, 0.95)
trimmed <- ps_trim(ps, method = "ps", lower = 0.1, upper = 0.9)

is_unit_trimmed(trimmed)

# Use to subset data to retained observations
kept <- !is_unit_trimmed(trimmed)
ps[kept]

Description

is_unit_truncated() returns a logical vector indicating which observations had their propensity scores modified by truncation. Use is_ps_truncated() to test whether an object has been truncated at all.

The answer comes from the truncation record, which is written for a fixed set of observations and can both be lost and outlive them. On a ps_trunc, vctrs::vec_slice() and c() drop it, and subassignment that grows the vector carries it across a length change; see ps_trunc() for the whole contract. On a psw vector built from truncated propensity scores, a subset drops it, while subassignment that grows the weights carries it across the length change.

is_unit_truncated() therefore checks that the record covers the object it is given, and raises an error of class propensity_missing_meta_error when it does not, or when an object marked as truncated carries no record at all, rather than name truncated units at stale positions. Query the ps_trunc object the record was written for instead.

That check counts observations, which a reordering does not change, so a ps_trunc or psw reordered through vctrs rather than through [, such as by vctrs::vec_slice(x, 5:1) or dplyr::arrange(), drops the positions instead, and is_unit_truncated() refuses the result. See ps_trunc() and psw for the whole contract.

Usage

is_unit_truncated(x)

Arguments

Value

A logical vector the same length as x (or number of rows for matrix input). TRUE marks observations whose values were winsorized.

See Also

ps_trunc(), is_ps_truncated(), ps_trunc_meta()

Examples

ps <- c(0.02, 0.3, 0.5, 0.7, 0.98)
ps_t <- ps_trunc(ps, method = "ps", lower = 0.05, upper = 0.95)
is_unit_truncated(ps_t)

Description

is_wt_truncated() returns TRUE for a psw vector whose own values were bounded, and is_unit_wt_truncated() returns which of its weights the bound moved.

Truncating the weights is recorded apart from truncating the propensity scores they were built from, which is_ps_truncated() and is_unit_truncated() report. The two are different operations, and a set of weights can carry either, both, or neither.

Usage

is_wt_truncated(x)

is_unit_wt_truncated(x)

Arguments

Details

The truncation leaves a record on the weights, the psw_trunc_meta attribute, holding the method, the bound as given and as applied, and the positions of the weights it moved. is_unit_wt_truncated() answers from those positions, so the record follows the rules psw describes for the records a modified propensity score leaves: it is kept through arithmetic and subassignment, re-indexed through the subscript of [, and loses its positions to any other slice and to a combine, which keeps the method and bounds. The wt_truncated flag describes the weights as a whole and is kept through all of those, so is_wt_truncated() keeps its answer where is_unit_wt_truncated() has none to give.

is_unit_wt_truncated() therefore checks that the record covers the weights it is given, and raises an error of class propensity_missing_meta_error when it does not, or when weights marked as truncated carry no record at all, rather than name truncated weights at stale positions. Weights that were not truncated have no record to read, and every weight is reported as untouched.

Value

• is_wt_truncated(): a single TRUE or FALSE. It is FALSE for anything that is not a psw vector.

• is_unit_wt_truncated(): a logical vector the same length as x, TRUE for each weight the bound moved.

See Also

psw for the weight vector class, and is_ps_truncated() for truncation of the propensity scores.

Examples

w <- psw(c(1.2, 0.8, 1.5), estimand = "ate")
is_wt_truncated(w)
is_unit_wt_truncated(w)

Description

joint_wt_models() records the two fitted treatment models a joint weight is built from, in the order the joint density factorizes: the first treatment's model, then the second treatment's model, which conditions on the first. is_joint_wt_models() reports whether an object is such a record.

The container exists to refuse. It checks that the pair describes a sequential factorization at all, which is the one thing about a joint weight that cannot be detected from the weight vector afterwards.

Usage

joint_wt_models(...)

is_joint_wt_models(x)

Arguments

Value

joint_wt_models() returns an S3 object of class joint_wt_models, a list with three fields:

names

The two treatment names, in the factorization's order.

models

The two fitted models, named for their treatments.

exposure_type

Each model's exposure type, named for its treatment: "binary", "categorical", or "continuous". The type follows the model's class, so a nnet::multinom() is recorded as "categorical" whatever its level count. One fit to two levels is a binary treatment model, and ipw() stacks it as one.

is_joint_wt_models() returns a single logical.

What is checked

Each argument is named for the treatment its model fits, and the name has to be that model's response: the factorization check looks for the first treatment by name in the second model's formula, so a model recorded under the wrong name would be checked for the wrong variable.

The second model must condition on the first treatment. The product of two marginal models, $f(A | L) f(E | L)$, is not the joint weight $f(A | L) f(E | A, L)$ for any data in which the second treatment depends on the first, and nothing downstream can tell the two apart. What is checked is whether the second model's right-hand side mentions the first treatment, so e ~ factor(a) * x1 and e ~ x1 + a:x1 both pass, as does the additive e ~ a + x1.

Passing is not the same as modeling the dependence well. An additive term satisfies the factorization and may still model it badly, and a badly modeled dependence biases the joint weight invisibly. Model the dependence on the first treatment flexibly, with an interaction or a spline, unless you have a reason to believe it is additive.

Neither may the two models condition on each other. A sequential factorization has a first factor that is marginal in the second treatment, so a pair that each read the other's treatment are not the two factors of any factorization, in either order.

See Also

wt_joint() for the product weight itself. These are used by ipw() methods for joint treatment models.

Examples

set.seed(4)
n <- 200
x1 <- rnorm(n)
a <- rbinom(n, 1, plogis(0.3 * x1))
e <- rbinom(n, 1, plogis(-0.2 + 0.5 * x1 - 0.8 * a))
dat <- data.frame(x1, a, e)

models <- joint_wt_models(
  a = glm(a ~ x1, data = dat, family = binomial()),
  e = glm(e ~ a * x1, data = dat, family = binomial())
)
models$names
models$exposure_type
is_joint_wt_models(models)

Description

numerator_model() returns the fitted model supplied to stabilize that estimated the numerator of a set of weights, and NULL for weights whose numerator was not estimated by one.

Usage

numerator_model(wt)

Arguments

Details

A model supplied to wt_ate()'s stabilize argument estimates the numerator of the weights, and what it estimates follows the exposure. A binomial() fit of a binary exposure reports the probability of the level each unit took; an nnet::multinom() fit of a categorical exposure reports one probability per level, and the numerator is the column named for the level each unit took; a model of a dose reports the conditional mean the density family is read at, and the numerator is that density. The model itself is recorded rather than the numerator it evaluates to, because ipw() rebuilds that numerator at every value of its parameter vector, which takes the model's design and its coefficients.

Where the record is kept differs by exposure type, and this accessor is the one place to read it either way. Weights for a continuous exposure are a ratio of densities and keep the model inside the density record density_meta() returns; weights for a binary or categorical exposure are a ratio of probabilities, leave no density record, and keep the model on its own. Both come back here.

The record describes the numerator rather than the units, so it holds no length of its own and survives subsetting, arithmetic, and anything else that changes the number of weights. Two sets of weights record the same numerator when their models write the same formula and were fit to the same coefficients; combining two that record different ones drops the record, with a warning of class propensity_metadata_conflict_warning. Weights stabilized without a model were divided by a numerator no model estimated, which is a different numerator again: combining those with weights stabilized on a model drops the model the same way, and the result is stabilized on nothing it can name.

A stabilization_score is a numerator the caller computed rather than one a model estimated, so weights stabilized on a score record no model, and stabilization_score() returns NULL for weights stabilized on a model.

Value

The fitted model supplied to stabilize, or NULL.

See Also

wt_ate() for the stabilize argument that writes this record, stabilization_score() for the numerator a caller supplies as a vector, and density_meta() for the rest of what a continuous exposure's weights record.

Examples

set.seed(1)
x <- rnorm(50)
v <- rbinom(50, 1, 0.5)
z <- rbinom(50, 1, plogis(0.4 * x - 0.7 * v))

ps_mod <- glm(z ~ x + v, family = binomial())
num_mod <- glm(z ~ v, family = binomial())

w <- wt_ate(ps_mod, stabilize = num_mod)
numerator_model(w)

# Weights stabilized on the marginal probability record no model.
numerator_model(wt_ate(ps_mod, stabilize = TRUE))

Description

se_method = "robust" reports the sandwich the weighted outcome model computes for itself, which treats the estimated weights as known. Printing such a result writes what print() writes for any ipw() result and then one line naming the method, so that the number beside each estimate is not read as one accounting for the propensity score model. See Standard errors as a diagnostic in ipw() for what the diagnostic does and does not cover.

Usage

## S3 method for class 'ipw_diagnostic_se'
print(x, ...)

## S3 method for class 'ipw_diagnostic_se'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'ipw_diagnostic_se'
tidy(x, ...)

Arguments

Value

x, invisibly.

as.data.frame() returns what it returns for any ipw() result, carrying an ipw_se_diagnostic attribute naming the method.

tidy() returns what it returns for any ipw() result, carrying an ipw_se_diagnostic attribute naming the method.

See Also

ipw() for the estimator and the standard error methods.

Examples

set.seed(2)
n <- 300
x <- rnorm(n)
z <- rbinom(n, 1, plogis(0.4 * x))
y <- rbinom(n, 1, plogis(-0.3 + 0.7 * z + 0.5 * x))
dat <- data.frame(x, z, y)

ps_mod <- glm(z ~ x, data = dat, family = binomial())
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = quasibinomial(), weights = wts)

print(ipw(ps_mod, outcome_mod, se_method = "robust"))

Description

ps_calibrate() adjusts estimated propensity scores so they better reflect true treatment probabilities. This can improve the accuracy of inverse probability weights derived from a misspecified propensity score model.

Usage

ps_calibrate(
  .propensity,
  .exposure,
  method = c("logistic", "isoreg"),
  smooth = TRUE,
  .focal_level = NULL,
  .reference_level = NULL,
  ...,
  .treated = NULL,
  .untreated = NULL,
  ps = lifecycle::deprecated()
)

Arguments

Details

Calibration is useful when the propensity score model is correctly specified in terms of variable selection but produces probabilities that are systematically too high or too low. Unlike ps_trim() and ps_trunc(), which handle extreme scores by removing or bounding them, calibration reshapes the entire distribution of scores.

Choosing a method:

• Use "logistic" (the default) as a first choice. It is stable and fast, and the smooth = TRUE option adds flexibility via a spline.

• Use "isoreg" when you suspect a non-smooth or irregular relationship between estimated and true probabilities and have a sufficiently large sample.

Fitted models: .propensity can be the fitted propensity score model instead of the scores it reports, and the exposure is then read off the model as well unless .exposure is supplied. The variable read is announced; options(propensity.quiet = TRUE) silences the announcement. A binomial stats::glm() and a two-level nnet::multinom() report the one score per unit calibration reads. A nnet::multinom() of three or more levels reports one column per level and no single score to calibrate, so it is refused with an error of class propensity_model_family_error; calibrate the columns one at a time against the exposure indicator each of them belongs to.

Missing values: A unit with a missing exposure tells the calibration model nothing, so it is dropped from the fit. What happens to that unit afterwards depends on the method. Logistic calibration, smoothed or not, fits a curve from propensity score to calibrated score and reads every unit off it, so a unit with a missing exposure but an observed score is still calibrated. Isotonic calibration fits separately within each exposure group and reads each unit from the fit for the group it belongs to; a unit with a missing exposure belongs to neither group and is returned as NA. Under both methods a missing propensity score yields a missing calibrated score.

The calibrated scores are returned as a ps_calib object, which can be passed directly to weight functions such as wt_ate().

Value

A ps_calib vector the same length as .propensity. The attribute ps_calib_meta stores calibration metadata (method and whether smoothing was applied). Use is_ps_calibrated() to test whether an object has been calibrated.

References

Platt, J. (1999). Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers, 61–74.

Zadrozny, B., & Elkan, C. (2002). Transforming classifier scores into accurate multiclass probability estimates. Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 694–699. doi:10.1145/775047.775151

See Also

is_ps_calibrated() to test for calibrated scores; ps_trim() and ps_trunc() for alternative approaches to extreme propensity scores; wt_ate() and other weight functions that accept ps_calib objects.

Examples

# Simulate data
set.seed(42)
ps <- runif(200)
exposure <- rbinom(200, 1, ps)

# Logistic calibration without smoothing (simple Platt scaling)
cal <- ps_calibrate(ps, exposure, smooth = FALSE)
cal

# Use calibrated scores to calculate weights
wt_ate(cal, exposure)

# Isotonic regression calibration
cal_iso <- ps_calibrate(ps, exposure, method = "isoreg")

if (rlang::is_installed("mgcv")) {
  # Logistic calibration with spline smoothing (default)
  cal_smooth <- ps_calibrate(ps, exposure)
}

# Calibrate the scores a fitted model reports, reading the exposure off it
x <- rnorm(200)
fit <- glm(exposure ~ x, family = binomial())
ps_calibrate(fit, smooth = FALSE)

if (rlang::is_installed("nnet")) {
  exposure_factor <- factor(exposure)
  multinomial_fit <- nnet::multinom(exposure_factor ~ x, trace = FALSE)
  ps_calibrate(multinomial_fit, smooth = FALSE)
}

Description

Re-estimates a propensity score model using only the observations retained after trimming. This is the recommended intermediate step between ps_trim() and weight calculation (e.g. wt_ate()):

ps_trim() -> ps_refit() -> ⁠wt_*()⁠

Trimming changes the target population by removing observations with extreme propensity scores. Refitting the model on the retained subset produces propensity scores that better reflect this population, improving both model fit and downstream weight estimation. Weight functions warn if a trimmed propensity score has not been refit.

Usage

ps_refit(trimmed_ps, model, .data = NULL, ...)

Arguments

Details

Composing with a subset

A subset in the original call has already chosen the sample the propensity scores are about, and the trimming record indexes that sample rather than every row the data carry. Refitting narrows that sample further, to the retained rows, so the original subset is dropped from the call rather than put to work a second time on rows it was never about. A subset passed through ... is an instruction of its own and is honored.

Which level a refit predicts

For a vector of scores, the refit predicts the probability of the level the scores describe: when ps_trim() was given a model with its first level named as focal, the refit reports one minus the model's prediction, and for scores supplied directly it reports the level the model predicts by default.

Arguments read from outside the formula

weights, offset, and na.action in the original call are re-evaluated against the retained rows. A weights or offset naming a column of the data the model was fit on is read from that column and follows the retained rows, whether the data are recovered from model or passed to .data. A vector held outside the data cannot follow them: it keeps the length it had and raises an error about differing variable lengths.

Scores predicted from a fit with na.action = na.exclude are padded back to the full length of the data, so they describe more observations than the fit read and the trimming record indexes a sample the model never analyzed. ps_refit() refuses such scores. Trim scores from a fit whose na.action drops those rows instead.

Value

A ps_trim object with re-estimated propensity scores for retained observations and NA for trimmed observations. Use is_refit() to confirm refitting was applied. Refitting replaces the values a calibrated score held with predictions from the refit model, so a score trimmed after ps_calibrate() is no longer calibrated and is_ps_calibrated() answers FALSE for the result.

For a trimmed dose model, the retained values are the refit model's conditional means, and the spread in the record is re-estimated from the retained residuals under the recorded density family, unless the trim was made at a .sigma the caller supplied, which is kept. The recorded threshold describes the cut that was made and is left as it was, so the refit model can place a retained unit's conditional density below it; nothing is trimmed again.

See Also

ps_trim() for the trimming step, is_refit() to check refit status, wt_ate() and other weight functions for the next step in the pipeline.

Examples

set.seed(2)
n <- 200
x <- rnorm(n)
z <- rbinom(n, 1, plogis(0.4 * x))

# fit a propensity score model
ps_model <- glm(z ~ x, family = binomial)
ps <- predict(ps_model, type = "response")

# trim -> refit -> weight pipeline
trimmed <- ps_trim(ps, lower = 0.1, upper = 0.9)
refit <- ps_refit(trimmed, ps_model)
wts <- wt_ate(refit, .exposure = z)

is_refit(refit)

Description

Every estimand this package supports targets a population whose covariate distribution is the study population reweighted by a tilting function $h$ of the propensity score. ps_tilt() evaluates $h$ at each observation's propensity score.

The tilt is the numerator of every weight: a weight is $h$ divided by the propensity score of the exposure level the unit actually received. It is also what standardizes a g-computation estimate to a target population, which for "atm", "ato", and "entropy" is the only route to the estimand, since those populations are not subsets of the rows and cannot be reached by filtering.

Usage

ps_tilt(
  .propensity,
  estimand,
  ...,
  .focal_level = NULL,
  ps = lifecycle::deprecated()
)

Arguments

Value

A plain double vector, unnamed, with one element per observation: the length of .propensity for the numeric method, and the number of rows of .propensity for the matrix and data frame methods.

Tilting functions

For a binary exposure with propensity score $e = P(Z = \text{focal} \mid X)$:

For a categorical exposure with propensity score vector $(e_1, \ldots, e_K)$ and focal level $f$: "ate" is $1$, "att" is $e_f$, "atu" is $1 - e_f$, "atm" is $\min_k e_k$, "ato" is $(\sum_k 1 / e_k)^{-1}$, and "entropy" is $-\sum_k e_k \log(e_k)$.

Stabilization and censoring weights are not tilts. Stabilization multiplies a weight by a marginal quantity that does not depend on the covariates, and wt_cens() reuses the "ate" formula, so neither has an entry here.

Standardizing model predictions

A g-computation estimate averages an outcome model's per-row counterfactual predictions, and where those predictions vary from row to row it is the weights of that average that standardize the estimate to a target population. The tilt is those weights, which is what the second example below supplies by hand. An average taken with equal weight standardizes a covariate-adjusted outcome model to the whole sample, so such a model targets the ATE however its own fitting weights were built: fit for a non-"ate" estimand and averaged that way, it reports the full-sample contrast rather than the one it was weighted for. Tools that average per-row model predictions, such as the avg_comparisons() function in the marginaleffects package, average with equal weight unless they are told otherwise, and the remedy is to hand them the tilt as the weights of the average: ps_tilt(ps, "att") for a binary exposure and ps_tilt(ps, "att", .focal_level = "b") for a categorical one.

The requirement is easy to miss because an outcome model saturated in the exposure hides it. Such a model predicts one value per exposure level, so every average of its predictions returns the same contrast, and there the estimand is settled by the weights the model was fit with rather than by the weights of the average: a model fit with wt_att() weights reports the ATT whether the average is taken with equal weight or with the tilt. The two averages come apart at the first covariate the outcome model adjusts for, where the predictions vary from row to row and the weights the average is taken under are what decides the population the estimate describes.

Propensity score range

Every propensity score in .propensity must lie strictly inside $(0, 1)$, the same requirement wt_ate() and the rest of the weight family impose before they tilt. The bound holds for a matrix or data frame .propensity entry by entry, and each row must sum to one on top of it. A fitted model that separates the exposure can return a probability of exactly zero or one; those scores have no weight to divide and are rejected here rather than tilted.

A missing propensity score gives a missing tilt under every estimand, "ate" included, so an observation whose propensity score is unknown never counts toward a tilted mean. A numeric .propensity propagates NA position by position, and a matrix .propensity gives NA for any row holding one: a probability vector with a missing entry is not one this can tilt on, whichever level the tilt reads.

Fitted models

.propensity can be the fitted propensity score model instead of the scores it reports. A binomial stats::glm() and a two-level nnet::multinom() report one score per unit and are tilted as a vector; a nnet::multinom() of three or more levels reports one column per level and is tilted as a matrix, with .focal_level naming the column "att" and "atu" read. Those are the shapes predict(fit, type = "response") and fitted(fit) give, and tilting a fit tilts exactly what tilting those values would. A tilt reads no exposure, so nothing but the scores is taken off the model.

Modified propensity scores

ps_tilt() takes plain propensity scores. A score modified by ps_trim(), ps_trunc(), or ps_calibrate() carries a class of its own and has no method here; pass the scores underneath it, with as.numeric(x) for a binary exposure or as.matrix(x) for a categorical one. Units ps_trim() set to NA stay NA through that extraction and take an NA tilt.

See Also

wt_ate() and the rest of the weight family, which divide the tilt by the propensity score of the received exposure level.

Examples

set.seed(1)
n <- 500
x <- rnorm(n)
z <- rbinom(n, 1, plogis(0.6 * x))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + x))
sim <- data.frame(x = x, z = z, y = y)
ps <- unname(
  predict(glm(z ~ x, data = sim, family = binomial()), type = "response")
)

# a weight is the tilt over the propensity of the received exposure level
received <- z * ps + (1 - z) * (1 - ps)
all.equal(
  as.numeric(wt_ato(ps, z, exposure_type = "binary")),
  ps_tilt(ps, "ato") / received
)

# tilted g-computation: standardize counterfactual predictions to the
# overlap population with an h-weighted mean, no weights in the outcome model
fit <- glm(y ~ z * x, data = sim, family = binomial())
m1 <- predict(fit, transform(sim, z = 1), type = "response")
m0 <- predict(fit, transform(sim, z = 0), type = "response")

h <- ps_tilt(ps, "ato")
weighted.mean(m1, h) - weighted.mean(m0, h)

# the same predictions standardized to everyone give the ATE instead
mean(m1) - mean(m0)

# the tilt of the scores a fitted model reports
ps_tilt(glm(z ~ x, data = sim, family = binomial()), "ato")

if (rlang::is_installed("nnet")) {
  sim$trt <- factor(sample(c("a", "b", "c"), n, replace = TRUE))
  multinomial_fit <- nnet::multinom(trt ~ x, data = sim, trace = FALSE)
  head(ps_tilt(multinomial_fit, "att", .focal_level = "b"))
}

Description

Trim observations with extreme propensity scores by replacing them with NA, effectively removing those units from downstream analyses. The returned object has the same length (or dimensions) as the input, with trimmed entries set to NA. After trimming, refit the propensity score model on the retained observations with ps_refit().

Usage

ps_trim(
  .propensity,
  method = c("ps", "adaptive", "pctl", "pref", "cr", "optimal", "density", "resid"),
  lower = NULL,
  upper = NULL,
  .exposure = NULL,
  .focal_level = NULL,
  .reference_level = NULL,
  ...,
  .sigma = NULL,
  .density = "normal",
  .treated = NULL,
  .untreated = NULL,
  ps = lifecycle::deprecated()
)

Arguments

Details

How trimming works

Trimming identifies observations with extreme (near 0 or 1) propensity scores and sets them to NA. These observations are excluded from subsequent weight calculations and effect estimation. The goal is to remove units that lack sufficient overlap between exposure groups, which would otherwise receive extreme weights and destabilize estimates.

Choosing a method

• Use "ps" when you have a specific threshold in mind or want a simple default.

• Use "adaptive" for a principled, data-driven cutoff that targets variance reduction.

• Use "pctl" to trim a fixed percentage of extreme values from each tail.

• Use "pref" when you want to restrict to the region of clinical equipoise based on the preference score.

• Use "cr" to restrict to the common support region where both exposure groups have observed propensity scores.

• Use "optimal" for multi-category (3+) exposures; this is the only data-driven method available for categorical treatments.

• Use "density" or "resid" for a continuous exposure, on its dose model.

For a binary or categorical exposure, trimming and refitting is the first recourse when the analysis can accept a change of estimand. When every unit must be kept, bound the scores with ps_trunc() instead, for a binary exposure with method = "adaptive", or bound the weights with wt_trunc().

Typical workflow

1. Fit a propensity score model

2. Apply ps_trim() to flag extreme values

3. Call ps_refit() to re-estimate propensity scores on the retained sample

4. Compute weights with wt_ate() or another weight function

Fitted models

.propensity can be the fitted propensity score model instead of the scores it reports. A binomial stats::glm() and a two-level nnet::multinom() are read as one score per unit; a nnet::multinom() of three or more levels is read as one column per level. Those are the shapes predict(fit, type = "response") and fitted(fit) give, and trimming a fit trims exactly what trimming those values would.

The methods that read an exposure ("pref", "cr", and every method on the categorical route) take it from the model when .exposure is not supplied, announcing the variable they read; options(propensity.quiet = TRUE) silences the announcement. An .exposure you supply is used instead, and a categorical model's columns are matched to its levels by name, so an exposure whose levels are ordered differently is still trimmed against the right column.

Trimming a dose model

For a continuous exposure, .propensity can be the model of the exposure's conditional mean: a stats::lm(), a glm of the gaussian() family (or another family whose variance is constant), a MASS::rlm(), or an mgcv::gam(), the models wt_ate() reads a dose from. A glm whose spread changes with its mean, such as poisson(), is refused with an error of class propensity_model_family_error. Only "density" and "resid" apply.

For a dose, the generalized propensity score is the conditional density $f(a \mid x)$ (Hirano and Imbens, 2004), and a unit whose observed dose has a small conditional density gets a large weight, as a unit with a propensity score near 0 does for a binary exposure. Both methods read the fitted conditional mean mu, the spread sigma of the residuals under the family in .density, the standardized residual z = (a - mu) / sigma, and the conditional density f = g(z) / sigma, where g is the family's standardized density. "density" keeps the units whose f is at least its lower quantile. "resid" keeps those whose ⁠|z|⁠ is at most upper; for the symmetric unimodal families it accepts, that is the floor g(upper) / sigma on f, which is the threshold it records. The trim is one-sided, at the low end of the density, since only a small density makes a large weight (Branson et al., 2024).

The retained values are the conditional means, and the trimmed ones are NA. A unit with a missing mean or dose takes no part in the trim. .focal_level, .reference_level, and their deprecated forms are refused with an error of class propensity_focal_level_error, since a dose has no levels.

A trim changes the estimand. The analysis describes the units whose observed dose was plausible under the model, not the full sample, and which units those are depends on the threshold. Branson et al. (2024) define the trimmed dose response at each dose among the units whose conditional density at that dose exceeds the threshold, a population that changes with the dose; this trim is its sample analogue at each unit's own observed dose rather than the same population. A dose-response curve fit to the trimmed sample is a curve for the units it keeps.

After trimming, ps_refit() refits the conditional mean on the retained rows and re-estimates the spread there under the recorded family, and wt_ate() builds the weights from the refit under the family and spread the record holds. ipw() refuses those weights, as it refuses every weight built from a trimmed score.

The other way to hold down extreme weights for a dose is wt_trunc(), which bounds the weights and keeps every unit; see Truncating weights or trimming the density there for the simulation results behind this summary. With a density family that fits the residuals, a density trim at lower = 0.01 did not help: at an effective sample size of 64% of the sample, its interval covered 0.843, against 0.933 for the untrimmed weights under the same fixed-weight interval. A loose bound on the weights left that coverage where it was and lowered the root mean squared error. With heavy-tailed residuals read through a normal density, or a residual spread the model left out, the trim covered 0.979 to 0.989, against 0.618 to 0.881 for the untrimmed weights and 0.644 to 0.906 for the bound, because it removed the units the misspecified density fit worst. Two cautions go with those numbers. The trim's intervals were computed on the trimmed sample, held the weights fixed, and were conservative, and its bias was measured against the effect in the full population, so part of that bias is the change of estimand. "resid" was not simulated; the results carry over to it only because a bound on the standardized residual is a floor on the conditional density.

A family that fits the residuals is the better repair where it applies. Under heavy tails, dens_t() with 4 degrees of freedom at its default scale had about a third of the root mean squared error of the bounded normal weights and an M-estimation interval that covered 0.926 to 0.939 with standard errors close to right, and it kept every unit. The trim's coverage was higher, but its interval was conservative. Check the family first; trim when the units at the low-density end are ones the analysis should not describe, and say so when reporting the estimand.

Object behavior

Arithmetic operations on ps_trim objects return plain numeric vectors, since transformed propensity scores (e.g., 1/ps) are no longer propensity scores. Trimmed values propagate as NA in calculations; use na.rm = TRUE where appropriate.

When combining ps_trim objects with c(), trimming parameters must match. Mismatched parameters trigger a warning and return a numeric vector.

Use ps_trim_meta() to inspect the trimming metadata, including the method, cutoffs, and which observations were retained or trimmed.

Missing values

A propensity score that arrives missing is not one this function removed, so it takes no part in the trimming record: it joins neither the retained nor the trimmed positions, and is_unit_trimmed() reports FALSE for it. The value propagates as NA into the result. For a matrix of categorical propensity scores, a row with a missing cell comes back exactly as it arrived, since there is no complete probability vector to place against a threshold.

The methods that read a cutoff off the scores read it off the scores they have. "adaptive", "pctl", "cr", and "optimal" all work their cutoffs out from the complete scores or rows, so the cutoff is the one the same call would produce with the missing observations dropped.

"pref" centers its preference scores on the proportion exposed across the whole sample, which is a fact about the exposure rather than about the scores, so a unit whose propensity score is missing still counts toward it and the cutoff is not the one the shorter call would produce. That unit's own preference score is missing, which leaves it outside both the retained and the trimmed positions like any other missing score.

A missing exposure is a different matter, and "pref" and "cr" refuse one with an error of class propensity_missing_value_error. Their cutoffs come from the exposure groups, and a unit that belongs to neither leaves them undefined. Remove or impute the missing exposure values first.

The trimming record

A ps_trim records which units were trimmed as positions among the observations it was written for, along with how many observations that was. Operations that hand this package the subscript re-index those positions onto the result: subsetting with [, sort(), rep(), and na.omit() all return a record written for what they return, and a subscript naming a position more than once reports that unit at every place it now holds. unique() does the same when it can, as described below.

Operations that change how many observations there are without supplying a subscript cannot re-index the record, and it is dropped rather than worked out from the values, since reading membership back from the NA pattern would report a propensity score that arrived missing as one this package removed. vctrs::vec_slice(), which is how filtering, joining, and grouped verbs in dplyr reach a column, is the usual route, and combining two or more vectors with c() is another, because concatenation appends one set of observations to another. The record is dropped without comment on every route: most of these length changes build vectors the caller never holds, such as the pieces a grouped verb slices a column into or the rows a tibble slices off to print, so a warning would mostly describe something that is not the result. The values, the class, and the method and its cutoffs are untouched.

A combine drops the positions even when it is handed a single vector. c() of one ps_trim returns it unchanged, record included, but vctrs::vec_c(x), dplyr::bind_rows(df), vctrs::vec_rbind(df), and an ungrouped dplyr::reframe() rebuild the column, so a later ps_refit() or is_unit_trimmed() on the result refuses it.

unique() keeps one element for each distinct value, or one row for each distinct row of a matrix of scores, and that element or row stands for every unit holding the same scores. A matrix comes back as a matrix of the same class with its column names. The record is re-indexed onto the result when all of the merged units share one status, and dropped otherwise: a trimmed score and one that arrived missing are both NA, and a trimmed row is NA throughout, so a vector or matrix holding both returns a single NA element or row that neither status describes.

as.data.frame() and tibble::as_tibble() turn a matrix of scores into a data frame whose columns are ps_trim vectors carrying the matrix's record, one row per unit, so the weight functions read the data frame as they read the matrix and a subset of its rows re-indexes the record. A data frame whose columns are not all trimmed with the same record is refused by the weight functions with an error of class propensity_matrix_type_error. A matrix without column names gives columns named V1, V2, and so on, which the weight functions then refuse because they name no exposure level, so name the columns after the exposure levels before converting.

A record can also outlive the observations it describes, because it travels by routes vctrs does not see: growing a ps_trim by subassignment carries it across the length change. is_unit_trimmed() and ps_refit() therefore check that the record covers the object they are given and raise an error of class propensity_missing_meta_error when it does not, rather than name trimmed units at stale positions.

That check compares how many observations the record was written for against how many the object holds, which a reordering does not change, so a route that could reorder the observations without saying how drops the positions instead, at any length: vctrs::vec_slice(), dplyr::arrange(), dplyr::filter(), and vctrs::vec_assign() and the helpers built on it return a ps_trim whose record keeps its method and cutoffs and names no units, and is_unit_trimmed() and ps_refit() refuse it. Subsetting with [, sort(), unique(), and rep() know where the units went and re-index the record, and ⁠[<-⁠ and ⁠is.na<-⁠ move no unit and keep it, so reorder with [, or put the propensity scores in the order you want before trimming them.

Casting a numeric vector into a ps_trim with vctrs::vec_cast() is a type operation and not a trimming. The result is described by the method and cutoffs of the target and records that none of the arriving values was trimmed, so it can hold scores outside those cutoffs, including 0 and 1. Call ps_trim() on the scores to trim them.

Value

A ps_trim object (a numeric vector with class "ps_trim", or a matrix with class "ps_trim_matrix"). Trimmed observations are NA. Metadata is stored in the "ps_trim_meta" attribute and can be accessed with ps_trim_meta(). Key fields include:

• method: the trimming method used

• keep_idx: integer indices of retained observations

• trimmed_idx: integer indices of trimmed (NA) observations

• n_obs: the number of observations those indices describe

• Method-specific fields such as cutoff (adaptive), q_lower/q_upper (pctl), cr_lower/cr_upper (cr), delta (categorical ps), or lambda (optimal)

• focal_inverted (vector scores only): TRUE when the scores are one minus the probability a fitted model reports, because the model was trimmed with its first level named as focal, and FALSE otherwise, including for every vector of scores supplied directly

A trim of a dose model records method, lower ("density", else NULL), upper ("resid", else NULL), threshold (the realized floor on the conditional density), sigma (the spread it was read at), sigma_kind ("pooled", "mle", or "supplied"), density (the density specification), keep_idx, trimmed_idx, and n_obs.

References

Branson, Z., Kennedy, E. H., Balakrishnan, S., & Wasserman, L. (2024). Causal effect estimation after propensity score trimming with continuous treatments. arXiv preprint arXiv:2309.00706.

Crump, R. K., Hotz, V. J., Imbens, G. W., & Mitnik, O. A. (2009). Dealing with limited overlap in estimation of average treatment effects. Biometrika, 96(1), 187–199.

Hirano, K., & Imbens, G. W. (2004). The propensity score with continuous treatments. In Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives (pp. 73–84).

Walker, A. M., Patrick, A. R., Lauer, M. S., et al. (2013). A tool for assessing the feasibility of comparative effectiveness research. Comparative Effectiveness Research, 3, 11–20.

Yang, S., Imbens, G. W., Cui, Z., Faries, D. E., & Kadziola, Z. (2016). Propensity score matching and subclassification in observational studies with multi-level treatments. Biometrics, 72(4), 1055–1065.

See Also

ps_trunc() for bounding (winsorizing) instead of discarding, ps_refit() to re-estimate propensity scores after trimming, ps_calibrate() for calibration-based adjustment, ps_trim_meta() to inspect trimming metadata, is_ps_trimmed() and is_unit_trimmed() for logical queries.

Examples

set.seed(2)
n <- 300
x <- rnorm(n)
z <- rbinom(n, 1, plogis(1.3 * x))
fit <- glm(z ~ x, family = binomial)
ps <- predict(fit, type = "response")

# Fixed threshold trimming (default)
trimmed <- ps_trim(ps, method = "ps", lower = 0.1, upper = 0.9)
trimmed

# How many observations were trimmed?
sum(is_unit_trimmed(trimmed))

# Data-driven adaptive trimming
ps_trim(ps, method = "adaptive")

# Quantile-based trimming at 5th and 95th percentiles
ps_trim(ps, method = "pctl")

# Refit after trimming, then compute weights
trimmed <- ps_trim(ps, method = "adaptive")
refitted <- ps_refit(trimmed, fit)
wt_ate(refitted, .exposure = z)

# Trim the scores a fitted model reports, reading the exposure off the model
ps_trim(fit, method = "cr")

# Trim a dose model on the scale of its conditional density
dose <- 1 + 0.5 * x + rt(n, df = 3)
dose_fit <- lm(dose ~ x)
ps_trim(dose_fit, method = "density", lower = 0.05)
ps_trim(dose_fit, method = "resid", upper = 3, .density = dens_t(4))

if (rlang::is_installed("nnet")) {
  trt <- factor(sample(c("a", "b", "c"), n, replace = TRUE))
  multinomial_fit <- nnet::multinom(trt ~ x, trace = FALSE)
  ps_trim(multinomial_fit, method = "optimal")
}

Description

ps_trim_meta() returns the metadata list attached to a ps_trim object by ps_trim().

Usage

ps_trim_meta(x)

Arguments

Value

A named list with elements:

method

Character string indicating the trimming method used.

keep_idx

Integer vector of retained observation indices.

trimmed_idx

Integer vector of trimmed observation indices.

n_obs

The number of observations those indices describe.

lower, upper

Numeric cutoffs, when applicable.

refit

Logical, TRUE if the model was refit via ps_refit().

focal_inverted

For a vector of scores, logical, TRUE if the scores are one minus the probability the fitted model reports, which ps_refit() reads to predict the same level.

Additional method-specific elements (e.g. cutoff, delta, lambda) may also be present.

keep_idx, trimmed_idx, and n_obs are absent when an operation that changed the number of observations dropped the record; see ps_trim().

See Also

ps_trim() for trimming propensity scores, is_ps_trimmed() and is_unit_trimmed() for predicate queries.

Examples

ps <- c(0.05, 0.3, 0.6, 0.95)
trimmed <- ps_trim(ps, method = "ps", lower = 0.1, upper = 0.9)

ps_trim_meta(trimmed)

Description

ps_trunc() bounds extreme propensity scores to fixed limits, replacing out-of-range values with the boundary value (a form of winsorizing). The result is a vector or matrix of the same length and dimensions as .propensity, with no observations removed. This contrasts with ps_trim(), which sets extreme values to NA (effectively removing those observations from analysis).

Usage

ps_trunc(
  .propensity,
  method = c("ps", "adaptive", "pctl", "cr"),
  lower = NULL,
  upper = NULL,
  .exposure = NULL,
  .focal_level = NULL,
  .reference_level = NULL,
  ...,
  .treated = NULL,
  .untreated = NULL,
  ps = lifecycle::deprecated()
)

Arguments

Details

Unlike ps_trim(), truncation preserves all observations. No NA values are introduced; out-of-range scores are replaced with the boundary value.

For binary exposures, each propensity score $e_i$ is bounded:

• If $e_i < l$, set $e_i = l$ (the lower bound).

• If $e_i > u$, set $e_i = u$ (the upper bound).

For categorical exposures, values below the threshold are set to the threshold and each row is renormalized to sum to 1.

The "ps" bounds

With method = "ps", both paths default to a floor of 0.1. For a binary exposure, each bound supplied must be a single number strictly between 0 and 1, and a bound supplied alone implies its mirror: lower alone bounds the scores at ⁠[lower, 1 - lower]⁠, and upper alone at ⁠[1 - upper, upper]⁠. With neither supplied, the bounds are ⁠[0.1, 0.9]⁠; with both supplied, both are used as written, so an asymmetric bound must be given in full. The mirror follows from the weights: an untreated unit's ATE weight is $1/(1 - e)$, so a floor on $e$ alone would bound the treated weights and leave the untreated ones unbounded. A lone lower of 0.5 or more, or a lone upper of 0.5 or less, meets or crosses its own mirror and is an error of class propensity_range_error. lower = 1/c gives the bounds method = "adaptive" computes for the same $c$.

For a categorical exposure, lower is the threshold $\delta$ (default 0.1) applied to every column, after which each row is renormalized, and upper is not read and is recorded as NA. With k exposure levels, a threshold of 1/k or larger cannot be met by every column of a row that sums to one and is an error, so an exposure with ten or more levels needs an explicit lower below 1/k.

The adaptive bound

method = "adaptive" bounds a binary propensity score at $[1/c, 1 - 1/c]$, where $c = \sqrt{n} \log(n) / 5$ is the weight bound of Gruber et al. (2022) and $n$ counts the scores that are present, as in wt_trunc(). At $n = 1000$, $c$ is 43.7 and the bounds are 0.023 and 0.977. The bound adapts to the sample size alone, loosening as $n$ grows, and leaves the estimand alone. ps_trim(method = "adaptive") instead adapts to the estimated scores, tightening under poor overlap, and changes the estimand to the population it keeps.

For unstabilized ATE weights the bound is the weight bound written on the score scale. A treated unit's weight is $1/e$, so the floor $1/c$ caps the treated weights at $c$; an untreated unit's weight is $1/(1 - e)$, so the ceiling $1 - 1/c$ caps the untreated weights at $c$. wt_ate() on the result therefore gives the weights wt_trunc(wt_ate(ps, .exposure), method = "adaptive") gives, with one difference, provided every unit's exposure is observed. (A unit with a present score and a missing exposure counts toward $n$ here but has no weight for wt_trunc() to count, so the two bounds then differ.) The difference is this: a treated unit above $1 - 1/c$ or an untreated unit below $1/c$, whose weight is already close to 1, is moved to that bound too, and its weight rises to $c / (c - 1)$. Every weight above $c / (c - 1)$ is the same under either route. A lower bound alone would not do the same work, since no floor on $e$ bounds an untreated unit's weight.

For stabilized ATE weights the equivalence does not carry over as it stands. Stabilizing multiplies each arm's weights by that arm's prevalence, so the treated weights are capped at $P(A = 1) c$ and the untreated weights at $P(A = 0) c$, which agree only when the two arms are equally common.

The bounds cross when $c < 2$, which holds for fewer than 15 scores, so "adaptive" needs at least 15 scores present and raises an error of class propensity_range_error otherwise. It is refused for categorical exposures with an error of class propensity_method_error: truncating a categorical score renormalizes each row, which moves scores the bound never reached, so the bound no longer caps the weights. A two-level nnet::multinom() fit is read as one binary score and accepted.

Arithmetic behavior: Arithmetic operations on ps_trunc objects return plain numeric vectors. Once propensity scores are transformed (e.g., into weights), the result is no longer a propensity score.

Combining behavior: Combining ps_trunc objects with c() requires matching truncation parameters. Mismatched parameters produce a warning and return a plain numeric vector.

Fitted models

.propensity can be the fitted propensity score model instead of the scores it reports. A binomial stats::glm() and a two-level nnet::multinom() are read as one score per unit; a nnet::multinom() of three or more levels is read as one column per level. Those are the shapes predict(fit, type = "response") and fitted(fit) give, and truncating a fit bounds exactly what bounding those values would. A model of a continuous exposure's conditional mean, such as a stats::lm() fit or a gaussian stats::glm(), has no propensity score to bound and raises an error of class propensity_model_family_error. A dose has two routes instead: trim the dose model with ps_trim(method = "density"), which sets aside the units whose dose is implausible under the model, or build the weights with wt_ate() and bound them with wt_trunc(), which keeps every unit. There is no density method here: a floor on a dose's conditional density is a bound on its unstabilized weights, the ones wt_ate(stabilize = FALSE) builds, which wt_trunc() already applies.

The methods that read an exposure ("cr", and every method on the categorical route) take it from the model when .exposure is not supplied, announcing the variable they read; options(propensity.quiet = TRUE) silences the announcement. An .exposure you supply is used instead, and a categorical model's columns are matched to its levels by name, so an exposure whose levels are ordered differently is still bounded against the right column.

Missing values

A propensity score that arrives missing is not one this function pinned to a bound, so it takes no part in the truncation record: it never joins truncated_idx, and is_unit_truncated() reports FALSE for it. The value stays missing in the result. For a matrix of categorical propensity scores, a row with a missing cell comes back missing throughout, because renormalizing a row divides by its sum and that sum is unknown.

The methods that read their bounds off the data read them off the data they have. "pctl" and "cr" work their bounds out from the complete scores, and the categorical "pctl" bounds come from the cells of the complete rows, so the bounds are the ones the same call would produce with the missing observations dropped.

A missing exposure is a different matter, and "cr" refuses one with an error of class propensity_missing_value_error. Its bounds come from the exposure groups, and a unit that belongs to neither leaves them undefined. Remove or impute the missing exposure values first.

The truncation record

A ps_trunc records which units were winsorized as positions among the observations it was written for, along with how many observations that was. Operations that hand this package the subscript re-index those positions onto the result: subsetting with [, sort(), and rep() all return a record written for what they return, and a subscript naming a position more than once reports that unit at every place it now holds. unique() does the same when it can, as described below.

Operations that change how many observations there are without supplying a subscript cannot re-index the record, and it is dropped rather than worked out from the values, since a score that arrived equal to a bound is indistinguishable from one this function pinned there. vctrs::vec_slice(), which is how filtering, joining, and grouped verbs in dplyr reach a column, is the usual route, and combining two or more vectors with c() is another, because concatenation appends one set of observations to another. The record is dropped without comment on every route, since most of these length changes build vectors the caller never holds, such as the pieces a grouped verb slices a column into. The values, the class, and the method and its bounds are untouched.

A combine drops the positions even when it is handed a single vector. c() of one ps_trunc returns it unchanged, record included, but vctrs::vec_c(x), dplyr::bind_rows(df), vctrs::vec_rbind(df), and an ungrouped dplyr::reframe() rebuild the column, so a later is_unit_truncated() on the result refuses it.

unique() keeps one element for each distinct value, or one row for each distinct row of a matrix of scores, and that element or row stands for every unit holding the same scores. A matrix comes back as a matrix of the same class with its column names. The record is re-indexed onto the result when all of the merged units share one status, and dropped otherwise: a score or row that arrived where truncation would put it and one truncation moved there merge into a single element or row that neither status describes.

as.data.frame() and tibble::as_tibble() turn a matrix of scores into a data frame whose columns are ps_trunc vectors carrying the matrix's record, one row per unit, so the weight functions read the data frame as they read the matrix and a subset of its rows re-indexes the record. A data frame whose columns are not all truncated with the same record is refused by the weight functions with an error of class propensity_matrix_type_error. A matrix without column names gives columns named V1, V2, and so on, which the weight functions then refuse because they name no exposure level, so name the columns after the exposure levels before converting.

A record can also outlive the observations it describes, because it travels by routes vctrs does not see: growing a ps_trunc by subassignment carries it across the length change. is_unit_truncated() therefore checks that the record covers the object it is given and raises an error of class propensity_missing_meta_error when it does not, rather than name truncated units at stale positions.

That check compares how many observations the record was written for against how many the object holds, which a reordering does not change, so a route that could reorder the observations without saying how drops the positions instead, at any length: vctrs::vec_slice(), dplyr::arrange(), dplyr::filter(), and vctrs::vec_assign() and the helpers built on it return a ps_trunc whose record keeps its method and bounds and names no units, and is_unit_truncated() refuses it. Subsetting with [, sort(), unique(), and rep() know where the units went and re-index the record, and ⁠[<-⁠ and ⁠is.na<-⁠ move no unit and keep it, so reorder with [, or put the propensity scores in the order you want before truncating them.

Casting a numeric vector into a ps_trunc with vctrs::vec_cast() is a type operation and not a truncation. The result is described by the method and bounds of the target and records that none of the arriving values was pinned to a bound, so it can hold scores outside those bounds, including 0 and 1. Call ps_trunc() on the scores to truncate them.

Value

A ps_trunc object (a numeric vector for binary exposures, or a matrix for categorical exposures). Use ps_trunc_meta() to inspect metadata including method, lower_bound, upper_bound, truncated_idx (positions of modified values), and n_obs (the number of observations those positions describe).

References

Crump, R. K., Hotz, V. J., Imbens, G. W., & Mitnik, O. A. (2009). Dealing with limited overlap in estimation of average treatment effects. Biometrika, 96(1), 187–199.

Gruber, S., Phillips, R. V., Lee, H., & van der Laan, M. J. (2022). Data-adaptive selection of the propensity score truncation level for inverse-probability-weighted and targeted maximum likelihood estimators of marginal point treatment effects. American Journal of Epidemiology, 191(9), 1640–1651.

Walker, A. M., Patrick, A. R., Lauer, M. S., et al. (2013). A tool for assessing the feasibility of comparative effectiveness research. Comparative Effectiveness Research, 3, 11–20.

See Also

ps_trim() for removing (rather than bounding) extreme values, ps_refit() for refitting the propensity model after trimming, is_ps_truncated(), is_unit_truncated(), ps_trunc_meta()

Examples

set.seed(2)
n <- 200
x <- rnorm(n)
z <- rbinom(n, 1, plogis(1.2 * x))
fit <- glm(z ~ x, family = binomial)
ps <- predict(fit, type = "response")

# Truncate to [0.1, 0.9]
ps_t <- ps_trunc(ps, method = "ps", lower = 0.1, upper = 0.9)
ps_t

# Bound at [1/c, 1 - 1/c], c = sqrt(n) log(n) / 5, which caps the
# unstabilized ATE weights at c
ps_trunc(ps, method = "adaptive")

# Truncate at the 1st and 99th percentiles
ps_trunc(ps, method = "pctl", lower = 0.01, upper = 0.99)

# Use truncated scores to calculate weights
wt_ate(ps_t, .exposure = z)

# Inspect which observations were truncated
is_unit_truncated(ps_t)

# Bound the scores a fitted model reports, reading the exposure off the model
ps_trunc(fit, method = "cr")

if (rlang::is_installed("nnet")) {
  trt <- factor(sample(c("a", "b", "c"), n, replace = TRUE))
  multinomial_fit <- nnet::multinom(trt ~ x, trace = FALSE)
  ps_trunc(multinomial_fit, method = "ps")
}

Description

Returns the metadata list attached to a ps_trunc object. The list includes fields such as method, lower_bound, upper_bound, and truncated_idx.

Usage

ps_trunc_meta(x)

Arguments

Value

A named list with truncation metadata, including:

• method – the truncation method used ("ps", "adaptive", "pctl", or "cr")

• lower_bound, upper_bound – the applied bounds

• truncated_idx – integer positions of values that were winsorized

• n_obs – the number of observations those positions describe

truncated_idx and n_obs are absent when an operation that changed the number of observations dropped the record; see ps_trunc().

See Also

ps_trunc(), is_ps_truncated(), is_unit_truncated()

Examples

ps <- c(0.02, 0.3, 0.5, 0.7, 0.98)
ps_t <- ps_trunc(ps, method = "ps", lower = 0.05, upper = 0.95)
ps_trunc_meta(ps_t)

Description

psw objects are numeric vectors that carry metadata about propensity score weights, including the target estimand and whether the underlying propensity scores were trimmed, truncated, or calibrated.

Most users will encounter psw objects as return values from wt_ate() and related weight functions. These constructor and helper functions are useful for inspecting weight objects or for package developers extending propensity.

Usage

new_psw(
  x = double(),
  estimand = NULL,
  stabilized = FALSE,
  trimmed = FALSE,
  truncated = FALSE,
  calibrated = FALSE,
  stabilization_score = NULL,
  wt_truncated = FALSE,
  ...
)

psw(
  x = double(),
  estimand = NULL,
  stabilized = FALSE,
  trimmed = FALSE,
  truncated = FALSE,
  calibrated = FALSE,
  stabilization_score = NULL,
  wt_truncated = FALSE
)

is_psw(x)

is_stabilized(wt)

stabilization_score(wt)

as_psw(x, estimand = NULL)

Arguments

Details

Constructors

• psw() is the user-facing constructor. It coerces x to double and validates inputs before creating the object.

• new_psw() is the low-level constructor intended for developers. It assumes x is already a double vector and performs minimal validation.

• as_psw() coerces an existing numeric vector to a psw object.

Queries

• is_psw() tests whether an object is a psw vector.

• is_stabilized() returns TRUE if the weights are stabilized.

• is_wt_truncated() returns TRUE if the weights themselves were truncated, which is recorded apart from whether the propensity scores they were built from were.

• stabilization_score() returns the user-supplied stabilization score, or NULL when none was recorded or when a per-observation score was dropped because an operation changed the length of the weights.

• numerator_model() returns the fitted model supplied to stabilize that estimated the numerator of the weights, or NULL when no model did.

The stabilization score

A stabilization_score is the multiplier the weights were stabilized on. It holds either a single value, which scales every weight, or one value per observation, which scales each weight by its own. Both forms are checked where the score is recorded: a score must be numeric, positive, and finite, and hold a length the weights can use, and one that does not is refused with an error of class propensity_stabilization_score_error. A score is recorded as a plain double vector, with its storage type normalized and any names dropped, so a score written 1L and one written 1 are the same score and combine without conflict.

The exposure records

Weights built by wt_ate() and its relatives also record the exposure they were built for: exposure_type names its type, and density_meta describes the ratio of densities that weights for a continuous exposure are. Both are read back with exposure_type() and density_meta(), and both describe the exposure rather than the units, so neither carries a length of its own.

A zero-length psw is a prototype: it records what a result built from it will carry rather than describing observations of its own, so a per-observation score on one is recorded as given and checked for length when the observations arrive. A cast whose data arrives at a length the score does not match drops the score rather than refusing the cast, since the score describes the prototype's observations and says nothing about the data being cast. A cast to no observations keeps it, having brought no observations for it to contradict.

psw objects also inherit the broader causal_wts class. The accessors that class carries, causalgenerics::is_causal_wt(), causalgenerics::estimand(), and ⁠estimand<-⁠, are re-exported by propensity and documented at causalgenerics::estimand().

Arithmetic and combining

Arithmetic operations on psw objects preserve the class and attributes, so operations like normalization (weights / sum(weights)) retain metadata.

An operation between two psw objects merges what each of them records. Two different estimands are pasted together, and an estimand only one operand names stands for the result; the result is stabilized only when both operands are, and it is marked as trimmed, truncated, calibrated, or weight-truncated when either operand is. The remaining attributes, the stabilization_score, the records left by a modified propensity score, the attributes describing a categorical exposure, and the exposure records, are carried by agreement: one only a single operand records carries, and one both record with the same value carries. One they record differently is dropped, since neither value describes the result, and a warning of class propensity_metadata_conflict_warning names it, once for each attribute dropped that way and whatever order the inputs were given in. The rule is applied per operation, so an attribute one operation drops for a disagreement can be carried again by a later operation whose operands agree. Two density records agree when they name the same family with the same parameters and the same numerator and residual spread, so weights built the same way in two calls agree even though each call builds its own copy of the function that evaluates the density. A numerator estimated by a fitted model is read the same way: two records agree when their models write the same formula and were fit to the same coefficients, and a model is not compared against the marginal density or against no numerator at all, both of which are disagreements.

A density_meta record describes a continuous exposure, so a result that drops exposure_type for a disagreement drops the density record with it, without a second warning: the record has nothing left to be about.

A stabilization_score is carried only when the result is stabilized, which takes both operands. A result that is not stabilized drops the score without comment.

Combining psw objects with c() preserves the class only when all metadata matches; mismatched metadata produces a warning and falls back to a plain numeric vector. The trimming, truncation, and weight truncation records are compared by what they say about the modification: weights built from scores trimmed or truncated differently, from a trim that was refit and one that was not, or truncated at different bounds describe different estimands and are combined only as numbers. Weights flagged as modified with no record agree with any record. Concatenation appends one set of observations to another, so the positions a record names would describe units from the other input; the result keeps each record without its positions, which is_refit(), the printed footer, and a later combine still read. The calibration record (ps_calib_meta) names the curve the scores were calibrated with rather than any position, the categorical attributes name exposure levels, and the exposure records describe the exposure rather than any unit, so all of them carry by the agreement rule above: two calibration records agree when they name the same method and smoothing.

c() of a single psw returns it unchanged, with every record. A combine through vctrs drops the positional records even when it is handed a single input: vctrs::vec_c(x), dplyr::bind_rows(df), vctrs::vec_rbind(df), and an ungrouped dplyr::reframe() all rebuild the column, so a later is_unit_trimmed(), is_unit_truncated(), or is_unit_wt_truncated() on the result refuses it.

Subsetting with [ preserves class and attributes for vector subscripts. Two kinds of attribute hold one value per observation. The records left by a trimmed or truncated propensity score (ps_trim_meta and ps_trunc_meta) or by truncating the weights themselves (psw_trunc_meta) name units by position, and [ is handed the subscript, so it re-indexes each record onto the result: rev(), sort(), x[order(x)], and a shorter subset all return records naming the units at their new positions. Any other operation through vctrs is not handed a subscript, so it keeps the positions only where every unit stays at its place, as in elementwise arithmetic, and otherwise drops them, at any length: vctrs::vec_slice(), dplyr::arrange(), dplyr::filter(), unique(), rep_len(), and vctrs::vec_c() of a single input all return weights whose records name no units. A dropped record keeps its method, its bounds, and whether the model was refit, so is_refit(), the printed footer, and the bound comparison of a later combine still read it. Subassignment with ⁠[<-⁠ and ⁠is.na<-⁠ moves no unit and keeps the records whole; ⁠[<-⁠ refuses a value whose records describe a different modification, and accepts one with no record. A ⁠[<-⁠ of a value modified the same way that writes a missing weight leaves that unit listed as retained in the record, as ⁠is.na<-⁠ on a psw does; ⁠is.na<-⁠ on a ps_trim score instead removes the unit from its retained set. vctrs' own assignment, vctrs::vec_assign(), reaches the same restore a slice does and cannot be told apart from one, so it drops the positions, and so do the helpers built on it or on a combine, such as tidyr::replace_na(), dplyr::coalesce(), dplyr::if_else(), and dplyr::case_when().

A stabilization_score with more than one value, on the weights or on a component of a wt_joint() product, is in the order of the units too. [ subsets it with the weights, elementwise arithmetic keeps it, and any other route through vctrs that brings back observations drops it with a warning of class propensity_stabilization_score_warning, since the score may no longer line up with the weights.

Dropping the stabilization_score warns, because the score was supplied by the user and the weights can be recomputed on the subset. Dropping a modification record is silent, because these records also travel by routes vctrs does not see, and a warning on the one route it controls would be neither complete nor about anything the user wrote. Subassignment that grows the vector carries a record across the length change, since ⁠[<-⁠ casts the replacement and then leaves base R to preserve the attributes; and model.frame() drops the NA-weighted rows from a weights column in C and re-attaches the original variable's attributes to the shortened result, so weights built on trimmed propensity scores come back out of every outcome model fit on them still carrying a record written for rows that are no longer there.

Honesty therefore lives at query time. is_unit_trimmed(), is_unit_truncated(), and is_unit_wt_truncated() answer by position, so each checks that the record covers the vector it is given and raises an error of class propensity_missing_meta_error when it does not, or when weights marked as modified carry no record at all, rather than name modified units at stale positions. is_refit() reads a single flag rather than a position, so it answers from any record present and refuses only when the record is absent entirely.

The result of any of these operations stays a psw and keeps every other attribute, including its stabilized, trimmed, truncated, and calibrated status, the calibration record, which names a calibration curve rather than any unit, the attributes describing a categorical exposure, which name the exposure levels rather than the units, and the exposure records, which describe the exposure rather than any unit, so all of them mean the same thing at any length.

Matrix or array subscripts intentionally drop the psw class and return a plain numeric vector via base R linear indexing; this is required so glm.fit()-style internal indexing works on psw-weighted GLMs. Summary functions (sum(), mean(), etc.) return plain numeric values.

Value

• new_psw(), psw(), as_psw(): A psw vector.

• is_psw(), is_stabilized(): A single logical value.

• stabilization_score(): A numeric value or vector, or NULL if none was recorded or a per-observation score was dropped.

See Also

wt_ate(), wt_att(), wt_atu(), wt_atm(), wt_ato() for calculating propensity score weights (which return psw objects).

ps_trim(), ps_trunc(), and ps_calibrate() for modifying propensity scores before weight calculation.

Examples

# Create psw objects directly
w <- psw(c(1.2, 0.8, 1.5), estimand = "ate")
w

# Query metadata
is_psw(w)
estimand(w)
is_stabilized(w)

# Coerce a plain numeric vector
as_psw(c(1.0, 2.0), estimand = "att")

# Arithmetic preserves the psw class
w / sum(w)

# Combining: compatible metadata is preserved
x <- psw(c(1.2, 0.8), estimand = "ate")
y <- psw(c(1.1, 0.9), estimand = "ate")
c(x, y)

# Combining: incompatible metadata warns and returns numeric
x <- psw(c(1.2, 0.8), estimand = "ate")
y <- psw(c(1.1, 0.9), estimand = "att")
c(x, y)

Description

tidy() returns the estimates of an ipw() result as a tibble using the column names broom conventions use. For a binary or categorical exposure there is one row per exposure level, under the term "mean", and then one row per effect measure per contrast of levels; for a continuous exposure there is one row per reported coefficient. The rows arrive in the order the result stores them and nothing is dropped.

Those columns are the ones the result's own coercion surface reports, so either reading is as.data.frame() read as a tibble rather than a second assembly of the same table. The values are the ones the result already holds, nothing being re-estimated: the confidence interval is the only thing rebuilt, and only when the requested conf.level differs from the level the result was fit at. The one thing the frame carries that the tibble does not is the covariance it attaches as an attribute, a tidied table being its columns.

conf.int, conf.level, and exponentiate are the arguments this method shares with that surface, so the surface validates them, in both readings and before either is assembled. A value it refuses is reported under the causalgenerics condition naming the argument at fault, which is what keeps one argument of one method from being well formed in one reading of a result and refused in the other.

A result reports its effects in one of two readings, and tidy() returns the one the result records unless effects names the other for the call. The marginal reading is the table of causal contrasts described above. The conditional reading is the outcome model's coefficient surface: one row per coefficient, with the standard errors of the block of the joint estimation that carries the uncertainty of having estimated the weights from the same data. The two readings return the same columns in the same order, with one exception: the contrast column belongs to the marginal reading. A binary and a categorical result both carry it there, naming the exposure level each mean row belongs to and the pair of levels each effect measure compares, and the conditional reading of either returns the same table one column narrower. A coefficient names neither a level nor a pair of them.

Usage

## S3 method for class 'ipw'
tidy(
  x,
  conf.int = FALSE,
  conf.level = 0.95,
  exponentiate = FALSE,
  ...,
  effects = NULL,
  parametric = NULL
)

Arguments

Value

A tibble with one row per estimate and the columns:

term

The effect measure, such as "mean" for a counterfactual mean, or "rd", "log(rr)", "log(or)", "diff", "slope", or "coef" for the measures built from them.

contrast

What the row is about: the exposure level a "mean" row belongs to, such as "0" or "b"; the pair of levels an effect measure compares, such as "1 vs 0" or "b vs a"; or the coefficient name on a surface whose rows are named after their coefficients. Present only where a row reports one.

estimate

The estimated effect.

std.error

The standard error of the estimate.

statistic

The z statistic, the estimate over its standard error.

p.value

The two-sided p-value of that statistic.

conf.low, conf.high

The interval bounds. Present only when conf.int is TRUE.

The conditional reading returns those columns in that order, holding one row per coefficient of the outcome model: term is the coefficient's name, estimate is the coefficient, std.error is the square root of the diagonal of the corrected covariance, statistic is the estimate over that standard error, and p.value is the two-sided normal p-value of the statistic. There is no contrast column, whatever the exposure, because a coefficient names neither an exposure level nor a pair of them, and the bounds are built at the level conf.level asks for, the stored ones belonging to the effects the marginal reading reports.

See Also

ipw() for the estimator, glance() for a one-row summary of the fit, augment() for its per-observation columns, as.data.frame() for the result's own columns, and as_marginal() and as_conditional() for the reading a result records. For results fitted to multiply imputed data, see the Multiple imputation section of ipw(), pool_ipw() for the pooling, and tidy() for what it returns.

Examples

set.seed(123)
n <- 200
x1 <- rnorm(n)
z <- rbinom(n, 1, plogis(0.5 * x1))
y <- rbinom(n, 1, plogis(-0.5 + 0.8 * z + 0.3 * x1))
dat <- data.frame(x1, z, y)

ps_mod <- glm(z ~ x1, data = dat, family = binomial())
wts <- wt_ate(ps_mod)
outcome_mod <- glm(y ~ z, data = dat, family = quasibinomial(), weights = wts)
result <- ipw(ps_mod, outcome_mod)

tidy(result)

# A 90% interval, with the ratio measures on their natural scale
tidy(result, conf.int = TRUE, conf.level = 0.9, exponentiate = TRUE)

# The outcome model's coefficients, with the standard errors the joint
# estimation of the weights and the outcome implies
tidy(result, conf.int = TRUE, effects = "conditional")

Description

tidy() returns the estimates of a pool_ipw() result as a tibble using the column names broom conventions use. It holds the rows the results themselves report: for a binary or categorical exposure one row per exposure level, under the term "mean", and then one row per effect measure per contrast of levels; for a continuous exposure one row per reported coefficient. The rows arrive in the order the pooled result stores them and nothing is dropped.

Those columns are the ones the pooled result's own coercion surface reports, so this method is that surface read as a tibble rather than a second assembly of the same table. The covariance of the effects the frame attaches is the one thing that does not travel: it is an attribute rather than a column, and a tidied table is its columns.

The pooling has already happened by the time this method sees the result. Nothing here re-pools and no argument changes an estimate: they say how the table already in the object is reported, which is which of the two readings it is, whether an interval comes with it, the level that interval is reported at, and the scale the ratio rows are put on.

A pooled result reports its effects in one of two readings, and tidy() returns the one the result records unless effects names the other for the call. pool_ipw() pools both readings of the results it is given, so naming one here reports a table the pooling already built rather than pooling again. The marginal reading is the table of pooled causal contrasts described above. The conditional reading is the outcome models' coefficient surface, pooled the same way: one row per coefficient, from the block of the joint estimation that carries the uncertainty of having estimated the weights from the same data. The two readings return the same columns in the same order, with one exception: the contrast column belongs to the marginal reading. A binary and a categorical pool both carry it there, naming the exposure level each mean row belongs to and the pair of levels each effect measure compares, and the conditional reading of either returns the same table one column narrower.

Usage

## S3 method for class 'ipw_pooled'
tidy(
  x,
  conf.int = FALSE,
  conf.level = NULL,
  exponentiate = FALSE,
  ...,
  effects = NULL
)

Arguments