Nothing
staggR
In staggR: Fit Difference-in-Differences Models with Staggered Interventions
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
# Check for sandwich package
if (!requireNamespace("sandwich", quietly = TRUE)) {
stop("Package 'sandwich' needed for this vignette to work. Please install it.")
}
staggR: Fit difference-in-differences models with staggered interventions
The staggR package fits linear difference-in-differences models in scenarios where intervention roll-outs are
staggered over time. The package implements a version of an approach proposed by Sun and Abraham @sun2021did to
estimate cohort- and time-since-treatment specific difference-in-differences parameters. This package provides an
alternative approach that proceeds in two steps:
-
Estimate a difference-in-differences regression with cohort- and (calendar) time-specific parameters;
-
Translate difference-in-differences estimates into average estimates for arbitrary time periods relative to the intervention (e.g., post-intervention period).
By separating these two steps, we provide an easy-to-implement, transparent, intuitive difference-in-differences
framework.
Background
Consider a scenario with $i=1,...,N$ units and $t=1,...,T$ periods. Units are exposed to the intervention at different
times. We group units having the intervention during the same time period into cohorts, which we denote by $c=1,...,C$.
To simplify the notation, we assume that units in cohort $c=1$ are exposed to treatment at time $t=1$, etc. Units
exposed to the intervention at time $t>T$ form the comparison group because they are not exposed to the intervention
during the study period.
Abraham and Sun's difference-in-differences specification interacts cohort indicators with time indicators relative to
intervention exposure start, and might be written as follows:
$$
y_{it} = \alpha_c D_c + \beta_t D_t + \left(\sum_c \sum_{l\neq -1} \gamma_{cl} D_c \cdot D^l_{it}\right) + \psi X_{it} + \varepsilon_{st},
$$
where $D_c$ are cohort indicators equal to one if unit $i$ is in cohort $c$ and zero otherwise, $D_t$ are time
indicators, and $D^l_{it}$ are indicators equal to one if unit $i$ is $l$ periods away from initial intervention
exposure, $X_{it}$ are covariates, and $\varepsilon_{st}$ is the error term. The last pre-intervention period $l=-1$
is omitted to ensure identification. The parameters $\gamma_{cl}$ measure difference-in-differences estimates for a
cohort $c$ and time since intervention $l$.
Our approach recognizes the following relationship between time since intervention $l$ and time $t$, which was also
noted by Sun and Abraham:
$$
l = t-c
$$
For instance, time $t=2$ for cohort $c=3$ is $l=-1$, i.e., the last period prior to intervention exposure,
because intervention exposure for cohort $c$ is $t=3$.
Using this insight, the difference-in-differences regression can alternatively be stated as follows:
$$
y_{it} = \alpha_c D_c + \beta_t D_t + \left(\sum_c \sum_{t\neq (c-1)} \gamma_{ct} D_c \cdot D_t\right) + \psi X_{it} + \varepsilon_{st}.
$$
Using this specification has two advantages:
- It clarifies that this staggered difference-in-differences method is a generalization of the basic
difference-in-differences design with two periods (a pre-intervention and a post-intervention period) and two cohorts (a
treatment group and a comparison group). In that setting, the indicator for the post-intervention period measures the
change between the pre- and post-intervention period for the comparison group, and the difference-in-differences
estimate measures the
differential change between the pre- and post-intervention period for the treatment group.
The above equation preserves this structure, while showing how it is a more general version of the
difference-in-differences design. Specifically, parameters $\beta_t$ measure changes over time for the comparison
group, which now can include multiple periods; and parameters $\gamma_{ct}$ measure differential changes, but for
multiple treatment cohorts instead of just one treatment group, and for multiple periods instead of just one.
- It can be easily implemented in R (or other statistical programs) because it does not require constructing
cohort-specific time-since-intervention indicators $D^l_{it}$. Instead, it uses cohort and time indicators and the
interaction between the two.
This packages provides functions to estimate such a difference-in-differences specification. It also provides functions
to translate difference-in-differences estimates into quantities that are useful for research purposes. For instance,
users can calculate average difference-in-differences estimates for each period relative to intervention exposure, or
they can calculate the average difference-in-differences estimates for the full post-intervention period. The
ave_coeff() function accomplishes this by calculating averages of user-specified sets of regression coefficients and
combining the corresponding standard errors using the estimated variance–covariance matrix, whose calculation method
users can also specify. This aggregation yields a single estimate that properly accounts for sampling variability and
the correlation among underlying event-study coefficients.
Example using simulated data
Explore data
Consider as a didactic example a policy intervention designed to reduce inpatient hospitalizations in 15 counties. This
longitudinal data set contains one row per individual-year. Each individual in the data set is identified by a globally
unique identifier (guid), and we have measures of the individuals' ages, sexes, and comorbidities, and an indicator
showing whether the individual was hospitalized during the current year. Each individual resides in a county, which
itself is grouped into cohorts based on intervention year. Finally, the yr column indicates the year of the current
observation.
wzxhzdk:1
The county is never exposed to the intervetnion if intervention_yr is NA. Among the 15 counties in the data set, 3
counties implemented the intervention in 2015; 2 counties implemented in 2016; 5 counties implemented in 2017; 1 county
implemented in 2018; and 4 counties did not implement the intervention at all during the study period.
with(hosp,
table(county, intervention_yr, useNA = "ifany"))
Counties are organized into cohorts by intervention year, i.e., counties with the same intervention year are assigned to
the same cohort.
county_cohorts <- unique(hosp[, c("county", "cohort", "intervention_yr")])
county_cohorts[order(county_cohorts$cohort), ]
The study period includes 11 years, from 2010 through 2020. Hospitalizations occurs in every county-year.
stats::xtabs(hospitalized ~ county + yr,
data = hosp)
Visualize hospitalization trends
To visualize trends in the proportion of individuals who were hospitalized over time, we can call the function
ts_plot(), which takes, at minimum, 3 arguments: a formula, in the format outcome ~ group + time_var, the source data set, and the name of the column that contains the time period during which each group implemented the intervention.
staggR::ts_plot(hospitalized ~ county + yr,
df = hosp,
intervention_var = "intervention_yr")
Equivalently, we can specify the outcome, group, and time variable using the parameters y, group, and time_var.
staggR::ts_plot(y = "hospitalized",
group = "county",
time_var = "yr",
df = hosp,
intervention_var = "intervention_yr")
While this visualization is helpful, it excludes counties that did not implement the intervention. This is because the
function ts_plot() omits any counties where any of y, group, time_var, (or our formula equivalent, hospitalized ~ county + yr) or the intervention_var are NA.
table(hosp$hospitalized, useNA = "always")
table(hosp$county, useNA = "always")
table(hosp$yr, useNA = "always")
table(hosp$intervention_yr, useNA = "always")
We can add these counties to the plot if we change NAs in intervention_yr to some other string value.
hosp2 <- hosp
hosp2$intervention_yr <- ifelse(is.na(hosp2$intervention_yr),
yes = "Comparison group",
no = hosp2$intervention_yr)
county_plot <- staggR::ts_plot(hospitalized ~ county + yr,
df = hosp2,
intervention_var = "intervention_yr")
county_plot
It might be easier to compare the panels in this plot if the X axes were all aligned at the point of intervention. We
can accomplish this by passing a time-since-intervention (tsi) object to ts_plot(). We create a tsi object by
passing the hosp2 data set to the id_tsi() function, and we identify the names of the columns that contain
intervention cohorts, time periods, and the time period during which each cohort implemented the intervention. The
returned tsi object contains a data frame, which includes columns cohort, time, intervention_time, and
tsi. Each row corresponds to a unique county-time period. The tsi column contains an integer that tells us how many
time periods the corresponding time period is away from the intervention. For example, Ashbrook County belongs to cohort
5, which means the intervention took place in 2015. If we examine the first row of the tsi object, we'll see that the
time period 2010 is 5 time periods prior to the intervention, so tsi has a value of -5. Values of tsi are negative
before the intervention, 0 during the year of intervention, and positive for the years following the intervention.
hosp_tsi <- staggR::id_tsi(df = hosp2,
cohort_var = "county",
time_var = "yr",
intervention_var = "intervention_yr")
head(hosp_tsi$data, 11)
Now we can pass hosp_tsi to ts_plot() to generate a set of plots with time since intervention on the X axis instead of years.
county_plot <- staggR::ts_plot(hospitalized ~ county + yr,
df = hosp2,
intervention_var = "intervention_yr",
tsi = hosp_tsi)
county_plot
We might like to customize this plot a little bit and use better labels for the X and Y axes. The function ts_plot()
returns a ggplot object, so we can append any ggplot2 functions to customize the plot.
library(ggplot2)
county_plot +
scale_x_continuous(name = "Years",
breaks = seq(-8, 5, by = 2)) +
scale_y_continuous(name = "Percent hospitalized",
breaks = seq(0, 1, by = 0.2),
labels = function(x) paste0(round(x * 100), "%")) +
theme_classic() +
theme(panel.grid.major.x = element_line(linewidth = 0.5))
Fit a difference-in-differences regression model
Now that we have some familiarity with the data, we can fit a staggered difference-in-differences model. The function
for estimating such a model is called sdid().
The first argument of the function is the regression formula. The order of terms in the regression formula is important:
sdid() will assume that the first term on the right-hand side of the formula is the variable that identifies cohorts,
and the second term is the variable that identifies the time period for each observation. All remaining terms on the
right-hand side of the formula are covariates that we would like to use to adjust our model. The dependent variable is
on the left-hand side of the formula.
Other important parameters include intervention_var, which specifies the name of the column that contains the time
period during which each cohort implemented the intervention (in our case, this is intervention_yr), and df, which specifies the data set.
The call for estimating the staggered difference-in-differences regression, adjusting for age, sex, and comorbitiy,
is as follows:
sdid_hosp <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb,
df = hosp,
intervention_var = "intervention_yr")
summary(sdid_hosp)
Although we supplied a formula to fit this model, it is important to recognize that this formula is specified according
to a convenience syntax, which assumes that we want to fit a Sun and Abraham-style staggered difference-in-differences
model, rather than the very simple model implied by our formula. If we would like to see the formula passed to lm() to
fit the staggered difference-in-differences model, we can examine the components of the sdid_hosp object.
names(sdid_hosp)
sdid_hosp$formula
Here we see both the supplied and fitted formulas. The fitted formula includes main effects for each cohort, time
periods, and interaction terms for each cohort-time period combination. You might notice that these terms do not exist
in our supplied data.frame; this is because sdid() calls prep_data() to create all necessary dummy variables. These
are needed because we have to be very explicit about which levels of each cohort and time period variable are used as
referents in the model, and these referent levels differ for main effects and for each set of interaction terms.
By default, sdid() assumes that the first observed levels of both the cohort column and the time-period column should
function as the referents for the main effects. This is the same as with lm() and most other modeling functions in
R. For the long list of cohort-time period interactions, however, sdid() assumes that the time period immediately
preceding the intervention time period (defined by intervention_var) should be the referent. In our example, that
means that the interaction terms for cohort_5 should omit 2014 as the referent time period, because that cohort's
intervention took place in 2015. Similarly, the interaction terms for cohort_6 should omit 2015 as the referent time
period, because that cohort's intervention took place in 2016. Indeed, we observe that the formula has interaction terms
cohort_5:yr_2010 through cohort_5:yr_2020, with cohort_5_2014 omitted for identification.
We can use the cohort_time_refs parameter to specify a different time period to function as the referent for
cohort-time interactions:
sdid_hosp_earlier_timeref <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb,
df = hosp,
cohort_time_refs = list(`5` = "2013",
`6` = "2014",
`7` = "2015",
`8` = "2016"),
intervention_var = "intervention_yr")
We can now proceed with the second, post-estimation step: translating difference-in-differences estimates into average estimates.
Calculating coefficient averages
Before we can calculate coefficient averages, we first need to identify the coefficients we want to combine. We do this by supplying the sdid object and either "pre" or "post" to the select_period() function to retrieve the coefficients for interaction terms corresponding to the pre-intervention or post-intervention periods, respectively.
# Identify the coefficients that compose the effect of the intervention during the
# post-intervention period
(post_coefs <- staggR::select_period(sdid = sdid_hosp,
period = "post"))
Now that we've identified the coefficients of interest, we will pass the sdid object and the selected coefficients to the ave_coeff() function, which calculates averages of any specified set of regression coefficients, along with the corresponding standard errors.
To demonstrate how this function works, we first calculate the average effect of the intervention during the
post-intervention period.
staggR::ave_coeff(sdid = sdid_hosp,
coefs = post_coefs)
The result shows that the intervention is associated with an increase in the probability of hospitalization.
Next, we can use the ave_coeff() function to examine whether the risk of hospitalization differs by treatment group during the pre-intervention period. For that, we pass "pre" to the period argument, as follows:
# Identify the coefficients that represent the difference between exposure groups
# during the pre-intervention period
(pre_coefs <- staggR::select_period(sdid = sdid_hosp,
period = "pre"))
staggR::ave_coeff(sdid = sdid_hosp,
coefs = pre_coefs)
Here we see risk of hospitalization during the pre-intervention period does not significantly differ by intervention status.
We can also examine whether there exists heterogeneity in difference-in-differences
estimates between various cohorts, using the cohorts parameter in the select_period() function. For instance, we could calculate the average of difference-in-differences estimates for the post-intervention period for cohorts 5 and 6.
# Identify the coefficients for the effect of the intervention on cohorts 5 and 6
# during the post-intervention period
(post_coefs_cohort5 <- staggR::select_period(sdid = sdid_hosp,
period = "post",
cohorts = c("5", "6")))
staggR::ave_coeff(sdid = sdid_hosp,
coefs = post_coefs_cohort5)
We can similarly calculate the average of difference-in-differences estimates for the post-intervention period for
cohorts 7 and 8.
# Identify the coefficients for the effect of the intervention on cohorts 7 and 8
# during the post-intervention period
(post_coefs_cohort8 <- staggR::select_period(sdid = sdid_hosp,
period = "post",
cohorts = c("7", "8")))
staggR::ave_coeff(sdid = sdid_hosp,
coefs = post_coefs_cohort8)
Comparing these averages shows that the association between the intervention and probability of
hospitalization for the last two cohorts is almost twice as large compared to the first two cohorts.
If we want to aggregate a set of coefficients other than those for the pre- or post-intervention periods, we select our coefficients using select_terms() instead of select_period(). We supply a named list containing values for cohorts and either times or tsi. The cohorts element should contain a character vector of cohorts; if we omit the cohorts element, the function will select all available cohorts. The times element contains a character vector of time periods as they are recorded in the time variable supplied to sdid(). If we would prefer to select time periods relative to the intervention period, we would instead use the tsi element, which contains a vector of integers representing the number of units of time relative to each cohort's intervention.
For example, we could calculate the added risk of hospitalization associated with the intervention across all cohorts for the year 2018.
# Identify the coefficients corresponding to all cohorts in 2018
(terms_2018 <- staggR::select_terms(sdid = sdid_hosp,
selection = list(times = "2018")))
Then we can pass the selected terms to the ave_coeff() function to aggregate the 2018 effect across all cohorts.
staggR::ave_coeff(sdid = sdid_hosp,
coefs = terms_2018)
If we are interested in the added risk of hospitalization associated with the intervention in the year 2018, but only for the first two cohorts (5 and 6), we could include a cohorts element in the selection list.
(terms_2018_cohorts56 <- staggR::select_terms(sdid = sdid_hosp,
selection = list(cohorts = c("5", "6"),
times = "2018")))
staggR::ave_coeff(sdid = sdid_hosp,
coefs = terms_2018_cohorts56)
Alternatively, if we already know exactly which coefficients we want to aggregate, we can simply list them in the coefs argument.
(terms_custom <- staggR::select_terms(sdid = sdid_hosp,
coefs = c("cohort_5:yr_2018", "cohort_6:yr_2018")))
staggR::ave_coeff(sdid = sdid_hosp,
coefs = terms_custom)
Calculating custom standard errors
Because outcomes within the same unit (in our case, county) are often correlated over time, assuming independent errors
could understate uncertainty. Clustering standard errors allows for correlation within clusters while maintaining
independence across them. Following Sun and Abraham @sun2021did, we can supply a custom variance–covariance function to
compute standard errors clustered at the county level.
# Fit SDID model with standard errors clustered at the county level
sdid_hosp <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb,
df = hosp,
intervention_var = "intervention_yr",
.vcov = sandwich::vcovCL,
cluster = hosp$county)
summary(sdid_hosp)
# Examine the association between the intervention and risk of hospitalization in
# the post-intervention period
staggR::ave_coeff(sdid = sdid_hosp,
coefs = staggR::select_period(sdid = sdid_hosp,
period = "post"))
Using aggregated data
Data sets containing staggered interventions can sometimes be very large, so it might be preferable to fit the model
using aggregated data. To this end, we created an alternate version of the data set that has been aggregated to the
county-year level.
data("hosp_agg", package = "staggR")
head(hosp_agg)
# This data set contains one row per county-year
nrow(hosp_agg); nrow(unique(hosp_agg[,c("yr", "county")]))
As before, we can use the function sdid() to fit the difference-in-differences regression model. The regression formula
now uses the variable representing the percent of hospitalized individuals in each county-year (pct_hospitalized) as
the dependent variable, and the covariates have been aggregated (i.e., mean_age, pct_fem, and pct_cmb). The
main difference is that we need to specify the weights argument for the number of observations contained within each
row of the data set, which is represented by the n_enr column.
sdid_hosp_agg <- staggR::sdid(pct_hospitalized ~ cohort + yr + mean_age + pct_fem + pct_cmb,
df = hosp_agg,
weights = "n_enr",
intervention_var = "intervention_yr",
# Cluster standard errors at the county level
.vcov = sandwich::vcovCL,
cluster = hosp_agg$county)
Aggregating means that we decoupled values of these covariates from outcomes within county-years, so coefficients are
slightly different but substantially similar compared to the regression using individual-level data. We can examine differences between models by summarizing them, as we would with a lm or glm model.
# Summarize the original model
summary(sdid_hosp)
# Summarize the model based on aggregated data
summary(sdid_hosp_agg)
If we would like to compare the two models' estimates of risk of hospitalization associated with the intervention during the post-intervention period, we can use ave_coeff() again.
# Summarize the post-intervention period from the original model
summary_sdid_hosp <-
staggR::ave_coeff(sdid = sdid_hosp,
coefs = staggR::select_period(sdid = sdid_hosp,
period = "post"))
# # Summarize the post-intervention period from the model fitted using aggregated data
summary_sdid_hosp_agg <-
staggR::ave_coeff(sdid = sdid_hosp_agg,
coefs = staggR::select_period(sdid = sdid_hosp_agg,
period = "post"))
# Combine the two summaries
summary_sdid_hosp$model <- "Individual-level data"
summary_sdid_hosp_agg$model <- "Aggregated data"
rbind(summary_sdid_hosp, summary_sdid_hosp_agg)
Adjusting for non-parallel trends
A key assumption in any difference-in-differences analysis is that the outcome trends among intervention and comparison groups are parallel. If examination of the pre-intervention trends does not support this assumption, we can adjust outcomes to account for trend divergence.
We start by fitting a staggered DiD model, using the original outcomes. In our example, we can borrow the model sdid_hosp from above.
# First examine the model we fit using the original outcome
sdid_hosp$formula$supplied
# Calculate de-trending adjustments
hosp_det <- staggR::detrend(sdid = sdid_hosp,
df = hosp)
# Then refit the same model, substituting the _detrended version of the outcome
sdid_hosp_det <- staggR::sdid(hospitalized_detrended ~ cohort + yr + age + sex + comorb,
df = hosp_det,
intervention_var = "intervention_yr")
# Organize and present results, comparing them to the non-trend-adjusted model
original_rslts <- staggR::ave_coeff(sdid = sdid_hosp,
coefs = staggR::select_period(sdid = sdid_hosp,
period = "post"))
original_rslts$model = "Original"
original_rslts <- original_rslts[, c("model", "est", "se", "pval", "lb", "ub", "n")]
detrend_rslts <- staggR::ave_coeff(sdid = sdid_hosp_det,
coefs = staggR::select_period(sdid = sdid_hosp_det,
period = "post"))
detrend_rslts$model = "Trend-adjusted"
detrend_rslts <- detrend_rslts[, c("model", "est", "se", "pval", "lb", "ub", "n")]
rbind(original_rslts, detrend_rslts)
References
Try the staggR package in your browser
Any scripts or data that you put into this service are public.
staggR documentation built on April 2, 2026, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) # Check for sandwich package if (!requireNamespace("sandwich", quietly = TRUE)) { stop("Package 'sandwich' needed for this vignette to work. Please install it.") }
staggR: Fit difference-in-differences models with staggered interventions
The staggR package fits linear difference-in-differences models in scenarios where intervention roll-outs are
staggered over time. The package implements a version of an approach proposed by Sun and Abraham @sun2021did to
estimate cohort- and time-since-treatment specific difference-in-differences parameters. This package provides an
alternative approach that proceeds in two steps:
-
Estimate a difference-in-differences regression with cohort- and (calendar) time-specific parameters;
-
Translate difference-in-differences estimates into average estimates for arbitrary time periods relative to the intervention (e.g., post-intervention period).
By separating these two steps, we provide an easy-to-implement, transparent, intuitive difference-in-differences framework.
Background
Consider a scenario with $i=1,...,N$ units and $t=1,...,T$ periods. Units are exposed to the intervention at different times. We group units having the intervention during the same time period into cohorts, which we denote by $c=1,...,C$. To simplify the notation, we assume that units in cohort $c=1$ are exposed to treatment at time $t=1$, etc. Units exposed to the intervention at time $t>T$ form the comparison group because they are not exposed to the intervention during the study period.
Abraham and Sun's difference-in-differences specification interacts cohort indicators with time indicators relative to intervention exposure start, and might be written as follows:
$$ y_{it} = \alpha_c D_c + \beta_t D_t + \left(\sum_c \sum_{l\neq -1} \gamma_{cl} D_c \cdot D^l_{it}\right) + \psi X_{it} + \varepsilon_{st}, $$
where $D_c$ are cohort indicators equal to one if unit $i$ is in cohort $c$ and zero otherwise, $D_t$ are time indicators, and $D^l_{it}$ are indicators equal to one if unit $i$ is $l$ periods away from initial intervention exposure, $X_{it}$ are covariates, and $\varepsilon_{st}$ is the error term. The last pre-intervention period $l=-1$ is omitted to ensure identification. The parameters $\gamma_{cl}$ measure difference-in-differences estimates for a cohort $c$ and time since intervention $l$.
Our approach recognizes the following relationship between time since intervention $l$ and time $t$, which was also noted by Sun and Abraham:
$$ l = t-c $$
For instance, time $t=2$ for cohort $c=3$ is $l=-1$, i.e., the last period prior to intervention exposure, because intervention exposure for cohort $c$ is $t=3$.
Using this insight, the difference-in-differences regression can alternatively be stated as follows:
$$ y_{it} = \alpha_c D_c + \beta_t D_t + \left(\sum_c \sum_{t\neq (c-1)} \gamma_{ct} D_c \cdot D_t\right) + \psi X_{it} + \varepsilon_{st}. $$
Using this specification has two advantages:
- It clarifies that this staggered difference-in-differences method is a generalization of the basic difference-in-differences design with two periods (a pre-intervention and a post-intervention period) and two cohorts (a treatment group and a comparison group). In that setting, the indicator for the post-intervention period measures the change between the pre- and post-intervention period for the comparison group, and the difference-in-differences estimate measures the differential change between the pre- and post-intervention period for the treatment group.
The above equation preserves this structure, while showing how it is a more general version of the difference-in-differences design. Specifically, parameters $\beta_t$ measure changes over time for the comparison group, which now can include multiple periods; and parameters $\gamma_{ct}$ measure differential changes, but for multiple treatment cohorts instead of just one treatment group, and for multiple periods instead of just one.
- It can be easily implemented in R (or other statistical programs) because it does not require constructing cohort-specific time-since-intervention indicators $D^l_{it}$. Instead, it uses cohort and time indicators and the interaction between the two.
This packages provides functions to estimate such a difference-in-differences specification. It also provides functions
to translate difference-in-differences estimates into quantities that are useful for research purposes. For instance,
users can calculate average difference-in-differences estimates for each period relative to intervention exposure, or
they can calculate the average difference-in-differences estimates for the full post-intervention period. The
ave_coeff() function accomplishes this by calculating averages of user-specified sets of regression coefficients and
combining the corresponding standard errors using the estimated variance–covariance matrix, whose calculation method
users can also specify. This aggregation yields a single estimate that properly accounts for sampling variability and
the correlation among underlying event-study coefficients.
Example using simulated data
Explore data
Consider as a didactic example a policy intervention designed to reduce inpatient hospitalizations in 15 counties. This
longitudinal data set contains one row per individual-year. Each individual in the data set is identified by a globally
unique identifier (guid), and we have measures of the individuals' ages, sexes, and comorbidities, and an indicator
showing whether the individual was hospitalized during the current year. Each individual resides in a county, which
itself is grouped into cohorts based on intervention year. Finally, the yr column indicates the year of the current
observation.
wzxhzdk:1
The county is never exposed to the intervetnion if intervention_yr is NA. Among the 15 counties in the data set, 3
counties implemented the intervention in 2015; 2 counties implemented in 2016; 5 counties implemented in 2017; 1 county
implemented in 2018; and 4 counties did not implement the intervention at all during the study period.
with(hosp, table(county, intervention_yr, useNA = "ifany"))
Counties are organized into cohorts by intervention year, i.e., counties with the same intervention year are assigned to the same cohort.
county_cohorts <- unique(hosp[, c("county", "cohort", "intervention_yr")]) county_cohorts[order(county_cohorts$cohort), ]
The study period includes 11 years, from 2010 through 2020. Hospitalizations occurs in every county-year.
stats::xtabs(hospitalized ~ county + yr, data = hosp)
Visualize hospitalization trends
To visualize trends in the proportion of individuals who were hospitalized over time, we can call the function
ts_plot(), which takes, at minimum, 3 arguments: a formula, in the format outcome ~ group + time_var, the source data set, and the name of the column that contains the time period during which each group implemented the intervention.
staggR::ts_plot(hospitalized ~ county + yr, df = hosp, intervention_var = "intervention_yr")
Equivalently, we can specify the outcome, group, and time variable using the parameters y, group, and time_var.
staggR::ts_plot(y = "hospitalized", group = "county", time_var = "yr", df = hosp, intervention_var = "intervention_yr")
While this visualization is helpful, it excludes counties that did not implement the intervention. This is because the
function ts_plot() omits any counties where any of y, group, time_var, (or our formula equivalent, hospitalized ~ county + yr) or the intervention_var are NA.
table(hosp$hospitalized, useNA = "always") table(hosp$county, useNA = "always") table(hosp$yr, useNA = "always") table(hosp$intervention_yr, useNA = "always")
We can add these counties to the plot if we change NAs in intervention_yr to some other string value.
hosp2 <- hosp hosp2$intervention_yr <- ifelse(is.na(hosp2$intervention_yr), yes = "Comparison group", no = hosp2$intervention_yr) county_plot <- staggR::ts_plot(hospitalized ~ county + yr, df = hosp2, intervention_var = "intervention_yr") county_plot
It might be easier to compare the panels in this plot if the X axes were all aligned at the point of intervention. We
can accomplish this by passing a time-since-intervention (tsi) object to ts_plot(). We create a tsi object by
passing the hosp2 data set to the id_tsi() function, and we identify the names of the columns that contain
intervention cohorts, time periods, and the time period during which each cohort implemented the intervention. The
returned tsi object contains a data frame, which includes columns cohort, time, intervention_time, and
tsi. Each row corresponds to a unique county-time period. The tsi column contains an integer that tells us how many
time periods the corresponding time period is away from the intervention. For example, Ashbrook County belongs to cohort
5, which means the intervention took place in 2015. If we examine the first row of the tsi object, we'll see that the
time period 2010 is 5 time periods prior to the intervention, so tsi has a value of -5. Values of tsi are negative
before the intervention, 0 during the year of intervention, and positive for the years following the intervention.
hosp_tsi <- staggR::id_tsi(df = hosp2, cohort_var = "county", time_var = "yr", intervention_var = "intervention_yr") head(hosp_tsi$data, 11)
Now we can pass hosp_tsi to ts_plot() to generate a set of plots with time since intervention on the X axis instead of years.
county_plot <- staggR::ts_plot(hospitalized ~ county + yr, df = hosp2, intervention_var = "intervention_yr", tsi = hosp_tsi) county_plot
We might like to customize this plot a little bit and use better labels for the X and Y axes. The function ts_plot()
returns a ggplot object, so we can append any ggplot2 functions to customize the plot.
library(ggplot2) county_plot + scale_x_continuous(name = "Years", breaks = seq(-8, 5, by = 2)) + scale_y_continuous(name = "Percent hospitalized", breaks = seq(0, 1, by = 0.2), labels = function(x) paste0(round(x * 100), "%")) + theme_classic() + theme(panel.grid.major.x = element_line(linewidth = 0.5))
Fit a difference-in-differences regression model
Now that we have some familiarity with the data, we can fit a staggered difference-in-differences model. The function
for estimating such a model is called sdid().
The first argument of the function is the regression formula. The order of terms in the regression formula is important:
sdid() will assume that the first term on the right-hand side of the formula is the variable that identifies cohorts,
and the second term is the variable that identifies the time period for each observation. All remaining terms on the
right-hand side of the formula are covariates that we would like to use to adjust our model. The dependent variable is
on the left-hand side of the formula.
Other important parameters include intervention_var, which specifies the name of the column that contains the time
period during which each cohort implemented the intervention (in our case, this is intervention_yr), and df, which specifies the data set.
The call for estimating the staggered difference-in-differences regression, adjusting for age, sex, and comorbitiy, is as follows:
sdid_hosp <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb, df = hosp, intervention_var = "intervention_yr") summary(sdid_hosp)
Although we supplied a formula to fit this model, it is important to recognize that this formula is specified according
to a convenience syntax, which assumes that we want to fit a Sun and Abraham-style staggered difference-in-differences
model, rather than the very simple model implied by our formula. If we would like to see the formula passed to lm() to
fit the staggered difference-in-differences model, we can examine the components of the sdid_hosp object.
names(sdid_hosp) sdid_hosp$formula
Here we see both the supplied and fitted formulas. The fitted formula includes main effects for each cohort, time
periods, and interaction terms for each cohort-time period combination. You might notice that these terms do not exist
in our supplied data.frame; this is because sdid() calls prep_data() to create all necessary dummy variables. These
are needed because we have to be very explicit about which levels of each cohort and time period variable are used as
referents in the model, and these referent levels differ for main effects and for each set of interaction terms.
By default, sdid() assumes that the first observed levels of both the cohort column and the time-period column should
function as the referents for the main effects. This is the same as with lm() and most other modeling functions in
R. For the long list of cohort-time period interactions, however, sdid() assumes that the time period immediately
preceding the intervention time period (defined by intervention_var) should be the referent. In our example, that
means that the interaction terms for cohort_5 should omit 2014 as the referent time period, because that cohort's
intervention took place in 2015. Similarly, the interaction terms for cohort_6 should omit 2015 as the referent time
period, because that cohort's intervention took place in 2016. Indeed, we observe that the formula has interaction terms
cohort_5:yr_2010 through cohort_5:yr_2020, with cohort_5_2014 omitted for identification.
We can use the cohort_time_refs parameter to specify a different time period to function as the referent for
cohort-time interactions:
sdid_hosp_earlier_timeref <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb, df = hosp, cohort_time_refs = list(`5` = "2013", `6` = "2014", `7` = "2015", `8` = "2016"), intervention_var = "intervention_yr")
We can now proceed with the second, post-estimation step: translating difference-in-differences estimates into average estimates.
Calculating coefficient averages
Before we can calculate coefficient averages, we first need to identify the coefficients we want to combine. We do this by supplying the sdid object and either "pre" or "post" to the select_period() function to retrieve the coefficients for interaction terms corresponding to the pre-intervention or post-intervention periods, respectively.
# Identify the coefficients that compose the effect of the intervention during the # post-intervention period (post_coefs <- staggR::select_period(sdid = sdid_hosp, period = "post"))
Now that we've identified the coefficients of interest, we will pass the sdid object and the selected coefficients to the ave_coeff() function, which calculates averages of any specified set of regression coefficients, along with the corresponding standard errors.
To demonstrate how this function works, we first calculate the average effect of the intervention during the post-intervention period.
staggR::ave_coeff(sdid = sdid_hosp, coefs = post_coefs)
The result shows that the intervention is associated with an increase in the probability of hospitalization.
Next, we can use the ave_coeff() function to examine whether the risk of hospitalization differs by treatment group during the pre-intervention period. For that, we pass "pre" to the period argument, as follows:
# Identify the coefficients that represent the difference between exposure groups # during the pre-intervention period (pre_coefs <- staggR::select_period(sdid = sdid_hosp, period = "pre")) staggR::ave_coeff(sdid = sdid_hosp, coefs = pre_coefs)
Here we see risk of hospitalization during the pre-intervention period does not significantly differ by intervention status.
We can also examine whether there exists heterogeneity in difference-in-differences
estimates between various cohorts, using the cohorts parameter in the select_period() function. For instance, we could calculate the average of difference-in-differences estimates for the post-intervention period for cohorts 5 and 6.
# Identify the coefficients for the effect of the intervention on cohorts 5 and 6 # during the post-intervention period (post_coefs_cohort5 <- staggR::select_period(sdid = sdid_hosp, period = "post", cohorts = c("5", "6"))) staggR::ave_coeff(sdid = sdid_hosp, coefs = post_coefs_cohort5)
We can similarly calculate the average of difference-in-differences estimates for the post-intervention period for cohorts 7 and 8.
# Identify the coefficients for the effect of the intervention on cohorts 7 and 8 # during the post-intervention period (post_coefs_cohort8 <- staggR::select_period(sdid = sdid_hosp, period = "post", cohorts = c("7", "8"))) staggR::ave_coeff(sdid = sdid_hosp, coefs = post_coefs_cohort8)
Comparing these averages shows that the association between the intervention and probability of hospitalization for the last two cohorts is almost twice as large compared to the first two cohorts.
If we want to aggregate a set of coefficients other than those for the pre- or post-intervention periods, we select our coefficients using select_terms() instead of select_period(). We supply a named list containing values for cohorts and either times or tsi. The cohorts element should contain a character vector of cohorts; if we omit the cohorts element, the function will select all available cohorts. The times element contains a character vector of time periods as they are recorded in the time variable supplied to sdid(). If we would prefer to select time periods relative to the intervention period, we would instead use the tsi element, which contains a vector of integers representing the number of units of time relative to each cohort's intervention.
For example, we could calculate the added risk of hospitalization associated with the intervention across all cohorts for the year 2018.
# Identify the coefficients corresponding to all cohorts in 2018 (terms_2018 <- staggR::select_terms(sdid = sdid_hosp, selection = list(times = "2018")))
Then we can pass the selected terms to the ave_coeff() function to aggregate the 2018 effect across all cohorts.
staggR::ave_coeff(sdid = sdid_hosp, coefs = terms_2018)
If we are interested in the added risk of hospitalization associated with the intervention in the year 2018, but only for the first two cohorts (5 and 6), we could include a cohorts element in the selection list.
(terms_2018_cohorts56 <- staggR::select_terms(sdid = sdid_hosp, selection = list(cohorts = c("5", "6"), times = "2018"))) staggR::ave_coeff(sdid = sdid_hosp, coefs = terms_2018_cohorts56)
Alternatively, if we already know exactly which coefficients we want to aggregate, we can simply list them in the coefs argument.
(terms_custom <- staggR::select_terms(sdid = sdid_hosp, coefs = c("cohort_5:yr_2018", "cohort_6:yr_2018"))) staggR::ave_coeff(sdid = sdid_hosp, coefs = terms_custom)
Calculating custom standard errors
Because outcomes within the same unit (in our case, county) are often correlated over time, assuming independent errors could understate uncertainty. Clustering standard errors allows for correlation within clusters while maintaining independence across them. Following Sun and Abraham @sun2021did, we can supply a custom variance–covariance function to compute standard errors clustered at the county level.
# Fit SDID model with standard errors clustered at the county level sdid_hosp <- staggR::sdid(hospitalized ~ cohort + yr + age + sex + comorb, df = hosp, intervention_var = "intervention_yr", .vcov = sandwich::vcovCL, cluster = hosp$county) summary(sdid_hosp) # Examine the association between the intervention and risk of hospitalization in # the post-intervention period staggR::ave_coeff(sdid = sdid_hosp, coefs = staggR::select_period(sdid = sdid_hosp, period = "post"))
Using aggregated data
Data sets containing staggered interventions can sometimes be very large, so it might be preferable to fit the model using aggregated data. To this end, we created an alternate version of the data set that has been aggregated to the county-year level.
data("hosp_agg", package = "staggR") head(hosp_agg) # This data set contains one row per county-year nrow(hosp_agg); nrow(unique(hosp_agg[,c("yr", "county")]))
As before, we can use the function sdid() to fit the difference-in-differences regression model. The regression formula
now uses the variable representing the percent of hospitalized individuals in each county-year (pct_hospitalized) as
the dependent variable, and the covariates have been aggregated (i.e., mean_age, pct_fem, and pct_cmb). The
main difference is that we need to specify the weights argument for the number of observations contained within each
row of the data set, which is represented by the n_enr column.
sdid_hosp_agg <- staggR::sdid(pct_hospitalized ~ cohort + yr + mean_age + pct_fem + pct_cmb, df = hosp_agg, weights = "n_enr", intervention_var = "intervention_yr", # Cluster standard errors at the county level .vcov = sandwich::vcovCL, cluster = hosp_agg$county)
Aggregating means that we decoupled values of these covariates from outcomes within county-years, so coefficients are
slightly different but substantially similar compared to the regression using individual-level data. We can examine differences between models by summarizing them, as we would with a lm or glm model.
# Summarize the original model summary(sdid_hosp) # Summarize the model based on aggregated data summary(sdid_hosp_agg)
If we would like to compare the two models' estimates of risk of hospitalization associated with the intervention during the post-intervention period, we can use ave_coeff() again.
# Summarize the post-intervention period from the original model summary_sdid_hosp <- staggR::ave_coeff(sdid = sdid_hosp, coefs = staggR::select_period(sdid = sdid_hosp, period = "post")) # # Summarize the post-intervention period from the model fitted using aggregated data summary_sdid_hosp_agg <- staggR::ave_coeff(sdid = sdid_hosp_agg, coefs = staggR::select_period(sdid = sdid_hosp_agg, period = "post")) # Combine the two summaries summary_sdid_hosp$model <- "Individual-level data" summary_sdid_hosp_agg$model <- "Aggregated data" rbind(summary_sdid_hosp, summary_sdid_hosp_agg)
Adjusting for non-parallel trends
A key assumption in any difference-in-differences analysis is that the outcome trends among intervention and comparison groups are parallel. If examination of the pre-intervention trends does not support this assumption, we can adjust outcomes to account for trend divergence.
We start by fitting a staggered DiD model, using the original outcomes. In our example, we can borrow the model sdid_hosp from above.
# First examine the model we fit using the original outcome sdid_hosp$formula$supplied # Calculate de-trending adjustments hosp_det <- staggR::detrend(sdid = sdid_hosp, df = hosp) # Then refit the same model, substituting the _detrended version of the outcome sdid_hosp_det <- staggR::sdid(hospitalized_detrended ~ cohort + yr + age + sex + comorb, df = hosp_det, intervention_var = "intervention_yr") # Organize and present results, comparing them to the non-trend-adjusted model original_rslts <- staggR::ave_coeff(sdid = sdid_hosp, coefs = staggR::select_period(sdid = sdid_hosp, period = "post")) original_rslts$model = "Original" original_rslts <- original_rslts[, c("model", "est", "se", "pval", "lb", "ub", "n")] detrend_rslts <- staggR::ave_coeff(sdid = sdid_hosp_det, coefs = staggR::select_period(sdid = sdid_hosp_det, period = "post")) detrend_rslts$model = "Trend-adjusted" detrend_rslts <- detrend_rslts[, c("model", "est", "se", "pval", "lb", "ub", "n")] rbind(original_rslts, detrend_rslts)
References
Try the staggR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.