Code
using CausalGraphs
using CausalEstimate
import CausalGraphs: identify, make_graph, to_mermaid, markov_pillow
import CausalEstimate: estimate, confint, pvalue
using DataFrames
using DelimitedFiles
using Downloads
using Logging
using Random4 Applied Graph Workflow: Smoking Cessation
The previous chapter used small graph examples. Here I apply the same graph-identification-estimation workflow to a real data set: the complete-case National Health and Nutrition Examination Survey I Epidemiologic Follow-up Study (NHEFS) data distributed through causaldata and mirrored by RDatasets.
The question is:
What is the effect of quitting smoking on weight change from 1971 to 1982?The treatment is qsmk, an indicator for quitting smoking between baseline and follow-up. The outcome is wt82_71, weight change over the follow-up period. This is observational data. The effect is not identified by the data alone. It is identified only after we state the adjustment assumptions in the graph.
4.1 Data
For this example, keep the treatment, outcome, and baseline covariates that are commonly used in NHEFS smoking-cessation examples.
Code
function read_rdatasets_csv(url, cols)
local_file = Downloads.download(url; timeout = 120)
x, header = readdlm(local_file, ',', Any, '\n'; header = true)
raw = DataFrame(x, Symbol.(vec(header)))
DataFrame((c => Float64.(raw[!, c]) for c in cols)...)
end
nhefs_url = "https://vincentarelbundock.github.io/Rdatasets/csv/causaldata/nhefs_complete.csv"
nhefs_cols = [
:qsmk, :wt82_71,
:sex, :race, :age, :school,
:smokeintensity, :smokeyrs,
:exercise, :active, :wt71,
]
nhefs = read_rdatasets_csv(nhefs_url, nhefs_cols);
(n = nrow(nhefs), variables = ncol(nhefs))(n = 1566, variables = 11)Code
preview_table(nhefs)5×11 DataFrame
| Row | qsmk | wt82_71 | sex | race | age | school | smokeintensity | smokeyrs | exercise | active | wt71 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 0.0 | -10.09 | 0.0 | 1.0 | 42.0 | 7.0 | 30.0 | 29.0 | 2.0 | 0.0 | 79.04 |
| 2 | 0.0 | 2.605 | 0.0 | 0.0 | 36.0 | 9.0 | 20.0 | 24.0 | 0.0 | 0.0 | 58.63 |
| 3 | 0.0 | 9.414 | 1.0 | 1.0 | 56.0 | 11.0 | 20.0 | 26.0 | 2.0 | 0.0 | 56.81 |
| 4 | 0.0 | 4.99 | 0.0 | 1.0 | 68.0 | 5.0 | 3.0 | 53.0 | 2.0 | 1.0 | 59.42 |
| 5 | 0.0 | 4.989 | 0.0 | 0.0 | 40.0 | 11.0 | 20.0 | 19.0 | 1.0 | 1.0 | 87.09 |
4.2 Adjustment DAG
Start with the grouped graph
X -> A -> Y
X -----> Ywhere X is the baseline covariate set, A is quitting smoking, and Y is later weight change. This graph says the baseline variables may affect both quitting and later weight change, and that there is no unmeasured common cause of quitting and later weight change after conditioning on those baseline variables.
Code
conceptual_graph = make_graph(
vertices = [:X, :A, :Y],
di_edges = [(:X, :A), (:X, :Y), (:A, :Y)],
);Code
println("```{mermaid}")
println(to_mermaid(conceptual_graph; direction = "LR"))
println("```")flowchart LR n1_X["X"] n2_A["A"] n3_Y["Y"] n1_X --> n2_A n1_X --> n3_Y n2_A --> n3_Y classDef fixed stroke:#111827,stroke-width:2px linkStyle 0 stroke:#2563eb,stroke-width:1.8px linkStyle 1 stroke:#2563eb,stroke-width:1.8px linkStyle 2 stroke:#2563eb,stroke-width:1.8px
For estimation, expand X into the observed NHEFS columns. The expanded graph is still an adjustment DAG: each baseline variable is allowed to point into quitting and into later weight change.
Code
nhefs_vertices = [
:sex, :race, :age, :school,
:smokeintensity, :smokeyrs,
:exercise, :active, :wt71,
:qsmk, :wt82_71,
]
baseline = setdiff(nhefs_vertices, [:qsmk, :wt82_71])
nhefs_edges = vcat(
[(w, :qsmk) for w in baseline],
[(w, :wt82_71) for w in baseline],
[(:qsmk, :wt82_71)],
)
nhefs_graph = make_graph(vertices = nhefs_vertices, di_edges = nhefs_edges);
(vertices = length(nhefs_vertices), directed_edges = length(nhefs_edges))(vertices = 11, directed_edges = 19)4.3 Identification
CausalGraphs.identify checks the graph before estimation. In this adjustment DAG, the quitting node is a-fixable. The Markov pillow is the sufficient adjustment set:
Code
nhefs_id = identify(nhefs_graph, :qsmk, :wt82_71)
DataFrame(
strategy = [string(nhefs_id.strategy)],
markov_pillow = [join(string.(markov_pillow(nhefs_graph, :qsmk; treatment = :qsmk)), ", ")],
)1×2 DataFrame
| Row | strategy | markov_pillow |
|---|---|---|
| String | String | |
| 1 | a_fixable | sex, race, age, school, smokeintensity, smokeyrs, exercise, active, wt71 |
The graph therefore identifies the causal mean by the adjustment formula
\[ E[Y(a)] = E_X\{E(Y \mid A=a, X)\}. \tag{4.1}\]
This formula is not a modeling claim by itself. It is the observed-data functional implied by the graph. Estimation still requires nuisance models for the outcome regression and treatment mechanism.
4.4 Estimation with CausalEstimate
The book keeps identification and estimation as separate steps. Here CausalGraphs supplies the graph result and CausalEstimate estimates the identified functional. Because the graph is a-fixable, GraphID routes to the AIPW average treatment effect machinery.
Code
Random.seed!(2026)
nhefs_result = estimate(
ATE(outcome = :wt82_71, treatment = :qsmk),
GraphID(graph = nhefs_graph),
AIPW(crossfit = 5),
nhefs,
)
effect_table("Quit smoking vs continued smoking", nhefs_result)1×5 DataFrame
| Row | estimand | estimate | lower_95 | upper_95 | se |
|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | |
| 1 | Quit smoking vs continued smoking | 3.452 | 1.986 | 4.918 | 0.748 |
Under the DAG, the estimate is the average effect of quitting smoking versus continuing smoking on weight change in the complete-case NHEFS population. The positive estimate is consistent with weight gain after quitting. The causal interpretation depends on the graph, especially the assumption that the measured baseline variables block the relevant backdoor paths.
Two settings deserve comment, because both move the number.
Folds. crossfit = 5 is the package default and this chapter’s seed makes the comparison reproducible. Dropping to crossfit = 2 gives 3.866 with a standard error of 1.331 — a 95% interval of \([1.26, 6.48]\), more than five kilograms wide. Five folds gives 3.452 (0.748) and ten gives 3.312 (0.737). Two folds fits each nuisance on 783 observations, and the resulting noise in the out-of-fold predictions inflates the point estimate and the influence-function variance at the same time. Use the default, or more.
Nuisance estimation. The companion R chapter fits the same adjustment set by full-sample parametric nuisances (a linear outcome model interacted with treatment, a logistic propensity) and reports 3.293 with a standard error of 0.498. The propensity scores in this data set stay inside \([0.05, 0.75]\), so there is no overlap problem to explain the gap: the wider interval here is the variance price of cross-fitting flexible learners rather than trusting a correctly specified parametric model. When the parametric model is right, it wins; the cross-fit estimator buys insurance against its being wrong.
4.5 Sensitivity Graph
Now encode a different assumption: quitting and later weight change share an unobserved common cause even after conditioning on the measured baseline variables. In an ADMG this is a bidirected edge, qsmk <-> wt82_71.
Code
sensitivity_graph = make_graph(
vertices = nhefs_vertices,
di_edges = nhefs_edges,
bi_edges = [(:qsmk, :wt82_71)],
)
sensitivity_id = identify(sensitivity_graph, :qsmk, :wt82_71);Code
println("```{mermaid}")
println(to_mermaid(sensitivity_graph; direction = "LR"))
println("```")flowchart LR n1_sex["sex"] n2_race["race"] n3_age["age"] n4_school["school"] n5_smokeintensity["smokeintensity"] n6_smokeyrs["smokeyrs"] n7_exercise["exercise"] n8_active["active"] n9_wt71["wt71"] n10_qsmk["qsmk"] n11_wt82_71["wt82_71"] n1_sex --> n10_qsmk n2_race --> n10_qsmk n3_age --> n10_qsmk n4_school --> n10_qsmk n5_smokeintensity --> n10_qsmk n6_smokeyrs --> n10_qsmk n7_exercise --> n10_qsmk n8_active --> n10_qsmk n9_wt71 --> n10_qsmk n1_sex --> n11_wt82_71 n2_race --> n11_wt82_71 n3_age --> n11_wt82_71 n4_school --> n11_wt82_71 n5_smokeintensity --> n11_wt82_71 n6_smokeyrs --> n11_wt82_71 n7_exercise --> n11_wt82_71 n8_active --> n11_wt82_71 n9_wt71 --> n11_wt82_71 n10_qsmk --> n11_wt82_71 n10_qsmk <--> n11_wt82_71 classDef fixed stroke:#111827,stroke-width:2px linkStyle 0 stroke:#2563eb,stroke-width:1.8px linkStyle 1 stroke:#2563eb,stroke-width:1.8px linkStyle 2 stroke:#2563eb,stroke-width:1.8px linkStyle 3 stroke:#2563eb,stroke-width:1.8px linkStyle 4 stroke:#2563eb,stroke-width:1.8px linkStyle 5 stroke:#2563eb,stroke-width:1.8px linkStyle 6 stroke:#2563eb,stroke-width:1.8px linkStyle 7 stroke:#2563eb,stroke-width:1.8px linkStyle 8 stroke:#2563eb,stroke-width:1.8px linkStyle 9 stroke:#2563eb,stroke-width:1.8px linkStyle 10 stroke:#2563eb,stroke-width:1.8px linkStyle 11 stroke:#2563eb,stroke-width:1.8px linkStyle 12 stroke:#2563eb,stroke-width:1.8px linkStyle 13 stroke:#2563eb,stroke-width:1.8px linkStyle 14 stroke:#2563eb,stroke-width:1.8px linkStyle 15 stroke:#2563eb,stroke-width:1.8px linkStyle 16 stroke:#2563eb,stroke-width:1.8px linkStyle 17 stroke:#2563eb,stroke-width:1.8px linkStyle 18 stroke:#2563eb,stroke-width:1.8px linkStyle 19 stroke:#dc2626,stroke-width:1.8px
Code
DataFrame(strategy = [string(sensitivity_id.strategy)])1×1 DataFrame
| Row | strategy |
|---|---|
| String | |
| 1 | not_identified |
The identifying conclusion changes. Under this graph, the observed NHEFS variables are not enough to identify the effect. This is why the graph has to come before the estimator. The same data can imply different causal answers under different assumptions.
4.6 Structured Baseline Graph
The simple adjustment graph is often the clearest starting point. A more detailed DAG can record baseline ordering: demographics may precede education, smoking history, activity, and baseline weight; those variables may then affect quitting and later weight change.
Code
structured_edges = [
(:age, :school), (:age, :smokeyrs), (:age, :wt71),
(:age, :active), (:age, :exercise), (:age, :qsmk), (:age, :wt82_71),
(:sex, :smokeintensity), (:sex, :smokeyrs), (:sex, :wt71),
(:sex, :qsmk), (:sex, :wt82_71),
(:race, :school), (:race, :smokeintensity),
(:race, :qsmk), (:race, :wt82_71),
(:school, :smokeintensity), (:school, :qsmk), (:school, :wt82_71),
(:smokeyrs, :smokeintensity), (:smokeyrs, :qsmk), (:smokeyrs, :wt82_71),
(:smokeintensity, :qsmk), (:smokeintensity, :wt82_71),
(:exercise, :wt71), (:exercise, :qsmk), (:exercise, :wt82_71),
(:active, :wt71), (:active, :qsmk), (:active, :wt82_71),
(:wt71, :qsmk), (:wt71, :wt82_71),
(:qsmk, :wt82_71),
]
structured_graph = make_graph(vertices = nhefs_vertices, di_edges = structured_edges);
structured_id = identify(structured_graph, :qsmk, :wt82_71)
DataFrame(
strategy = [string(structured_id.strategy)],
same_adjustment_set = [
Set(markov_pillow(structured_graph, :qsmk; treatment = :qsmk)) == Set(baseline)
],
)1×2 DataFrame
| Row | strategy | same_adjustment_set |
|---|---|---|
| String | Bool | |
| 1 | a_fixable | true |
The check returns true, and the reason is worth spelling out, because it is not an accident of this data set. Every arrow we added runs between baseline variables, and every baseline variable is pre-treatment and in the adjustment set. A back-door path from qsmk to wt82_71 has to leave qsmk through an incoming arrow, travel among the baseline variables, and then enter wt82_71 from one of them. That last baseline variable has an arrow pointing into the outcome, so on this path it is a non-collider, and we are conditioning on it. The path is blocked whatever the arrows upstream of it look like. The measured parents of qsmk are therefore still the full baseline, and the identifying functional is unchanged.
That does not make the detail useless. It records what we believe about the baseline ordering, which a reader may want to see. But each arrow is a causal assumption that can be wrong, and here none of them buys anything for this estimand. The grouped adjustment graph asserts less and identifies the same thing.