Visual Demonstration of Residual Confounding

Question

To drive home the point that dichotomizing a continuous confounder creates problems of residual confounding in regression analyses, I'm wondering if there are any available visual demonstrations (not DAGS) for the average biomedical researcher (amongst whom I include myself).

Answer 1

Confounding is about bias. There are at least two ways I can think of in which dichotomization can bias a regression:

• Model Bias. Dichotomizing a continuous predictor transforms the conditional mean from a continuous (possibly smooth, depending on the functional form) function to a discontinuous constant function. It would be good to ask one's self if we truly believe the phenomenon to display this sort of behaviour (I tend to think most things in medicine do not behave like this). Hence, dichotomizing biases our estimates towards the set of functions which are discontinuous and constant. That's bad, as we will see.

• Bias in Estimates of Association. Regression Analysis, so far as it is from observational data and not part of some causal framework, is all about associations. The coefficients of the model $\hat{\beta}_j$ are the observed association between $x_j$ and the outcome. The estimated association can differ from the true association in expectation. This is the sort of standard statistical bias you might see in many senior undergrad or graduate level texts

$$ \operatorname{Bias} = E(\hat{\beta_j} - \beta_j) $$

with the estimate being called unbiased if this expectation is 0.

In what follows, I demonstrate that dichomotization leads to both kinds of bias in a very simple model.

---

Let's simulate some data from the following data generating process

$$ y \vert x \sim \mathcal{N}(\mu(x), 1) $$

$$ \mu(x) = 2x + 1$$

Let's treat $x$ as if it were continuous, when in reality I will bucket it into bins of length 0.1 for means of exposition. Recall linear model fit via OLS has a few key properties, one being that the expectation of the residual should be 0 and independent of the predictor.

Let's generate data, fit a model, and get the residuals say 1000 times. Let's take the expectation of the residual across each $x$ value. The expectation should be 0 for each $x$ value. We can demonstrate that with some R code (which I will include at the end of the answer).

I've gone ahead and included approximate 95% confidence intervals as well. Not every residual CI covers 0, but this is to be expected (in fact, 5% should fail to do so), but this is largely in line with what we would expect. Let's do this again, but now we will dichotomize $x$ at the sample median (perhaps the most charitable point to split at so that both groups have approximately the same number of observations).

The effect of dichotomization is that the expectation for the residual is no longer 0 and is dependent on the covariate. This is residual confounding, as I understand it. In particular, this is the result of model bias. Our chosen class of functions (discontinuous constant functions) can not properly accommodate our data, leading to a relationship between the predictor and the residuals. This is the same sort of diagnostic one might perform to see if there are any non-linear relationships in the predictor (the difference being here that I have aggregated residuals over 1000 simulations, and in practice one usually has 1 dataset). After all, a prediction is an estimate, and it seems that many estimates in this model are biased (even within groups determined by the dichotomization).

Not only are the residuals not consistent with the modelling assumptions, the estimate of the observed association is also biased. If the assumptions of OLS are met, then the estimate of $\beta$ (the slope, or really any parameter in the model) should be unbiased. This means that, under repeated simulation, the expectation of the difference between estimate and truth should be 0. We can demonstrate that this is not the case when the data are dichotomized. In fact, for this problem, the estimated relationship is about 4 units higher (which makes sense after you think about the structure of the problem).

The difference is even more striking when you plot the difference between estimated and truth on the same axis.

Note that the bias can be fixed by correctly adjusting for the distance between mean of points which belong to each group (shown below) but this does not improve the model bias (the systematic tendency for the model to under/over estimate given a particular value of the predictor). Additionally, if you're going to correct the dichotmomized estimate in this way, you might as well just estimate the slope correctly. This is just getting the right answer via the wrong way and calling the entire approach valid.

---

What can we conclude from this answer? Somewhat tongue in cheek, we might conclude dichotmization bad. More sincerely, dichotmization needlessly biases our class of possible conditional expectations towards a class of functions which, at least in my opinion, is rarely observed in practice. Indeed, if I could opine slighly longer, dichotmization is not so much about statistical efficiency more than it is cognitive economy; people dichotmize because it is easier to think about 2 numbers than a continuum of them.

Second, the observed relationship is biased (in this case away from the null, making the relationship seem more extreme than it actually is, but I could just as easily create an example where the bias is towards the null). This is a particularly pernicious form of bias, because it can go largely undetected without critical analysis of the model, unlike model bias. I will leave you to wonder about how this bias can effect care, should the regression analysis be used in some sort of evidence based investigation.

I should add that the simplicity of the model is not a limitation. These sorts of phenomenon can appear in multiple regressions, and in different ways, in GLMs (like logistic regression).

# Predictor buckted at bins of width 0.1
x = seq(-3.0, 3.0, 0.1)
N = length(x)

# Simulate the data 1000 times and return a residual for each predictor.
r = replicate(1000,{
  y = 2*x + 1 + rnorm(N, 0, 1)
  fit = lm(y~x)
  resid(fit)
  })


e = rowMeans(r)
s = apply(r, 1, function(x) sd(x)/sqrt(length(x)))
plot(x, e, pch=20, ylim = c(-0.15, 0.15), main='Expectation of Residuals Under No Dichotimization', ylab=expression(y-hat(y)))
arrows(x0=x, y0=e-2*s, x1=x, y1=e+2*s, code=3, angle=90, length=0, col="black", lwd=2)
abline(h=0, col='dark grey')

r2 = replicate(1000,{
  y = 2*x + 1 + rnorm(N, 0, 1)
  z = as.numeric(x>0)
  fit = lm(y~z)
  resid(fit)
})


e = rowMeans(r2)
s = apply(r2, 1, function(x) sd(x)/sqrt(length(x)))
plot(x, e, pch=20,, main='Expectation of Residuals Under Dichotomization', ylab=expression(y-hat(y)))
arrows(x0=x, y0=e-2*s, x1=x, y1=e+2*s, code=3, angle=90, length=0, col="black", lwd=2)
abline(h=0, col='dark grey')




no_dichot = replicate(1000, {
  x = seq(-3.0, 3.0, 0.1)
  N = length(x)
  true_slope=2
  y =true_slope*x + 1 + rnorm(N, 0, 1)
  fit = lm(y~x)
  estimated_slope = coef(fit)[2]
  estimated_slope -true_slope
  
})

dichot = replicate(1000, {
  x = seq(-3.0, 3.0, 0.1)
  N = length(x)
  true_slope=2
  z = as.numeric(x>0)
  y =true_slope*x + 1 + rnorm(N, 0, 1)
  fit = lm(y~z)
  estimated_slope = coef(fit)[2]
  estimated_slope -true_slope
  
})


hist(no_dichot,
     col=rgb(1, 0, 0, 0.5),
     breaks = seq(-2, 6, 0.05),
     xlab=expression(hat(beta)-beta),
     main='Difference between estimated relationship and true relationship (1000 simulations)',
     xlim = c(-2, 6))

hist(dichot,
     col=rgb(0, 0, 1, 0.5),
     breaks = seq(-2, 6, 0.05),
     add=T)

legend("topright", c("Continuous", "Dichotomize"), col=c(rgb(1, 0, 0, 0.5), rgb(0, 0, 1, 0.5)), lwd=10)