十四、回归的推断

十四、回归的推断

原文：Inference for Regression

译者：飞龙

协议：CC BY-NC-SA 4.0

自豪地采用谷歌翻译

到目前为止，我们对变量之间关系的分析纯粹是描述性的。我们知道如何找到穿过散点图的最佳直线来绘制。在所有直线中它的估计的均方误差最小，从这个角度来看，这条线是最好的。

但是，如果我们的数据是更大总体的样本呢？如果我们在样本中发现了两个变量之间的线性关系，那么对于总体也是如此嘛？它会是完全一样的线性关系吗？我们可以预测一个不在我们样本中的新的个体的响应变量吗？

如果我们认为，散点图反映了被绘制的两个变量之间的基本关系，但是并没有完全规定这种关系，那么就会出现这样的推理和预测问题。例如，出生体重与孕期的散点图，显示了我们样本中两个变量之间的精确关系；但是我们可能想知道，对于抽样总体中的所有新生儿或实际中的一般新生儿，这样的关系是否是真实的，或者说几乎是正确的。

一如既往，推断思维起始于仔细检查数据的假设。一组假设被称为模型。大致线性的散点图中的一组随机性的假设称为回归模型。

回归模型

简而言之，这样的模型认为，两个变量之间的底层关系是完全线性的；这条直线是我们想要识别的信号。但是，我们无法清楚地看到这条线。我们看到的是分散在这条线上的点。在每一点上，信号都被随机噪声污染。因此，我们的推断目标是将信号从噪声中分离出来。

更详细地说，回归模型规定了，散点图中的点是随机生成的，如下所示。

x和y之间的关系是完全线性的。我们看不到这个“真实直线”，但它是存在的。
散点图通过将线上的点垂直移动，或上或下来创建，如下所示：
对于每个x，找到真实直线上的相应点（即信号），然后生成噪声或误差。
误差从误差总体中带放回随机抽取，总体是均值为 0 的正态分布。
创建一个点，横坐标为x，纵坐标为“x处的真实高度加上误差”。
最后，从散点图中删除真正的线，只显示创建的点。

基于这个散点图，我们应该如何估计真实直线？我们可以使其穿过散点图的最佳直线是回归线。所以回归线是真实直线的自然估计。

下面的模拟显示了回归直线与真实直线的距离。第一个面板显示如何从真实直线生成散点图。第二个显示我们看到的散点图。第三个显示穿过散点图的回归线。第四个显示回归线和真实直线。

为了运行模拟，请使用三个参数调用draw_and_compare函数：真实直线的斜率，真实直线的截距以及样本量。

运行模拟几次，用不同的斜率和截距，以及不同的样本量。因为所有的点都是根据模型生成的，所以如果样本量适中，你会看到回归线是真实直线的一个良好估计。

# The true line,
# the points created,
# and our estimate of the true line.
# Arguments: true slope, true intercept, number of points
draw_and_compare(4, -5, 10)

实际上，我们当然不会看到真实直线。模拟结果表明，如果回归模型看起来合理，并且如果我们拥有大型样本，那么回归线就是真实直线的一个良好近似。

真实斜率的推断

我们的模拟表明，如果回归模型成立，并且样本量很大，则回归线很可能接近真实直线。这使我们能够估计真实直线的斜率。

我们将使用我们熟悉的母亲和她们的新生儿的样本，来开发估计真实直线的斜率的方法。首先，我们来看看我们是否相信，回归模型是一系列适当假设，用于描述出生体重和孕期之间的关系。

correlation(baby, 'Gestational Days', 'Birth Weight')
0.40754279338885108

总的来说，散点图看起来相当均匀地分布在这条线上，尽管一些点分布在主云形的周围。相关系数为 0.4，回归线斜率为正。

这是否反映真实直线斜率为正的事实？为了回答这个问题，让我们看看我们能否估计真实斜率。我们当然有了一个估计：我们的回归线斜率。这大约是 0.47 盎司每天。

slope(baby, 'Gestational Days', 'Birth Weight')
0.46655687694921522

但是如果散点图出现的方式不同，回归线会有所不同，可能会有不同的斜率。我们如何计算，斜率可能有多么不同？

我们需要点的另一个样本，以便我们可以绘制回归线穿过新的散点图，并找出其斜率。但另一个样本从哪里得到呢？

你猜对了 - 我们将自举我们的原始样本。这会给我们自举的散点图，通过它我们可以绘制回归线。

自举散点图

我们可以通过对原始样本带放回地随机抽样，来模拟新样本，它的次数与原始样本量相同。这些新样本中的每一个都会给我们一个散点图。我们将这个称为自举散点图，简而言之，我们将调用整个过程来自举散点图。

这里是来自样本的原始散点图，以及自举重采样过程的四个复制品。请注意，重采样散点图通常比原始图稀疏一点。这是因为一些原始的点没有在样本中被选中。

估计真实斜率

我们可以多次自举散点图，并绘制穿过每个自举图的回归线。每条线都有一个斜率。我们可以简单收集所有的斜率并绘制经验直方图。回想一下，在默认情况下，sample方法带放回地随机抽取，次数与表中的行数相同。也就是说，sample默认生成一个自举样本。

slopes = make_array()
for i in np.arange(5000):
    bootstrap_sample = baby.sample()
    bootstrap_slope = slope(bootstrap_sample, 'Gestational Days', 'Birth Weight')
    slopes = np.append(slopes, bootstrap_slope)
Table().with_column('Bootstrap Slopes', slopes).hist(bins=20)

然后，我们可以使用percentile方法，为真实直线的斜率构建约 95% 置信区间。置信区间从 5000 个自举斜率的第 2.5 百分位数，延伸到第 97.5 百分位数。

left = percentile(2.5, slopes)
right = percentile(97.5, slopes)
left, right
(0.38209399211893086, 0.56014757838023777)

用于自举斜率的函数

让我们收集我们估计斜率的方法的所有步骤，并定义函数bootstrap_slope来执行它们。它的参数是表的名称，预测变量和响应变量的标签，以及自举复制品的所需数量。在每个复制品中，该函数自举原始散点图并计算所得回归线的斜率。然后绘制所有生成的斜率的直方图，并打印由斜率的“中间 95%”组成的区间。

def bootstrap_slope(table, x, y, repetitions):
    # For each repetition:
    # Bootstrap the scatter, get the slope of the regression line,
    # augment the list of generated slopes
    slopes = make_array()
    for i in np.arange(repetitions):
        bootstrap_sample = table.sample()
        bootstrap_slope = slope(bootstrap_sample, x, y)
        slopes = np.append(slopes, bootstrap_slope)
    # Find the endpoints of the 95% confidence interval for the true slope
    left = percentile(2.5, slopes)
    right = percentile(97.5, slopes)
    # Slope of the regression line from the original sample
    observed_slope = slope(table, x, y)
    # Display results
    Table().with_column('Bootstrap Slopes', slopes).hist(bins=20)
    plots.plot(make_array(left, right), make_array(0, 0), color='yellow', lw=8);
    print('Slope of regression line:', observed_slope)
    print('Approximate 95%-confidence interval for the true slope:')
    print(left, right)

当响应变量为出生体重，预测变量为孕期时，我们调用bootstrap_slope来找到真实斜率的置信区间，我们得到了一个区间，非常接近我们之前获得的东西：大约 0.38 到 0.56 盎司每天。

bootstrap_slope(baby, 'Gestational Days', 'Birth Weight', 5000)
Slope of regression line: 0.466556876949
Approximate 95%-confidence interval for the true slope:
0.378663152966 0.555005146304

现在我们有一个函数，可以自动完成估计在回归模型中展示斜率的过程，我们也可以在其他变量上使用它。

例如，我们来看看出生体重与母亲身高的关系。更高的女性往往有更重的婴儿吗？

回归模型似乎是合理的，基于散点图，但相关性不高。这只有大约 0.2。

scatter_fit(baby, 'Maternal Height', 'Birth Weight')

correlation(baby, 'Maternal Height', 'Birth Weight')
0.20370417718968034

像之前一样，我们使用bootstrap_slope来估计回归模型中真实直线的斜率。

bootstrap_slope(baby, 'Maternal Height', 'Birth Weight', 5000)
Slope of regression line: 1.47801935193
Approximate 95%-confidence interval for the true slope:
1.0403083964 1.91576886223

真实斜率的 95% 置信区间，从约 1 延伸到约 1.9 盎司每英寸。

真实斜率可能为 0 嘛？

假设我们相信我们的数据遵循回归模型，并且我们拟合回归线来估计真实直线。如果回归线不完全是平的，几乎总是如此，我们将观察到散点图中的一些线性关联。

但是，如果这种观察是假的呢？换句话说，如果真实直线是平的 - 也就是说，这两个变量之间没有线性关系 - 我们观察到的联系，只是由于从样本中产生点的随机性。

这是一个模拟，说明了为什么会出现这个问题。我们将再次调用draw_and_compare函数，这次要求真实斜率为 0。我们的目标是，观察我们的回归线是否显示不为 0 的斜率。

请记住函数draw_and_compare的参数是真实直线的斜率和截距，以及要生成的点的数量。

draw_and_compare(0, 10, 25)

运行模拟几次，每次保持真实直线的斜率为 0 。你会注意到，虽然真实直线的斜率为 0，但回归线的斜率通常不为 0。回归线有时会向上倾斜，有时会向下倾斜，每次都给我们错误的印象，即这两个变量是相关的。

为了确定我们所看到的斜率是否真实，我们想测试以下假设：

原假设。真实直线的斜率是 0。

备选假设。真实直线的斜率不是 0。

我们很有条件来实现它。由于我们可以为真实斜率构建一个 95% 的置信区间，我们所要做的就是看区间是否包含 0。

如果没有，那么我们可以拒绝原假设（P 值为 5% 的截断值）。

如果真实斜率的置信区间确实包含 0，那么我们没有足够的证据来拒绝原假设。也许我们看到的斜率是假的。

我们在一个例子中使用这个方法。假设我们试图根据母亲的年龄来估计新生儿的出生体重。根据样本，根据母亲年龄估计出生体重的回归线的斜率为正，约为 0.08 盎司每年。

slope(baby, 'Maternal Age', 'Birth Weight')
0.085007669415825132

虽然斜率为正，但是很小。回归线非常接近平的，这就产生了一个问题，真实直线是否是平的。

scatter_fit(baby, 'Maternal Age', 'Birth Weight')

我们可以使用bootstrap_slope来估计真实直线的斜率。计算表明，真实斜率的约 95% 的自举置信区间左端为负，右端为正 - 换句话说，区间包含 0。

bootstrap_slope(baby, 'Maternal Age', 'Birth Weight', 5000)
Slope of regression line: 0.0850076694158
Approximate 95%-confidence interval for the true slope:
-0.104335243815 0.272791852339

因为区间包含 0，所以我们不能拒绝原假设，母亲年龄与新生儿出生体重之间的真实线性关系的斜率为 0。基于此分析，使用母亲年龄作为预测变量，基于回归模型预测出生体重是不明智的。

预测区间

回归的主要用途之一是对新个体进行预测，这个个体不是我们原始样本的一部分，但是与样本个体相似。在模型的语言中，我们想要估计新值x的y。

我们的估计是真实直线在x处的高度。当然，我们不知道真实直线。我们使用我们的样本点的回归线来代替。

给定值x的拟合值，是基于x值的y的回归估计。换句话说，给定值x的拟合值就是回归线在x处的高度。

假设我们试图根据孕期天数来预测新生儿的出生体重。我们在前面的章节中看到，这些数据非常适合回归模型，真实直线的斜率的 95% 置信区间不包含 0。因此，我们的预测似乎是合理的。

下图显示了预测位于回归线上的位置。红线是x = 300。

红线与回归线的相交点的高度是孕期天数 300 的拟合值。

函数fitted_value计算这个高度。像函数的相关性，斜率和截距一样，它的参数是表的名称和x和y的列标签。但是它也需要第四个参数，即x的值，在这个值上进行估算。

def fitted_value(table, x, y, given_x):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * given_x  + b

孕期天数 300 的拟合值约为 129.2 盎司。换句话说，对于孕期为 300 天的孕妇，我们估计的新生儿体重约为 129.2 盎司。

fit_300 = fitted_value(baby, 'Gestational Days', 'Birth Weight', 300)
fit_300
129.2129241703143

预测的可变性

我们已经开发了一种方法，使用我们样本中的数据，根据孕期天数预测新生儿的体重。但作为数据科学家，我们知道样本可能有所不同。如果样本不同，回归线也会不一样，我们的预测也是。为了看看我们的预测有多好，我们必须了解预测的可变性。

为此，我们必须生成新的样本。我们可以像上一节那样，通过自举散点图来实现。然后，我们为每个散点图的复制品拟合回归线，并根据每一行进行预测。下图显示了 10 条这样的线，以及孕期天数 300 对应的出生体重预测。

lines

slope	intercept	prediction at x=300
0.503931	-21.6998	129.479
0.53227	-29.5647	130.116
0.518771	-25.363	130.268
0.430556	-1.06812	128.099
0.470229	-11.7611	129.308
0.48713	-16.5314	129.608
0.51241	-23.2954	130.428
0.52473	-27.2053	130.214
0.409943	5.22652	128.21
0.468065	-11.6967	128.723

每一行的预测都不相同。下表显示了 10 条线的斜率、截距以及预测。

自举预测区间

如果我们增加重采样过程的重复次数，我们可以生成预测的经验直方图。这将允许我们创建预测区间，使用为斜率创建自举置信区间时的相同的百分比方法。

让我们定义一个名为bootstrap_prediction的函数来实现。该函数有五个参数：

表的名称
预测变量和响应变量的列标签
用于预测的x的值
所需的自举重复次数

在每次重复中，函数将自举原始散点图，并基于x的指定值查找y的预测值。具体来说，它调用我们在本节前面定义的函数fitted_value，来寻找指定x处的拟合值。

最后，绘制所有预测值的经验直方图，并打印由预测值的“中间 95%”组成的区间。它还打印基于穿过原始散点图的回归线的预测值。

# Bootstrap prediction of variable y at new_x
# Data contained in table; prediction by regression of y based on x
# repetitions = number of bootstrap replications of the original scatter plot
def bootstrap_prediction(table, x, y, new_x, repetitions):
    # For each repetition:
    # Bootstrap the scatter; 
    # get the regression prediction at new_x; 
    # augment the predictions list
    predictions = make_array()
    for i in np.arange(repetitions):
        bootstrap_sample = table.sample()
        bootstrap_prediction = fitted_value(bootstrap_sample, x, y, new_x)
        predictions = np.append(predictions, bootstrap_prediction)
    # Find the ends of the approximate 95% prediction interval
    left = percentile(2.5, predictions)
    right = percentile(97.5, predictions)
    # Prediction based on original sample
    original = fitted_value(table, x, y, new_x)
    # Display results
    Table().with_column('Prediction', predictions).hist(bins=20)
    plots.xlabel('predictions at x='+str(new_x))
    plots.plot(make_array(left, right), make_array(0, 0), color='yellow', lw=8);
    print('Height of regression line at x='+str(new_x)+':', original)
    print('Approximate 95%-confidence interval:')
    print(left, right)
bootstrap_prediction(baby, 'Gestational Days', 'Birth Weight', 300, 5000)
Height of regression line at x=300: 129.21292417
Approximate 95%-confidence interval:
127.300774171 131.361729528

上图显示了基于 5000 次重复的自举过程，孕期天数 300 的预测出生体重的自举经验直方图。经验分布大致是正泰的。

我们已经通过预测的“中间 95%”，即预测的第 2.5 百分位数到第 97.5 百分位数的区间，构建了分数的约 95% 的预测区间。区间范围从大约 127 到大约 131。基于原始样本的预测是大约 129，接近区间的中心。

改变预测变量的值的效果

下图显示了孕期天数 285 的 5,000 次自举预测的直方图。基于原始样本的预测是约 122 盎司，区间范围从约 121 盎司到约 123 盎司。

bootstrap_prediction(baby, 'Gestational Days', 'Birth Weight', 285, 5000)
Height of regression line at x=285: 122.214571016
Approximate 95%-confidence interval:
121.177089926 123.291373304

请注意，这个区间比孕妇天数 300 的预测区间更窄。让我们来调查其原因。

孕妇天数均值约为 279 天：

np.mean(baby.column('Gestational Days'))
279.10136286201021

所以 285 比 300 更接近分布的中心。通常，基于自举样本的回归线，在预测变量的分布中心附近彼此更接近。因此，所有的预测值也更接近。这解释了预测区间的宽度更窄。

你可以在下面的图中看到这一点，它显示了 10 个自举复制品中每一个的x = 285和x = 300的预测值。通常情况下，直线在x = 300处比x = 285处相距更远，因此x = 300的预测更加可变。

注意事项

我们在本章中进行的所有预测和测试，都假设回归模型是成立的。具体来说，这些方法假设，散点图中的点由直线上的点产生，然后通过添加随机正态噪声将它们推离直线。

如果散点图看起来不像那样，那么模型可能不适用于数据。如果模型不成立，那么假设模型为真的计算是无效的。

因此，在开始基于模型进行预测，或者对模型参数进行假设检验之前，我们首先要确定回归模型是否适用于我们的数据。一个简单的方法就是，按照我们在本节所做的操作，即绘制两个变量的散点图，看看它看起来是否大致线性，并均匀分布在一条线上。我们还应该使用残差图，执行我们在前一节中开发的诊断。