Substitution analyses

Coding examples for “what if we replaced X with Y” - dietary substitution models on a linked cohort

Published

July 21, 2026

Under development. These are coding examples, not a methods lecture - see the reference at the bottom for the theory, assumptions and interpretation.

Substitution analysis answers a question ordinary regression cannot: what happens if we replace one unit of A with one unit of B, while the total stays the same? A and B can be foods (A = red meat, B = fish) or drugs (A = NSAID, B = paracetamol).

It makes sense whenever the exposure is a simultaneous composition - the parts coexist and are bound by a whole, so you cannot add more of one without taking away from another (foods within total energy, or minutes within a 24-hour day). Drug classes are the exception: switching happens over time, not as a simultaneous composition, so the code below does not apply to them - see the takeaway under the table.

The code examples below use dietary examples, but the same substitution idea carries over to other compositional exposures. What changes is the constraint (“total”) and how you build the exposure - and for drugs the study design takes over entirely (see the takeaway below the table):

	Dietary substitution	Physical activity / time-use	Drug-class switching
Components	foods (g / servings)	minutes in activity types	drug classes (DDDs, fills, days on drug)
“Total” held constant	total energy intake	total time (1440 min/day)	total treatment time / exposure
Exposure timing	usually baseline questionnaire	baseline or repeated	time-varying (LMDB fills -> episodes)
Natural design	cross-sectional + logistic/Cox	regression / time-use models	active-comparator new-user / target trial

Takeaway. The substitution regression transfers cleanly to a simultaneous composition held at a total - foods within total energy, activity minutes within a day. Drugs are different: drug-class switching is a time-varying treatment change, not a simultaneous composition, so the code above does not apply (it would ignore when people switch - immortal time, treatment-confounder feedback). There the substitution question (“the effect of using drug B instead of A”) is answered by the design: an active-comparator, new-user study, ideally a target-trial emulation, where the B-vs-A contrast in the outcome model already is the substitution effect. The partition/leave-one-out code only carries over if you genuinely have a drug composition - simultaneous shares of several classes within total treatment - which is uncommon.

Compositional data (CoDA). Physical activity over a day is a closed composition (sleep + sedentary + light + vigorous = 1440 min), where you use formal compositional data analysis (CoDA) with log-ratios (e.g. the isometric log-ratio, ilr; the compositions R package). Diet is not strictly closed - total energy varies between people - so nutritional epidemiology uses the energy-adjusted substitution models below rather than formal ilr-CoDA.

This sits at the edge of the guide’s scope - read the disclaimers.

DST holds no diet (or physical-activity) data. These exposures come from a linked cohort or survey (e.g. a dietary questionnaire), joined to register outcomes and covariates in Phase 12. The code below runs on that assembled analysis dataset - so it opens an .rds/.dta file, not a register (readRDS, not read_register).
These are coding examples, not a methods lecture. For the theory, assumptions and interpretation, see Ibsen et al. (reference below).
The estimates are model-based. They assume you adjusted for the right total (here total energy), no residual confounding, linearity, and no substantial measurement error in the diet variables.

The idea, in one place

A substitution model estimates the effect of swapping one component for another while holding the total constant. Following Ibsen et al. 2021, there are two equivalent families:

Leave-one-out (nonspecified) model: put the food of interest in the model together with total energy and a composite total of the food group. Its coefficient is the effect of a higher intake of that food and a concomitantly lower intake of the omitted food in the group.
Partition (specified) model: put all component foods in the model plus total energy. The substitution of food A for food B is then the difference of their coefficients (b_A − b_B), whose variance is var(b_A) + var(b_B) − 2·cov(b_A, b_B).

Both routes target the same estimand - a swap holding the total fixed. The partition model makes the “for what?” explicit; done right, the two give the same answer.

Setup

library(dplyr) # mutate, %>%
# library(haven) # only if the linked cohort ships Stata files -> read_dta()

# The linked analysis dataset from Phase 12: diet (from the cohort/survey) plus
# outcome and covariates (from registers), one row per person.
df <- readRDS("path/to/analysis_dataset.rds")
# df <- haven::read_dta("path/to/analysis_dataset.dta")  # if it is a Stata .dta

Method 1 - partition (difference of coefficients)

Put every component food in the model, then read the substitution off as the difference between two coefficients. This is the R equivalent of Stata’s lincom.

library(dplyr)

# 1. Put the component foods on a common substitution unit (here per 100 g/day).
df <- df %>%
  mutate(
    redmeat100 = redmeat_g / 100, # red meat, per 100 g/day
    fish100 = fish_g / 100, # fish
    poultry100 = poultry_g / 100, # poultry
    mixedmeat100 = mixedmeat_g / 100 # mixed meat
  )

# 2. Partition model: ALL component foods + total energy in the same model.
#    Adjust for your confounders (age, sex, ...); logistic here, but for a
#    survival outcome swap glm() for survival::coxph() - the logic is identical.
partition <- glm(
  diabetes ~ redmeat100 +
    fish100 +
    poultry100 +
    mixedmeat100 +
    total_energy +
    age +
    sex, # + your confounder set
  family = binomial(link = "logit"),
  data = df
)

# 3. Substitution of FISH for RED MEAT = difference of their coefficients.
#    Var(fish - meat) = Var(fish) + Var(meat) - 2*Cov(fish, meat)
b <- coef(partition)
vc <- vcov(partition)

diff_beta <- b["fish100"] - b["redmeat100"] # log-OR of the swap
se_diff <- sqrt(
  vc["fish100", "fish100"] +
    vc["redmeat100", "redmeat100"] -
    2 * vc["fish100", "redmeat100"]
)

# 4. Odds ratio + 95% CI: replacing 100 g/day red meat with 100 g/day fish.
ci_logit <- diff_beta + c(-1, 1) * qnorm(0.975) * se_diff
cbind(OR = exp(diff_beta), LCL = exp(ci_logit[1]), UCL = exp(ci_logit[2]))

Method 2 - leave-one-out

Build a composite total of the food group and include every food except the one you substitute for. That omitted food (here red meat) becomes the implicit reference, so a food’s coefficient is its effect relative to red meat.

library(dplyr)

# 1. Per-100 g foods + a composite total of the whole group (all four foods).
df <- df %>%
  mutate(
    fish100 = fish_g / 100,
    poultry100 = poultry_g / 100,
    mixedmeat100 = mixedmeat_g / 100,
    totalmeat_g = redmeat_g + poultry_g + mixedmeat_g + fish_g, # group total
    meat100 = totalmeat_g / 100
  )

# 2. Leave-one-out model: composite total + every individual food EXCEPT the food
#    you substitute for (red meat is omitted). Total energy still adjusted for.
leaveoneout <- glm(
  diabetes ~ fish100 +
    poultry100 +
    mixedmeat100 +
    meat100 +
    total_energy +
    age +
    sex, # + your confounder set
  family = binomial(link = "logit"),
  data = df
)

# 3. fish100 = substituting 100 g/day fish for 100 g/day of the omitted food
#    (red meat) - the same swap as Method 1, so the estimates should agree.
OR <- exp(coef(leaveoneout)["fish100"])
CI <- exp(confint(leaveoneout)["fish100", ])
cbind(OR = OR, LCL = CI[1], UCL = CI[2])

Many swaps at once, and survival outcomes

Both extensions build on the partition model from Method 1: the substitution of A for B is always b_A - b_B, so one fitted model can answer many swaps.

A reusable helper turns any fitted model (glm or coxph) into a swap estimate with a 95% CI:

library(purrr) # map()
library(dplyr) # tibble(), bind_rows()

# Effect of swapping food `a` FOR food `b`, from a fitted partition model.
# Returns an odds ratio (glm) or a hazard ratio (coxph) with a 95% CI.
swap_effect <- function(model, a, b) {
  est <- coef(model)[a] - coef(model)[b] # difference of coefficients
  vc <- vcov(model)
  se <- sqrt(vc[a, a] + vc[b, b] - 2 * vc[a, b])
  tibble(
    swap = paste(a, "for", b),
    est = exp(est), # OR (glm) or HR (coxph)
    lcl = exp(est - qnorm(0.975) * se),
    ucl = exp(est + qnorm(0.975) * se)
  )
}

# Several foods swapped for red meat, from the ONE partition model (Method 1):
pairs <- list(
  c("fish100", "redmeat100"),
  c("poultry100", "redmeat100"),
  c("mixedmeat100", "redmeat100")
)
map(pairs, ~ swap_effect(partition, .x[1], .x[2])) %>% bind_rows()

Survival outcome (Cox). Fit the same partition model with coxph(), and the helper returns hazard ratios instead of odds ratios:

library(survival) # Surv(), coxph()

partition_cox <- coxph(
  Surv(followup_years, event) ~ redmeat100 +
    fish100 +
    poultry100 +
    mixedmeat100 +
    total_energy +
    age +
    sex, # + your confounder set
  data = df,
  ties = "breslow" # Breslow handling of tied event times
)

swap_effect(partition_cox, "fish100", "redmeat100") # HR: 100 g/day fish for red meat

Interpreting and refining

Substitution unit. The estimate is per unit swapped (here 100 g/day). Pick a unit that is meaningful and realistic, and always report both the unit and the pairing (“100 g/day red meat replaced by 100 g/day fish”).
What total-energy adjustment does. Holding total energy fixed is what turns “more fish” into “more fish instead of something else” rather than “more fish on top of the current diet”.
Non-linearity. A linear term assumes a constant swap effect across the whole range. If that is implausible, model the foods with splines (splines::ns(), see Regression) and evaluate the swap at chosen values.

The idea, in one place

Setup

Method 1 - partition (difference of coefficients)

Method 2 - leave-one-out

Many swaps at once, and survival outcomes

Interpreting and refining

See also