Introduction to propertee
propertee Authors
2025-01-17
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 or standardization for experiments and observational studies, with the option of using a separately fitted prognostic regression model for covariance adjustment. Propertee calls for a self-standing specification of a study’s design and treatment allocations, using these to create weights needed to align treatment and control samples for estimation of treatment effects averaged across a common standard population, with standard errors that appropriately reflect the study design. For covariance adjustment 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:
- Encode the study’s design, including the unit of assignment,
treatment status of each unit of assignment, and block membership
information as appropriate, in a
StudySpecification
object. This is accomplished with theobs_spec()
,rct_spec()
orrd_spec()
functions. - Optionally, fit a covariate adjustment model. This is done with any
of a number of modeling functions external to
propertee
, for examplelm()
orglm()
from thestats
package, orlmrob
from therobustbase
package. - Use
lm()
to contrast the separately reweighted treatment and control groups, incorporating optional covariance adjustments as anoffset
tolm()
. TheStudySpecification
is unobtrusively retained, as is error propagation from the optional covariance adjustment, for use in subsequent standard error calculations. This is done using thelmitt()
function, which either takes and enriches a fittedlm
or, in the more fully supported usage, internally callslm()
to itself create a directly adjusted contrast of treatment conditions.
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
We take the random assignment scheme to have been either complete or
Bernoulli randomization of students, separately by school.1 Accordingly,
experimental blocks are given by 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
We need a unique identifier for the unit by which treatment was
allocated. STAR treatments were assigned to students, as opposed to
classrooms or other aggregates of students; accordingly we need a unique
identifier per student. The studentid
variable fills this
role.
A Basic Example
Defining the StudySpecification
The three _spec
functions (rct_spec()
,
obj_spec()
, and rd_spec()
) operate similarly.
The first argument is the most important, and encodes all the
specification 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 specifications. The aliasuoa()
can be used in its place. -
block()
: The identifies the variable(s) which contain block information. Optional. -
forcing()
: In regression discontinuity specifications (rd_spec()
), this identifies the variable(s) which contain forcing information.
To define a StudySpecification
in our example:
spec <- obs_spec(starkbinary ~ unit_of_assignment(studentid) + block(schoolidk),
data = STARdata, na.fail = FALSE)
summary(spec)
#> 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 additional variables be needed to identify the unit of
assignment, block, or forcing, they can be included. Had the variable
identifying the years in which participants entered the study been
available, for example, we might have identified the experimental blocks
as combinations of school and year of study entry, rather than as the
school alone. To do this we would pass
block(entryyear, schoolidk)
, rather than
block(schoolidk)
, in the above.
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
StudySpecification
.
Note that the data set does not need to be the same
data set which generated the StudySpecification
; 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 StudySpecification
,
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, specification = spec)
#> The StudySpecification 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, specification = spec)
#>
#> 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, specification = spec)
#> The StudySpecification 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, specification = spec)
#>
#> 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.93e-06 ***
#> starkbinary.TRUE_ethnicityasian -9.939 21.021 -0.473 0.6363
#> starkbinary.TRUE_ethnicityhispanic 35.667 18.634 1.914 0.0557 .
#> starkbinary.TRUE_ethnicityamindian 19.000 NA NA NA
#> starkbinary.TRUE_ethnicityother 13.333 26.856 0.496 0.6196
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Std. Error calculated via type "CR0"
Including specification weights
Study specification 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, specification = spec, weights = "ate")
#> starkbinary.TRUE
#> 6.116683
lmitt(readk ~ 1, data = STARdata, specification = spec, 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, specification = spec,
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(spec),
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(spec),
#> 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, specification = spec, absorb = TRUE)
#> starkbinary.TRUE
#> 5.14933