ovbias() function

To use the ovbias() function, the user will need to run first collect the parameters from the short, intermediate and auxiliary regressions, using the collect_par() function and then feed it into the ovbias() function, along with the following:

• deltalow: lower limit of delta (e.g. 0.01)
• deltahigh: upper limit of delta (e.g. 0.99)
• Rhigh: upper limit of Rmax (e.g. 0.61)
• e: step size in defining the grid (e.g. 0.01)

In using the collect_par() function, the user will need to specify the following:

• data: the name of the data frame (e.g. NLSY_IQ)
• outcome: name of the outcome variable in double quotes (e.g. “iq_std”)
• treatment: name of the treatment variable in double quotes (e.g. “BF_months”)
• control: names of additional regressors to include in the intermediate regression, supplied as a vector (e.g. c(“age”,“sex”,“income”,“motherAge”,“motherEDU”,“mom_married”,“race”))
• other_regressors: names of regressors in the short regression, other than the treatment variable, supplied as a vector (e.g. c(“sex”,“age”))

The output of the ovbias() function is a list of three elements.

• Data: A data frame containing the bias (bias) and bias-adjusted treatment effect (bstar) for each point on the grid.
• bias_Distribution: Quantiles (2.5,5.0,50,95,97.5) of the empirical distribution of bias.
• bstar_Distribution: Quantiles (2.5,5.0,50,95,97.5) of the empirical distribution of the BATE (bstar).

ovbias_par() function

To use the ovbias_par() function, the user needs to specify the following:

• data: the name of the data frame (e.g. NLSY_IQ)
• outcome: name of the outcome variable in double quotes (e.g. “iq_std”)
• treatment: name of the treatment variable in double quotes (e.g. “BF_months”)
• control: names of additional regressors to include in the intermediate regression, supplied as a vector (e.g. c(“age”,“sex”,“income”,“motherAge”,“motherEDU”,“mom_married”,“race”))
• other_regressors: names of regressors in the short regression, other than the treatment variable, supplied as a vector (e.g. c(“sex”,“age”))
• deltalow: lower limit of delta (e.g. 0.01)
• deltahigh: upper limit of delta (e.g. 0.99)
• Rhigh: upper limit of Rmax (e.g. 0.61)
• e: step size in defining the grid (e.g. 0.01)

The output of the ovbias_par() function is identically same with the output of the ovbias() function.

ovbias_lm() function

To use the ovbias_lm() function, the user needs to specify three lm objects that capture the short, intermediate and auxiliary regressions:

• lm_shrt: lm object for the short regression
• lm_int: lm object for the intermediate regression
• lm_aux: lm object for the auxiliary regression
• deltalow: lower limit of delta (e.g. 0.01)
• deltahigh: upper limit of delta (e.g. 0.99)
• Rhigh: upper limit of Rmax (e.g. 0.61)
• e: step size in defining the grid (e.g. 0.01)

The output of the ovbias_lm() function is identically same with the output of the ovbias() function.

Additional Functions

Using the output from ovbias(), ovbias_par() or ovbias_lm(), users can construct various plots to visualize the results:

• cplotbias(): contour plot of the bias over the bounded box; the output of this function is a plot object;
• dplotbate(): histogram and density plot of BATE; the output of this function is a plot object;

The methodology proposed in Basu (2022) is slightly different from, and a critique of, Oster (2019). Hence, it might be useful to compare the results of the two methodologies. The methodology proposed in Oster (2019) is implemented via these functions:

• osterbds(): identified sets according to Oster’s methodology; the output of this function is a data frame;
• osterdelstar(): the value of δ^* for a chosen value of R_max; the output of this function is a data frame;
• delfplot(): a plot of the graph of the function, δ = f(R_max); the output of this function is a plot object.

Setting Up

Let us load the data set.

data("NLSY_IQ")

The data set has two .RData objects: NLSY_IQ (to be used for the analysis of maternal behavior on child IQ) and NLSY_BW (to be used for the analysis of maternal behavior on child birthweight).

In this example, we will analyse the effect of maternal behavior on child IQ scores. Let us start out by seeing the names of the variables in the NLSY_IQ data set.

names(NLSY_IQ)
#>  [1] "iq_std"             "BF_months"          "mom_drink_preg_all"
#>  [4] "lbw_preterm"        "age"                "female"            
#>  [7] "black"              "motherAge"          "motherEDU"         
#> [10] "mom_married"        "income"             "sex"               
#> [13] "race"

For use in the example below, let us set age and race as factor variables

NLSY_IQ$age <- factor(NLSY_IQ$age)
NLSY_IQ$race <- factor(NLSY_IQ$race)

Let us use the vtable package to see the summary statistics of the variables in our data set.

library(vtable)
#> Loading required package: kableExtra
vtable::st(NLSY_IQ)

Setting Up the Analysis

We will work with an example where the effect of months of breastfeeding on children’s IQ score is analyzed. Thus, here, the outcome variable is a child’s IQ score and the treatment variable is the months of breastfeeding by the mother. Additional control variables are: sex and age of the child, and the mother’s age, the mother’s years of formal education, whether the mother is married and the race of the mother. For further details of this example, see section 4.2 in Oster (2019). For ease of reference, let us note that we are working with the model reported in the first row of Table 3 in Oster (2019) and the first block of 4 rows in Table 2 in Basu (2022).

Using the names of variables in the data set, we have: * short regression: iq_std ~ BF_months + sex + age * intermediate regression: iq_std ~ BF_months + sex + age + income + motherAge + motherEDU + mom_married + race.

Let us use the collect_par() function to collect parameters from the short, intermediate and auxiliary regressions. It is important to note that other_parameters option in this function should refer to a subset of control. If the researcher fails to ensure this, the collect_par() function will throw an error.

parameters <- bate::collect_par(data=NLSY_IQ,
            outcome="iq_std",
            treatment="BF_months",
            control=c("age","sex","income","motherAge","motherEDU","mom_married","race"),
            other_regressors = c("sex","age"))

Let us see the parameters by looking at the object parameters that we used to store the output of the collect_par() function.

(parameters)
#>                beta0         R0  betatilde   Rtilde    sigmay   sigmax     taux
#> BF_months 0.04447926 0.04465201 0.01740748 0.255621 0.9900242 4.629618 18.99883

Our next task is to choose the dimensions of the bounded box over which we want the bias computation to be carried out. It is here that the researcher needs to use domain knowledge to set limits over which δ and R_max can run.

# Upper bound of Rmax
Rhigh <- 0.61
# Lower bound of delta
deltalow <- 0.01
# Upper bound of delta
deltahigh <- 0.99
# step size to construct grid
e <- 0.01

Note that while setting the dimensions of the bounded box, we have not chosen a value for Rlow (lower bound for R_max). This is because Rlow is chosen by default to be equal to parameters$Rtilde (the R-squared in the intermediate regression).

Let us see the division of the bounded box into the URR (unique real root) and NURR (nonunique real root) regions using the urrplot() function.

bate::urrplot(parameters = parameters, deltalow = deltalow, 
              deltahigh = deltahigh, Rlow = parameters$Rtilde,
              Rhigh = 0.61, e=0.01)

Using ovbias()

Now we can use the ovbias() function to compute the empirical distribution of omitted variable bias and BATE. Note that this step make take a few minutes, depending on the dimensions of the box and the size of e, to complete itself. As the function works, it will print a message in the console informing the user of the progress of the computation. Here, I have suppressed these messages.

OVB <- bate::ovbias(
  parameters = parameters,
  deltalow=deltalow, 
  deltahigh=deltahigh,
  Rhigh=Rhigh, 
  e=e)

Once the computation is completed, we can see the quantiles of omitted variable bias

(OVB$bias_Distribution)
#>  2.5%    5%   50%   95% 97.5% 
#> 0.000 0.000 0.009 0.034 0.039

and quantiles of the BATE (computed over the bounded box we chose above).

(OVB$bstar_Distribution)
#>   2.5%     5%    50%    95%  97.5% 
#> -0.021 -0.017  0.009  0.017  0.017

We can create a contour plot of BATE over the bounded box using the cplotbias() function.

cplotbias(OVB$Data)

We can create the histogram and density plot of the β^*, the bias-adjusted treatment effect using the dplotbate() function.

dplotbate(OVB$Data)
#> Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(density)` instead.
#> ℹ The deprecated feature was likely used in the bate package.
#>   Please report the issue to the authors.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Comparing Our Results with Oster (2019)

Let us compare our results with the methods proposed in Oster (2019). The methodology proposed in Oster (2019) relies on computing two things: (a) identified sets, and (b) δ^*.

Let us compute the identified sets.

bate::osterbds(parameters = parameters, Rmax=0.61)
#>       Discriminant          Interval1          Interval2 
#> "20.3210250174806"  "[0.0174,0.3752]" "[-0.0337,0.0174]"

The output from the above function contains three things: the value of the discriminant of the quadratic equation that is solved to compute the identified sets, and the two identified sets. The identified set is the interval formed by β̃ (treatment effect in the intermediate regression) and β^* = β̃ − ν, where ν is a root of a quadratic equation.

It has been shown in Proposition 5 in Basu (2022) that the discriminant of this quadratic will always be positive. Hence, there will be two real roots of the quadratic. Hence, there will be two identified sets, instead of a unique identified set. That is why we see two identified sets in the result above.

Let us also compute δ^*, the value of δ that is associated with a given value of Rmax (in this case Rmax=0.61) such that the treatment effect is zero.

bate::osterdelstar(parameters = parameters, Rmax=0.61)
#>          deltastar      discontinuity              slope 
#> "0.36619461884058"            "FALSE"         "Negative"

In addition to the value of δ^*, the output has two other things:

• discontinuity: this tells us whether the interval formed by R̃ and 1 contains a point of discontinuity; if it is FALSE, then the interval does not contain the point of discontinuity; if it is TRUE, then the interval contains the point of discontinuity and the analysis of Oster (2019) should be avoided; for more details see Section 5.2 in Basu (2022).

• slope: this gives the slope of the graph of the function δ = f(R_max); the slope can be either positive or negative and helps in interpreting the meaning of δ^*.

The value of δ^* printed above is just one value picked up from the function δ = f(R_max), for the specific choice of R_max. To see the graph of this function, we can use the function delfplot().

bate::delfplot(parameters = parameters)

The red vertical dashed line in the plot identifies the value of R_max that corresponds to δ = 1. In this example, this value of R_max is 0.38. This means that if a researcher uses any value of R_max that is greater than 0.38, she will get a value of δ^* that is less than 1, and if she uses a value of R_max that is less than 0.38, she will get a value of δ^* that is greater than 1.

Using the ovbias_par() function

We could have carried out the same analysis using the ovbias_par() function.

OVB.par <- ovbias_par(
  data=NLSY_IQ,
  outcome="iq_std",
  treatment="BF_months",
  control=c("age","sex","income","motherAge","motherEDU","mom_married","race"),
  other_regressors = c("sex","age"),
  deltalow=deltalow, 
  deltahigh=deltahigh,
  Rhigh=Rhigh, 
  e=e)

We can now see the quantiles of omitted variable bias

(OVB.par$bias_Distribution)
#>  2.5%    5%   50%   95% 97.5% 
#> 0.000 0.000 0.009 0.034 0.039

and quantiles of the BATE (computed over the bounded box we chose above).

(OVB.par$bstar_Distribution)
#>   2.5%     5%    50%    95%  97.5% 
#> -0.021 -0.017  0.009  0.017  0.017

Using the ovbias_lm() function

We could have also carried out the same analysis using the ovbias_lm() function. To use this function, we need to estimate the short regression

reg_col1 <- lm(
  iq_std ~ BF_months + factor(age) + sex,
  data = NLSY_IQ
)

and the intermediate regression

reg_col2 <- lm(
  iq_std ~ BF_months + factor(age) + sex +
    income + motherAge + motherEDU + mom_married +
    factor(race),
  data = NLSY_IQ
)

and the auxiliary regression

reg_aux <- lm(
  BF_months ~ factor(age) + sex +
    income + motherAge + motherEDU + mom_married +
    factor(race),
  data = NLSY_IQ
)

and then call the ovbias_lm() function and provide the three lm objects created above

OVB.lm <- ovbias_lm(
  lm_shrt = reg_col1,
  lm_int = reg_col2,
  lm_aux = reg_aux,
  deltalow=deltalow, 
  deltahigh=deltahigh,
  Rhigh=Rhigh, 
  e=e)

We can now see the quantiles of omitted variable bias

(OVB.lm$bias_Distribution)
#>  2.5%    5%   50%   95% 97.5% 
#> 0.000 0.000 0.009 0.034 0.039

and quantiles of the BATE (computed over the bounded box we chose above).

(OVB.lm$bstar_Distribution)
#>   2.5%     5%    50%    95%  97.5% 
#> -0.021 -0.017  0.009  0.017  0.017