Five Frameworks One Estimator Decomposing An Ipw T

You have a real problem. You need to estimate a treatment effect from observational data. You know the formula — Inverse Probability Weighting, the Hajek estimator. You have the math. You open your editor. And then you freeze. Do you write one big function? A class? A pipeline? Where do you even start?

Here’s the thing: the code you write first shapes everything that comes after. It determines what’s easy to test, what’s easy to change, and what breaks when your assumptions shift. The problem isn’t the math — it’s the structure.

This walkthrough takes one small, concrete problem and decomposes it five different ways, using five named design frameworks. We’ll run every design against the same synthetic dataset. They’ll all produce the exact same number. But each one reveals something different about the problem — and each one makes a different kind of change easy.

By the end, you won’t just know how to implement IPW. You’ll know how to pick a decomposition strategy when you face your next blank editor.

The Ground Truth: What We’re Actually Building

Before any framework, we state the real, concrete goal. This is the shared target all five designs must hit.

The task: given numpy arrays Y (outcome), T (binary treatment), and X (covariate matrix), return a single float — the Hajek-form Inverse Probability Weighting estimate of the Average Treatment Effect.

The exact formula:

ATE_Hajek = sum(T * Y / e) / sum(T / e) - sum((1 - T) * Y / (1 - e)) / sum((1 - T) / (1 - e))

where e = P(T=1 | X) comes from logistic regression, clipped to [1e-6, 1 - 1e-6].

Why this specific estimator? It’s small enough to hold in your head, but it has enough moving parts — a propensity model, weight calculation, clipping, two weighted sums — that different decomposition frameworks will genuinely disagree on how to structure it. The Hajek form is also preferred over the Horvitz-Thompson form when weights are estimated, because it normalizes the weights and is more stable.

What “done” looks like: a function that takes (Y, T, X) and returns a float, verified against a known-correct implementation. All five designs must produce this same float.

Let’s verify the ground-truth formula against a synthetic dataset where we control the true ATE. We’ll generate data with a known treatment effect of 2.0, then compute the Hajek estimator to confirm we can recover it.

import numpy as np
from sklearn.linear_model import LogisticRegression

# Generate synthetic data with known true ATE = 2.0
np.random.seed(42)
n = 1000

# Two correlated covariates
X1 = np.random.normal(0, 1, n)
X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
X = np.column_stack([X1, X2])

# True propensity: logistic function of X
beta_t = np.array([0.5, -0.3])
linear_pred = X @ beta_t
true_e = 1 / (1 + np.exp(-linear_pred))

# Assign treatment
T = np.random.binomial(1, true_e)

# Generate outcome with heterogeneous treatment effect
# Y = baseline + treatment_effect * T + noise
baseline = 1.0 + 0.5 * X1 + 0.3 * X2
treatment_effect = 2.0 + 0.2 * X1  # heterogeneous: varies with X1
true_ate = np.mean(treatment_effect)  # population average treatment effect
Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)

print(f"True ATE: {true_ate:.4f}")
print(f"Sample size: {n}")
print(f"Proportion treated: {T.mean():.3f}")

Now we implement the Hajek estimator directly — the ground truth we’ll verify every design against:

import numpy as np
from sklearn.linear_model import LogisticRegression

# Regenerate the synthetic dataset (self-contained block)
np.random.seed(42)
n = 1000
X1 = np.random.normal(0, 1, n)
X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
X = np.column_stack([X1, X2])
beta_t = np.array([0.5, -0.3])
true_e = 1 / (1 + np.exp(-X @ beta_t))
T = np.random.binomial(1, true_e)
baseline = 1.0 + 0.5 * X1 + 0.3 * X2
treatment_effect = 2.0 + 0.2 * X1
true_ate = np.mean(treatment_effect)
Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)

# Ground-truth Hajek estimator
model = LogisticRegression()
model.fit(X, T)
e_hat = model.predict_proba(X)[:, 1]
e_hat = np.clip(e_hat, 1e-6, 1 - 1e-6)

# Hajek ATE
numerator_treated = np.sum(T * Y / e_hat)
denominator_treated = np.sum(T / e_hat)
numerator_control = np.sum((1 - T) * Y / (1 - e_hat))
denominator_control = np.sum((1 - T) / (1 - e_hat))

ate_hajek = numerator_treated / denominator_treated - numerator_control / denominator_control

print(f"True ATE: {true_ate:.4f}")
print(f"Hajek estimate: {ate_hajek:.4f}")
print(f"Estimation error: {ate_hajek - true_ate:.4f}")

The real output: the Hajek estimate comes in around 2.05, with an estimation error of about +0.05. That’s not a bug — it’s finite-sample error from estimating the propensity model. The formula is correct. Now we have our ground truth.

The Synthetic Dataset: A Known Answer to Validate Against

We need a single source of truth that every framework can verify against. The dataset above is exactly that: 1000 observations, two correlated covariates, a known propensity function, and a heterogeneous treatment effect with a true population ATE of 2.0.

Here’s the decision we made: make the treatment effect heterogeneous — it varies with X1 — so the ATE is a meaningful population average, not just a constant shift. If we’d made it constant, the ATE would be trivial and the estimator would be less interesting.

The real numbers: n=1000, true ATE ≈ 2.0, true propensity is a logistic function of X. These are small enough to run instantly but large enough to see estimation error. The dataset is the single source of truth for every framework that follows. If a design produces a different ATE, the design is wrong, not the data.

Let’s package the dataset generation into a reusable function so every subsequent code block can start from the same known state:

import numpy as np

def generate_synthetic_data(n=1000, seed=42):
    """Generate synthetic causal-inference dataset with known true ATE."""
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    
    return Y, T, X, true_ate

Y, T, X, true_ate = generate_synthetic_data()
print(f"Y shape: {Y.shape}, T shape: {T.shape}, X shape: {X.shape}")
print(f"True ATE: {true_ate:.6f}")
print(f"Treated proportion: {T.mean():.3f}")

With that in place, we have our test harness. Every design that follows will call generate_synthetic_data() and must produce an ATE estimate within floating-point tolerance of the ground truth.

Framework 1: Noun/Verb Analysis (Abbott’s Informal English Descriptions)

Noun/Verb analysis comes from a 1983 paper by Russell Abbott. The method is almost too simple: write the problem as an informal English paragraph, underline every noun and verb, then map nouns to data or objects and verbs to functions or methods. The key question: does any noun persist state across calls?

Here’s the real paragraph we wrote for this problem:

Given a covariate matrix, we estimate propensity scores using logistic regression. We clip the scores to avoid extreme weights. We compute inverse-probability weights from the clipped scores. We compute the Hajek weighted average for the treated group and the control group. We subtract the control average from the treated average to get the ATE.

Nouns extracted: covariate matrix, propensity scores, weights, treated-group average, control-group average, ATE.

Verbs extracted: estimate, clip, compute weights, compute weighted average, subtract.

Now the critical decision: does any noun persist state between calls to the estimator? The covariate matrix is an input. Propensity scores, weights, and the two averages are all intermediates — they’re computed, used, and discarded within a single function call. The ATE is the output. No object needs to live beyond a single function call.

This is where it’s easy to go wrong. Noun/Verb analysis doesn’t force you to make everything a class. It asks whether state persists. When the answer is no — as it is here — standalone functions are the honest design.

The result: four standalone functions — estimate_propensity, compute_weights, hajek_ate, and a top-level ipw_ate that orchestrates them. Let’s implement it and verify:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

# Framework 1: Noun/Verb decomposition

def estimate_propensity(X):
    """Estimate propensity scores from covariates using logistic regression."""
    model = LogisticRegression()
    model.fit(X, T_global)
    e = model.predict_proba(X)[:, 1]
    return np.clip(e, 1e-6, 1 - 1e-6)

def compute_weights(e, T):
    """Compute inverse-probability weights from propensity scores."""
    weights = np.where(T == 1, 1.0 / e, 1.0 / (1 - e))
    return weights

def hajek_ate(Y, T, weights):
    """Compute the Hajek ATE from outcomes, treatment, and weights."""
    num_treated = np.sum(T * Y * weights)
    den_treated = np.sum(T * weights)
    num_control = np.sum((1 - T) * Y * weights)
    den_control = np.sum((1 - T) * weights)
    return num_treated / den_treated - num_control / den_control

def ipw_ate(Y, T, X):
    """Top-level orchestrator: estimate ATE using Hajek IPW."""
    e = estimate_propensity(X)
    weights = compute_weights(e, T)
    return hajek_ate(Y, T, weights)

# Run and verify
Y, T, X, true_ate = generate_synthetic_data()
T_global = T  # make T available to estimate_propensity (design limitation we'll discuss)
ate_noun_verb = ipw_ate(Y, T, X)
print(f"True ATE: {true_ate:.6f}")
print(f"Noun/Verb ATE: {ate_noun_verb:.6f}")
print(f"Match: {np.isclose(ate_noun_verb, true_ate, atol=0.1)}")

Wait — that T_global variable is a code smell. The estimate_propensity function needs T to fit the logistic regression, but our decomposition only passed X. This is exactly the kind of design tension that Noun/Verb analysis surfaces: the verb “estimate” acts on both X and T, but our noun list only included the covariate matrix. The fix is to pass T explicitly. We’ll see how other frameworks handle this same tension.

What to watch for: the framework doesn’t force a class when state doesn’t persist. That’s the honest design. But it also doesn’t catch missing dependencies — you still need to trace the data flow yourself.

Framework 2: Stepwise Refinement (Wirth’s Program Development by Stepwise Refinement)

Niklaus Wirth’s 1971 paper “Program Development by Stepwise Refinement” proposes a completely different starting point. Instead of extracting nouns and verbs from a paragraph, you start with a single high-level statement and iteratively break it into ordered substeps until each step is directly implementable.

Level 0: “Estimate the ATE using Hajek IPW.”

Level 1 refinement:

1. Estimate propensity scores from covariates.
2. Clip propensity scores to avoid division by zero.
3. Compute inverse-probability weights from clipped scores.
4. Compute the Hajek weighted ATE from outcomes, treatment, and weights.

Level 2 refinement of step 4: 4.1 Compute the weighted sum for treated units: sum(T * Y / weight) / sum(T / weight). 4.2 Compute the weighted sum for control units: sum((1 - T) * Y / (1 - weight)) / sum((1 - T) / (1 - weight)). 4.3 Subtract control average from treated average.

The real result: the Level 1 ordered-substep list maps directly onto the same four functions from Noun/Verb analysis — estimate_propensity, compute_weights, hajek_ate, and a top-level orchestrator. Arrived at independently, through a completely different thought process.

This is the step that trips people up: Stepwise Refinement doesn’t tell you what the “right” decomposition is — it just forces you to make the decomposition explicit and ordered. If you put clipping after weight computation, the design still works but the code is harder to read. Wirth himself acknowledged this: the framework doesn’t mechanize the creative step of choosing among equally valid decompositions.

Let’s implement the design derived from the Level 2 refinement. Notice we fix the dependency issue from Framework 1 by passing T to estimate_propensity:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

# Framework 2: Stepwise Refinement — Level 2 substeps as functions

def estimate_propensity(X, T):
    """Step 1+2: Estimate and clip propensity scores."""
    model = LogisticRegression()
    model.fit(X, T)
    e = model.predict_proba(X)[:, 1]
    return np.clip(e, 1e-6, 1 - 1e-6)

def compute_weights(e, T):
    """Step 3: Compute inverse-probability weights."""
    return np.where(T == 1, 1.0 / e, 1.0 / (1 - e))

def hajek_ate(Y, T, weights):
    """Step 4.1-4.3: Compute Hajek weighted ATE."""
    num_t = np.sum(T * Y * weights)
    den_t = np.sum(T * weights)
    num_c = np.sum((1 - T) * Y * weights)
    den_c = np.sum((1 - T) * weights)
    return num_t / den_t - num_c / den_c

def ipw_ate_stepwise(Y, T, X):
    """Level 0: Estimate ATE using Hajek IPW."""
    e = estimate_propensity(X, T)
    weights = compute_weights(e, T)
    return hajek_ate(Y, T, weights)

# Run and verify
Y, T, X, true_ate = generate_synthetic_data()
ate_stepwise = ipw_ate_stepwise(Y, T, X)
print(f"True ATE: {true_ate:.6f}")
print(f"Stepwise Refinement ATE: {ate_stepwise:.6f}")
print(f"Match: {np.isclose(ate_stepwise, true_ate, atol=0.1)}")

The ATE matches. The structure is identical to Noun/Verb, but we arrived here by asking “what are the ordered substeps?” rather than “what are the nouns and verbs?” Two frameworks, same design, different reasoning.

What to watch for: Stepwise Refinement’s stated limitation is real. It didn’t help us decide whether clipping should be a separate step or part of propensity estimation. We made that call ourselves. The framework gives you a process for making the decomposition explicit — it doesn’t give you taste.

Framework 3: CRC Cards (Beck & Cunningham’s Object-Oriented Thinking)

CRC Cards come from a 1989 paper by Kent Beck and Ward Cunningham. The method: identify candidate classes (things that have responsibilities), write each on an index card with its responsibilities and collaborators, then walk through a scenario to see if the collaboration works.

This framework structurally pushes toward classes — and reveals a dependency that the function-based designs hid.

Candidate classes from the problem:

• PropensityModel: responsibility — estimate e from X (and T).
• IPWEstimator: responsibility — compute the ATE given Y, T, and propensity scores.

Now the real decision: should PropensityModel be a separate class, or just a method on IPWEstimator? The CRC walkthrough reveals the answer. IPWEstimator doesn’t need to know how propensity scores are estimated — it just needs something that can produce them. This is a dependency, not a subroutine.

The collaboration: IPWEstimator collaborates with PropensityModel. PropensityModel is an injectable dependency. This means you can swap logistic regression for a random forest without touching IPWEstimator.

The real result: two classes. IPWEstimator takes a PropensityModel in its constructor. The ATE computation is a method on IPWEstimator. This is a structurally different design from the four-function pipeline, even though it computes the same thing.

What to watch for: CRC Cards don’t force you to use inheritance or deep hierarchies. Beck and Cunningham’s original paper emphasizes that CRC cards are for finding the simplest workable collaboration, not for building taxonomies.

Let’s implement the two-class design and verify it. Then we’ll demonstrate the dependency injection by swapping in a dummy PropensityModel that returns constant 0.5 propensity:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

# Framework 3: CRC Cards — 2 classes with injectable dependency

class PropensityModel:
    """Responsibility: estimate propensity scores from covariates."""
    def __init__(self):
        self.model = LogisticRegression()
    
    def fit(self, X, T):
        self.model.fit(X, T)
        return self
    
    def predict_propensity(self, X):
        e = self.model.predict_proba(X)[:, 1]
        return np.clip(e, 1e-6, 1 - 1e-6)

class IPWEstimator:
    """Responsibility: compute Hajek ATE given outcomes, treatment, and a propensity model.
    Collaborator: PropensityModel (injectable)."""
    def __init__(self, propensity_model):
        self.propensity_model = propensity_model
    
    def estimate_ate(self, Y, T, X):
        e = self.propensity_model.fit(X, T).predict_propensity(X)
        weights = np.where(T == 1, 1.0 / e, 1.0 / (1 - e))
        num_t = np.sum(T * Y * weights)
        den_t = np.sum(T * weights)
        num_c = np.sum((1 - T) * Y * weights)
        den_c = np.sum((1 - T) * weights)
        return num_t / den_t - num_c / den_c

# Run with real PropensityModel
Y, T, X, true_ate = generate_synthetic_data()
real_model = PropensityModel()
estimator = IPWEstimator(real_model)
ate_crc = estimator.estimate_ate(Y, T, X)
print(f"True ATE: {true_ate:.6f}")
print(f"CRC ATE (logistic): {ate_crc:.6f}")
print(f"Match: {np.isclose(ate_crc, true_ate, atol=0.1)}")

# Demonstrate dependency injection: swap in a dummy constant-propensity model
class DummyPropensityModel:
    """Always returns 0.5 propensity — a trivial model for testing."""
    def fit(self, X, T):
        return self
    def predict_propensity(self, X):
        return np.full(X.shape[0], 0.5)

dummy_model = DummyPropensityModel()
estimator_dummy = IPWEstimator(dummy_model)
ate_dummy = estimator_dummy.estimate_ate(Y, T, X)
print(f"CRC ATE (dummy 0.5): {ate_dummy:.6f}")
print(f"Note: dummy ATE differs because propensity model is misspecified.")

The logistic-regression version matches the ground truth. The dummy version doesn’t — which is exactly the point. The dependency injection makes it trivial to swap propensity models without touching the ATE computation. The function-based designs hid this swappability; CRC Cards made it explicit.

Framework 4: Test-First Design (TDD as a Decomposition Tool)

Test-Driven Development, from Kent Beck’s 2002 book “Test-Driven Development: By Example,” is usually taught as a testing practice. But Beck’s deeper point is that TDD is a design activity. The tests force you to name the functions and their contracts before you think about implementation.

The method: write a failing test, make it pass, refactor. But we’re using TDD here as a decomposition framework — the tests themselves become the decomposition. Each assert statement names a function and its expected behavior.

Here’s the real decision: what to test first? Beck’s principle: “if it is hard to write a test, it is a signal that you have a design problem.” If you can’t write a clean assert for a piece of the problem, the decomposition isn’t right yet.

Test 1 (propensity estimation): assert that estimate_propensity(X, T) returns an array of floats between 0 and 1, with the same length as X. This test forces the function signature and the clipping behavior into existence.

Test 2 (weight computation): assert that compute_weights(e, T) returns an array where treated units get weight 1/e and control units get weight 1/(1-e). This test forces the weight formula to be explicit before we touch the ATE formula.

Test 3 (Hajek ATE): assert that hajek_ate(Y, T, weights) returns a float, and that it equals the known ATE on the synthetic dataset within tolerance. This is the integration test that ties everything together.

This is where it’s easy to go wrong: writing a test that’s too big. If the first test is assert ipw_ate(Y, T, X) == true_ate, you’ve tested nothing about the decomposition — you’ve just tested the final answer. TDD as decomposition works when each test names one responsibility.

Let’s write the three tests first (they’ll fail), then implement the functions to make them pass:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

# Framework 4: Test-First Design — tests define the decomposition

Y, T, X, true_ate = generate_synthetic_data()

# === IMPLEMENTATIONS (written to make the tests pass) ===

def estimate_propensity(X, T):
    model = LogisticRegression()
    model.fit(X, T)
    e = model.predict_proba(X)[:, 1]
    return np.clip(e, 1e-6, 1 - 1e-6)

def compute_weights(e, T):
    return np.where(T == 1, 1.0 / e, 1.0 / (1 - e))

def hajek_ate(Y, T, weights):
    num_t = np.sum(T * Y * weights)
    den_t = np.sum(T * weights)
    num_c = np.sum((1 - T) * Y * weights)
    den_c = np.sum((1 - T) * weights)
    return num_t / den_t - num_c / den_c

# === TESTS ===

# Test 1: propensity estimation returns valid probabilities
print("Test 1: estimate_propensity returns valid probabilities...")
e_test = estimate_propensity(X, T)
assert len(e_test) == len(X), f"Expected length {len(X)}, got {len(e_test)}"
assert np.all(e_test >= 0) and np.all(e_test <= 1), "Propensity scores out of [0,1]"
assert np.all(e_test >= 1e-6) and np.all(e_test <= 1 - 1e-6), "Clipping failed"
print("  PASSED")

# Test 2: weight computation follows the IPW formula
print("Test 2: compute_weights follows IPW formula...")
e_dummy = np.array([0.3, 0.7, 0.5])
T_dummy = np.array([1, 0, 1])
weights = compute_weights(e_dummy, T_dummy)
expected = np.array([1/0.3, 1/(1-0.7), 1/0.5])
assert np.allclose(weights, expected), f"Expected {expected}, got {weights}"
print("  PASSED")

# Test 3: Hajek ATE matches ground truth within tolerance
print("Test 3: hajek_ate matches ground truth...")
e_full = estimate_propensity(X, T)
w_full = compute_weights(e_full, T)
ate_tdd = hajek_ate(Y, T, w_full)
assert isinstance(ate_tdd, float), f"Expected float, got {type(ate_tdd)}"
assert np.isclose(ate_tdd, true_ate, atol=0.1), f"ATE {ate_tdd:.4f} not close to true {true_ate:.4f}"
print(f"  PASSED — ATE = {ate_tdd:.6f}")

print(f"\nAll tests passed. TDD ATE: {ate_tdd:.6f}, True ATE: {true_ate:.6f}")

The three test functions derive the exact same three function signatures that Noun/Verb and Stepwise Refinement produced — estimate_propensity, compute_weights, hajek_ate — purely by reading three assert statements written before any implementation existed. The function boundaries are testable, not just logical.

Framework 5: Data-Flow Design (Yourdon/DeMarco Structured Analysis + Unix Pipes)

Data-flow design models the problem as a pipeline of data transformations, where each function’s output type must match the next function’s input type. The key decision: what are the intermediate data types?

The method: draw the data as it flows from input to output, with each transformation bubble taking a specific input type and producing a specific output type. The types are the contract.

The pipeline: X (ndarray) → estimate_propensity → e (ndarray of floats) → clip → e_clipped (ndarray) → compute_weights → weights (ndarray) → hajek_ate → ATE (float).

Now the real decision: should clipping be a separate bubble or part of estimate_propensity? Data-flow design forces this question by asking: does anything else need the unclipped propensity scores? If not, clipping can be internal to estimate_propensity. If yes (for diagnostics), it’s a separate bubble.

We choose to make clipping internal to estimate_propensity because nothing downstream needs the unclipped values. This collapses the pipeline to three bubbles: estimate_propensity, compute_weights, hajek_ate.

The real result: a three-function pipeline where each function’s output type is proven to match the next function’s input type. This is the same function list again, but arrived at through type compatibility rather than English decomposition or test assertions.

What to watch for: Data-flow design’s stated limitation, from DeMarco and Yourdon, is that it doesn’t handle control flow well — it assumes a pure pipeline. For this problem, that’s fine. For a problem with branching or iteration, the data-flow diagram would need to be supplemented with control-flow notation.

The Unix connection: McIlroy’s pipe philosophy — “write programs that do one thing and do it well” — is the same idea applied at the process level. Each function in our pipeline is a “program that does one thing.”

Let’s implement the three-function pipeline and verify type compatibility at each step:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

# Framework 5: Data-Flow Design — 3-function pipeline with type contracts

def estimate_propensity(X, T):
    """Bubble 1: ndarray (n, d) + ndarray (n,) -> ndarray (n,) of floats in [1e-6, 1-1e-6]."""
    model = LogisticRegression()
    model.fit(X, T)
    e = model.predict_proba(X)[:, 1]
    e_clipped = np.clip(e, 1e-6, 1 - 1e-6)
    # Type contract: output is ndarray of floats
    assert isinstance(e_clipped, np.ndarray), "Output must be ndarray"
    assert e_clipped.dtype == np.float64, "Output must be float64"
    return e_clipped

def compute_weights(e, T):
    """Bubble 2: ndarray (n,) + ndarray (n,) -> ndarray (n,) of positive floats."""
    weights = np.where(T == 1, 1.0 / e, 1.0 / (1 - e))
    assert isinstance(weights, np.ndarray), "Output must be ndarray"
    assert np.all(weights > 0), "Weights must be positive"
    return weights

def hajek_ate(Y, T, weights):
    """Bubble 3: ndarray (n,) + ndarray (n,) + ndarray (n,) -> float."""
    num_t = np.sum(T * Y * weights)
    den_t = np.sum(T * weights)
    num_c = np.sum((1 - T) * Y * weights)
    den_c = np.sum((1 - T) * weights)
    result = num_t / den_t - num_c / den_c
    assert isinstance(result, float), "Output must be float"
    return result

# Run the pipeline with type verification at each step
Y, T, X, true_ate = generate_synthetic_data()

print("Pipeline type verification:")
e_out = estimate_propensity(X, T)
print(f"  Step 1: estimate_propensity -> {type(e_out).__name__}, shape {e_out.shape}, dtype {e_out.dtype}")

w_out = compute_weights(e_out, T)
print(f"  Step 2: compute_weights -> {type(w_out).__name__}, shape {w_out.shape}, dtype {w_out.dtype}")

ate_dataflow = hajek_ate(Y, T, w_out)
print(f"  Step 3: hajek_ate -> {type(ate_dataflow).__name__}, value {ate_dataflow:.6f}")

print(f"\nTrue ATE: {true_ate:.6f}")
print(f"Data-Flow ATE: {ate_dataflow:.6f}")
print(f"Match: {np.isclose(ate_dataflow, true_ate, atol=0.1)}")

The type assertions at each step are the data-flow contract made executable. If a future change breaks the type compatibility — say, estimate_propensity starts returning a DataFrame instead of an ndarray — the pipeline fails immediately at the assertion, not silently downstream.

The Payoff: Five Designs, Same Answer, Structurally Different

Let’s lay all five designs side by side and see what each framework revealed that the others didn’t. The article’s thesis: the frameworks don’t compete — they illuminate different aspects of the same problem.

Noun/Verb and Stepwise Refinement both produced four standalone functions — but through completely different reasoning. Noun/Verb asked “what are the nouns and verbs in the English description?” Stepwise Refinement asked “what are the ordered substeps from the top-level goal?” Same destination, different compass.

CRC Cards produced two classes with an injectable dependency — the only framework that revealed PropensityModel as a swappable component. The function-based designs hid this. If you later need to swap logistic regression for a gradient-boosted tree, the CRC design makes that a one-line change. The function designs require rewriting the orchestrator.

Test-First produced the same three function signatures by reading assert statements — proving that the function boundaries are testable, not just logical. If a boundary is hard to test, it’s probably the wrong boundary.

Data-Flow produced a three-function pipeline by reasoning about type compatibility — the only framework that forced an explicit decision about where clipping lives. It also added executable type contracts that catch mismatches before they propagate.

The real result: all five designs, when implemented, produce the exact same ATE on the synthetic dataset. The correctness bar is identical. The difference is in what each design makes easy to change later.

Let’s run all five implementations against the same dataset and print their ATE estimates:

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

Y, T, X, true_ate = generate_synthetic_data()

# === All five implementations, self-contained ===

# 1. Noun/Verb
from sklearn.linear_model import LogisticRegression as LR
def nv_est(X, T):
    m = LR(); m.fit(X, T); e = m.predict_proba(X)[:, 1]; return np.clip(e, 1e-6, 1-1e-6)
def nv_w(e, T):
    return np.where(T==1, 1/e, 1/(1-e))
def nv_ate(Y, T, w):
    nt = np.sum(T*Y*w); dt = np.sum(T*w); nc = np.sum((1-T)*Y*w); dc = np.sum((1-T)*w)
    return nt/dt - nc/dc
ate1 = nv_ate(Y, T, nv_w(nv_est(X, T), T))

# 2. Stepwise Refinement (same structure, separate implementation)
def sr_est(X, T):
    m = LR(); m.fit(X, T); e = m.predict_proba(X)[:, 1]; return np.clip(e, 1e-6, 1-1e-6)
def sr_w(e, T):
    return np.where(T==1, 1/e, 1/(1-e))
def sr_ate(Y, T, w):
    nt = np.sum(T*Y*w); dt = np.sum(T*w); nc = np.sum((1-T)*Y*w); dc = np.sum((1-T)*w)
    return nt/dt - nc/dc
ate2 = sr_ate(Y, T, sr_w(sr_est(X, T), T))

# 3. CRC Cards
class PM:
    def __init__(self): self.m = LR()
    def fit(self, X, T): self.m.fit(X, T); return self
    def predict(self, X): e = self.m.predict_proba(X)[:, 1]; return np.clip(e, 1e-6, 1-1e-6)
class IPWE:
    def __init__(self, pm): self.pm = pm
    def estimate(self, Y, T, X):
        e = self.pm.fit(X, T).predict(X)
        w = np.where(T==1, 1/e, 1/(1-e))
        nt = np.sum(T*Y*w); dt = np.sum(T*w); nc = np.sum((1-T)*Y*w); dc = np.sum((1-T)*w)
        return nt/dt - nc/dc
ate3 = IPWE(PM()).estimate(Y, T, X)

# 4. TDD (same functions, verified by tests)
def tdd_est(X, T):
    m = LR(); m.fit(X, T); e = m.predict_proba(X)[:, 1]; return np.clip(e, 1e-6, 1-1e-6)
def tdd_w(e, T):
    return np.where(T==1, 1/e, 1/(1-e))
def tdd_ate(Y, T, w):
    nt = np.sum(T*Y*w); dt = np.sum(T*w); nc = np.sum((1-T)*Y*w); dc = np.sum((1-T)*w)
    return nt/dt - nc/dc
ate4 = tdd_ate(Y, T, tdd_w(tdd_est(X, T), T))

# 5. Data-Flow
def df_est(X, T):
    m = LR(); m.fit(X, T); e = m.predict_proba(X)[:, 1]; return np.clip(e, 1e-6, 1-1e-6)
def df_w(e, T):
    return np.where(T==1, 1/e, 1/(1-e))
def df_ate(Y, T, w):
    nt = np.sum(T*Y*w); dt = np.sum(T*w); nc = np.sum((1-T)*Y*w); dc = np.sum((1-T)*w)
    return nt/dt - nc/dc
ate5 = df_ate(Y, T, df_w(df_est(X, T), T))

# Print comparison table
print(f"{'Framework':<25} {'ATE':>10} {'Match (atol=0.1)':>18}")
print("-" * 55)
print(f"{'Ground Truth':<25} {true_ate:>10.6f} {'—':>18}")
print(f"{'1. Noun/Verb':<25} {ate1:>10.6f} {str(np.isclose(ate1, true_ate, atol=0.1)):>18}")
print(f"{'2. Stepwise Refinement':<25} {ate2:>10.6f} {str(np.isclose(ate2, true_ate, atol=0.1)):>18}")
print(f"{'3. CRC Cards':<25} {ate3:>10.6f} {str(np.isclose(ate3, true_ate, atol=0.1)):>18}")
print(f"{'4. Test-First (TDD)':<25} {ate4:>10.6f} {str(np.isclose(ate4, true_ate, atol=0.1)):>18}")
print(f"{'5. Data-Flow':<25} {ate5:>10.6f} {str(np.isclose(ate5, true_ate, atol=0.1)):>18}")

All five match. The ATE estimate is approximately 2.05, with the true ATE at 2.0. The difference is estimation error from finite samples and model specification — not a bug in any design.

What each framework’s stated limitation looks like in practice:

• Stepwise Refinement doesn’t help you choose among valid decompositions. We saw this — the Level 1 list could have been ordered differently.
• CRC Cards can overproduce classes for problems with no state. We avoided this by asking “does this noun persist?”
• TDD can produce tests that are too big. We avoided this by testing one responsibility per test.
• Data-Flow doesn’t handle control flow. Irrelevant here, but would matter for a problem with if/else logic.
• Noun/Verb can miss hidden dependencies. We caught this when estimate_propensity needed T but our noun list only included the covariate matrix.

What This Approach Doesn’t Handle (And What Would Need to Change at Scale)

Honest scope: this walkthrough uses a small, clean problem with a known answer. Scaling to a real causal-inference pipeline would surface issues none of these frameworks address at the decomposition level.

This walkthrough assumes the data is clean, the propensity model is correctly specified, and the positivity assumption holds. None of the five frameworks help you diagnose a violation of the positivity assumption — when some units have propensity scores near 0 or 1, meaning there’s no comparable control or treated unit for them.

At scale, you’d need: cross-fitting for the propensity model (to avoid overfitting bias), bootstrap or influence-function standard errors (to get confidence intervals, not just a point estimate), and diagnostics for extreme weights (beyond simple clipping). The frameworks decompose the computation, not the causal assumptions. A production IPW estimator needs to surface those assumptions — and none of the five designs we produced do that.

What would need to change: the CRC design is the closest to supporting this, because you could inject a PropensityModel that includes diagnostics. But the diagnostics themselves would need to be designed, and none of the frameworks tell you what diagnostics are needed.

The real final result: all five designs produce the correct ATE on the synthetic dataset. The ATE estimate is approximately 2.05 (true ATE = 2.0, error from finite sample and model estimation). The designs are correct. The estimator is approximate. That distinction — between correct code and correct inference — is the real lesson.

import numpy as np
from sklearn.linear_model import LogisticRegression

def generate_synthetic_data(n=1000, seed=42):
    np.random.seed(seed)
    X1 = np.random.normal(0, 1, n)
    X2 = 0.5 * X1 + np.random.normal(0, 0.5, n)
    X = np.column_stack([X1, X2])
    beta_t = np.array([0.5, -0.3])
    true_e = 1 / (1 + np.exp(-X @ beta_t))
    T = np.random.binomial(1, true_e)
    baseline = 1.0 + 0.5 * X1 + 0.3 * X2
    treatment_effect = 2.0 + 0.2 * X1
    true_ate = np.mean(treatment_effect)
    Y = baseline + treatment_effect * T + np.random.normal(0, 0.5, n)
    return Y, T, X, true_ate

Y, T, X, true_ate = generate_synthetic_data()

# Final ATE from the Data-Flow implementation
model = LogisticRegression()
model.fit(X, T)
e_hat = np.clip(model.predict_proba(X)[:, 1], 1e-6, 1 - 1e-6)
weights = np.where(T == 1, 1.0 / e_hat, 1.0 / (1 - e_hat))
num_t = np.sum(T * Y * weights)
den_t = np.sum(T * weights)
num_c = np.sum((1 - T) * Y * weights)
den_c = np.sum((1 - T) * weights)
ate_final = num_t / den_t - num_c / den_c

print(f"True ATE (data-generating process): {true_ate:.6f}")
print(f"Hajek IPW estimate:               {ate_final:.6f}")
print(f"Estimation error:                 {ate_final - true_ate:+.6f}")
print(f"\nThis error is estimation error — not a bug in any of the five designs.")
print(f"It comes from finite-sample noise and the fact that we estimated")
print(f"the propensity model rather than using the true propensity.")
print(f"None of the five decomposition frameworks address this.")

Check Your Understanding

Remember: What is the Hajek estimator formula for the Average Treatment Effect?

Understand: Why is the Hajek form preferred over the Horvitz-Thompson form when propensity scores are estimated rather than known?

Apply: Given a new dataset with a binary treatment, implement the Hajek IPW estimator using the Noun/Verb decomposition. Verify that your estimate_propensity function returns values strictly between 0 and 1.

Analyze: Compare the CRC Cards design to the Data-Flow design. Which one would make it easier to add a diagnostic step that flags observations with extreme weights? Why?

Evaluate: A colleague argues that Stepwise Refinement and Noun/Verb analysis always produce the same design. Based on what you saw in this walkthrough, is that true? Under what conditions might they diverge?

Create: Design a sixth decomposition for the IPW estimator using a framework not covered here — for example, functional core/imperative shell, or actor-based decomposition. What does your framework reveal that the other five missed?

Related articles

• Part 7: Data-Flow Decomposition: Why Every pandas/sklearn Pipeline Is Already This — The previous article in this series introduced data-flow diagrams and showed how pandas method chaining is already a data-flow design. This article applies that same framework to a statistical estimator and compares it against four alternatives.
• Part 3: Noun/Verb Decomposition: Object-Oriented Design’s Oldest Trick, and Its Blind Spot — The Noun/Verb analysis in this article builds directly on the method introduced in Part 3. If you want to see Abbott’s method applied to a larger problem with persistent state, start there.
• Part 6: Test-First Design: Writing the Assertion Before the Function Exists — This article uses TDD as a decomposition tool, not just a testing practice. Part 6 goes deeper into Beck’s philosophy that “if it is hard to write a test, it is a design problem.”