4  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 the causaldata package.

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

We keep the treatment qsmk, the outcome wt82_71 (weight change in kg), and the nine baseline covariates commonly used in NHEFS smoking-cessation examples: sex, race, age, school (years of schooling), smokeintensity (cigarettes per day in 1971), smokeyrs (years of smoking), exercise and active (ordinal activity codes), and wt71 (baseline weight in kg). All are measured at baseline, before quitting.

nhefs_cols <- c(
  "qsmk", "wt82_71",
  "sex", "race", "age", "school",
  "smokeintensity", "smokeyrs",
  "exercise", "active", "wt71"
)
nhefs <- causaldata::nhefs_complete |>
  dplyr::select(all_of(nhefs_cols)) |>
  dplyr::mutate(across(everything(), as.numeric))

c(n = nrow(nhefs), variables = ncol(nhefs))

        n variables 
     1566        11 

The complete-case sample has 1,566 people and these 11 columns. The first five rows:

head(nhefs, 5) |>
  dplyr::mutate(across(everything(), \(x) signif(x, 4))) |>
  kable()

qsmk	wt82_71	sex	race	age	school	smokeintensity	smokeyrs	exercise	active	wt71
0	-10.090	1	2	42	7	30	29	3	1	79.04
0	2.605	1	1	36	9	20	24	1	1	58.63
0	9.414	2	2	56	11	20	26	3	1	56.81
0	4.990	1	2	68	5	3	53	3	2	59.42
0	4.989	1	1	40	11	20	19	2	2	87.09

Weight change is already in kilograms and can be negative, so the estimand is in kilograms too.

4.2 Adjustment DAG

The graph is written and queried with dagitty (Textor et al. 2016), applying Pearl’s backdoor criterion (1995).

Start with the grouped graph

X -> A -> Y
X -----> Y

where 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.

conceptual_graph <- dagitty("dag {
  X -> A
  X -> Y
  A -> Y
}")
coordinates(conceptual_graph) <- list(
  x = c(X = 1, A = 0, Y = 2),
  y = c(X = 1, A = 0, Y = 0)
)

ggdag(conceptual_graph) + theme_dag_blank()

Figure 4.1: Conceptual adjustment DAG: baseline X confounds A → Y.

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. We build it and count what it contains.

baseline_vars <- setdiff(nhefs_cols, c("qsmk", "wt82_71"))

build_grouped_graph <- function(baseline) {
  edges_X_to_A <- paste0(baseline, " -> qsmk")
  edges_X_to_Y <- paste0(baseline, " -> wt82_71")
  body <- paste(c(edges_X_to_A, edges_X_to_Y, "qsmk -> wt82_71"),
                collapse = "\n  ")
  dagitty(paste0("dag {\n  ", body, "\n}"))
}

nhefs_graph <- build_grouped_graph(baseline_vars)
c(vertices = length(names(nhefs_graph)),
  directed_edges = length(edges(nhefs_graph)$v))

      vertices directed_edges 
            11             19 

Eleven vertices and 19 directed edges: two out of each of the nine baseline variables, plus qsmk -> wt82_71.

4.3 Identification

dagitty::adjustmentSets checks the graph before estimation. In this adjustment DAG, the minimal sufficient adjustment set is the full baseline:

adj_sets <- adjustmentSets(nhefs_graph, exposure = "qsmk",
                           outcome = "wt82_71", type = "minimal")
adj_sets

{ active, age, exercise, race, school, sex, smokeintensity, smokeyrs,
  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 AIPW

The estimator is that of Robins, Rotnitzky and Zhao (1994). The NHEFS analysis it mirrors is Hernán and Robins (2020).

The book keeps identification and estimation as separate steps. The graph supplies the adjustment set; the estimator uses that set.

aipw_backdoor <- function(data, treatment, outcome, adjust) {
  fmla_y <- as.formula(paste0(outcome, " ~ ", treatment, " * (",
                              paste(adjust, collapse = " + "), ")"))
  fmla_a <- as.formula(paste0(treatment, " ~ ",
                              paste(adjust, collapse = " + ")))

  mu_fit <- lm(fmla_y, data = data)
  pi_fit <- glm(fmla_a, data = data, family = binomial())

  d1 <- data; d1[[treatment]] <- 1
  d0 <- data; d0[[treatment]] <- 0
  mu1 <- predict(mu_fit, newdata = d1)
  mu0 <- predict(mu_fit, newdata = d0)
  pi1 <- predict(pi_fit, newdata = data, type = "response")
  # Truncate (winsorize) the propensity score away from 0/1. AIPW divides by
  # pi1 and (1-pi1), so an estimated propensity near the boundary produces
  # huge inverse-probability weights and blows up the estimator's variance.
  # Bounding at [0.01, 0.99] trades a small amount of bias (we no longer use
  # the true estimated propensity for the most extreme units) for much
  # lower variance -- a standard and usually favorable trade-off, but worth
  # flagging explicitly since it is a bias-variance choice, not a free lunch.
  pi1 <- pmin(pmax(pi1, 0.01), 0.99)

  w  <- data[[treatment]]; y <- data[[outcome]]
  if1 <- w     * (y - mu1) / pi1       + mu1
  if0 <- (1-w) * (y - mu0) / (1 - pi1) + mu0
  psi  <- if1 - if0

  est <- mean(psi); se <- sd(psi) / sqrt(length(psi))
  ci  <- est + c(-1, 1) * qnorm(0.975) * se
  data.frame(estimand = "Quit smoking vs continued smoking",
             estimate = signif(est, 4),
             lower_95 = signif(ci[1], 4),
             upper_95 = signif(ci[2], 4),
             se       = signif(se, 4))
}

aipw_backdoor(nhefs, "qsmk", "wt82_71", adjust = baseline_vars) |> kable()

estimand	estimate	lower_95	upper_95	se
Quit smoking vs continued smoking	3.293	2.317	4.269	0.4981

Quitters gain 3.29 kg more than continuing smokers, with a standard error of 0.50 and a 95% interval of \([2.32, 4.27]\). Under the DAG this is the average effect of quitting versus continuing on weight change in the complete-case NHEFS population. The causal interpretation depends on the graph, especially the assumption that the measured baseline variables block the relevant backdoor paths.

The sd(psi) / sqrt(n) standard error is the AIPW influence-function SE and is valid here because the nuisances are parametric logistic/linear fits. If you replace them with flexible machine learners, use cross-fitting and an influence-function (or bootstrap) variance instead; the plug-in SE on full-sample ML fits is not generally valid.

4.5 Sensitivity Graph

The sensitivity parameters are those of Cinelli and Hazlett (2020).

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. We rebuild the same graph with that one edge added and count the sufficient adjustment sets it admits.

build_sensitivity_graph <- function(baseline) {
  edges_X_to_A <- paste0(baseline, " -> qsmk")
  edges_X_to_Y <- paste0(baseline, " -> wt82_71")
  body <- paste(c(edges_X_to_A, edges_X_to_Y,
                  "qsmk -> wt82_71", "qsmk <-> wt82_71"),
                collapse = "\n  ")
  dagitty(paste0("dag {\n  ", body, "\n}"))
}

sensitivity_graph <- build_sensitivity_graph(baseline_vars)
sens_adj <- adjustmentSets(sensitivity_graph, exposure = "qsmk",
                           outcome = "wt82_71", type = "minimal")
cat("Number of sufficient adjustment sets:", length(sens_adj), "\n")

Number of sufficient adjustment sets: 0 

The count is zero. One edge changes the identifying conclusion: under this graph the observed NHEFS variables are not enough to identify the effect, and no adjustment set exists. The data are the same as in the previous section, and the 3.29 kg is still computable, but nothing in this graph justifies reading it as a causal effect. This is why the graph has to come before the estimator.

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. We write that graph out, ask for its minimal adjustment set, and check whether the set equals the full baseline we used before.

structured_graph <- dagitty("dag {
  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
}")

struct_adj <- adjustmentSets(structured_graph, exposure = "qsmk",
                             outcome = "wt82_71", type = "minimal")
struct_adj

{ active, age, exercise, race, school, sex, smokeintensity, smokeyrs,
  wt71 }

data.frame(
  matches_full_baseline = setequal(unlist(struct_adj[[1]]), baseline_vars)
) |> kable()

matches_full_baseline
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 minimal sufficient set is therefore still the full baseline, the identifying functional is unchanged, and the estimate is still 3.29 kg.

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.

Cinelli, Carlos, and Chad Hazlett. 2020. “Making Sense of Sensitivity: Extending Omitted Variable Bias.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 82 (1): 39–67.

Hernán, Miguel A., and James M. Robins. 2020. Causal Inference: What If. Chapman & Hall/CRC.

Pearl, Judea. 1995. “Causal Diagrams for Empirical Research.” Biometrika 82 (4): 669–88.

Robins, James M., Andrea Rotnitzky, and Lue Ping Zhao. 1994. “Estimation of Regression Coefficients When Some Regressors Are Not Always Observed.” Journal of the American Statistical Association 89 (427): 846–66.

Textor, Johannes, Benito van der Zander, Mark S. Gilthorpe, Maciej Liśkiewicz, and George T. H. Ellison. 2016. “Robust Causal Inference Using Directed Acyclic Graphs: The R Package dagitty.” International Journal of Epidemiology 45 (6): 1887–94.