Richardson et al. (2023): Bespoke Instrumental Variables
The following is a replication of Richardson et al. (2023), which considered the use of Bespoke Instrumental Variable Analysis (BSIV) for estimation of the effect of a treatment under one-sided nonadherence. The empirical example in that paper comes from the Obstetrics and Periodontal Therapy Study. This data can be obtained from https://www.causeweb.org/tshs/category/dataset/.
As in the paper, we consider the effect of nonsurgical periodontal treatment during pregnancy on the fraction of gingival sites bleeding on probing at 29–32 weeks’ gestation. Besides the proposed BSIV approach, the authors considered others as well.
Setup
[22]:
import numpy as np
import pandas as pd
import delicatessen as deli
from delicatessen import MEstimator
from delicatessen.estimating_equations import ee_regression, ee_2sls, ee_gestimation_snmm_iv
from delicatessen.sandwich import delta_method
print("Versions")
print("NumPy: ", np.__version__)
print("pandas: ", pd.__version__)
print("Delicatessen: ", deli.__version__)
Versions
NumPy: 2.3.5
pandas: 2.3.3
Delicatessen: 4.2
[6]:
d = pd.read_csv("data/obp_study.csv")
d['C'] = 1
d.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 640 entries, 0 to 639
Data columns (total 7 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 Unnamed: 0 640 non-null int64
1 R 640 non-null int64
2 A 640 non-null int64
3 Y 640 non-null float64
4 L1 640 non-null float64
5 L2 640 non-null float64
6 C 640 non-null int64
dtypes: float64(3), int64(4)
memory usage: 35.1 KB
Intent-to-Treat
As a starting point, we can estimate the Intent-to-Treat (ITT) effect. As this was a randomized trial, we can do this by simply fitting a linear regression model for the outcome (Y) given the assignment (R). The following code does this process
[7]:
def psi_itt(theta):
return ee_regression(theta=theta, X=d[['C', 'R']], y=d['Y'], model='linear')
[8]:
estr = MEstimator(psi_itt, init=[0, 0])
estr.estimate()
[9]:
estr.print_results(subset=[1, ])
==============================================================
Estimation Method: M-estimator
--------------------------------------------------------------
No. Observations: 640 | No. Parameters: 2
Solving algorithm: lm | Max Iterations: 5000
Solving tolerance: 1e-09 | Allow P-Inverse: 1
Derivative Method: approx | Deriv Approx: 1e-09
Small N Correction: None | Distribution: Z-stat
==============================================================
Theta StdErr Z-score LCL UCL P-value S-value
--------------------------------------------------------------
-0.24 0.02 -14.65 -0.27 -0.21 0.00 159.11
==============================================================
These results match those reported in the paper (Table 2; -0.24, 95% CI: -0.27, -0.21). As there is non-adherence, this estimate corresponds to the effect of being assigning periodontal treatment, not the effect of treatment itself.
Instrumental Variable
An alternative approach is to leverage the randomization as an instrumental variable for adherence. From this, we can estimate the local average causal effect (also known as the complier average causal effect). In the special case where the control group cannot access the treatment (for example a novel treatment that is not available outside of the trial), the local average causal effect is equal to the effect of treatment on the treated randomized to treatment.
The following code applies 2-stage least squares regression to estimate the local average causal effect.
[13]:
def psi_2sls(theta):
return ee_2sls(theta=theta, y=d['Y'], A=d['A'], Z=d[['R', ]],
W=d[['C', ]])
[14]:
estr = MEstimator(psi_2sls, init=[0, 0, 0, 0])
estr.estimate()
[16]:
estr.print_results(subset=[0, ])
==============================================================
Estimation Method: M-estimator
--------------------------------------------------------------
No. Observations: 640 | No. Parameters: 4
Solving algorithm: lm | Max Iterations: 5000
Solving tolerance: 1e-09 | Allow P-Inverse: 1
Derivative Method: approx | Deriv Approx: 1e-09
Small N Correction: None | Distribution: Z-stat
==============================================================
Theta StdErr Z-score LCL UCL P-value S-value
--------------------------------------------------------------
-0.48 0.04 -11.17 -0.57 -0.40 0.00 93.79
==============================================================
These results match the book (with a slight difference in the confidence intervals due to different choices of variance estimators).
One of the key assumptions for the previous IV approach is the so-called exclusion restriction assumption. This assumption stipulates that the randomization has no effect on the outcome except through the treatment. However, this assumption might be suspect, particularly when the treatment cannot be double blinded. In this example, the participants could not be blinded to periodontal treatment. Therefore, whether the exclusion restriction assumption is true is suspect. This is the problem that the BSIV approach seeks to address.
Bespoke Instrumental Variable
The BSIV approach operates by decomposes the effect of assignment and the effect of treatment into separate parameters. Briefly, this approach operates by using the random assignment and lack of access to the treatment to turn a measured covariate (related to the outcome) into a “bespoke instrumental variable”. See the paper for further details and assumptions.
[19]:
y = np.asarray(d['Y'])
r = np.asarray(d['R'])
a = np.asarray(d['A'])
L = np.asarray(d[['C', 'L1', 'L2']])
[17]:
def psi_bsiv_2sl2(theta):
# Parameters
beta = theta[0:2] # Parameters of interest
alpha = theta[2:2+3] # Nuisance parameters for E[Y | L]
gamma = theta[2+3:] # Nuisance parameters for E[A | L]
# Control arm model (E[Y | L])
ee_control = ee_regression(alpha, L, y, # Regression for Y given L
model='linear') # ... via linear model
ee_control = ee_control * (1 - r) # Restricting to control arm for fit
yhat = np.dot(L, alpha) # Predicted values for Y
# Received treatment model (E[A | L, R=1])
ee_receipt = ee_regression(gamma, L, a, # Regression for A given L
model='linear') # ... via linear model
ee_receipt = ee_receipt * r # Restricting to non-control arm for fit
ahat = np.dot(L, gamma) # Predicted values for A
# Controlled direct effect
X = np.asarray([np.ones(ahat.shape), ahat]).T # Stacking a new design matrix together
ee_cde = ee_regression(beta, X=X, y=y, # Regression for Y given A-hat
model='linear', # ... via linear model
offset=yhat) # ... with Y-hat offset term
ee_cde = ee_cde * r # Restricting to non-control arm for fit
# Returning stacked estimating functions
return np.vstack([ee_cde, ee_control, ee_receipt])
[20]:
estr = MEstimator(psi_bsiv_2sl2, init=[0., ]*8)
estr.estimate()
[21]:
estr.print_results(subset=[0, 1])
==============================================================
Estimation Method: M-estimator
--------------------------------------------------------------
No. Observations: 640 | No. Parameters: 8
Solving algorithm: lm | Max Iterations: 5000
Solving tolerance: 1e-09 | Allow P-Inverse: 1
Derivative Method: approx | Deriv Approx: 1e-09
Small N Correction: None | Distribution: Z-stat
==============================================================
Theta StdErr Z-score LCL UCL P-value S-value
--------------------------------------------------------------
-0.33 0.05 -6.67 -0.43 -0.23 0.00 35.17
0.18 0.09 2.01 0.00 0.36 0.04 4.49
==============================================================
Here, the first parameter corresponds to the effect of assignment to treatment versus control. So, these results indicate there is an effect of assignment and thus the exclusion restriction assumption is suspect (provided that the BSIV assumptions hold). The second parameter corresponds to the effect of receiving treatment among persons who were assigned to treatment. These results are similar to those reported in the paper. However, note that Richardson et al. use a g-estimator for estimation rather than two-stage least squares. These results do match those reported in Zivich (2025).
We can also examine the joint effect of assignment and treatment from this model by simply adding the two parameters together. While that can be done in the estimating functions themselves, it can also be computed external using the delta_method function
[27]:
def joint_effect(theta):
return [theta[0] + theta[1], ]
[29]:
var = delta_method(theta=estr.theta[:2], g=joint_effect, covariance=estr.variance[:2, :2])
[35]:
joint = joint_effect(theta=estr.theta)
ci = joint[0] - 1.96*np.sqrt(var), joint[0] + 1.96*np.sqrt(var)
print(np.round(joint, decimals=2))
print(np.round(ci, decimals=2))
[-0.15]
[-0.24 -0.06]
These are similar to those reported in Richardson et al. The advantage of the approach taken here (as described in Zivich 2025) is that we can avoid having to apply the bootstrap and rely on the empirical sandwich variance estimator instead.
References
Richardson DB, Dukes O, & Tchetgen Tchetgen EJ. (2023). Estimating the effect of a treatment when there is nonadherence in a trial. American Journal of Epidemiology, 192(10), 1772-1780.
Zivich PN. (2025). Re:”Estimating the effect of a treatment when there is nonadherence in a trial”. American Journal of Epidemiology, 194(2), 552-553.