混淆是关于偏见的。我能想到的至少有两种方法可以使二分法使回归产生偏差：

• 模型偏差。对连续预测变量进行二分法将条件均值从连续（可能是平滑的，取决于函数形式）函数转换为不连续的常数函数。如果我们真的相信这种现象会表现出这种行为，最好问问自己（我倾向于认为医学中的大多数事情都不是这样）。因此，二分法使我们的估计偏向于不连续且恒定的函数集。这很糟糕，我们将看到。

• 关联估计中的偏差。回归分析，只要它来自观察数据而不是某些因果框架的一部分，都是关于关联的。模型的系数$\hat{\beta}_j$是观察到的关联$x_j$和结果。估计的关联可能与期望的真实关联不同。这是您可能在许多高年级本科生或研究生水平的文本中看到的那种标准统计偏差

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

如果该期望为 0 ，则该估计称为无偏估计。

在接下来的内容中，我将在一个非常简单的模型中证明二分法会导致两种偏差。

---

让我们从以下数据生成过程中模拟一些数据

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

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

让我们对待$x$就好像它是连续的，而实际上我会将它放入长度为 0.1 的箱中以作为展示的手段。通过 OLS 进行召回线性模型拟合具有一些关键属性，一个是残差的期望应该为 0 并且独立于预测变量。

让我们生成数据，拟合模型，然后得到 1000 次残差。让我们取每个残差的期望$x$价值。每个人的期望应该是 0$x$价值。我们可以用一些 R 代码来证明这一点（我将在答案末尾包含这些代码）。

我已经继续并包括了大约 95% 的置信区间。并非每个剩余 CI 都覆盖 0，但这是意料之中的（事实上，5% 应该没有这样做），但这在很大程度上符合我们的预期。让我们再做一次，但现在我们将二分法$x$在样本中位数（也许是最慈善的分裂点，这样两组的观察次数大致相同）。

二分法的效果是残差的期望不再是 0，而是依赖于协变量。据我了解，这是残留的混淆。特别是，这是模型偏差的结果。我们选择的函数类别（不连续常数函数）不能正确容纳我们的数据，导致预测变量和残差之间存在关系。这是为了查看预测变量中是否存在任何非线性关系而可能执行的相同类型的诊断（不同之处在于我已经汇总了 1000 次模拟的残差，而实际上通常有 1 个数据集）。毕竟，预测是一个估计，而且这个模型中的许多估计似乎是有偏差的（即使在由二分法确定的组内）。

不仅残差与建模假设不一致，观察到的关联的估计也有偏差。如果满足 OLS 的假设，则估计$\beta$（斜率，或者模型中的任何参数）应该是无偏的。这意味着，在重复模拟的情况下，估计值与真实值之间的差异的期望值应该为 0。我们可以证明，当数据被二分法时，情况并非如此。事实上，对于这个问题，估计的关系大约要高 4 个单位（在考虑问题的结构后才有意义）。

当您在同一轴上绘制估计值和真实值之间的差异时，差异会更加显着。

请注意，可以通过正确调整属于每个组的点的平均值之间的距离来修复偏差（如下所示），但这并不能改善模型偏差（模型在给定特定值的情况下低估/高估的系统趋势的预测器）。此外，如果您要以这种方式更正二分法估计，您不妨正确估计斜率。这只是通过错误的方式得到正确的答案，并称整个方法有效。

---

我们可以从这个答案中得出什么结论？有点开玩笑，我们可能会得出二分法不好的结论。更真诚地，二分法不必要地使我们可能的条件期望类别偏向一类功能，至少在我看来，这在实践中很少被观察到。确实，如果我可以稍微长一点，二分法与其说是统计效率，不如说是认知经济。人们之所以二分法，是因为考虑 2 个数字比考虑它们的连续统一体更容易。

其次，观察到的关系是有偏见的（在这种情况下，远离零，使关系看起来比实际更极端，但我可以很容易地创建一个偏向零的例子）。这是一种特别有害的偏差形式，因为与模型偏差不同，如果不对模型进行批判性分析，它可能在很大程度上未被发现。如果回归分析用于某种基于证据的调查，我会让你想知道这种偏见如何影响护理。

我应该补充一点，模型的简单性不是限制。这些现象可以出现在多重回归中，并以不同的方式出现在 GLM（如逻辑回归）中。

# 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)