Have you ever been confused about data and evidence? Frankly, I find myself confused in data many times. The real goal isn’t to eliminate confusion entirely but to reduce it and gain some sense of certainty.
Can’t you see this is a land of confusion? This is the world we live in. —Land of confusion (song by Genesis)
One such thing that can be truly confusing in medical research is confounding.
But what is confounding? Confounding bias occurs when a third variable is associated with a risk factor (or treatment) and also influences the health outcome, leading to a wrong result. This can either create a false relationship or hide a real one.
I make up an example. Say, doctors have higher cancer rates. The confounding factor may be smoking if doctors tend to smoke more, and smoking causes cancer. So, no, being a doctor per se does not cause cancer at all.
After a literature search, I found two excellent examples from published medical studies where confounding occurred.
In the 1990s, observational studies suggested that hormone replacement therapy (HRT) reduced the risk of cardiovascular disease (CVD) in women [1]. However, this association was most likely false. Subsequent randomized controlled trials and a meta-analysis from 2002 found no evidence that HRT reduced the risk of CVD [2]. The apparent benefits were likely influenced by confounders such as socioeconomic status and education. Lower socioeconomic status is a strong risk factor for CVD, coronary artery disease, and many other adverse health outcomes. Additionally, women using HRT tend to have higher socioeconomic status, which may partially explain their better outcomes in earlier studies.
A 2016 observational study compared tracheal intubation to bag-valve-mask ventilation during pediatric cardiac arrest [3]. The study found that survival was lower in those children intubated compared with those not intubated. The authors concluded that the study does not support early tracheal intubation for pediatric in-hospital cardiac arrest. However, the study suffered from confounding by indication [4]. It was likely that children with a more severe condition and worse overall prognosis had a greater probability of being intubated. Consequently, the treatment group showed a higher risk of death.
The two examples show that confounding can confuse doctors in making the right choice based on solid evidence. Needless to say, this can potentially harm patients.
There are several ways to address confounding in medical research. One is study design, including randomized controlled trials (RCTs), which are highly effective and provide strong evidence.
Say, we are interested in the therapeutic effect of hypothermic machine perfusion (HMP) in kidney transplantation. We could perform an RCT and randomize the kidneys into a control group (standard cold storage) and an intervention group (HMP).
However, an RCT is not always possible. For example, if we want to assess the risk of donation after brain death (DBD) versus donation after circulatory death (DCD) kidney transplants, it is not possible to randomly assign kidneys to either group. Similarly, just as we cannot randomly assign individuals to be smokers or non-smokers, certain real-world conditions make randomization impossible.
There are also other approaches to tackle confounding, including matching and propensity score methods, which pair participants with similar clinical characteristics to reduce confounding [5,6]. It is important to carefully select the matching variables, choose an appropriate matching method, and perform a diagnostic assessment of the matching process. Additionally, stratified analysis helps control for confounding by analyzing data within specific subgroups.
In the remainder of this post, I will focus on another key approach: multivariable regression modeling. This method allows researchers to account for multiple confounders simultaneously by including them as covariates in a regression model, helping to isolate the true relationship between the exposure and outcome of interest.
A model is a lie that helps you see the truth. —Howard E. Skipper
I’m a big fan of model-based approaches [7,8]. Models provide the most comprehensive description of data. Unlike traditional statistical hypothesis tests, which focus primarily on detecting differences, models emphasize estimation, making results more meaningful for patients and clinicians.
One major advantage of models is their ability to adjust for multiple factors simultaneously, providing more accurate effect estimates. While researchers often focus on multiplicity adjustments for p-values, they rarely apply similar adjustments to effect estimates. Models, however, allow for these adjustments because the effect of interest is fully conditional on covariate values. In other words, model-based methods are much more powerful and reliable than analyzing variables in isolation.
For instance, consider a scenario in kidney transplantation, where we want to assess the impact of the deceased donor’s age increasing from 50 to 65 while holding other risk factors constant. A simple subgroup comparison—such as analyzing younger versus older donors doesn’t provide a reasonable estimate of this specific effect. A model, however, allows us to quantify the risk associated with this change while controlling for other variables as well.
Models also generalize statistical tests. For example, a t-test and ANOVA can be performed using multiple linear regression. Similarly, the Wilcoxon test and Kruskal-Wallis test correspond to a proportional odds model, while the log-rank test is a special case of the Cox proportional hazards model.
Ultimately, a model-based approach is the most accurate description of your data. However, it cannot definitively establish causality without careful attention to the research question, causal estimand, study design, and underlying assumptions [9].
In the end, you have the choice: get confused by analyzing variables individually or let a statistical model confuse you all at once.
The title of this post is a song by Genesis; the preview image of this post is a photo by Katie Moum on Unsplash.
I performed a simulation study with an example from the heart transplant waiting list. Again, this is not real data; I made this up. All data were generated from the following assumptions.
set.seed(2025)
# Install the swt package
# remotes::install_github("Swisstransplant/swt")
library(swt)
library(rms)
library(ggplot2)
N = 2000
n1 = 0.25*N # female
n2 = 0.75*N # male
heart = data.frame(weight = rep(NA, N),
sex = c(rep("female", n1), rep("male", n2))
)
heart$weight[heart$sex == "female"] = rnorm(n1, mean = 65, sd = 10)
heart$weight[heart$sex == "male"] = rnorm(n2, mean = 80, sd = 10)Next, I generated the probabilities for a heart transplant as an inverse function of body weight. As body weight increases, the likelihood of receiving a matching heart transplant decreases. Research supports this inverse relationship between body weight and transplant probability: larger individuals require greater cardiac output, and a smaller donor heart may struggle to meet these demands, potentially leading to poor circulation and organ failure [10].
I sampled from a binomial distribution to determine whether each patient received a transplant or not based on their individual probability. In my simulation, the probability of a transplant depended solely on body weight and not on sex. Precisely, the probability was constant at 0.80 until 65 kg; then it decreased as a quadratic function.
# function with constant prob until 65, then declines
prob_fun <- function(x) {
prob = ifelse(x <= 65, 0.8, 0.8 - (x - 65)^2 / 2200)
return(pmax(prob, 0))
}
heart$P = sapply(heart$weight, FUN = prob_fun)
heart$P = heart$P + rnorm(N, 0, sd = 0.1) # add some noise
heart$P[heart$P > 1] = 1 # truncate probabilities
heart$P[heart$P < 0] = 0
heart$event = rbinom(N, size = 1, prob = heart$P)
heart$event = heart$event == 1 # binary conversionThe simulated data for 10 randomly picked patients are shown in Table 1. Figure 1 and Table 2 display descriptives of the data.
set.seed(1986)
idx = sample(nrow(heart), 7)
heart[idx,]| weight | sex | P | event | |
|---|---|---|---|---|
| 840 | 72.26087 | male | 0.8648124 | TRUE |
| 1434 | 93.79369 | male | 0.3070354 | FALSE |
| 1083 | 67.71200 | male | 0.7036489 | TRUE |
| 1219 | 93.54775 | male | 0.5773796 | TRUE |
| 406 | 73.69477 | female | 0.7522511 | TRUE |
| 142 | 45.84125 | female | 0.8791807 | TRUE |
| 737 | 83.48071 | male | 0.6698299 | TRUE |
par(mfrow = c(1, 3))
hist(heart$weight, ylab = "Count", xlab = "Body weight (kg)", col = "#818C70", main=NULL,
cex.lab = 1.3, cex.axis = 1.3)
hist(heart$P, ylab = "Count", xlab = "Transplant probability", col = "#BF505A", main=NULL,
cex.lab = 1.3, cex.axis = 1.3)
#plot(heart$weight, sapply(heart$weight, FUN = prob_fun))
plot(heart$weight, heart$P,
xlab = "Body weight (kg)", ylab = "Transplant probability", col = "#6F87A6",
cex.lab = 1.3, cex.axis = 1.3)
tab = data.frame(N = nrow(heart),
weight = swt::mean_sd(heart$weight),
transplant = swt::freq_perc(heart$event)
)
heart.males = subset(heart, subset = sex == "male")
heart.females = subset(heart, subset = sex == "female")
tab.m = data.frame(N = nrow(heart.males),
weight = swt::mean_sd(heart.males$weight),
transplant = swt::freq_perc(heart.males$event)
)
tab.f = data.frame(N = nrow(heart.females),
weight = swt::mean_sd(heart.females$weight),
transplant = swt::freq_perc(heart.females$event)
)
tab = rbind(tab, tab.f, tab.m)
tab = as.data.frame(t(tab))
colnames(tab) = c("Overall", "Female", "Male")
rownames(tab) = c("Sample size N", "Body weight in kg, mean (SD)",
"Transplant, count (%)")
tab| Overall | Female | Male | |
|---|---|---|---|
| Sample size N | 2000 | 500 | 1500 |
| Body weight in kg, mean (SD) | 76.0 (12.0) | 65.0 (9.9) | 79.7 (10.2) |
| Transplant, count (%) | 1374 (68.7) | 384 (76.8) | 990 (66.0) |
As expected, females had an average weight of approximately 65 kg; males around 80 kg. Among females, about 77% received a transplant, compared to 66% of males.
I tested whether sex was independent of a heart transplant using Pearson’s \(\chi^2\)-test for count data. I knew this was a bad idea, but since looking at variables individually is common in medical research, doing it one more time didn’t hurt.
M = table(heart$event, heart$sex)
nontpx.fmt = sprintf("%d (%.1f)", M[1,], M[1,]/marginSums(M, margin = 2)*100)
tpx.fmt = sprintf("%d (%.1f)", M[2,], M[2,]/marginSums(M, margin = 2)*100)
tab.fmt = rbind(nontpx.fmt, tpx.fmt)
tab = as.data.frame(tab.fmt)
rownames(tab) = c("no transplant", "transplant")
colnames(tab) = c("female", "male")
tab| female | male | |
|---|---|---|
| no transplant | 116 (23.2) | 510 (34.0) |
| transplant | 384 (76.8) | 990 (66.0) |
stats = chisq.test(M, correct = FALSE)I found that more females received a transplant than expected (384 versus 343.5), while fewer males received a transplant than expected (990 versus 1030.5). This provided strong evidence that sex was related to receiving a heart transplant (p < 0.001).
But what was the issue? This simple analysis failed to capture the whole picture. Weight was related to both sex and the likelihood of receiving a heart transplant (Figure 2). Females typically have lower body weight, and heavier individuals are less likely to receive a heart transplant. Therefore, the body weight acts as a confounder. If I were unaware of body weight, I might incorrectly conclude that being female increases the likelihood of receiving a heart transplant. However, this was not the case and would only mislead us, as a heavy female would have a low chance of a transplant, while a lightweight male would have a higher chance.
flowchart LR A(Sex) --> B(Likelihood of a heart transplant) A <--> C(Body weight) C --> B
Next, I analyzed the data all at once using a logistic regression model. The results are shown in Table 4. Figure 3 shows the relationship between weight and the chance of a heart transplant.
dd = datadist(heart)
options(datadist='dd')
fit = lrm(event ~ sex + rcs(weight, 3), data = heart)
swt::tidy_rmsfit(fit)| Interquartile difference | Odds ratio (95%-CI) | Chi-Square | d.f. | p-value | |
|---|---|---|---|---|---|
| sex | – | – | 0.0138 | 1 | 0.91 |
| sex female | – | 1.02 (from 0.76 to 1.36) | – | – | – |
| weight | 16.12 (68.07–84.19) | 0.49 (from 0.42 to 0.57) | 119 | 2 | < 0.001 *** |
| weight nonlinear | – | – | 33.5 | 1 | < 0.001 *** |
| TOTAL | – | – | 139 | 3 | < 0.001 *** |
av = anova(fit)
ggplot(Predict(fit, weight), anova = av, colfill = "#6F87A6", adj.subtitle = FALSE) +
xlab("Body weight (kg)") +
theme_minimal() 
In conclusion, I found that body weight had a strong relationship with the chance of a heart transplant. Increasing weight from 68 kg to 84 kg was associated with roughly a 50% reduction in the chance of a heart transplant. No evidence was found for sex.
A model is often the best description of the data. Regression models are part of university-level statistics courses across life sciences. Thus, there is really no excuse for testing variables individually instead of using a statistical model.
@misc{schwab2025,
author = {Schwab, Simon},
title = {Land of Confusion},
date = {2025},
url = {https://www.statsyup.org/posts/confusion/},
langid = {en}
}