Introduction to propertee
propertee Authors
2024-10-29
Source:vignettes/intro-to-propertee.Rmd
intro-to-propertee.Rmd
Main Features
The propertee package, Prognostic Regression Offsets with Propagation of ERrors (for Treatment Effect Estimation), facilitates direct adjustment for experiments and observational studies with design-informed standard errors and flexible options for covariance adjustment. It uses explicit specification of study design to provide probability of assignment weights and standard errors that appropriately reflect the design. For covariance adjustment of its Hajek and (one-way) fixed effects estimates, it enables offsetting the outcome against predictions from a dedicated covariance model, with standard error calculations propagating error as appropriate from the covariance model.
The main workflow consists of two main steps and one optional step:
- Generate a
Design
object which encodes the study design, including the unit of assignment, treatment status of each unit of assignment, and optionally block information. This is accomplished with theobs_design()
,rct_design()
orrd_design()
functions. - Optionally, fit a covariate adjustment model.
- Fit a model to estimate a treatment effect, accounting for the
design information, and optionally the covariate adjustment. This is
done via the
lmitt()
function.
Example Data
The example dataset comes from the state of Tennessee’s Student-Teacher Achievement Ratio (STAR) experiment. Students were randomly assigned to three possible classroom conditions: small (13 to 17 students per teacher), regular class (22 to 25 students per teacher), and regular-with-aide class (22 to 25 students with a full-time teacher’s aide).
For simplicity for this first example, we will examine a single binary treatment - “small” classrooms versus “regular” and “regular+aide” classrooms.
STARdata$starkbinary <- STARdata$stark == "small"
table(STARdata$starkbinary)
#>
#> FALSE TRUE
#> 4425 1900
After this basic example, we will see how propertee
makes it easy to handle non-binary treatment variables by introducing
dichotomy
s.
The outcome of interest is a reading score at the end of kindergarten.
summary(STARdata$readk)
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> 315.0 414.0 433.0 436.7 453.0 627.0 5809
The students were blocked into schools via the schoolidk
variable:
length(unique(STARdata$schoolidk))
#> [1] 80
head(table(STARdata$schoolidk))
#>
#> 1 2 3 4 5 6
#> 74 54 100 60 61 54
Students were the assigned units, so we need a unique identifier per student. If this does not currently exist, it can easily be generated:
A Basic Example
Defining the Design
The three _design
functions (rct_design()
,
obj_design()
, and rd_design()
) operate
similarly. The first argument is the most important, and encodes all the
design information through the use of a formula. The left-hand side of
the formula identifies the treatment variable. The right-hand side of
the formula consists of the following potential pieces of
information:
-
unit_of_assignment()
: This identifies the variable(s) which indicate the units of assignment. This is required for all designs. The aliasuoa()
can be used in its place. -
block()
: The identifies the variable(s) which contain block information. Optional. -
forcing()
: In regression discontinuity designs (rd_design()
), this identifies the variable(s) which contain forcing information.
To define a Design
in our example:
des <- obs_design(starkbinary ~ unit_of_assignment(studentid) + block(schoolidk),
data = STARdata, na.fail = FALSE)
summary(des)
#> Observational Study
#>
#> Structure Variables
#> --------- ---------
#> Treatment starkbinary
#> Unit of Assignment studentid
#> Block schoolidk
#>
#> Number of units per Treatment group:
#> Txt Grp Num Units
#> FALSE 4425
#> TRUE 1900
Should more than one variable be needed to identify the unit of
assignment, block, or forcing, they can be included. For example,
perhaps schoolidk
may be unique within district, but
potentially not unique across districts. Then we’d use something like
block(districtid, schoolidk)
in the _design
function.
Estimating the treatment effect
The main function for estimating treatment effects is the
lmitt()
function. It takes in three main required
arguments:
- A formula specifying the outcome and the desired treatment effect.
- The data set containing the outcome information.
- A
Design
.
Note that the data set does not need to be the same
data set which generated the Design
; it does however need
to include the same variables to identify the units of assignment. (If
the variable names differ, the by=
argument can be used to
link them, though we recommend renaming to reduce the likelihood of
issues.)
For example, you may have one dataset containing school-level
information, and a separate dataset containing student-level
information. Assume school is the unit of assignment. While you could of
course merge those two data-sets, propertee can instead
use the school-level data to define the Design
, and the
student-level data to estimate the treatment effect.
The formula entering lmitt()
can take on one of two
forms:
y ~ 1
will estimate the main treatment effect on outcome variable
y
, and
y ~ x
will estimate subgroup specific treatment effects for each level of
x
for the outcome y
, if x
is
categorical. For continuous x
, a main effect and a
treatment-x
interaction effect is estimated.
Therefore, to estimate the treatment effect in our example, we can run:
te <- lmitt(readk ~ 1, data = STARdata, design = des)
#> The Design object contains block-level information, but it is not used in this model. Block information is used when weights are defined via `ate()` or `ett()` or if the `absorb=TRUE` argument is passed.
summary(te)
#>
#> Call:
#> lmitt(readk ~ 1, data = STARdata, design = des)
#>
#> Treatment Effects :
#> Estimate Std. Error t value Pr(>|t|)
#> starkbinary.TRUE 5.4632 0.9207 5.934 3.13e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Std. Error calculated via type "CR0"
The data includes ethnicity; we can estimate subgroup effects by ethnicity:
te_s <- lmitt(readk ~ ethnicity, data = STARdata, design = des)
#> The Design object contains block-level information, but it is not used in this model. Block information is used when weights are defined via `ate()` or `ett()` or if the `absorb=TRUE` argument is passed.
summary(te_s)
#> Warning: The following subgroups do not have sufficient degrees of freedom for
#> standard error estimates and will be returned as NA: ethnicityamindian
#>
#> Call:
#> lmitt(readk ~ ethnicity, data = STARdata, design = des)
#>
#> Treatment Effects :
#> Estimate Std. Error t value Pr(>|t|)
#> starkbinary.TRUE_ethnicitycauc 4.717 1.147 4.112 3.97e-05 ***
#> starkbinary.TRUE_ethnicityafam 6.607 1.468 4.500 6.94e-06 ***
#> starkbinary.TRUE_ethnicityasian -9.939 21.021 -0.473 0.6364
#> starkbinary.TRUE_ethnicityhispanic 35.667 18.635 1.914 0.0557 .
#> starkbinary.TRUE_ethnicityamindian 19.000 NA NA NA
#> starkbinary.TRUE_ethnicityother 13.333 26.857 0.496 0.6196
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Std. Error calculated via type "CR0"
Including design weights
Design weights can be easily included in this estimation. propertee supports average treatment effect (ATE) and effect of the treatment on the treated (ETT) weights.
To include one of the weights, simply include the
weights = "ate"
or weights = "ett"
argument to
lmitt()
:
lmitt(readk ~ 1, data = STARdata, design = des, weights = "ate")
#> starkbinary.TRUE
#> 6.116683
lmitt(readk ~ 1, data = STARdata, design = des, weights = "ett")
#> starkbinary.TRUE
#> 5.650771
Internally, these call the ate()
or ett()
functions which can be used directly.
When included inside lmitt()
, you do not need to specify
any additional arguments to ate()
or ett()
,
enabling easy functions of weights. For example if you had some other
weight variable, say wgt
, you could include
weights = wgt*ate()
in the lmitt()
call.
Covariance Adjustment models
By itself, lmitt()
does not allow for other covariates;
e.g. something like lmitt(y ~ 1 + control_var,...
will
fail. To adjust for covariates, a separate covariate model should be
fit. Any model which supports a predict()
function should
work.
camod <- lm(readk ~ gender + birth + lunchk, data = STARdata)
The cov_adj()
function can be used to process the
covariance adjustment model and produce the required values; and its
output can be passed as an offset=
.
lmitt(readk ~ 1, data = STARdata, design = des,
weights = "ate", offset = cov_adj(camod))
#> Warning in validityMethod(object): Some covariance adjustments are NA; be
#> careful of dropping these observations when fitting the ITT effect model
#> starkbinary.TRUE
#> 6.050889
Similarly to the weight functions, cov_adj()
attempts to
locate the correct arguments (in this case, mainly the
data=
argument) to use in the model command; while
cov_adj()
does fall back to using the data which is in the
covariance model, its safer to use the newdata=
argument if
calling cov_adj()
outside of the model.
Also, similarly to weights, cov_adj()
can be used in
normal modeling commands as well.
lm(readk ~ starkbinary, data = STARdata, weights = ate(des),
offset = cov_adj(camod))
#> Warning in validityMethod(object): Some covariance adjustments are NA; be
#> careful of dropping these observations when fitting the ITT effect model
#>
#> Call:
#> lm(formula = readk ~ starkbinary, data = STARdata, weights = ate(des),
#> offset = cov_adj(camod))
#>
#> Coefficients:
#> (Intercept) starkbinaryTRUE
#> -1.684 6.051
Absorbing Blocks
If fixed effects for blocks are desired, which can be absorbed away
to avoid estimating, the absorb=TRUE
argument can be
passed.
lmitt(readk ~ 1, data = STARdata, design = des, absorb = TRUE)
#> starkbinary.TRUE
#> 5.14933