Example: Smoking and Homocysteine Levels
In this vignette, we demonstrate the use of the natstrat package with a motivating example of the effect of smoking on homocysteine levels using data from NHANES 2005-2006 that is provided in the package.
library(natstrat)
#> Loading required package: caret
#> Loading required package: ggplot2
#> Loading required package: lattice
data('nh0506')Creating strata
Our first step is to choose factors for our strata. We choose to stratify by age and sex, cutting age into three categories:
age_cat <- cut(nh0506$age,
breaks = c(19, 39, 50, 85),
labels = c('< 40 years', '40 - 50 years', '> 50 years'))
strata <- age_cat : nh0506$sexWe can now see how many treated and control units are in each of our strata. Our treated units are composed of daily smokers, whereas control units have never smoked. We look at the ratio of controls to treated units available in each stratum:
Selecting a ratio of controls to treated
Natural strata are defined by having a fixed ratio of controls to treated units within each stratum, so we now select the fixed ratio of controls we will strive for. As the minimum ratio of controls to treated units available is 2.3, we set a fixed ratio of 1 for our natural strata. Note that we can select a noninteger ratio, such as 1.5 controls for every treated unit. However, if we select too large a ratio compared to the minimum available ratio, it will become difficult to balance covariates in this small dataset.
If we had multiple control or treated groups, as in nh0506_3groups, we could also use natural strata, specifying a ratio for each group relative to the primary treated group.
Generating covariate balance constraints
Our data contains 6 covariates, sex, age, race, education, poverty ratio, and bmi, that we’d like to be balanced between the control and treated groups. Sex is already balanced through the definition of the strata. While the strata force balance of our coarse age groups, we still seek to balance age more finely. We thus seek to balance age and our additional four covariates both across and within strata.
To do this, we generate balance constraints using the function generate_constraints() in the package. We must specify balance_formulas, z, and data.
balance_formulas: Here we specify a list of formulas that demonstrate the ways in which we want to balance. These formulas are extensions toRformulas, which follow many of the same conventions but have added functionality. Specifically, we include multiple terms on the left hand side, specifying the covariates we’d like to balance. On the right hand side, we explicitly specify1which indicates to balance these covariates across the whole population and is included implicitly unless- 1or+ 0are specified. We also include on the right hand side ourstratavariable, indicating that we’d like to balance the covariates within each level ofstrataas well. To learn more about how to specify balance formulas, skip to the Balance formula syntax section below.z: Here we specify the 0/1 vector indicating control/treated status.data: Here we specify thedata.framethat contains some or all of the covariates included in ourbalance_formulas.
There are many additional arguments which can be adjusted to control the desired covariate balance. These can be found in the Additional generate_constraints() arguments section below. Here, we leave all additional arguments as the default.
constraints <- generate_constraints(
balance_formulas = list(age + race + education + povertyr + bmi ~ 1 + strata),
z = nh0506$z,
data = nh0506)
names(constraints)
#> [1] "X" "importances"
dim(constraints$X)
#> [1] 2928 105
length(constraints$importances)
#> [1] 105We see that the return of generate_constraints() is a list with two elements.
The first element, X, has a row for each unit and a column for each constraint. We see there are 105 constraints generated, which is much more than the 5 covariates entered on the left hand side of balance_formulas. This is because those 5 covariates imply 15 variables to balance: age, 5 levels of race, 5 levels of education, poverty ratio, missing poverty ratio, bmi, and missing bmi. Each of these 15 variables is balanced across the population and within each of the 6 strata. This leads to \(15 \cdot (1 + 6) = 105\) balance constraints.
The second element, importances, contains an entry for each of the 105 constraints, stating how important it is to balance that constraint. This plays a bigger role as you do more fine tuning as illustrated in the Balance formula syntax and Additional generate_constraints() arguments sections.
Optimizing the control group
We now run the linear program and randomized rounding that will choose the best controls to balance our covariates. We input 5 arguments:
z: Treatment vectorX: Constraint matriximportances: Named vector of the importance of balancing each of the constraints inXst: Strata vectorratio: Desired ratio of controls
Additional arguments can be found in in the Additional optimize_controls() arguments section. One notable argument is solver, which is set by default to ‘Rglpk’. The other solver implemented is ‘gurobi’, and is highly recommended for use on larger data sets if available to you.
results <- optimize_controls(z = nh0506$z,
X = constraints$X,
importances = constraints$importances,
st = strata,
ratio = 1)
names(results)
#> [1] "objective" "objective_wo_importances"
#> [3] "eps" "eps_star"
#> [5] "importances" "weight_star"
#> [7] "selected" "selected_star"
#> [9] "pr" "pr_star"
#> [11] "rrdetails" "lpdetails"We see the return of optimize_controls() is a list. The most useful element is selected, which is a boolean vector stating whether each unit should be included in the analysis or not. This includes both treated and control units. The rest of the elements in the list are useful in analyzing the performance of the linear program and randomized rounding solutions, and are described in the help page for the function.
This package also supports optimizing for multiple comparisons simultaneously; see the Additional optimize_controls() arguments section to learn more.
Checking balance achieved
We now see how well we have balanced our covariates of interest, specified in the cov_data dataframe below. The check_balance function below will generate standardized differences for all variables in the specified X dataframe, whether or not they were included in the original optimization. The selected argument simply takes the selected element of the results of the optimization in the last step, and the plot argument is a boolean stating whether to generate ggplot objects containing love plots demonstrating balance.
cov_data <- nh0506[, c('sex', 'age', 'race', 'education', 'povertyr', 'bmi')]
stand_diffs <- check_balance(z = nh0506$z,
X = cov_data,
st = strata,
selected = results$selected,
plot = TRUE)
#> Excellent balance achieved across strata (<.1 standardized difference in all covariates).
#> The mean standardized difference across strata for the covariates is 0.017,
#> the median is 0.004,
#> and the maximum is 0.072.
#> Taking a weighted average of standardized differences within strata,
#> we have the mean covariate imbalance is 0.024,
#> the median is 0.009,
#> and the maximum is 0.085.
#> Within strata, we see imbalances (>=.2 standardized difference in some covariates in some strata).
#> The mean standardized difference within strata for the covariates is 0.024,
#> the median is 0,
#> and the maximum is 0.268.We see below that there is excellent balance across strata as all standardized differences in means after refining the control group fall below the first dashed line at 0.1. We also see that when we take the weighted average of the standardized differences in means within strata, where the average is weighted by stratum size, these values also fall below 0.1.
stand_diffs$plot_pair
We also take a look within each of the 6 strata, and find only two cases of slight imbalance where the standardized differences in means fall outside the second dashed line at 0.2: the ‘Mexican American’ category for males less than 40 years old, and the ‘9-11th grade’ category for females 40-50 years old.
stand_diffs$plot_strata
#> $`< 40 years:Male`
#>
#> $`< 40 years:Female`
#>
#> $`40 - 50 years:Male`
#>
#> $`40 - 50 years:Female`
#>
#> $`> 50 years:Male`
#>
#> $`> 50 years:Female`
As the imbalances are not too severe, we will accept this balance. However, there are many ways we could fine-tune the balance by shifting the arguments specified in the Further details sections.
Outcome analysis
Once you have balanced treated and control groups formed by taking the units with TRUE values in results$selected, you can use these to analyze your outcomes. Recall that since the individuals in these groups are not weighted and there is covariate balance both within and across strata, it is very easy to make plots and conduct analyses.
For example, we can make boxplots of the log of the homocysteine levels in the two groups both within each stratum and across the strata:
nh0506$group <- factor(nh0506$z, levels = c(0, 1),
labels = c("Never smoker", "Daily smoker"))
nh0506$strata <- strata
ggplot(nh0506[results$selected, ], aes(x = group, y = log(homocysteine))) +
geom_boxplot() +
facet_wrap(~ strata)
ggplot(nh0506[results$selected, ], aes(x = group, y = log(homocysteine))) +
geom_boxplot()
We can also calculate the average treatment effect on the treated (ATT) in terms of log of homocysteine, both within the strata and across:
att <- mean(log(nh0506$homocysteine[results$selected & nh0506$z == 1])) -
mean(log(nh0506$homocysteine[results$selected & nh0506$z == 0]))
cat(paste0("\nThe ATT across the strata is ", round(att, 3), " on the log scale,
which converted back to a regular scale is ", round(exp(att), 3), " umol/L."))
#>
#> The ATT across the strata is 0.127 on the log scale,
#> which converted back to a regular scale is 1.135 umol/L.
for (st in levels(nh0506$strata)) {
att <- mean(log(nh0506$homocysteine[results$selected & nh0506$z == 1 & nh0506$strata == st])) -
mean(log(nh0506$homocysteine[results$selected & nh0506$z == 0 & nh0506$strata == st]))
cat(paste0("\nThe ATT for stratum ", st, " is ", round(att, 3), " on the log scale,
which converted back to a regular scale is ", round(exp(att), 3), " umol/L."))
}
#>
#> The ATT for stratum < 40 years:Male is 0.078 on the log scale,
#> which converted back to a regular scale is 1.081 umol/L.
#> The ATT for stratum < 40 years:Female is 0.237 on the log scale,
#> which converted back to a regular scale is 1.268 umol/L.
#> The ATT for stratum 40 - 50 years:Male is 0.045 on the log scale,
#> which converted back to a regular scale is 1.046 umol/L.
#> The ATT for stratum 40 - 50 years:Female is 0.213 on the log scale,
#> which converted back to a regular scale is 1.237 umol/L.
#> The ATT for stratum > 50 years:Male is 0.077 on the log scale,
#> which converted back to a regular scale is 1.08 umol/L.
#> The ATT for stratum > 50 years:Female is 0.162 on the log scale,
#> which converted back to a regular scale is 1.175 umol/L.We see a similar effect of smoking increasing homocysteine levels across all the strata, but the effect is greater for females of all ages than for males.