Hazards on the Causal Path

Three Ways to Say ‘Because’ in Survival Analysis

Nathaniel Forde

Data Science @ Personio

and Open Source Contributor @ PyMC

2026-06-14

Preliminaries

Who am I?

• I’m a data scientist at Personio
  • Bayesian statistician,
  • Reformed philosopher and logician.
• Website: https://nathanielf.github.io/

Link to my website

Code or it didn’t Happen

The full worked example can be found on my blog

“Causation is not a relation between events frozen in amber.

It is the process by which one thing makes another happen.”

The Pitch

Data scientists use causal inference every day.

But “causal inference” is not one thing. There are at least three schools, and they disagree about mediation.

This talk: build a model that works for all of them. Along the way, dissolve a philosophical puzzle.

Agenda

• Three ways to say “because” — and where they disagree

• Building a generative process model in PyMC

• Dissolving the cross-world counterfactual objection

Act I: Three Ways to Say “Because”

The Pre-emption Puzzle

Two assassins are hired. Assassin A poisons the coffee. Assassin B waits with a sniper rifle as backup.

A’s poison works. B never fires. The target dies.

Counterfactual test: “Would the target have survived if A hadn’t acted?”

No — B would have fired. So A fails the counterfactual test.

But A clearly caused the death.

This is the classic pre-emption case from the philosophy of causation (Lewis 1973). The audience doesn’t need to memorise the details — just the punchline: asking “would the outcome have been different?” is not always enough. Sometimes you need to know which process actually did the work. This isn’t an academic curiosity. In mediation analysis, the same issue arises: the standard counterfactual framework denies that certain mediation questions are even well-defined.

School 1: Potential Outcomes

The Framework

• Rubin, Holland (1986): “No causation without manipulation”
• Cause = difference between potential outcomes: \(Y(1) - Y(0)\)
• Clean, experiment-centric
• The gold standard for average treatment effects

The Limitation

• Mediation requires \(Y(1, M(0))\)
• “What would your outcome be under treatment, if your mediator took the value it would have had under control?”
• This nests two worlds
• Many PO theorists reject this as ill-defined

The potential outcomes framework is enormously successful for total effects. But it has a principled objection to mediation: the Natural Direct Effect requires imagining a person who is simultaneously treated (X=1) but whose mediator behaves as if untreated (M(0)). Robins, Richardson and others have argued this is incoherent because no single experiment could produce this state. If your framework says mediation is ill-defined, you cannot decompose effects into direct and indirect components.

School 2: Structural Causal Models

The Framework

• Pearl (2001): DAGs + structural equations + do-calculus
• Allows cross-world counterfactuals via structural equations
• \(\text{NDE} = E[Y(1, M(0))] - E[Y(0, M(0))]\) is well-defined
• Mediation is a first-class citizen

The Cost

• Requires commitment to structural equations as real
• The cross-world counterfactual remains philosophically contentious
• “Your structural equations had better be right”

Pearl’s framework solves the mediation problem by embracing structural equations. If you believe the equations describe the actual data-generating process, then cross-world counterfactuals are just computations over the model. The NDE formula falls out naturally. But critics (Dawid, Robins in some moods) object that this puts too much weight on untestable structural assumptions. We have a stalemate between two interventionist frameworks.

School 3: Process / Aalen

Aalen (1987, 1989): Causation as unfolding temporal process.

Effects are increments to risk over time, not static comparisons between frozen worlds.

The causal story is told by the trajectory, not the endpoint.

Time is the medium, not just a container.

This is where survival analysis lives naturally. And it sidesteps the cross-world problem entirely.

Aalen’s view is less well-known in the ML/data science community but standard in biostatistics and survival analysis. The key insight: causation is not a comparison between two frozen worlds — it is the accumulation of influence over time. A treatment does not “cause” the outcome at a single moment; it continuously shapes the hazard. This process-oriented view handles pre-emption naturally: the cause is the process that actually transmitted the influence. And it handles mediation naturally: the indirect effect is the product of two time-varying path coefficients, observable as the process unfolds.

The Stalemate

	Potential Outcomes	Pearl SCM	Process (Aalen)
Cause is…	Difference between worlds	Structural mechanism	Unfolding process
Mediation?	Ill-defined (cross-world)	Well-defined (structural eqns)	Naturally decomposable
Strength	Clean identification	Full causal calculus	Temporal dynamics
Weakness	No mediation	Structural faith	No marginal estimand?

What if one model could give you all three rows of the last column?

The table oversimplifies, but the pattern is real. Each school has a genuine strength and a genuine limitation. The punchline of this talk: dynamic path analysis in PyMC gives you the process decomposition (Aalen), the structural flexibility to compute NDE/NIE (Pearl), and a marginal survival contrast that even a potential outcomes purist can accept (Rubin). The trick is that the generative model does the heavy lifting.

Back to Earth: A Concrete Problem

A patient starts an exercise programme.

It reduces cardiovascular risk directly. Good.

But it also causes joint inflammation. The inflammation increases risk over time.

At the six-month endpoint, the treatment looks modestly protective.

You have no idea it’s being undermined.

Three schools, three different things they can tell you about this situation.

Same motivating example, but now the audience has a reason to care. The potential outcomes framework can tell you the total effect but refuses to decompose it. Pearl’s SCM can decompose it but requires commitment to structural equations at a single time point. Aalen’s process model shows you the decomposition unfolding in time — and that is what we will build.

What the Decomposition Reveals

Three panels. Three stories your endpoint missed.

Left: the direct effect is protective throughout. Centre: the indirect effect is inert for 100 days, then spikes. Right: the total effect is attenuated. This is the Aalen/process view made visible. Neither a potential outcomes ATE nor a single-timepoint Pearl NDE can give you this temporal decomposition.

The Structure: A Sequence of DAGs

The arrows don’t just exist. They strengthen and weaken over time.

This is Aalen’s insight. The DAG is not static. The causal structure evolves. Treatment X influences mediator M, and both influence the hazard. But the weights on every arrow are functions of time. This is what it means to take the process view seriously.

Act II: Building the Process Model

Why a Generative Model?

The process view demands a model that generates trajectories, not just estimates coefficients.

A partial likelihood (Cox) conditions out time. We need time in.

A full Poisson likelihood gives us: explicit baseline hazard, time-varying coefficients, and the ability to simulate counterfactual worlds.

Three technical moves: discretise, smooth, bridge.

This is the transition from philosophy to engineering. The process view of causation requires a specific technical capability: you must be able to simulate what would have happened under alternative interventions. You cannot simulate from a partial likelihood. This is not about fanciness — it is about what the philosophical commitment demands of the model.

Discretise + Smooth

Discretise time into bins

# Each row: who was at risk, how long, what happened
df["dt"] = df["stop"] - df["start"]
bin_idx = df["bin_idx"].values

Time is no longer background noise. It’s a covariate.

Smooth with B-splines

# Random walk prior on spline coefficients
beta_raw = pm.Normal("beta_raw", 0, 1,
                     dims=('splines', 'tv'))
coef_alpha = pt.cumsum(beta_raw * sigma_smooth,
                       axis=0)
alpha_1_t = pt.dot(basis, coef_alpha[:, 1])

The cumulative sum is a random walk. The data decides how much adjacent coefficients differ.

Discretisation gives us a time index so coefficients can vary. B-splines with a random walk prior enforce smoothness — adjacent time bins are probably similar. The sigma_smooth parameter controls wiggliness. We test this with LOO cross-validation across different knot counts. Details are in the blog post.

The Poisson Bridge

\[d_{it} \sim \text{Poisson}(\lambda(t)), \quad \lambda(t) = Y_{it} \cdot \Delta t \cdot f(\text{linear predictor})\]

Binary events are counts capped at 1.

The baseline hazard becomes an explicit parameter, not a nuisance to be conditioned away.

This is what makes the model generative. We can sample from it. We can intervene on it.

The Poisson likelihood cares about intensities — how fast is risk accumulating in each bin. This gives us a full generative model we can simulate from. A partial likelihood cannot do this. The ability to simulate is what converts a process model into a counterfactual reasoning engine.

The Joint Model

Mediator equation:

\[M_t = \beta(t) \cdot X + \rho \cdot M_{t-1} + \epsilon\]

The mediator has inertia. Its dynamics are modelled, not assumed away.

Hazard equation:

\[\lambda(t) = f\big(\alpha_0(t) + \alpha_1(t) X + \alpha_2(t) M_{t-1}\big)\]

Direct and indirect paths, both time-varying.

These are estimated simultaneously. One likelihood, one posterior.

This is a structural equation model in Pearl’s sense, with time-varying coefficients in Aalen’s sense.

The joint estimation is important. Information flows between the mediator model and the hazard model. If the mediator equation fits well, the hazard model gets a cleaner signal about M’s causal contribution. You’re fitting a system, not two regressions. Note the dual identity: this is simultaneously a Pearl-style SCM (structural equations, identifiable paths) and an Aalen-style process model (time-varying coefficients, additive hazard).

The Softplus Link

\[f(x) = \ln(1 + \exp(x))\]

Ensures the hazard is always positive. Smooth. Differentiable. MCMC-friendly.

This is not Aalen’s identity link. The path decomposition is approximate on the hazard scale.

That’s fine. The decomposition is diagnostic. The causal estimand lives elsewhere.

A common objection: “the decomposition doesn’t hold exactly under a nonlinear link.” Correct. And it doesn’t matter, because the coefficients are not the final causal estimand. They show you the mechanism — the process view. The causal quantity comes from g-computation, where the link function has already been applied — the interventionist view. Two views, one model.

Diagnostics: Does the Model Fit?

Left: mediator model fits. Right: event model fits. The generative model is internally coherent.

The PPC is the one diagnostic that matters for the talk: the model can reproduce the data it was trained on, for both equations simultaneously. LOO comparisons across spline counts and forest plots showing stability of the direct effect are in the blog post.

Act III: Dissolving the Puzzle

The Path Decomposition

These time-varying coefficients show which channel dominates when:

• Direct effect: protective throughout
• Indirect effect: inert, then harmful after day 100
• Total effect: attenuated

This is the model’s internal mechanics.

The structure of mediation, not the magnitude.

But mechanism is not yet an estimand. The potential outcomes school is unimpressed.

The path decomposition is the Aalen/process contribution. It tells you which channel dominates when. But a potential outcomes theorist would say: “These are just regression coefficients. Show me a well-defined causal contrast.” Fair enough. That is what g-computation gives us next.

G-Computation: Two Counterfactual Worlds

# World 1: everyone treated
pm.set_data({'trt': np.ones(n)})
idata_trt1 = pm.sample_posterior_predictive(
    idata, var_names=['Lambda', 'obs_m', 'obs_event'])

# World 0: nobody treated
pm.set_data({'trt': np.zeros(n)})
idata_trt0 = pm.sample_posterior_predictive(
    idata, var_names=['Lambda', 'obs_m', 'obs_event'])

Two counterfactual worlds. Same population. Same model.

The mediator evolves naturally under each intervention.

When you set treatment to 1 for everyone, the mediator equation generates what M would have been under treatment. When you set it to 0, M evolves under control. These counterfactual mediator trajectories feed into the hazard equation. The survival curves reflect the full mediated process. This is g-computation: simulate, don’t derive.

The Survival Contrast

The causal estimand:

\[\bar{S}_{\text{treated}}(t) - \bar{S}_{\text{control}}(t)\]

Computed by:

1. Individual hazard trajectories for each posterior draw
2. Transform to survival: \(S_i(t) = \exp(-\Lambda_i(t))\)
3. Average across individuals
4. Compare

Defined interventionally.

Estimated processually.

Full posterior uncertainty.

This marginal survival contrast is a perfectly valid causal estimand in the potential outcomes framework. It is a comparison of two well-defined interventions on the same population. No cross-world counterfactuals required. The potential outcomes school can accept this number. But we can go further.

The NDE requires \(Y(1, M(0))\). This nests two worlds.

. . .

The potential outcomes school says: “Ill-defined.”

. . .

Pearl says: “Well-defined, via structural equations.”

. . .

The process model says: “I just simulated it.”

This is the central claim of the talk. When you have a generative model of the process — mediator dynamics, hazard dynamics, time-varying coefficients — the NDE is not a philosophical commitment to cross-world counterfactuals. It is a computational operation. You run the model with X=1 but feed in the mediator trajectories from the X=0 simulation. The philosophical objection dissolves because the model gives you a concrete mechanism for what “M(0) under Y(1)” means: the mediator trajectories the model generates when treatment is absent, plugged into the hazard equation with treatment present. No metaphysics required — just simulation from a generative model.

NDE from the Process Model

# Mediator trajectories under control (the "natural" M(0))
m_nat = idata_trt0.posterior_predictive['obs_m']

# Hazard under treatment, but with control mediator
eta_nde = (alpha_0(t) + alpha_1(t) * 1
           + alpha_2(t) * m_nat)

# Transform through the link and simulate
hazard_nde = dt * log1pexp(eta_nde)

The cross-world counterfactual is just plugging one simulation’s mediator into another simulation’s hazard.

~25% risk reduction via direct path alone. The mediator claws back ~8%.

The blog post computes this in full and gets roughly 25% risk reduction for the NDE versus 17% for the total effect. The key philosophical point: this NDE computation does not require you to believe in cross-world counterfactuals as metaphysical entities. It requires you to believe your generative model is a reasonable approximation of the data-generating process. That is a much weaker commitment — and one that every Bayesian modeller already makes.

The same model generates both counterfactual worlds.

. . .

Modeling choices affect both equally.

. . .

What survives in the difference is the causal signal.

The softplus link, the hazard extrapolation for censored subjects, the spline flexibility — all modeling choices that could distort absolute survival curves. But because the same model generates both worlds for the same population, the distortions cancel in the contrast. Process gives you the mechanism. Intervention gives you the estimand. The Bayesian implementation is what makes both available simultaneously.

What the Practitioner Does

The treatment reduces events by ~17% overall.

But: the direct effect gives ~25% reduction. The mediator claws back ~8%.

The next intervention targets the mediator after day 100.

Without the decomposition, you’d report “modest benefit” and stop.

With it, you know exactly where to push next.

This is the actionable output. The marginal contrast tells you the treatment is worth doing. The decomposition tells you it could be worth twice as much if you manage the mediator. That is information an endpoint analysis cannot provide, and information that requires both the process view and the interventionist view working together.

When to Use What

If you need total effects and have a clean experiment:

Potential outcomes + Cox is fine.

If you need mediation but effects are time-constant:

Pearl’s mediation formulas work.

If effects change over time and you suspect a mediator whose influence evolves:

Dynamic path analysis gives you the process and the estimands.

How to tell: plot Schoenfeld residuals. If they trend, you have time-varying effects. If you also have a candidate mediator, you have a dynamic path problem.

What We Built

1. A Bayesian dynamic path analysis in PyMC. A generative model of the causal process, with full posterior uncertainty.

2. G-computation that converts process to estimand — including NDE/NIE — without requiring philosophical commitment to cross-world counterfactuals.

3. A practical resolution of a real philosophical debate. The generative model dissolves the objection by making the counterfactual concrete.

The process tells you why.

The counterfactual tells you how much.

The generative model tells you both.

Thank You

Resources

• Blog post: nathanielf.github.io
• PyMC: pymc.io
• Fosen et al. (2006), “Dynamic path analysis,” Lifetime Data Analysis
• Aalen (1989), “A linear regression model for the analysis of life times,” Statistics in Medicine
• Holland (1986), “Statistics and Causal Inference,” JASA
• Pearl (2001), “Direct and Indirect Effects,” UAI
• Robins & Richardson (2010), “Alternative graphical causal models”