Causal effects from a fitted model

Estimates the average or conditional effect of a treatment by g-computation over the posterior draws. estimate_effect() predicts every unit under every level of the treatment, contrasts those potential outcomes within each draw, and averages over the units the estimand asks for.

Usage

estimate_effect(
  object,
  treat = NULL,
  estimand = "ATE",
  comparison = "difference",
  by = NULL,
  newdata = NULL,
  level = 0.95,
  interval = "eti",
  focal = NULL,
  type = "response"
)

# S3 method for class 'bartisan_effect'
print(x, digits = 3L, contrasts = NULL, potential_outcomes = TRUE, ...)

# S3 method for class 'bartisan_effect'
plot(x, marginal = TRUE, ...)

Arguments

object: a <bartisan_fit> object; the output of a call to bcf() or bartisan().
treat: string; the name of the treatment variable. A fit from bcf() carries its own and needs none, so this is for a fit from bartisan(), where nothing marks one predictor as the treatment.
estimand: string; which units the effect is averaged over. "ATE" (the default) uses all of them, and "ATT" and "ATC" those in the focal group, which by default is the treatment's second level for the former and its first for the latter. With a binary treatment, that makes them the treated and the untreated without focal being named. "CATE" does not average at all and returns the conditional effect at each unit's covariates, which is not that unit's own individual effect; see Details. Abbreviations and lowercase spellings are allowed.
comparison: string; how the two potential outcomes are contrasted. Allowable options include "difference" (the default), "ratio", "lnratio", "or", and "lnor". The last two are available only when the response is a probability.
by: optional; a one-sided formula or a variable name naming a grouping variable, in which case the effect is averaged within each of its levels rather than over the whole sample. A formula is evaluated with stats::model.frame() rather than read for the names it mentions, so it can define a grouping the data has no column for, as in by = ~ age > 50 or by = ~ interaction(sex, region). It must give exactly one grouping variable, and the term as written names the column it occupies in the output.
newdata: optional; a data frame of units to average over. Default is the data the model was fit to, so that the average runs over the covariates of the fitted sample; see Details for what that average is.
level: numeric; the level of the credible interval. Default is .95.
interval: string; "eti" (the default) for an equal-tailed interval from the quantiles of the draws, or "hpdi" for the highest posterior density interval, which is the shortest interval containing level of the posterior mass.
focal: optional; the treatment level whose units the effect is averaged over for "ATT" and "ATC". With a binary treatment it defaults to the second level for "ATT" and the first for "ATC", which makes them the effect among the treated and among the untreated. With a treatment of more than two levels it has no default and must be supplied, and then "ATT" and "ATC" differ only in the level named.
type: string; the prediction scale the effect is computed on, passed to predict.bartisan_fit(). Default is "response", which is the scale on which an average of unit-level differences is the marginal effect. See Details before changing it.
x: a <bartisan_effect> object; the output of a call to estimate_effect().
digits: integer; the number of significant digits to print.
contrasts: string; which contrasts to display when the treatment has more than two levels. "focal" (the default when a focal group is known) shows only the contrasts involving it, "all" shows every pairwise contrast, and a character vector of contrast labels shows those. All of them are computed either way; this only decides what is printed.
potential_outcomes: logical; whether to print the average response under each treatment level below the contrasts, those being what the contrasts were computed from. Default is TRUE. They are in the result's "potential_outcomes" attribute either way.
...: ignored.
marginal: logical; for plot() on conditional effects (i.e., estimand = "CATE"), whether to draw the marginal effect beside the units. Default is TRUE. Ignored by the other plots, which draw no marginal effect beside their estimates.

Value

A <bartisan_effect> object, a data frame with one row per reported effect and columns contrast, estimate, lower, and upper, plus unit when estimand = "CATE" and the by variable's name when by is used. The posterior draws of every reported quantity are kept in the "draws" attribute, and the marginal mean of each potential outcome in "potential_outcomes", so a caller can re-contrast or re-summarize without refitting.

print() returns its input invisibly, and prints the potential outcomes below the contrasts unless potential_outcomes = FALSE. plot() returns a ggplot2 object.

Details

The Estimand

The ATE, ATT and ATC are averages of the conditional effect over the covariates of the units averaged over, computed within each posterior draw. This is sometimes called the mixed average treatment effect. The covariates are treated as fixed, so the interval reflects uncertainty about the outcome model but not the further variation a population average would carry, and the units' observed outcomes are not conditioned on, as they would be for a sample average of individual effects. Likewise, estimand = "CATE" reports the expected effect at each unit's covariates rather than the unit's own effect, which depends on how its two potential outcomes are associated and is not identified. See vignette("causal") for the distinction between these estimands.

Setting `type`

A varying coefficient is a contrast on the link scale: on a binomial("logit") fit, coef(object) is a per-unit difference in log odds, and the average of those is the average conditional log odds ratio rather than the marginal one a treatment question usually asks for. The default is therefore type = "response", where every unit's contrast is on the scale the response is measured on and averaging them gives the marginal effect.

Setting type = "link" is the right choice for looking at how the effect varies, since that is the scale the forest models it on, and the wrong one for reporting an average.

Setting `comparison`

For "ratio", "lnratio", "or" and "lnor", the potential outcomes are averaged over units first and contrasted afterward, which gives the marginal ratio. With estimand = "CATE" there is no averaging to do, so a ratio reported there is a conditional ratio and the print() method says so.

For the same reason, "or" and "lnor" are not each other's exp() and log(). Each summarizes the posterior of the quantity it names, and a posterior mean does not survive a nonlinear transformation: "or" reports the mean of the odds ratio and "lnor" the mean of its logarithm, which exponentiates to something smaller. Report whichever scale the interval should be symmetric on, which for a ratio is usually the log. The same holds of "ratio" against "lnratio".

Setting `focal`

With a binary treatment and no focal, which level is the treated one is worked out from the treatment's values by the rules WeightIt::weightit() uses, so that the same variable is read the same way by both packages. In order: a logical treatment is treated at TRUE; a numeric one is untreated at 0; a character or factor one whose levels parse as numbers is untreated at 0, or else treated at the largest; and otherwise the conventional names are matched, "t", "tr", "treat", "treated" and "exposed" against "c", "co", "ctrl", "control" and "unexposed". When none of those applies the second level is taken as the treated one and a message says so, that being the case where the guess can be wrong. "ATC" takes the other level as its focal group, so both estimands average over the group focal names and only the default differs.

Multi-category Treatments

There is nothing to work out from the values, so the focal argument is required for "ATT" and "ATC", and those two then name the same estimand: the effect among the units in the level named.

Every pairwise contrast is computed. Which ones are shown is a display choice, made by contrasts in the print() method: with a focal group the default is to show only the contrasts involving it, matching how WeightIt::weightit() and cobalt::bal.tab() use focal.

The effect of a continuous treatment is a slope or a dose-response curve rather than a contrast of levels, and the adrftools package has tools for summarizing and visualizing one.

Examples

data("rhc")
set.seed(123)

# Note: we use small sampler settings that would not be
# appropriate for analysis; this is just for the
# examples
fit <- bcf(death ~ age + sex + meanbp + aps,
           treat = ~ rhc, data = rhc, num_trees = 10,
           num_burn = 50, num_draws = 50,
           propensity_args = list(num_trees = 10,
                                  num_burn = 50,
                                  num_draws = 50),
           verbose = FALSE)
#> ℹ Using `family = binomial()`.
#> ℹ Set `family` explicitly to silence this message.

# The risk difference, averaged over everyone
estimate_effect(fit)
#> Average treatment effect (difference)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> 
#>     contrast estimate   lower  upper    n
#>  Y[1] - Y[0]   0.0159 -0.0126 0.0419 1500
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.649 0.627 0.679
#>      Y[1]    0.665 0.638 0.688
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".

# Among the treated, and as a risk ratio rather than a
# difference
estimate_effect(fit, estimand = "ATT",
                comparison = "ratio")
#> Average treatment effect on the treated (ratio)
#> 
#> Treatment: `rhc`
#> Averaged over the 565 units in group "1"
#> 
#>     contrast estimate lower upper   n
#>  Y[1] / Y[0]     1.03 0.979  1.08 565
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.675 0.641 0.703
#>      Y[1]    0.694 0.667 0.719
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".

# The effect at each unit's covariates, ordered, with
# the marginal effect beside them
cate <- estimate_effect(fit, estimand = "CATE")
plot(cate)


# By a subgroup, which is where a forest plot earns its keep
estimate_effect(fit, by = ~ sex)
#> Average treatment effect (difference)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> Within levels of `sex`
#> 
#>     sex    contrast estimate    lower  upper   n
#>  female Y[1] - Y[0]   0.0167 -0.00959 0.0411 676
#>    male Y[1] - Y[0]   0.0152 -0.01520 0.0419 824
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.649 0.627 0.679
#>      Y[1]    0.665 0.638 0.688
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".