Ross et al. (2024): Introduction to M-estimation

The following is a replication of the cases described in Ross et al. (2024). The original paper provides a tutorial on M-estimation and provides several introductory examples. Examples are provided in the context of regression, standardization, and measurement error. Here, we recreate the cases described in the paper. For finer details, see the original publication.

Setup

import numpy as np
import pandas as pd

import delicatessen
from delicatessen import MEstimator
from delicatessen.estimating_equations import ee_regression, ee_rogan_gladen
from delicatessen.utilities import inverse_logit

print("Versions")
print("NumPy:        ", np.__version__)
print("pandas:       ", pd.__version__)
print("Delicatessen: ", delicatessen.__version__)

Versions
NumPy:         1.25.2
pandas:        1.4.1
Delicatessen:  2.1

The following data is used for the first and second cases.

# From Table 1
d = pd.DataFrame()
d['X'] = [0, 0, 0, 0, 1, 1, 1, 1]            # X values
d['W'] = [0, 0, 1, 1, 0, 0, 1, 1]            # W values
d['Y'] = [0, 1, 0, 1, 0, 1, 0, 1]            # Y values
d['n'] = [496, 74, 113, 25, 85, 15, 15, 3]   # Counts
d['intercept'] = 1                           # Intercept term (always 1)

# Expanding rows by n
d = pd.DataFrame(np.repeat(d.values,         # Converting tabled data
                           d['n'], axis=0),  # ... by replicating counts
                 columns=d.columns)          # ... into rows for each X,W,Y
d = d[['intercept', 'X', 'W', 'Y']].copy()   # Dropping extra rows

Example 1: Logistic Regression

For the first example, we fit a logistic regression model for the variable \(Y\) given \(X,W\).

# Extracting arrays for easier coding later on
X = np.asarray(d[['intercept', 'X', 'W']])   # Design matrix for regression
y = np.asarray(d['Y'])                       # Outcome in regression

For estimation, we can use the built-in function for regression

def psi(theta):
    return ee_regression(theta=theta,
                         y=y, X=X,
                         model='logistic')

estr = MEstimator(psi, init=[0, 0, 0,])
estr.estimate()

# Formatting results into a nice table
result = pd.DataFrame()
result['Param'] = ['beta_0', 'beta_1', 'beta_2']
result['Coef'] = estr.theta
ci = estr.confidence_intervals()
result['LCL'] = ci[:, 0]
result['UCL'] = ci[:, 1]
result.round(2)

	Param	Coef	LCL	UCL
0	beta_0	-1.89	-2.13	-1.66
1	beta_1	0.12	-0.43	0.67
2	beta_2	0.36	-0.11	0.83

These results match those reported in the paper. They also match the results if you were to fit this using statsmodels GLM instead.

Example 2: Estimating the marginal risk difference

In this example, we build on the prior example. Using the logistic model, we can generate predictions and use those predictions to estimate the marginal risk difference by \(X\) via g-computation. Alternatively, we can use a logistic model to estimate propensity score and then estimate the marginal risk difference with inverse probability weighting.

First, let’s apply g-computation

# Copies of data with policies applied
d1 = d.copy()
d1['X'] = 1
X1 = np.asarray(d1[['intercept', 'X', 'W']])
d0 = d.copy()
d0['X'] = 0
X0 = np.asarray(d0[['intercept', 'X', 'W']])

def psi(theta):
    # Dividing parameters into corresponding parts and labels from slides
    beta = theta[0:3]              # Logistic model coefficients
    mu0, mu1 = theta[3], theta[4]  # Causal risks
    delta1 = theta[5]              # Causal contrast

    # Logistic regression model for outcome
    ee_logit = ee_regression(theta=beta,
                             y=y, X=X,
                             model='logistic')

    # Transforming logistic model coefficients into causal parameters
    y0_hat = inverse_logit(np.dot(X0, beta))  # Prediction under a=0
    y1_hat = inverse_logit(np.dot(X1, beta))  # Prediction under a=1

    # Estimating function for causal risk under a=1
    ee_r1 = y1_hat - mu1             # Simple mean

    # Estimating function for causal risk under a=0
    ee_r0 = y0_hat - mu0             # Simple mean

    # Estimating function for causal risk difference
    ee_rd = np.ones(d.shape[0])*((mu1 - mu0) - delta1)

    # Returning stacked estimating functions in order of parameters
    return np.vstack([ee_logit,   # EF of logistic model
                      ee_r0,      # EF of causal risk a=0
                      ee_r1,      # EF of causal risk a=1
                      ee_rd])     # EF of causal contrast

estr = MEstimator(psi, init=[0, 0, 0, 0.5, 0.5, 0])
estr.estimate()

# Formatting results into a nice table
result = pd.DataFrame()
result['Param'] = ['beta_0', 'beta_1', 'beta_2', 'mu_0', 'mu_1', 'delta']
result['Coef'] = estr.theta
ci = estr.confidence_intervals()
result['LCL'] = ci[:, 0]
result['UCL'] = ci[:, 1]
result.round(2)

	Param	Coef	LCL	UCL
0	beta_0	-1.89	-2.13	-1.66
1	beta_1	0.12	-0.43	0.67
2	beta_2	0.36	-0.11	0.83
3	mu_0	0.14	0.11	0.17
4	mu_1	0.15	0.09	0.22
5	delta	0.01	-0.06	0.09

These results match those reported in the publication (some differences in rounding by the chosen root-finding algorithm).

Now consider estimating the marginal risk difference using inverse probability weighting. The following code implements this estimator by-hand

a = np.asarray(d['X'])
W = np.asarray(d[['intercept', 'W']])

def psi(theta):
    # Dividing parameters into corresponding parts and labels from slides
    alpha = theta[0:2]              # Logistic model coefficients
    mu0, mu1 = theta[2], theta[3]   # Causal risks
    delta1 = theta[4]               # Causal contrast

    # Logistic regression model for propensity score
    ee_logit = ee_regression(theta=alpha,       # Regression model
                             y=a,               # ... for exposure
                             X=W,               # ... given confounders
                             model='logistic')  # ... logistic model

    # Transforming logistic model coefficients into causal parameters
    pscore = inverse_logit(np.dot(W, alpha))    # Propensity score
    wt = d['X']/pscore + (1-d['X'])/(1-pscore)  # Corresponding weights

    # Estimating function for causal risk under a=1
    ee_r1 = d['X']*d['Y']*wt - mu1         # Weighted conditional mean

    # Estimating function for causal risk under a=0
    ee_r0 = (1-d['X'])*d['Y']*wt - mu0     # Weighted conditional mean

    # Estimating function for causal risk difference
    ee_rd = np.ones(d.shape[0])*((mu1 - mu0) - delta1)

    # Returning stacked estimating functions in order of parameters
    return np.vstack([ee_logit,   # EF of logistic model
                      ee_r0,      # EF of causal risk a=0
                      ee_r1,      # EF of causal risk a=1
                      ee_rd])     # EF of causal contrast

estr = MEstimator(psi, init=[0, 0, 0.5, 0.5, 0])
estr.estimate()

# Formatting results into a nice table
result = pd.DataFrame()
result['Param'] = ['alpha_0', 'alpha_1', 'mu_0', 'mu_1', 'delta']
result['Coef'] = estr.theta
ci = estr.confidence_intervals()
result['LCL'] = ci[:, 0]
result['UCL'] = ci[:, 1]
result.round(2)

	Param	Coef	LCL	UCL
0	alpha_0	-1.74	-1.95	-1.53
1	alpha_1	-0.30	-0.83	0.24
2	mu_0	0.14	0.11	0.17
3	mu_1	0.15	0.09	0.22
4	delta	0.01	-0.06	0.08

Again, these results match those reported in the publication

Example 3: Outcome misclassification

For the final example

# Loading in data for the fusion example
d = pd.DataFrame()
d['R'] = [1, 1, 0, 0, 0, 0]           # R or population indicator
d['Y'] = [0, 0, 1, 1, 0, 0]           # True outcome
d['W'] = [1, 0, 1, 0, 1, 0]           # Measured outcome
d['n'] = [680, 270, 204, 38, 18, 71]  # Counts
d['intercept'] = 1                    # Intercept is always 1

# Expanding out data
d = pd.DataFrame(np.repeat(d.values, d['n'], axis=0),  # Expanding compact data frame
                 columns=d.columns)                    # ... keeping column names
d = d[['intercept', 'R', 'W', 'Y']].copy()             # Dropping the n column

# Converting to arrays to simplify process
r = np.asarray(d['R'])
w = np.asarray(d['W'])
y = np.asarray(d['Y'])

def psi(theta):
    return ee_rogan_gladen(theta=theta,    # Parameter of interest
                           y=d['Y'],       # ... correct measure
                           y_star=d['W'],  # ... mismeasure
                           r=d['R'])       # ... sample indicator

estr = MEstimator(psi, init=[0.75, 0.75, 0.75, 0.75])
estr.estimate()

# Formatting results into a nice table
result = pd.DataFrame()
result['Param'] = ['Corrected', 'Mismeasured', 'Sensitivity', 'Specificity']
result['Coef'] = estr.theta
ci = estr.confidence_intervals()
result['LCL'] = ci[:, 0]
result['UCL'] = ci[:, 1]
result.round(2)

	Param	Coef	LCL	UCL
0	Corrected	0.80	0.72	0.88
1	Mismeasured	0.72	0.69	0.74
2	Sensitivity	0.84	0.80	0.89
3	Specificity	0.80	0.71	0.88

These results match Ross et al.

Conclusions

Here, we replicated the introduction to M-estimation tutorial provided in Ross et al. The basic ideas illustrated in that paper and replicated here can be extended in numerous directions.

References

Ross RK, Zivich PN, Stringer JS, & Cole SR. (2024). M-estimation for common epidemiological measures: introduction and applied examples. International Journal of Epidemiology, 53(2), dyae030.