Start with smoking, income, and hypertension

Imagine a fixed follow-up study:

• A = smoking during the exposure window
• M = low income, measured before that exposure window
• Y = incident hypertension by the end of follow-up

Before any model, read the table in two ways: smoking’s effect within income groups, or the joint pattern across all four cells.

Question 1: effect modification

One exposure: smoking. Income defines the groups. Ask whether smoking’s effect on hypertension is the same in each income stratum.

• Additive scale: yes, smoking’s effect differs (+5 vs +15 points)
• Multiplicative scale: no, the smoking RR is 2 in both income groups

In R: read smoking’s effect within each income group

toy <- data.frame(
  income = c("Higher", "Low"),
  risk_nonsmoke = c(0.05, 0.15),
  risk_smoke    = c(0.10, 0.30)
)

toy$RD <- toy$risk_smoke - toy$risk_nonsmoke
toy$RR <- toy$risk_smoke / toy$risk_nonsmoke
toy
#>   income risk_nonsmoke risk_smoke   RD RR
#> 1 Higher          0.05        0.1 0.05  2
#> 2    Low          0.15        0.3 0.15  2

Question 2: interaction

Now make a stronger move: treat smoking and low income as two exposures and use the higher-income, non-smoking cell as the reference.

• RR for smoking alone = 10% / 5% = 2
• RR for low income alone = 15% / 5% = 3
• RR for both = 30% / 5% = 6

Multiplicative: 6 = 2 × 3 → no multiplicative interaction

Additive interaction contrast (IC): compare the smoking risk difference across income groups:

\[ IC=(R_{11}-R_{01})-(R_{10}-R_{00}) =R_{11}-R_{10}-R_{01}+R_{00}. \]

Here, IC = (30 − 15) − (10 − 5) = +10 percentage points. The additive null is IC = 0.

In R: interaction on two scales

r00 <- 0.05  # higher income, not smoking
r10 <- 0.10  # higher income, smoking
r01 <- 0.15  # low income, not smoking
r11 <- 0.30  # low income, smoking

rr10 <- r10 / r00
rr01 <- r01 / r00
rr11 <- r11 / r00

mult_ratio <- rr11 / (rr10 * rr01)
IC <- r11 - r10 - r01 + r00
round(c(mult_ratio = mult_ratio, IC = IC), 2)
#> mult_ratio         IC
#>        1.0        0.1

Same table, different question

Confounding is a third idea

A five-step workflow

1. Choose A and M. Give a substantive reason for M.
2. Check that A cannot cause M.
3. Count people and outcomes in A × M before modelling.
4. Estimate A’s effect within each M on a named scale.
5. Report the two effects, their contrast and uncertainty, and what you adjusted for.

The modifier must survive the “could A change it?” test

What do we adjust for?

Effect modification of A by M

We need to identify the effect of A on Y within levels of M.

• Control the confounding of A → Y appropriate to that question
• We do not need to identify the causal effect of M itself
• But conditioning on M still has to be causally defensible

Interaction of A and M

To interpret both as exposures, we need to identify the joint effects of A and M — a stronger requirement.

Name the scale before saying “no modification”

Week 4’s made-up table, now read for effect modification.

• RR: modified by age (1.60 vs 2.50)
• OR: not modified (4.00 vs 4.00)
• RD: not modified (0.30 vs 0.30)

Two scales are enough for the core lecture

The MASLD question we will try to answer

• A: high waist fat, WHtR >= 0.6
• M: obesity
• Y: all-cause death
• L: age, sex, race/ethnicity, smoking, sedentary time

This lecture’s analysis: the same 5,911 rows as Week 4, adjusted for L, unweighted (weights are Week 6). Slides that use other rows say so.

Why the MASLD example becomes causally difficult

The paper’s phenotype analysis is reasonable for risk stratification. Extra problems appear when we reuse the same data for a causal effect-modification question:

• BMI/obesity and WHtR/waist were measured at the same examination; the data do not establish which came first
• both are related measures of adiposity, so treating obesity as clearly upstream of waist fat is a substantive assumption
• the cohort is restricted to people with MASLD; MASLD ascertainment depends partly on BMI and waist-related information
• one important cell is sparse: the paper had 5 deaths in Group I; our complete-case analysis has 4

Causal audit: can obesity really be the modifier?

What the effect-modification interpretation assumes

For the obese/not-obese strata to be interpretable as modifiers of waist fat’s effect, we need:

• a time-ordering story in which obesity is established before the waist-fat exposure of interest
• no important unmeasured common cause of obesity and waist fat that becomes problematic when we stratify on obesity
• L to be adequate for confounding of waist fat → death within each obesity stratum
• conditioning on MASLD not to induce important selection bias in the within-stratum waist-fat comparison

Count first

• Group I anchors the obese waist-fat comparison II vs I and several phenotype contrasts; it contains only 4 deaths
• Among 2,432 obese women, only 57 have lower waist fat: limited support for that exposure level

In R: count before you model

table(group = d$group, died = d$dead)
#>      died
#> group    0    1
#>   I    351    4
#>   II  3903  362
#>   III  665   67
#>   IV   458  101

with(subset(d, obese == 1 & female == 1),
     table(high_waist = central))
#> high_waist
#>    0    1
#>   57 2375

One product-term model, two stratum-specific HRs

fit <- coxph(Surv(time_yr, dead) ~ obese * central +
               age + female + race + smk_current + sed, data = d)

• Among the not obese: HR for high vs lower waist fat = 1.09 (0.79–1.50)
• Among the obese: HR = 2.99 (1.11–8.09)
• Ratio of HRs = 2.74 (0.97–7.71) (Wald 95% CI)
• Tests for allowing the HR to differ: Wald p = 0.06; likelihood-ratio p = 0.03

\[ \exp(\beta_{A\times M}) =\frac{ HR_\text{waist, obese} }{ HR_\text{waist, not obese} } =\frac{2.99}{1.09}=2.74. \]

In R: recover both HRs and test the product term

b <- coef(fit)

round(exp(b[["central"]]), 2)                      # not obese: IV vs III
#> [1] 1.09

round(exp(b[["central"]] + b[["obese:central"]]), 2)   # obese: II vs I
#> [1] 2.99

fit_main <- coxph(Surv(time_yr, dead) ~ obese + central +
                    age + female + race + smk_current + sed, data = d)
round(anova(fit_main, fit)$`Pr(>|Chi|)`[2], 3)     # likelihood-ratio p
#> [1] 0.029

Why compute standardised six-year risks?

This is our re-analysis, not a result reported by the paper. The paper reports Kaplan-Meier curves and Cox hazard ratios. We add standardised risks because effect modification can also be asked on an absolute, additive scale.

For a survival outcome we must choose a time horizon. At 6 years, within each obesity stratum:

1. take the same people with their observed age, sex, race/ethnicity, smoking, and sedentary time
2. set everyone’s waist-fat exposure to lower, predict risk by 6 years, and average
3. set those same people’s waist-fat exposure to high, predict again, and average
4. subtract the two average risks

Additive scale: the standardised six-year results

Difference of risk differences = +3.4 percentage points
Bootstrap 95% CI: −0.1 to +6.6 points

In R: use the same people twice, changing only waist fat

# fit is the Cox model from the preceding slides
risk6 <- function(people, waist) {
  nd <- people                 # same people, same L values
  nd$central <- waist          # change only high-waist exposure
  s6 <- summary(survfit(fit, newdata = nd), times = 6)$surv
  mean(1 - s6)                 # average predicted 6-year risk
}

ob    <- subset(d, obese == 1)
notob <- subset(d, obese == 0)

r <- c(
  obese_low     = risk6(ob,    0),
  obese_high    = risk6(ob,    1),
  notobese_low  = risk6(notob, 0),
  notobese_high = risk6(notob, 1)
)
round(100 * r, 1)
#>     obese_low    obese_high  notobese_low notobese_high
#>           2.3           6.4           9.1           9.8

In R: turn the four predictions into risk differences

rd <- 100 * c(
  obese     = r[["obese_high"]]    - r[["obese_low"]],
  not_obese = r[["notobese_high"]] - r[["notobese_low"]]
)
round(rd, 1)
#>     obese not_obese
#>       4.1       0.7

round(rd[["obese"]] - rd[["not_obese"]], 1)   # difference of risk differences
#> [1] 3.4

Put both scales on one page

The paper’s IV vs I comparison asks a different question

Our effect-modification question uses the two within-row comparisons: II vs I and IV vs III. The paper’s IV vs I contrast is diagonal across the 2 × 2 table.

If obesity is a second exposure, the estimand changes

Now deliberately treat both variables as exposures, keeping the coding consistent:

• A = high waist fat (1 = high, 0 = lower)
• M = obesity (1 = obese, 0 = not obese)
• each person now has four joint potential outcomes, \(Y(a,m)\)

The interaction contrasts ask whether the waist-fat effect changes when obesity is also changed:

• Multiplicative: ratio of waist-fat HRs = 2.74 (0.97–7.71)
• Additive: \(IC=(R_{11}-R_{01})-(R_{10}-R_{00})\); using one common standardised population, IC = +3.9 points

In R: interaction with the same A and M coding

# Multiplicative interaction from the Cox product term
b <- coef(fit)
round(exp(b[["obese:central"]]), 2)   # HR(waist | obese) / HR(waist | not obese)
#> [1] 2.74

In R: the additive interaction contrast

For the additive joint-exposure contrast, put all four joint states in the same target population:

# A = high waist, M = obesity
# R00 = Group III; R10 = IV; R01 = I; R11 = II
risk6_all <- function(waist, obese) {   # everyone, set to one joint state
  nd <- d; nd$central <- waist; nd$obese <- obese
  mean(1 - summary(survfit(fit, newdata = nd), times = 6)$surv)
}
R <- c(R00 = risk6_all(0, 0), R10 = risk6_all(1, 0),
       R01 = risk6_all(0, 1), R11 = risk6_all(1, 1))
round(100 * R, 1)
#> R00 R10 R01 R11
#> 6.9 7.5 2.5 7.0

IC <- 100 * (R[["R11"]] - R[["R01"]] - R[["R10"]] + R[["R00"]])
round(IC, 1)   # percentage points
#> [1] 3.9

The right survival tool depends on the scale

• Do not turn a censored survival outcome into “ever died” and fit ordinary logistic regression
• A fixed-follow-up binary outcome is different; then a binary regression model may be appropriate

A dramatic subgroup claim can be built on almost no information

These are crude diagonal phenotype comparisons (Group IV vs Group I), fitted separately by sex. They are not sex-specific waist-fat effects.

Crude, unweighted, on the 6,048-row locked domain (dom in the code), so Group I has 375 people here, not the 355 of the 5,911.

In R: reproduce the crude IV-vs-I HRs

dom$group <- relevel(factor(dom$group), ref = "I")

fit_men <- coxph(Surv(time_yr, dead) ~ group,
                 data = subset(dom, female == 0))
fit_women <- coxph(Surv(time_yr, dead) ~ group,
                   data = subset(dom, female == 1))

round(summary(fit_men)$conf.int[
  "groupIV", c("exp(coef)", "lower .95", "upper .95")], 2)
#> exp(coef) lower .95 upper .95
#>     36.14      8.86    147.48

round(summary(fit_women)$conf.int[
  "groupIV", c("exp(coef)", "lower .95", "upper .95")], 2)
#> exp(coef) lower .95 upper .95
#>      4.11      0.99     17.12

In R: now count the reference events

with(subset(dom, group == "I"), table(female, died = dead))
#>       died
#> female   0   1
#>      0 315   2
#>      1  56   2

Ask the subgroup question cleanly

For the locked waist-fat contrast IV vs III, with sex as the modifier (all unweighted):

Additive six-year RDs are also imprecise: about +0.8 points in men and +0.7 in women (P3, Model 1).

1,296: P3’s paper-replication refits (Groups III + IV complete on the paper’s Model 1). 1,291 (168 deaths): Groups III + IV of our 5,911. CIs use Wald; p-values use likelihood-ratio tests.

In R: one test for sex modification on the locked contrast

s34 <- subset(d, group %in% c("III", "IV"))
s34$IV <- as.integer(s34$group == "IV")

f0 <- coxph(Surv(time_yr, dead) ~ IV + female +
              age + race + smk_current + sed, data = s34)
f1 <- coxph(Surv(time_yr, dead) ~ IV * female +
              age + race + smk_current + sed, data = s34)

round(summary(f1)$conf.int["IV:female", c(1, 3, 4)], 2)   # ratio of HRs
#> exp(coef) lower .95 upper .95
#>      1.13      0.48      2.66

round(anova(f0, f1)$`Pr(>|Chi|)`[2], 2)   # likelihood-ratio p
#> [1] 0.78

Honest subgroup reporting

1. Pre-specify a plausible modifier and explain why A cannot cause it
2. State the question: effect modification or joint interaction
3. Count people and outcomes in every A × M cell
4. Report A’s effect within each M on a named scale
5. Give the contrast and interval, not just a p-value
6. Do not invent a biological explanation for a noisy subgroup difference

The algorithm to take into your own paper

Key takeaways

• Effect modification: A’s effect differs across levels of M; for the usual causal subgroup interpretation, A should not determine M
• Interaction: the joint effect of A and M departs from additivity or multiplicativity
• The scale matters: the same data can modify one measure and not another
• The question determines confounding control; a product term does not
• With survival outcomes: Cox product term for the multiplicative scale; standardised risks at a fixed time for an additive scale
• In MASLD, the point estimates suggest more waist-fat harm among the obese, but the conclusion is fragile because the obese lower-waist cell anchoring II vs I has only 4 deaths

Appendix: optional — one model, two parameterisations

fit_joint <- coxph(Surv(time_yr, dead) ~ group +
                     age + female + race + smk_current + sed, data = d)
fit_prod  <- coxph(Surv(time_yr, dead) ~ obese * central +
                     age + female + race + smk_current + sed, data = d)
all.equal(fit_joint$loglik, fit_prod$loglik)
#> [1] TRUE

round(exp(coef(fit_prod))[c("obese", "central", "obese:central")], 2)
#>         obese       central obese:central
#>          0.34          1.09          2.74

The four-level factor and product-term parameterisations give the same fitted model.

Appendix: RERI, AP, S, and IC

For a four-cell joint-exposure table with a clearly defined (00) state:

• RERI = RR₁₁ − RR₁₀ − RR₀₁ + 1; null = 0
• AP = RERI / RR₁₁; null = 0
• S = (RR₁₁ − 1) / [(RR₁₀ − 1) + (RR₀₁ − 1)]; null = 1
• IC = R₁₁ − R₁₀ − R₀₁ + R₀₀; null = 0

In the core MASLD slides we keep A = high waist fat and M = obesity throughout. With the four joint states standardised to one common population, the six-year risk-scale IC = +3.9 points.

Appendix: why diabetes is a risky modifier here

tab <- table(diabetes = d$t2dm, high_waist = d$central)
round(100 * prop.table(tab, 2))
#>         high_waist
#> diabetes  0  1
#>        0 87 72
#>        1 13 28

Diabetes is more common among people with high waist fat (28% vs 13%). That does not prove waist fat causes diabetes, but it makes a same-visit diabetes stratum a poor choice for the simple pre-exposure-modifier interpretation.

Appendix: what splitting by diabetes produces

f_no  <- coxph(Surv(time_yr, dead) ~ central + obese + age + female +
                 race + smk_current + sed, data = subset(d, t2dm == 0))
f_yes <- coxph(Surv(time_yr, dead) ~ central + obese + age + female +
                 race + smk_current + sed, data = subset(d, t2dm == 1))

The waist-fat HR is about 1.23 without diabetes and 0.93 with diabetes, but neither number automatically has a clean causal subgroup interpretation if waist fat can affect diabetes.

Where to read more

References (1/2)

References (2/2)

Appendix: the two questions as dagitty code

Paste into the Model code box at dagitty.net. With Waist as the only exposure, dagitty returns {L, Obese}. Mark Obese as an exposure too, and it reports that no adjustment set exists — until you remove the latent tag from U.

dag {
  Waist [exposure]  Death [outcome]  U [latent]
  L -> Obese  L -> Waist  L -> Death
  U -> Obese  U -> Death
  Obese -> Waist  Obese -> Death  Waist -> Death
}

L = age, sex, race/ethnicity, smoking, sedentary time. U = serious illness and disease that came before obesity (Week 4’s Ill and Before). Then test the assumptions: add G [latent] G -> Obese G -> Waist, or FLI [adjusted] Waist -> FLI Waist -> TG TG -> FLI TG -> Death, and watch every adjustment set disappear.