tnkstat3e-merge-2
统计思维（程序员的概率统计）第三版（三）
原文：
allendowney.github.io/ThinkStats/index.html译者：飞龙
协议：CC BY-NC-SA 4.0
估计
假设你居住在一个有 10,000 人的小镇上，你想预测即将到来的选举中谁会获胜。从理论上讲，你可以询问镇上的每个人他们打算投给谁，如果他们都诚实地回答，你就可以做出可靠的预测。
但即使在小镇上，调查整个种群可能也不太实际。幸运的是，这并不必要。如果你调查了人群的随机子集，你可以使用样本来推断人群的投票偏好。这个过程——使用样本对人群进行推断——被称为统计推断。
统计推断包括估计，这是本章的主题，以及假设检验，这是下一章的主题。
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("Downloaded " + local)
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")
try:
import empiricaldist
except ImportError:
%pip install empiricaldist
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from thinkstats import decorate
``` ## 企鹅称重
假设你是一名在南极洲的研究员，正在研究当地的企鹅种群。你的一个任务是监控企鹅在整个年度中的平均体重变化。由于不可能对环境中的每只企鹅进行称重，因此你的计划是每周随机抽取 10 只企鹅进行称重，并使用这些样本来估计整个种群的平均值——这被称为**种群均值**。
你可以使用许多方法来使用样本估计种群均值，但我们将考虑其中的两种：样本均值和样本中位数。它们都是合理的选择，但让我们看看哪个更好——并思考一下我们所说的“更好”是什么意思。
为了演示目的，我们假设企鹅的体重服从一个已知均值和标准差的正态分布，我将用`mu`和`sigma`表示，并赋予千克单位。
```py
mu = 3.7
sigma = 0.46
这些值是正态分布的参数，这意味着它们指定了特定的分布。给定这些参数，我们可以使用 NumPy 来模拟抽样过程并生成任何大小的样本。例如，这里是一个假设的 10 个重量的样本。
# Seed the random number generator so we get the same results every time
np.random.seed(1)
sample = np.random.normal(mu, sigma, size=10)
sample
array([4.44719887, 3.41859205, 3.45704099, 3.20643443, 4.09808751,
2.6412922 , 4.50261341, 3.34984483, 3.84675798, 3.58528963])
这里是样本的平均值和中位数。
np.mean(sample), np.median(sample)
(np.float64(3.6553151902291945), np.float64(3.521165310619601))
平均值和中位数差异足够大，以至于我们应该思考哪个是更好的估计。为了找出答案，我们将使用以下函数生成给定大小n的假设样本。
def make_sample(n):
return np.random.normal(mu, sigma, size=n)
作为第一次实验，让我们看看样本均值和样本中位数随着样本大小的增加是如何表现的。我们将使用 NumPy 函数logspace来创建一个从 10 到 100,000 的对数尺度上的ns范围。
ns = np.logspace(1, 5).astype(int)
我们可以使用列表推导来生成每个n值的假设样本，计算平均值，并收集结果：
means = [np.mean(make_sample(n)) for n in ns]
我们也会对中位数做同样的处理。
medians = [np.median(make_sample(n)) for n in ns]
用来估计总体属性的统计量，如样本均值或中位数，被称为估计量。
以下图示展示了随着样本量的增加，这些估计量是如何表现的。水平线显示了总体中的实际均值。
plt.axhline(mu, color="gray", lw=1, alpha=0.5)
plt.plot(ns, means, "--", label="mean")
plt.plot(ns, medians, alpha=0.5, label="median")
decorate(xlabel="Sample size", xscale="log", ylabel="Estimate")
对于这两个估计量，随着样本量的增加，估计值会收敛到实际值。这证明了它们是一致的，这是良好估计量应具备的一个属性。基于这个属性，均值和中位数看起来同样好。
在前面的图中，你可能会注意到估计值有时过高，有时过低——看起来变化大致对称于真实值。这表明另一个实验：如果我们收集许多相同大小的样本并计算许多估计值，估计值的平均值是多少？
下面的循环通过生成 10,001 个 10 只企鹅的样本并计算每个样本的均值来模拟这种场景。
means = [np.mean(make_sample(n=10)) for i in range(10001)]
np.mean(means)
np.float64(3.70034508492869)
均值平均值接近我们用来生成样本的实际均值：3.7 千克。
下面的循环模拟了相同的场景，但这次它计算每个样本的中位数。
medians = [np.median(make_sample(n=10)) for i in range(10001)]
np.mean(medians)
np.float64(3.701214089907223)
这些假设的中位数平均值也非常接近实际总体均值。
这些结果表明，样本均值和中位数是无偏估计量，这意味着它们在平均上是正确的。在不同的语境中，“偏差”这个词有不同的含义，这可能会引起混淆。在这个语境中，“无偏”意味着估计的平均值是实际值。
到目前为止，我们已经表明这两个估计量都是一致且无偏的，但仍然不清楚哪个更好。让我们再做一个实验：让我们看看哪个估计量更准确。在不同的语境中，“准确”这个词也有不同的含义——作为一种量化方法，让我们考虑均方误差（MSE）。以下函数计算估计值与实际值之间的差异，并返回这些误差平方的平均值。
def mse(estimates, actual):
"""Mean squared error of a sequence of estimates."""
errors = np.asarray(estimates) - actual
return np.mean(errors**2)
注意，我们只能计算均方误差如果我们知道实际值。在实践中，我们通常不知道——毕竟，如果我们知道实际值，我们就不需要估计它了。但在我们的实验中，我们知道实际总体均值是 3.7 千克，因此我们可以用它来计算样本均值的均方误差。
mse(means, mu)
np.float64(0.020871984891289382)
如果我们有大小为 10 的样本，并使用样本均值来估计总体均值，平均平方误差大约为 0.021 千克平方。现在来看样本中位数的均方误差。
mse(medians, mu)
np.float64(0.029022273128644173)
如果我们使用样本中位数来估计总体均值，平均平方误差大约是 0.029 千克平方。在这个例子中，样本均值比样本中位数更好；而且一般来说，如果数据来自正态分布，它就是总体均值的最优无偏估计量，因为它最小化了均方误差（MSE）。
最小化 MSE 是一个估计量应该具备的良好属性，但 MSE 并不总是总结误差的最佳方式。一方面，它很难解释。在这个例子中，MSE 的单位是千克平方，因此很难说这意味着什么。
一个解决方案是使用均方误差的平方根，称为“均方根误差”（RMSE）。另一个选项是使用误差绝对值的平均值，称为“平均绝对误差”（MAE）。以下函数计算估计序列的 MAE。
def mae(estimates, actual):
"""Mean absolute error of a sequence of estimates."""
errors = np.asarray(estimates) - actual
return np.mean(np.abs(errors))
这是样本均值的 MAE。
mae(means, mu)
np.float64(0.11540433749505272)
以及样本中位数。
mae(medians, mu)
np.float64(0.13654429774596036)
平均来说，我们预计样本均值偏离大约 0.115 千克，样本中位数偏离 0.137 千克。所以样本均值可能是更好的选择，至少在这个例子中。
稳健性
现在让我们考虑一个不同的场景。假设有 2%的时间，当你试图称量企鹅时，它意外地按下了秤上的单位按钮，重量被记录为磅而不是千克。假设这个错误没有被注意到，它会在样本中引入一个异常值。
以下函数模拟了这种情况，将 2%的重量乘以每千克 2.2 磅的转换系数。
def make_sample_with_errors(n):
sample = np.random.normal(mu, sigma, size=n)
factor = np.random.choice([1, 2.2], p=[0.98, 0.02], size=n)
return sample * factor
为了看到这对此分布有什么影响，我们将生成一个大样本。
sample = make_sample_with_errors(n=1000)
为了绘制样本的分布，我们将使用核密度估计（KDE）和第六章中的Pdf对象。
from scipy.stats import gaussian_kde
from thinkstats import Pdf
kde = gaussian_kde(sample)
domain = 0, 10
pdf = Pdf(kde, domain)
pdf.plot(label='estimated density')
decorate(xlabel="Penguin weight (kg)", ylabel="Density")
除了接近 3.7 千克的峰值外，测量误差还引入了一个接近 8 千克的第二个峰值。
现在，让我们重复之前的实验，模拟许多大小为 10 的样本，计算每个样本的均值，然后计算样本均值的平均值。
means = [np.mean(make_sample_with_errors(n=10)) for i in range(10001)]
np.mean(means)
np.float64(3.786352945690677)
测量误差导致样本均值平均高于 3.7 千克。
现在让我们用样本中位数进行相同的实验。
medians = [np.median(make_sample_with_errors(n=10)) for i in range(10001)]
np.mean(medians)
np.float64(3.7121869836715353)
样本中位数的平均值也高于 3.7 千克，但它偏离的程度并不大。如果我们比较估计的 MSE，我们会看到样本中位数更准确。
mse(means, mu), mse(medians, mu)
(np.float64(0.06853430354724438), np.float64(0.031164467796883758))
如果测量实际上来自正态分布，样本均值最小化 MSE，但这种情况违反了该假设，因此样本均值不最小化 MSE。样本中位数对异常值不太敏感，因此它更少有偏差，其 MSE 更小。能够很好地处理异常值（以及类似的假设违反）的估计量被称为稳健的。
估计方差
作为另一个例子，假设我们想要估计企鹅体重的方差。在第一章中，我们看到了计算样本方差有两种方法。我承诺稍后解释差异——现在就是时候了。
有两种计算样本方差的方法的原因是其中一种是对总体方差的偏估计，另一种是无偏估计。以下函数计算的是偏估计，即平方偏差之和除以n。
def biased_var(xs):
# Compute variance with n in the denominator
n = len(xs)
deviations = xs - np.mean(xs)
return np.sum(deviations**2) / n
为了测试它，我们将模拟许多大小为 10 的样本，计算每个样本的偏方差，然后计算这些方差的平均值。
biased_vars = [biased_var(make_sample(n=10)) for i in range(10001)]
np.mean(biased_vars)
np.float64(0.19049277659404473)
结果大约是 0.19，但在这个例子中，我们知道实际总体方差大约是 0.21，因此这种样本方差的版本平均来说太低——这证实了它是偏的。
actual_var = sigma**2
actual_var
0.2116
以下函数计算的是无偏估计，即平方偏差之和除以n-1。
def unbiased_var(xs):
# Compute variance with n-1 in the denominator
n = len(xs)
deviations = xs - np.mean(xs)
return np.sum(deviations**2) / (n - 1)
我们可以通过生成许多样本并计算每个样本的无偏方差来测试它。
unbiased_vars = [unbiased_var(make_sample(n=10)) for i in range(10001)]
np.mean(unbiased_vars)
np.float64(0.21159109492300626)
无偏样本方差的平均值非常接近实际值——如果它是无偏的，我们预期会是这样的。
当样本大小为 10 时，偏估计和无偏估计之间的差异大约是 10%，这可能是不可忽视的。当样本大小为 100 时，差异仅为 1%，这足够小，在实际情况中可能不会产生影响。
n = 10
1 - (n - 1) / n
0.09999999999999998
n = 100
1 - (n - 1) / n
0.010000000000000009
``` ## 抽样分布
到目前为止，我们一直在处理模拟数据，假设企鹅体重是从具有已知参数的正态分布中抽取的。现在让我们看看真实数据会发生什么。
在 2007 年至 2010 年期间，南极帕默站的研究人员测量并称重了来自当地种群中的 342 只企鹅。他们收集的数据是免费提供的——下载说明在本书的笔记本中。
以下单元格从由 Allison Horst 创建的存储库中下载数据。
Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R package version 0.1.0\. [`allisonhorst.github.io/palmerpenguins/`](https://allisonhorst.github.io/palmerpenguins/). doi: 10.5281/zenodo.3960218.
数据是作为导致这篇论文的研究的一部分收集的：Gorman KB, Williams TD, Fraser WR (2014).南极企鹅（属 Pygoscelis）群落中的生态性别二态性和环境变异性。PLoS ONE 9(3):e90081\. [`doi.org/10.1371/journal.pone.0090081`](https://doi.org/10.1371/journal.pone.0090081)
```py
download(
"https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)
我们可以使用 Pandas 读取数据。
penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])
penguins.shape
(342, 17)
数据集包括三种企鹅物种。
penguins["Species"].value_counts()
Species
Adelie Penguin (Pygoscelis adeliae) 151
Gentoo penguin (Pygoscelis papua) 123
Chinstrap penguin (Pygoscelis antarctica) 68
Name: count, dtype: int64
对于第一个例子，我们将只选择帝企鹅。
chinstrap = penguins.query('Species.str.startswith("Chinstrap")')
我们将使用这个函数来绘制估计的 PDF。
def plot_kde(sample, name="estimated density", **options):
kde = gaussian_kde(sample)
m, s = np.mean(sample), np.std(sample)
plt.axvline(m, color="gray", ls=":")
domain = m - 4 * s, m + 4 * s
pdf = Pdf(kde, domain, name)
pdf.plot(**options)
这里是帝企鹅体重（千克）的分布情况。垂直虚线表示样本均值。
weights = chinstrap["Body Mass (g)"] / 1000
plot_kde(weights, "weights")
decorate(xlabel="Penguin weight (kg)", ylabel="Density")
样本均值约为 3.7 千克。
sample_mean = np.mean(weights)
sample_mean
np.float64(3.733088235294118)
如果你被要求估计总体均值，3.7 千克是一个合理的选择——但这个估计有多精确呢？
回答那个问题的方法之一是计算均值的抽样分布，它显示了估计均值从一个样本到另一个样本的变化程度。如果我们知道总体中的实际均值和标准差，我们可以模拟抽样过程并计算抽样分布。但如果我们知道实际总体均值，我们就不需要估计它了！
幸运的是，有一种简单的方法可以近似抽样分布，称为重采样。其核心思想是使用样本来构建总体模型，然后使用该模型来模拟抽样过程。
更具体地说，我们将使用参数重采样，这意味着我们将使用样本来估计总体的参数，然后使用理论分布来生成新的样本。
以下函数使用正态分布实现了这个过程。请注意，新的样本大小与原始样本相同。
def resample(sample):
# Generate a sample from a normal distribution
m, s = np.mean(sample), np.std(sample)
return np.random.normal(m, s, len(sample))
此循环使用resample生成许多样本并计算每个样本的均值。
sample_means = [np.mean(resample(weights)) for i in range(1001)]
以下图显示了这些样本均值的分布。
plot_kde(sample_means, "sample means")
decorate(xlabel="Sample mean of weight (kg)", ylabel="Density")
此结果近似了样本均值的抽样分布。它显示了如果我们收集许多相同大小的样本，我们预计样本均值会有多大的变化——假设我们的总体模型是准确的。
非正式地，我们可以看到，如果收集另一个相同大小的样本，样本均值可能低至 3.55，或者高至 3.9。
标准误
为了量化抽样分布的宽度，一个选项是计算其标准差——结果被称为标准误。
standard_error = np.std(sample_means)
standard_error
np.float64(0.04626531069684985)
在这种情况下，标准误约为 0.045 千克——因此如果我们收集许多样本，我们预计样本均值平均变化约为 0.045 千克。
人们经常混淆标准误和标准差。记住：
-
标准差量化了测量中的变化。
-
标准误量化了估计的精确度。
在这个数据集中，帝企鹅重量的标准差约为 0.38 千克。
np.std(weights)
np.float64(0.3814986213564681)
平均重量的标准误约为 0.046 千克。
np.std(sample_means)
np.float64(0.04626531069684985)
标准差告诉你企鹅体重差异有多大。标准误告诉你估计有多精确。它们是针对不同问题的答案。
然而，它们之间存在一种关系。如果我们知道标准差和样本大小，我们可以这样近似均值的标准误：
def approximate_standard_error(sample):
n = len(sample)
return np.std(sample) / np.sqrt(n)
approximate_standard_error(weights)
np.float64(0.046263503290595163)
此结果接近我们通过重采样得到的结果。
置信区间
另一种总结抽样分布的方法是计算置信区间。例如，90%的置信区间包含了抽样分布中的 90%的值，我们可以通过计算第 5 百分位数和第 95 百分位数来找到。以下是帝企鹅平均体重的 90%置信区间。
ci90 = np.percentile(sample_means, [5, 95])
ci90
array([3.6576334 , 3.80737506])
解释置信区间时，人们往往会说有 90%的概率，真实值落在 90%的置信区间内。在这个例子中，我们会说有 90%的概率，帝企鹅的总体均值在 3.66 至 3.81 公斤之间。
在一种称为频率主义的严格概率哲学下，这种解释是不被允许的，在许多统计学书籍中，你会被告知这种解释是错误的。
在我看来，这个禁令过于严格了。在合理的概率哲学下，置信区间意味着人们期望它意味着的东西：有 90%的概率，真实值落在 90%的置信区间内。
然而，置信区间仅量化了由于抽样引起的变异性——也就是说，只测量了部分人口。抽样分布没有考虑到其他错误来源，特别是抽样偏差和测量误差，这些是下一节的主题。
错误来源
假设你想要知道的不是南极企鹅的平均体重，而是你所在城市的女性的平均体重。你无法随机选择一个代表性的女性样本并对其进行称重。
一个简单的替代方案是“电话抽样”——也就是说，你可以从电话簿中随机选择号码，打电话并要求与成年女性通话，询问她的体重。但电话抽样存在明显的问题。
例如，样本仅限于电话号码被列出的人，因此排除了没有电话的人（他们可能比平均水平更穷）和电话号码未列出的人（他们可能比平均水平更富）。此外，如果你在白天打电话回家电话，你不太可能抽样到有工作的人。而且，如果你只抽样接电话的人，你不太可能抽样到共享电话线的人。
如果像收入、就业和家庭规模这样的因素与体重相关——而且这种可能性是合理的——那么你的调查结果会受到这样或那样的影响。这个问题被称为抽样偏差，因为它抽样过程的一个属性。
这个抽样过程也容易受到自我选择的影响，这是一种抽样偏差。有些人会拒绝回答问题，如果拒绝的倾向与体重相关，那么这会影响结果。
最后，如果你询问人们他们的体重，而不是实际称重，结果可能不准确。即使是有帮助的受访者，如果他们对自己的实际体重感到不舒服，也可能向上或向下取整。并不是所有受访者都有帮助。这些不准确是测量误差的例子。
当你报告一个估计量时，通过报告标准误差或置信区间来量化由于采样引起的变异性是有用的。但请记住，这种变异性只是错误的一个来源，而且通常不是最大的来源。
术语表
-
总体均值（population mean）：整个总体中某个量的真实均值，与从子集计算出的样本均值相对。
-
参数（parameter）：一组分布中指定特定分布的值之一——例如，正态分布的参数是均值和标准差。
-
估计量（estimator）：从样本中计算出的统计量，用于估计总体参数。
-
一致（consistent）：如果估计量随着样本大小的增加而收敛到参数的实际值，则估计量是一致的。
-
无偏（unbiased）：如果对于特定的样本大小，样本估计值的平均值等于参数的实际值，则估计量是无偏的。
-
均方误差（MSE）：估计量准确性的度量——它是估计值与真实参数值之间平均平方差的平均值，假设真实值是已知的。
-
稳健（robust）：如果估计量即使在数据集中包含异常值或错误，或者不完全遵循理论分布时仍保持准确，则估计量是稳健的。
-
重采样（resampling）：通过模拟采样过程来近似估计量的采样分布的方法。
-
参数重采样（parametric resampling）：一种重采样方法，它从样本数据中估计总体参数，然后使用理论分布来模拟采样过程。
-
采样分布（sampling distribution）：从同一总体中可能的样本中统计量的分布。
-
标准误差（standard error）：采样分布的标准差，它量化了由于随机采样（但不是测量误差或非代表性采样）导致的估计量的变异性。
-
置信区间（confidence interval）：包含采样分布中最可能值的区间。
-
采样偏差（sampling bias）：在收集样本的方式中存在的缺陷，使得样本无法代表总体。
-
测量误差（measurement error）：数据观察、测量或记录中的不准确。
练习
练习 8.1
重采样方法的一个优点是它们很容易扩展到其他统计量。在本章中，我们计算了企鹅重量的样本均值，然后使用重采样来近似均值的采样分布。现在让我们对标准差做同样的处理。
计算帝企鹅重量的样本标准差。然后使用resample来近似标准差的抽样分布。使用抽样分布来计算估计的标准误差和 90%的置信区间。
练习 8.2
行为风险因素监测系统（BRFSS）数据集包括美国成年人的自我报告身高和体重样本。使用这些数据来估计成年男性的平均身高。使用重采样来近似抽样分布并计算 90%的置信区间。
由于样本量非常大，置信区间非常小，这意味着随机抽样的变异性很小。但其他误差来源可能更大——你认为还有哪些误差来源会影响结果？
以下单元格下载数据，将其读入DataFrame，并选择男性受访者的身高。
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/CDBRFS08.ASC.gz")
from thinkstats import read_brfss
brfss = read_brfss()
male = brfss.query("sex == 1")
heights = male["htm3"]
heights.describe()
count 154407.000000
mean 178.066221
std 7.723563
min 61.000000
25% 173.000000
50% 178.000000
75% 183.000000
max 236.000000
Name: htm3, dtype: float64
练习 8.3
在足球和曲棍球等比赛中，进球之间的时间往往遵循指数分布（如我们在第六章中看到的章节 6）。假设我们观察到了进球之间的时间样本。如果我们假设样本来自指数分布，我们如何估计分布的实际均值？我们可能会考虑使用样本均值或样本中位数。让我们看看它们是否是一致的无偏估计量。对于实验，我们假设实际进球之间的平均时间是 10 分钟。
actual_mean = 10
以下函数从一个具有给定均值和样本大小的指数分布中生成样本。
def make_exponential(n):
return np.random.exponential(actual_mean, size=n)
使用此功能生成一系列不同大小的样本，并计算每个样本的平均值。随着n的增加，样本均值是否会收敛到实际均值？
接下来，生成一系列不同大小的样本，并计算每个样本的中位数。样本中位数是否会收敛到实际中位数？
这是给定均值的指数分布的实际中位数。
actual_median = np.log(2) * actual_mean
actual_median
np.float64(6.931471805599453)
接下来，生成许多大小为 10 的样本，并检查样本均值是否是总体均值的无偏估计量。
最后，检查样本中位数是否是总体中位数的无偏估计量。
练习 8.4
在本章中，我们测试了一个有偏方差估计量，并表明它实际上是有偏的。我们还表明无偏估计量是无偏的。现在让我们尝试标准差。
为了估计总体标准差，我们可以计算方差的有偏或无偏估计值的平方根，如下所示：
def biased_std(sample):
# Square root of the biased estimator of variance
var = biased_var(sample)
return np.sqrt(var)
def unbiased_std(sample):
# Square root of the unbiased estimator of variance
var = unbiased_var(sample)
return np.sqrt(var)
使用make_sample从均值为 3.7 和标准差为 0.46 的正态分布中计算许多大小为 10 的样本。检查这些样本中的任何一个是否是标准差的无偏估计量。
# Here's an example using `make_sample`
mu, sigma = 3.7, 0.46
make_sample(n=10)
array([4.5279695 , 3.75698359, 4.09347143, 3.56308034, 3.17123233,
4.40734952, 3.70858308, 4.15706704, 4.06716703, 3.7203591 ])
练习 8.5
这个练习基于德国坦克问题，这是第二次世界大战期间美国驻伦敦大使馆经济战部门实际分析的一个简化版本。
假设你是一名盟军间谍，你的任务是估计德国建造了多少辆坦克。作为数据，你有从k辆被捕获的坦克中恢复的序列号。
如果我们假设德国有N辆坦克，编号从 1 到N，并且这个范围内的所有坦克被捕获的可能性是相等的，我们可以这样估计N：
def estimate_tanks(sample):
m = np.max(sample)
k = len(sample)
return m + (m - k) / k
以一个例子来说明，假设N是 122。
N = 122
tanks = np.arange(1, N + 1)
我们可以使用以下函数来生成k辆坦克的随机样本。
def sample_tanks(k):
return np.random.choice(tanks, replace=False, size=k)
这里有一个例子。
np.random.seed(17)
sample = sample_tanks(5)
sample
array([74, 71, 95, 10, 17])
下面是基于这个样本的估计。
estimate_tanks(sample)
np.float64(113.0)
检查这个估计器是否有偏差。
关于这个问题的更多内容，请参阅[这个维基百科页面][en.wikipedia.org/wiki/German_tank_problem]以及 Ruggles 和 Brodie 的《二战期间经济情报的经验方法》，美国统计学会杂志，1947 年 3 月，可在这里找到。
关于这个估计器如何工作的解释，你可能喜欢这个视频。
练习 8.6
在几个体育项目中——尤其是篮球——许多球员和球迷相信一种称为“热手”的现象，这表明连续命中几个投篮的球员更有可能命中下一个，而连续失球的球员更有可能失球。
一篇著名的论文提出了一种测试“热手”现象是否真实或是一种错觉的方法，即通过观察职业篮球比赛中的连续命中和失球序列。对于每位球员，作者计算了投篮的整体概率和连续三次命中后的投篮条件概率。对于九位球员中的八位，他们发现连续三次命中后投篮的概率是更低的。基于这一点和其他结果，他们得出结论，没有“连续投篮结果之间存在正相关性的证据”。几十年来，许多人相信“热手”已经被驳倒了。
然而，这个结论是基于一个统计错误，至少部分如此。2018 年的一篇论文显示，第一篇论文中使用的统计量——连续三次命中后的投篮概率——是有偏差的。即使每次投篮的概率正好是 0.5，实际上没有相关性，连续三次命中后投篮的概率也是小于 0.5的。
这并不明显为什么这是真的，这也是为什么错误长期未被察觉的原因，我这里不会尝试解释它。但我们可以使用本章中的方法来检查它。我们将使用以下函数来生成概率为 0.5 且无相关性的 0 和 1 的序列。
def make_hits_and_misses(n):
# Generate a random sequence of 0s and 1s
return np.random.choice([0, 1], size=n)
在本章的笔记本中，我提供了一个函数，该函数可以找到所有三个连续击中（1s）的子序列，并返回序列中跟随的元素。
import numpy as np
def get_successors(seq, target_sum=3):
"""Get the successors of each subsequence that sums to a target value.
Parameters:
seq (array-like): Sequence of 1s and 0s.
target_sum (int): The target sum of the subsequence. Default is 3.
Returns:
np.ndarray: Array of successors to subsequences that sum to `target_sum`.
"""
# Check if the input sequence is too short
if len(seq) < 3:
return np.array([])
# Compute the sum of each subsequence of length 3
kernel = [1, 1, 1]
corr = np.correlate(seq, kernel, mode="valid")
# Find the indices where the subsequence sums to the target value
indices = np.nonzero(corr == target_sum)[0]
# Remove cases where the subsequence is at the end of the sequence
indices = indices[indices < len(seq) - 3]
# Find the successors of each valid subsequence
successors = seq[indices + 3] if len(indices) > 0 else np.array([])
return successors
生成大量长度为 100 的序列，并对每个序列，找出每个紧随三个击中的投篮。计算这些投篮中击中的百分比。提示：如果序列不包含三个连续的击中，该函数将返回一个空序列，因此你的代码必须处理这种情况。
如果你多次运行这个模拟，平均击中百分比是多少？随着序列长度的增加或减少，这个结果如何变化？
著名的论文是 Gilovich, T., Vallone, R., & Tversky, A. (1985). The hot hand in basketball: On the misperception of random sequences. 认知心理学，17(3)，295-314。
展示统计错误的论文是 Miller, J. B., & Sanjurjo, A. (2018). Surprised by the hot hand fallacy? A truth in the law of small numbers. 计量经济学，86(6)，2019-2047。
第一篇论文在此处可用。第二篇在此处可用。关于这个主题的概述和错误解释，你可能喜欢这个视频。
Think Stats: Exploratory Data Analysis in Python, 3rd Edition
版权所有 2024 艾伦·B·唐尼
代码许可：MIT License
文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
假设检验
在这本书中我们探索的数据集中，我们看到了人群——以及企鹅——变量之间的相关性，以及回归线的斜率。这样的结果被称为观察到的效应，因为它们出现在样本中，与通常无法直接观察到的总体中的实际效应形成对比。当我们看到明显的效应时，我们应该考虑它是否可能存在于更大的总体中，或者它是否可能只是由于样本中的偶然性而出现。
有几种方式来表述这个问题，包括费舍尔零假设检验、奈曼-皮尔逊决策理论和贝叶斯假设检验。我这里展示的是这些方法在实践中常用的混合体。
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("Downloaded " + local)
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")
try:
import empiricaldist
except ImportError:
%pip install empiricaldist
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from thinkstats import decorate
抛硬币
我们将从 David MacKay 的书籍《信息论、推理和学习算法》中的一个例子开始，进行一个简单的示例。
当欧元硬币在 2002 年引入时，一位好奇的硬币爱好者将一枚比利时 1 欧元硬币翻转 250 次，并记录下正面朝上 140 次，反面朝上 110 次。如果硬币是完全平衡的，我们预计只有 125 个正面，所以这些数据表明硬币是有偏的。另一方面，我们并不期望每次都能得到正好 125 个正面，所以硬币实际上可能是公平的，而与期望值的偏差可能是由于偶然。为了看看这是否合理，我们可以进行假设检验。
我们将使用以下函数来计算如果硬币是公平的，观察到的数字与期望数字之间的绝对差异。
n = 250
p = 0.5
def abs_deviation(heads):
expected = n * p
return np.abs(heads - expected)
在观察数据中，这个偏差是 15。
heads = 140
tails = 110
observed_stat = abs_deviation(heads)
observed_stat
np.float64(15.0)
如果硬币实际上是公平的，我们可以通过生成一个随机字符串序列来模拟抛硬币实验——要么是 'H' 或 'T'，概率相等——并计算 'H' 出现的次数。
def simulate_flips():
flips = np.random.choice(["H", "T"], size=n)
heads = np.sum(flips == "H")
return heads
每次我们调用这个函数时，我们都会得到一个模拟实验的结果。
np.random.seed(1)
simulate_flips()
np.int64(119)
下面的循环多次模拟实验，计算每个实验的偏差，并使用列表推导式将结果收集到一个列表中。
simulated_stats = [abs_deviation(simulate_flips()) for i in range(10001)]
结果是从假设硬币是公平的偏差分布中抽取的样本。以下是这些值的分布图。
from empiricaldist import Pmf
pmf_effects = Pmf.from_seq(simulated_stats)
pmf_effects.bar()
decorate(xlabel="Absolute deviation", ylabel="PMF")
接近 0 的值是最常见的；大于 10 的值较少见。记住观察数据的偏差是 15，我们可以看到这种程度的偏差是罕见的，但并非不可能。在这个例子中，模拟结果等于或超过 15 的大约 7.1%的时间。
(np.array(simulated_stats) >= 15).mean() * 100
np.float64(7.079292070792921)
因此，如果硬币是公平的，我们预计会偶然出现我们看到的偏差的大约 7.1%的时间。
我们可以得出结论，这种规模的效果并不常见，但即使硬币是公平的，这当然也不是不可能的。基于这个实验，我们不能排除硬币是公平的可能性。
这个例子展示了统计假设检验的逻辑。
-
我们从一个观察开始，250 次投掷中出现了 140 次正面，并提出了一个假设，即硬币是偏斜的——也就是说，正面的概率与 50%不同。
-
我们选择了一个测试统计量，它量化了观察到的效应的大小。在这个例子中，测试统计量是预期结果与实际观察到的绝对偏差。
-
我们定义了一个零假设，这是一个基于观察到的效应是由于偶然性这一假设的模型。在这个例子中，零假设是硬币是公平的。
-
接下来，我们计算了一个p 值，这是在零假设为真的情况下观察到观察到的效应的概率。在这个例子中，p 值是偏差达到 15 或更大的概率。
最后一步是解释结果。如果 p 值很小，我们得出结论，这种效应不太可能偶然发生。如果它很大，我们得出结论，这种效应可能合理地由偶然性解释。如果它位于中间，就像这个例子一样，我们可以说这种效应不太可能偶然发生，但我们不能排除这种可能性。
所有假设检验都基于这些元素——一个测试统计量、一个零假设和一个 p 值。
测试均值差异
在 NSFG 数据中，我们看到了初生儿的平均孕期略长于其他婴儿。现在让我们看看这种差异是否可能是偶然的。
下面的单元格下载数据并安装statadict，这是我们读取数据所需的。
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz")
try:
import statadict
except ImportError:
%pip install statadict
函数get_nsfg_groups读取数据，选择活产，并将活产分为初生儿和其他婴儿。
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
现在我们可以选择两组的孕期，以周为单位。
data = firsts["prglngth"].values, others["prglngth"].values
下面的函数接受数据作为两个序列的元组，并计算均值的绝对差异。
def abs_diff_means(data):
group1, group2 = data
diff = np.mean(group1) - np.mean(group2)
return np.abs(diff)
在初生儿和其他婴儿之间，观察到的孕期差异是 0.078 周。
observed_diff = abs_diff_means(data)
observed_diff
np.float64(0.07803726677754952)
因此，我们将要测试的假设是初生儿的孕期是否通常更长。零假设是两组的孕期实际上相同，而明显的差异是由于偶然性造成的。如果两组的孕期相同，我们可以将两组合并成一个单一的池子。为了模拟实验，我们可以使用 NumPy 函数shuffle将合并的值随机排序，然后使用切片索引选择与原始大小相同的两组。
def simulate_groups(data):
group1, group2 = data
n, m = len(group1), len(group2)
pool = np.hstack(data)
np.random.shuffle(pool)
return pool[:n], pool[-m:]
每次我们调用这个函数时，它都会返回一个序列的元组，我们可以将其传递给abs_diff_means。
abs_diff_means(simulate_groups(data))
np.float64(0.031193045602279312)
以下循环模拟实验多次，并计算每个模拟数据集的均值差异。
simulated_diffs = [abs_diff_means(simulate_groups(data)) for i in range(1001)]
为了可视化结果，我们将使用以下函数，它接受模拟结果的样本并创建一个近似其分布的Pmf对象。
from scipy.stats import gaussian_kde
from empiricaldist import Pmf
def make_pmf(sample, low, high):
kde = gaussian_kde(sample)
qs = np.linspace(low, high, 201)
ps = kde(qs)
return Pmf(ps, qs)
我们还将使用这个函数，它填充分布的尾部。
from thinkstats import underride
def fill_tail(pmf, observed, side, **options):
"""Fill the area under a PMF, right or left of an observed value."""
options = underride(options, alpha=0.3)
if side == "right":
condition = pmf.qs >= observed
elif side == "left":
condition = pmf.qs <= observed
series = pmf[condition]
plt.fill_between(series.index, 0, series, **options)
这是模拟结果的分布情况。阴影区域表示在零假设下均值差异超过观察到的差异的情况。这个区域的面积就是 p 值。
pmf = make_pmf(simulated_diffs, 0, 0.2)
pmf.plot()
fill_tail(pmf, observed_diff, "right")
decorate(xlabel="Absolute difference in means (weeks)", ylabel="Density")
以下函数计算 p 值，即模拟值中大于或等于观察值的比例。
def compute_p_value(simulated, observed):
"""Fraction of simulated values as big or bigger than the observed value."""
return (np.asarray(simulated) >= observed).mean()
在这个例子中，p 值大约是 18%，这意味着有 0.078 周这样的差异可能是偶然发生的。
compute_p_value(simulated_diffs, observed_diff)
np.float64(0.1838161838161838)
根据这个结果，我们不能确定怀孕长度对于第一个孩子来说通常会更长——这个数据集中的差异可能是偶然发生的。
注意，我们在两个假设检验的例子中都看到了相同的元素：一个检验统计量、一个零假设和零假设的模型。在这个例子中，检验统计量是均值之间的绝对差异。零假设是怀孕长度的分布实际上在两组中是相同的。我们通过将两组数据合并到一个池中，打乱池中的数据，并将其分成与原始大小相同的两组来模拟零假设。这个过程被称为排列，也就是打乱的意思。
这种假设检验的计算方法使得将不同统计量结合起来进行测试变得容易。
其他检验统计量
我们可能会想知道，对于第一个孩子来说，怀孕长度不仅仅是更长，也许还更不稳定。为了测试这个假设，我们可以使用两组标准差之间的绝对差异作为检验统计量。以下函数计算这个检验统计量。
def abs_diff_stds(data):
group1, group2 = data
diff = np.std(group1) - np.std(group2)
return np.abs(diff)
在 NSFG 数据集中，标准差之间的差异大约是 0.18。
observed_diff = abs_diff_stds(data)
observed_diff
np.float64(0.17600895913991677)
为了看看这个差异是否可能是偶然发生的，我们可以再次使用排列。以下循环模拟零假设多次，并计算每个模拟数据集的标准差差异。
simulated_diffs = [abs_diff_stds(simulate_groups(data)) for i in range(1001)]
这是结果分布的情况。同样，阴影区域表示在零假设下检验统计量超过观察到的差异的情况。
pmf = make_pmf(simulated_diffs, 0, 0.5)
pmf.plot()
fill_tail(pmf, observed_diff, "right")
decorate(xlabel="Absolute difference in std (weeks)", ylabel="Density")
我们可以通过计算结果中大于或等于观察到的差异的比例来估计这个区域的面积。
compute_p_value(simulated_diffs, observed_diff)
np.float64(0.17082917082917082)
p 值大约是 0.17，所以即使两个组相同，我们也能看到这么大的差异是有可能的。总之，我们不能确定孕妇的怀孕长度通常对第一个孩子来说更具有变异性——我们在这个数据集中看到的不同可能是偶然的。
测试相关性
我们可以使用相同的框架来测试相关性。例如，在 NSFG 数据集中，出生体重和母亲年龄之间存在相关性——年龄较大的母亲平均有更重的婴儿。但这种效果可能是偶然的吗？
为了找出答案，我们将首先准备数据。从活产中，我们将选择母亲年龄和出生体重都已知的情况。
valid = live.dropna(subset=["agepreg", "totalwgt_lb"])
valid.shape
(9038, 244)
然后，我们将选择相关的列。
ages = valid["agepreg"]
birthweights = valid["totalwgt_lb"]
以下函数接受一个包含xs和ys的元组，并计算相关性的幅度，无论是正还是负。
def abs_correlation(data):
xs, ys = data
corr = np.corrcoef(xs, ys)[0, 1]
return np.abs(corr)
在 NSFG 数据集中，相关性约为 0.07。
data = ages, birthweights
observed_corr = abs_correlation(data)
observed_corr
np.float64(0.0688339703541091)
原假设是母亲年龄和出生体重之间没有相关性。通过打乱观察到的值，我们可以模拟一个年龄和出生体重分布相同，但变量无关的世界。
以下函数接受一个包含xs和ys的元组，打乱xs并返回一个包含打乱后的xs和原始ys的元组。如果我们将ys打乱，或者两者都打乱，它也会工作。
def permute(data):
xs, ys = data
new_xs = xs.values.copy()
np.random.shuffle(new_xs)
return new_xs, ys
打乱值的相关性通常接近 0。
abs_correlation(permute(data))
np.float64(0.0019269515502894237)
以下循环生成许多打乱的数据集，并计算每个数据集的相关性。
simulated_corrs = [abs_correlation(permute(data)) for i in range(1001)]
下面是结果分布的图示。垂直虚线表示观察到的相关性。
pmf = make_pmf(simulated_corrs, 0, 0.07)
pmf.plot()
plt.axvline(observed_corr, color="gray", ls=":")
decorate(xlabel="Absolute value of correlation", ylabel="Density")
我们可以看到观察到的相关性位于分布的尾部，曲线下没有可见的区域。如果我们尝试计算 p 值，结果是 0，这表明在模拟中，打乱数据中的相关性没有超过观察到的值。
compute_p_value(simulated_corrs, observed_corr)
np.float64(0.0)
根据这个计算，我们可以得出结论，p 值可能小于 1000 分之一，但实际上并不是零。打乱数据的相关性超过观察值的可能性不大——但并非不可能。
当 p 值很小，传统上小于 0.05 时，我们可以说结果是统计显著的。但这种解释 p 值的方式一直存在问题，并且它正在逐渐被更少地使用。
一个问题是传统的阈值是任意的，并不适用于所有应用。另一个问题是这种对“显著”的使用具有误导性，因为它暗示这种效果在实践中很重要。母亲年龄与出生体重的相关性是一个很好的例子——它在统计上显著，但如此之小，以至于并不重要。
另一种方法是定性解释 p 值。
-
如果 p 值很大，那么观察到的效果可能只是偶然发生的。
-
比例测试
作为最后的例子，让我们考虑一个需要仔细思考测试统计量选择的情况。假设你经营一家赌场，你怀疑一位顾客正在使用一个不诚实的骰子——也就是说，一个被修改过，使得其中一个面比其他面更有可能出现的骰子。你逮捕了这位被怀疑的作弊者并没收了骰子，但现在你必须证明这个骰子是不诚实的。你掷了 60 次骰子，并记录了从 1 到 6 每个结果的频率。这里是在Hist对象中的结果。
from empiricaldist import Hist
qs = np.arange(1, 7)
freqs = [8, 9, 19, 5, 8, 11]
observed = Hist(freqs, qs)
observed.index.name = "outcome"
observed
| 频率 | |
|---|---|
| 5 | 8 |
| 观察数据的卡方统计量是 11.6。这个数字本身并没有什么意义，但我们可以将其与模拟掷骰子的结果进行比较。以下循环生成许多模拟数据集，并为每个数据集计算卡方统计量。 | |
| 如果 p 值很小，我们通常可以排除效果是偶然发生的可能性——但我们应该记住，它仍然可能是由于非代表性抽样或测量误差。 | |
| 2 | 9 |
| 3 | 19 |
| 4 | 5 |
| 1 | 8 |
| 6 | 11 |
平均来说，你期望每个值出现 10 次。在这个数据集中，值 3 比预期出现得更频繁，而值 4 出现得更少。但这些差异可能是偶然发生的吗？
为了测试这个假设，我们首先将计算每个结果的预期频率。
num_rolls = observed.sum()
outcomes = observed.qs
expected = Hist(num_rolls / 6, outcomes)
以下函数接受观察到的和预期的频率，并计算绝对差异的总和。
def total_abs_deviation(observed):
return np.sum(np.abs(observed - expected))
在观察到的数据集中，这个测试统计量是 20。
observed_dev = total_abs_deviation(observed)
observed_dev
np.float64(20.0)
以下函数接受观察到的数据，模拟掷一个公平的骰子相同次数，并返回一个包含模拟频率的Hist对象。
def simulate_dice(observed):
num_rolls = np.sum(observed)
rolls = np.random.choice(observed.qs, num_rolls, replace=True)
hist = Hist.from_seq(rolls)
return hist
以下循环多次模拟实验，并计算每个的绝对偏差总和。
simulated_devs = [total_abs_deviation(simulate_dice(observed)) for i in range(1001)]
这里是测试统计量在零假设下的分布情况。注意，总数总是偶数，因为每次一个结果比预期出现得更频繁时，另一个结果就必须出现得更少。
pmf_devs = Pmf.from_seq(simulated_devs)
pmf_devs.bar()
decorate(xlabel="Total absolute deviation", ylabel="PMF")
我们可以看到，总偏差为 20 并不罕见——p 值约为 13%，这意味着我们无法确定骰子是不诚实的。
compute_p_value(simulated_devs, observed_dev)
np.float64(0.13086913086913088)
但我们选择的测试统计量并不是唯一的选择。对于这样的问题，使用卡方统计量会更传统，我们可以这样计算。
def chi_squared_stat(observed):
diffs = (observed - expected) ** 2
return np.sum(diffs / expected)
将偏差平方（而不是取绝对值）会给大的偏差更多的权重。通过除以expected来标准化偏差——尽管在这种情况下它对结果没有影响，因为预期的频率都是相等的。
observed_chi2 = chi_squared_stat(observed)
observed_chi2
np.float64(11.6)
| 结果 | |
simulated_chi2 = [chi_squared_stat(simulate_dice(observed)) for i in range(1001)]
下面是这种检验统计量在零假设下的分布情况。阴影区域显示了超过观察值的结果。
pmf = make_pmf(simulated_chi2, 0, 20)
pmf.plot()
fill_tail(pmf, observed_chi2, "right")
decorate(xlabel="Chi-Squared statistic", ylabel="Density")
再次，阴影区域的面积是 p 值。
compute_p_value(simulated_chi2, observed_chi2)
np.float64(0.04495504495504495)
使用卡方统计量的 p 值大约为 0.04，这比我们使用总偏差得到的 0.13 小得多。如果我们认真对待 5%的阈值，我们会认为这个效应具有统计学意义。但考虑到两个测试的结果，我会说结论并不明确。我不会排除骰子是弯曲的可能性，但我也不会宣判这个被指控的作弊者有罪。
这个例子演示了一个重要的观点：p 值取决于检验统计量的选择和零假设模型，有时这些选择决定了效应是否具有统计学意义。
术语表
-
假设检验：一组用于检查观察到的效应是否可能是随机抽样的结果的方法。
-
检验统计量：在假设检验中用于量化观察到的效应大小的统计量。
-
零假设：基于样本中观察到的效应在总体中不存在的假设的系统模型。
-
排列：通过随机重新排列数据集来模拟零假设的方法。
-
p 值：在零假设下，效应与观察到的效应一样大的概率。
-
具有统计学意义：如果 p 值小于选定的阈值，通常为 5%，则效应具有统计学意义。在大数据集中，即使观察到的效应在实际中太小而不重要，也可能具有统计学意义。
练习
练习 9.1
让我们尝试使用第八章中的企鹅数据来进行假设检验。下载数据的说明在本书的笔记本中。
下面的单元格从由 Allison Horst 创建的存储库中下载数据。
Horst AM, Hill AP, Gorman KB (2020). palmerpenguins：帕默群岛（南极）企鹅数据。R 包版本 0.1.0。allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218。
这些数据是作为导致本文的研究的一部分收集的：Gorman KB, Williams TD, Fraser WR (2014).南极企鹅（属 Pygoscelis）群落中的生态性两性差异和环境可变性。PLOS ONE 9(3):e90081。doi.org/10.1371/journal.pone.0090081
download(
"https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)
下面是如何读取数据和选择帝企鹅的示例。
penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])
chinstrap = penguins.query('Species.str.startswith("Chinstrap")')
chinstrap.shape
(68, 17)
以下是我们可以如何提取雄性和雌性企鹅的重量（千克）。
male = chinstrap.query("Sex == 'MALE'")
weights_male = male["Body Mass (g)"] / 1000
weights_male.mean()
np.float64(3.9389705882352937)
female = chinstrap.query("Sex == 'FEMALE'")
weights_female = female["Body Mass (g)"] / 1000
weights_female.mean()
np.float64(3.5272058823529413)
使用abs_diff_means和simulate_groups在两组权重分布相同的零假设下生成大量模拟数据集，并计算每个数据集的均值差异。将模拟结果与观察到的差异进行比较，并计算 p 值。两组之间明显的差异是否可能是偶然造成的？
练习 9.2
使用前一个练习中的企鹅数据，我们可以提取雌性企鹅的喙顶深度和长度（喙顶是喙的顶部脊）。
data = female["Culmen Depth (mm)"], female["Culmen Length (mm)"]
这些变量之间的相关系数约为 0.26。
observed_corr = abs_correlation(data)
observed_corr
np.float64(0.2563170802728449)
让我们看看这种相关性是否可能偶然发生，即使实际上测量之间没有相关性。使用permute生成数据的多种排列，并使用abs_correlation计算每个排列的相关性。在零假设下绘制相关性的分布，并计算观察到的相关性的 p 值。你如何解释这个结果？
Think Stats: Exploratory Data Analysis in Python, 3rd Edition
版权所有 2024 艾伦·B·唐尼
代码许可：MIT License
文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 国际
最小二乘法
本章和下一章介绍了将模型拟合到数据中的概念。在这个上下文中，模型由变量之间关系的数学描述（例如一条直线）和随机变异性描述（例如正态分布）组成。
当我们说一个模型拟合数据时，我们通常意味着它最小化了误差，即模型和数据之间的距离。我们将从最广泛使用的一种拟合模型的方法开始，即最小化平方误差之和，这被称为最小二乘拟合。
我们还将从一次只处理两个变量的模型开始。下一章将介绍可以处理超过两个变量的模型。
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("Downloaded " + local)
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")
try:
import empiricaldist
except ImportError:
%pip install empiricaldist
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from thinkstats import decorate
最小二乘拟合
作为第一个例子，让我们回到第八章的场景章节 8。假设你是南极的研究员，研究当地的企鹅种群。作为数据收集的一部分，你捕获了一组企鹅，测量并称重它们——然后未受伤害地释放它们。
如你很快就会了解到，让企鹅在秤上停留足够长的时间以获得准确的测量可能很困难。假设对于一些企鹅，我们有鳍长和喙大小的测量数据，但没有体重。让我们看看我们是否可以使用其他测量数据来填补缺失的数据——这个过程被称为插补。
我们将首先探索重量和测量之间的关系，使用 2007 年至 2010 年间在南极帕默站收集的数据。他们收集的数据是免费提供的——下载说明在本书的笔记本中。
下一个单元将下载由 Allison Horst 创建的存储库中的数据。
Horst AM, Hill AP, Gorman KB (2020)。palmerpenguins：帕默群岛（南极）企鹅数据。R 包版本 0.1.0。allisonhorst.github.io/palmerpenguins/。doi: 10.5281/zenodo.3960218。
这些数据是作为导致这篇论文的研究的一部分收集的：Gorman KB, Williams TD, Fraser WR (2014)。南极企鹅（属 Pygoscelis）群落中的生态性别二态性和环境变异性。PLoS ONE 9(3):e90081。doi.org/10.1371/journal.pone.0090081
download(
"https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)
我们可以使用read_csv来读取数据。
penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])
penguins.shape
(342, 17)
数据集包括对 151 只阿德利企鹅的测量数据。我们可以使用query来选择包含这些数据的行。
adelie = penguins.query('Species.str.startswith("Adelie")')
len(adelie)
151
现在假设我们知道一只阿德利企鹅的鳍长——让我们看看我们如何预测它的体重。首先，我们将从DataFrame中选择这些列。
xvar = "Flipper Length (mm)"
yvar = "Body Mass (g)"
flipper_length = adelie[xvar]
body_mass = adelie[yvar]
这里有一个散点图显示了这些数量之间的关系。
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
decorate(xlabel=xvar, ylabel=yvar)
看起来它们是相关的——我们可以通过计算相关系数来量化这种关系的强度。
np.corrcoef(flipper_length, body_mass)[0, 1]
np.float64(0.4682016942179394)
相关系数约为 0.47，这意味着鳍较长的企鹅往往更重。这很有用，因为它意味着如果我们知道企鹅的鳍长，我们就可以更准确地猜测它的体重——但仅仅相关系数本身并不能告诉我们如何做出这些猜测。为此，我们需要选择一条最佳拟合线。
定义“最佳”线的方法有很多，但对于此类数据，一个常见的选择是线性最小二乘法拟合，这是最小化均方误差（MSE）的直线。
SciPy 提供了一个名为linregress的函数，用于计算最小二乘拟合。这个名字是线性回归的简称，这是此类模型的另一种说法。linregress的参数是x值和y值，顺序如下。
from scipy.stats import linregress
result = linregress(flipper_length, body_mass)
result
LinregressResult(slope=np.float64(32.83168975115009), intercept=np.float64(-2535.8368022002514), rvalue=np.float64(0.46820169421793933), pvalue=np.float64(1.3432645947790051e-09), stderr=np.float64(5.076138407990821), intercept_stderr=np.float64(964.7984274994059))
结果是一个包含拟合线的斜率和截距的LinregressResult对象，以及我们很快就会解包的其他信息。斜率约为 32.8，这意味着每增加 1 毫米的鳍长度，体重就会增加 32.8 克。
截距为-2535 克，这可能看起来有些不合逻辑，因为测量的重量不能是负数。如果我们使用斜率和截距来评估平均鳍长度的拟合线，这可能更有意义。
x = flipper_length.mean()
y = result.intercept + result.slope * x
x, y
(np.float64(189.95364238410596), np.float64(3700.662251655629))
对于鳍长平均约为 190 毫米的企鹅，预期的体重约为 3700 克。
以下函数从linregress的结果和一个xs序列中找到每个x值的拟合线上的点。
def predict(result, xs):
ys = result.intercept + result.slope * xs
return ys
predict这个名字在这里可能看起来有些奇怪——在自然语言中，预测通常与未来发生的事情有关，但在回归的上下文中，拟合线上的点也被称为预测。
我们可以使用predict来计算一系列鳍长度的线上的点。
fit_xs = np.linspace(np.min(flipper_length), np.max(flipper_length))
fit_ys = predict(result, fit_xs)
这是拟合线以及数据的散点图。
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel=xvar, ylabel=yvar)
如预期的那样，拟合线穿过数据的中心并遵循趋势。一些预测是准确的——但许多数据点离线很远。为了了解预测有多好（或有多坏），我们可以计算预测误差，即每个点与线的垂直距离。以下函数计算这些误差，这些误差也称为残差。
def compute_residuals(result, xs, ys):
fit_ys = predict(result, xs)
return ys - fit_ys
这里是体重作为鳍长函数的残差。
residuals = compute_residuals(result, flipper_length, body_mass)
例如，我们可以查看数据集中第一只企鹅的结果。
x = flipper_length[0]
y = predict(result, x)
x, y
(np.float64(181.0), np.float64(3406.699042757914))
被选中的企鹅鳍长为 181 毫米，预测的体重为 3407 克。现在让我们看看实际的体重是多少。
body_mass[0], residuals[0]
(np.float64(3750.0), np.float64(343.30095724208604))
这只企鹅的实际重量是 3750 克，残差（减去预测值后）是 343 克。
平方残差的平均值是预测的均方误差（MSE）。
mse = np.mean(residuals**2)
mse
np.float64(163098.85902884745)
仅凭这个数字本身并没有太多意义。我们可以通过计算决定系数来使它更有意义。
决定系数
假设你想猜测一只企鹅的重量。如果你知道它的鳍长，你可以使用最小二乘法来提供你的猜测，MSE 量化了你的猜测的平均准确性。
但如果你不知道鳍长，你会猜什么？结果证明，猜测均值是最好的策略，因为它最小化了 MSE。如果我们总是猜测均值，预测误差就是与均值的偏差。
deviations = body_mass - np.mean(body_mass)
MSE 是均方偏差。
np.mean(deviations**2)
np.float64(208890.28989956583)
你可能还记得，均方偏差是方差。
np.var(body_mass)
np.float64(208890.28989956583)
因此，如果我们总是猜测均值，可以将质量的方差视为 MSE，如果我们使用回归线，则将残差的方差视为 MSE。如果我们计算这些方差的比率并从 1 中减去它，结果表示如果我们使用鳍长来提供猜测信息，MSE 减少了多少。
以下函数计算这个值，技术上称为 决定系数，但由于它表示为 (R²)，大多数人称之为“R 平方”。
def coefficient_of_determination(ys, residuals):
return 1 - np.var(residuals) / np.var(ys)
在这个例子中，(R²) 大约是 0.22，这意味着拟合线将 MSE 减少了 22%。
R2 = coefficient_of_determination(body_mass, residuals)
R2
np.float64(0.21921282646854912)
结果表明，决定系数 (R²) 和相关系数 (r) 之间存在关系。根据符号，你可能已经猜到了，(r² = R²)。
我们可以通过计算 (R²) 的平方根来证明这一点。
r = np.sqrt(R2)
r
np.float64(0.4682016942179397)
并将其与我们之前计算的相关性进行比较。
corr = np.corrcoef(flipper_length, body_mass)[0, 1]
corr
np.float64(0.4682016942179394)
它们除了由于浮点近似产生的小差异外，都是相同的。
linregress 函数也计算这个值，并将其作为 RegressionResult 对象的属性返回。
result.rvalue
np.float64(0.46820169421793933)
决定系数和相关性传达了大部分相同的信息，但它们的解释不同：
-
相关系数量化了关系强度，其范围从 -1 到 1。
-
(R²) 量化了拟合线减少 MSE 的能力。
此外，(R²) 总是正的，所以它不指示相关性是正还是负。
最小化 MSE
之前我说过，最小二乘法是使均方误差（MSE）最小化的直线。我们不会证明这一点，但我们可以通过向截距和斜率添加小的随机值来测试它，并检查 MSE 是否变得更糟。
intercept = result.intercept + np.random.normal(0, 1)
slope = result.slope + np.random.normal(0, 1)
要运行测试，我们需要创建一个具有 intercept 和 slope 属性的对象——我们将使用 types 模块提供的 SimpleNamespace 对象。
from types import SimpleNamespace
fake_result = SimpleNamespace(intercept=intercept, slope=slope)
fake_result
namespace(intercept=np.float64(-2533.193389187175),
slope=np.float64(32.15393529728294))
我们可以将此对象传递给 compute_residuals 并使用残差来计算 MSE。
fake_residuals = compute_residuals(fake_result, flipper_length, body_mass)
fake_mse = np.mean(fake_residuals**2)
如果我们将结果与最小二乘线的均方误差（MSE）进行比较，它总是更差。
mse, fake_mse, fake_mse > mse
(np.float64(163098.85902884745), np.float64(179019.20812685465), np.True_)
最小化均方误差很好，但“最佳”的定义不止这一个。一个替代方案是最小化误差的绝对值。另一个替代方案是最小化每个点到拟合线的最短距离，这被称为“总误差”。在某些情况下，过高估计可能比过低估计更好（或更差）。在这种情况下，你可能想为每个残差计算一个成本函数，并最小化总成本。
但是，最小二乘法拟合比这些替代方案更广泛地被使用，主要是因为它计算效率高。以下函数展示了如何做到这一点。
def least_squares(xs, ys):
xbar = np.mean(xs)
ybar = np.mean(ys)
xdev = xs - xbar
ydev = ys - ybar
slope = np.sum(xdev * ydev) / np.sum(xdev**2)
intercept = ybar - slope * xbar
return intercept, slope
为了测试这个函数，我们将再次使用翻车鱼长度和体重。
intercept, slope = least_squares(flipper_length, body_mass)
intercept, slope
(np.float64(-2535.8368022002524), np.float64(32.831689751150094))
我们可以确认我们得到的结果与linregress得到的结果相同。
np.allclose([intercept, slope], [result.intercept, result.slope])
True
当计算效率比选择最适合当前问题的方法更重要时，最小化均方误差是有意义的。但情况已不再如此，因此考虑是否应该最小化平方残差是有价值的。
估计
slope（斜率）和intercept（截距）参数是基于样本的估计。像其他估计一样，它们容易受到非代表性抽样、测量误差和随机抽样引起的变异性的影响。通常，很难量化非代表性抽样和测量误差的影响。而量化随机抽样的影响则更容易一些。
实现这一目标的一种方法是称为自助法的重新抽样方法：我们将样本视为整个总体，并从观察到的数据中用放回抽样的方式抽取新的样本。以下函数接受一个DataFrame，使用sample方法重新抽样行，并返回一个新的DataFrame。
def resample(df):
n = len(df)
return df.sample(n, replace=True)
以下函数接受一个DataFrame，找到最小二乘法拟合，并返回拟合线的斜率。
def estimate_slope(df):
xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
result = linregress(xs, ys)
return result.slope
我们可以使用这些函数生成许多模拟数据集，并对每个数据集计算斜率。
# Seed the random number generator so we get the same results every time
np.random.seed(1)
resampled_slopes = [estimate_slope(resample(adelie)) for i in range(1001)]
结果是从斜率抽样分布中的一个样本。它看起来是这样的。
from thinkstats import plot_kde
plot_kde(resampled_slopes)
decorate(xlabel="Slope of the fitted line (g / mm)", ylabel="Density")
我们可以使用percentile来计算 90%置信区间。
ci90 = np.percentile(resampled_slopes, [5, 95])
print(result.slope, ci90)
32.83168975115009 [25.39604591 40.21054526]
因此，我们可以报告估计的斜率为 33 克/毫米，90%置信区间为[25, 40]克/毫米。
估计的标准误差是抽样分布的标准差。
stderr = np.std(resampled_slopes)
stderr
np.float64(4.570238986584832)
从linregress得到的RegressionResult对象提供了一个基于分布形状假设的标准误差近似值。
result.stderr
np.float64(5.076138407990821)
通过重新抽样计算的标准误差略小，但在实际应用中这种差异可能并不重要。
可视化不确定性
每次我们对数据集进行重采样时，都会得到不同的拟合线。为了了解这些线之间的变异性，一个选项是遍历它们并将它们全部绘制出来。以下函数接受一个重采样的 DataFrame，计算最小二乘拟合，并为一系列 xs 生成预测值。
def fit_line(df, fit_xs):
xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
result = linregress(xs, ys)
fit_ys = predict(result, fit_xs)
return fit_ys
这里是我们将使用的 xs 序列。
xs = adelie["Flipper Length (mm)"]
fit_xs = np.linspace(np.min(xs), np.max(xs))
下面是拟合线以及数据的散点图的样子。
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
for i in range(101):
fit_ys = fit_line(resample(adelie), fit_xs)
plt.plot(fit_xs, fit_ys, color="C1", alpha=0.05)
decorate(xlabel=xvar, ylabel=yvar)
在中间附近，拟合线彼此靠近 - 在极端情况下，它们相距较远。
另一种表示拟合线变异性的方法是绘制每个预测值的 90%置信区间。我们可以通过收集拟合线作为一个数组列表来实现这一点。
fitted_ys = [fit_line(resample(adelie), fit_xs) for i in range(1001)]
我们可以将这个数组列表视为一个二维数组，其中每一行对应一条拟合线，每一列对应于每个 xs。
我们可以使用 percentile 函数，并设置 axis=0 参数来找到与每个 xs 对应的 ys 的第 5、第 50 和第 95 百分位数。
low, median, high = np.percentile(fitted_ys, [5, 50, 95], axis=0)
现在我们将使用 fill_between 来绘制位于第 5 百分位数和第 95 百分位数之间的区域，这代表了 90%的置信区间，同时还包括每列的中值和数据的散点图。
plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
plt.fill_between(fit_xs, low, high, color="C1", lw=0, alpha=0.2)
plt.plot(fit_xs, median, color="C1")
decorate(xlabel=xvar, ylabel=yvar)
这是表示由于随机抽样而导致的拟合线变异性的一种我最喜欢的方式。
转换
在拟合数据线之前，有时对其中一个或两个变量进行转换是有用的，例如通过计算值的平方、平方根或对数。为了演示，我们将使用行为风险因素监测系统（BRFSS）中的身高和体重数据，这些数据在第五章中描述。
下一个单元格将下载 BRFSS 数据。
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/CDBRFS08.ASC.gz")
我们可以像这样加载 BRFSS 数据。
from thinkstats import read_brfss
brfss = read_brfss()
接下来，我们将找到包含有效数据的行，并选择包含身高和体重的列。
valid = brfss.dropna(subset=["htm3", "wtkg2"])
heights, weights = valid["htm3"], valid["wtkg2"]
我们可以使用 linregress 函数来计算最小二乘拟合的斜率和截距。
result_brfss = linregress(heights, weights)
result_brfss.intercept, result_brfss.slope
(np.float64(-82.65926054409877), np.float64(0.957074585033226))
斜率大约为 0.96，这意味着平均而言，每增加 1 厘米对应增加近 1 公斤。我们可以再次使用 predict 函数来生成一系列 xs 的预测值。
fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss, fit_xs)
在我们绘制数据的散点图之前，对身高和体重进行抖动是有用的。
from thinkstats import jitter
jittered_heights = jitter(heights, 2)
jittered_weights = jitter(weights, 1.5)
我们将使用身高和体重的平均值和标准差来选择 (x) 轴的界限。
m, s = heights.mean(), heights.std()
xlim = m - 4 * s, m + 4 * s
ylim = 0, 200
这里展示的是经过抖动处理的数据的散点图以及拟合线。
plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Height (cm)", ylabel="Weight (kg)", xlim=xlim, ylim=ylim)
拟合线没有穿过散点图中最密集的部分。这是因为重量并不遵循正态分布。正如我们在第五章中看到的，成年人的体重倾向于遵循对数正态分布，这种分布偏向于较大的值——这些值将拟合线向上拉。
另一个令人担忧的原因是残差的分布，看起来像这样。
residuals = compute_residuals(result_brfss, heights, weights)
from thinkstats import make_pmf
pmf_kde = make_pmf(residuals, -60, 120)
pmf_kde.plot()
decorate(xlabel="Residual (kg)", ylabel="Density")
残差的分布向右偏斜。单从这一点来看，这并不一定是问题，但它表明最小二乘拟合并没有正确地描述这些变量之间的关系。
如果重量遵循对数正态分布，它们的对数遵循正态分布。那么让我们看看如果我们拟合一个线来表示重量对高度的对数函数会发生什么。
log_weights = np.log10(weights)
result_brfss2 = linregress(heights, log_weights)
result_brfss2.intercept, result_brfss2.slope
(np.float64(0.9930804163932876), np.float64(0.005281454169417777))
因为我们变换了一个变量，所以斜率和截距更难解释。但我们可以使用 predict 来计算拟合线。
fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss2, fit_xs)
然后将其与变换数据的散点图一起绘制。
jittered_log_weights = jitter(log_weights, 1.5)
plt.scatter(jittered_heights, jittered_log_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Height (cm)", ylabel="Weight (log10 kg)", xlim=xlim)
拟合线穿过图表中最密集的部分，实际值在直线上下延伸大约相同的距离——因此残差的分布大致是对称的。
residuals = compute_residuals(result_brfss2, heights, log_weights)
pmf_kde = make_pmf(residuals, -0.6, 0.6)
pmf_kde.plot()
decorate(xlabel="Residual (kg)", ylabel="Density")
散点图的外观和残差的分布表明，高度和对数变换后的重量之间的关系很好地被拟合线描述。如果我们比较两个回归的 (r) 值，我们会看到高度与对数变换后的重量的相关性略高。
result_brfss.rvalue, result_brfss2.rvalue
(np.float64(0.5087364789734582), np.float64(0.5317282605983435))
这意味着 (R²) 值也会稍微高一些。
result_brfss.rvalue**2, result_brfss2.rvalue**2
(np.float64(0.2588128050383119), np.float64(0.28273494311893993))
如果我们用高度来猜测重量，如果我们处理的是对数变换后的重量，那么猜测会更好一些。
然而，变换数据使得模型的参数更难解释——在呈现结果之前反转变换可能会有所帮助。例如，以 10 为底的对数的逆是 10 为底的指数。这是逆变换后的拟合线，以及未变换的数据。
plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, 10**fit_ys, color="C1")
decorate(xlabel="Height (cm)", ylabel="Weight (kg)", xlim=xlim, ylim=ylim)
与变换后的数据相比是直线的拟合线，与原始数据相比是曲线的。
术语表
-
模型：在回归的上下文中，模型是变量之间关系的数学描述——例如一条直线——以及随机变化的描述——例如正态分布。
-
插补：一种估计和填补数据集中缺失值的过程。
-
最佳拟合线：一条（或曲线）最好地描述变量之间关系的线（或曲线），通过某种“最佳”的定义。
-
线性回归：一种寻找最佳拟合线的方法。
-
预测：最佳拟合线上的一个点——在回归的上下文中，它不一定是关于未来的声明。
-
残差：观察值与最佳拟合线预测值之间的差异。
-
线性最小二乘拟合：一条使残差平方和最小的线。
-
决定系数：一个统计量，表示为(R²)，通常读作“R 平方”，它量化了模型拟合数据的好坏。
-
自助重采样：一种通过将样本视为总体并以与原始样本相同的大小进行有放回抽样的重采样方法。
练习
练习 10.1
在本章中，我们计算了企鹅体重作为鳍长度的函数的最小二乘拟合。数据集中还有两个其他测量值可以考虑：喙长和喙深（喙是喙顶部的脊）。
计算体重作为喙长函数的最小二乘拟合。绘制这些变量的散点图并绘制拟合线。
根据RegressionResult对象的rvalue属性，这些变量的相关系数是多少？决定系数是多少？哪个是体重更好的预测指标，喙长还是鳍长？
xvar = "Culmen Length (mm)"
yvar = "Body Mass (g)"
culmen_length = adelie[xvar]
body_mass = adelie[yvar]
练习 10.2
在本章中，我们使用重采样来近似拟合线的斜率的抽样分布。我们可以以相同的方式近似截距的抽样分布：
-
编写一个名为
estimate_intercept的函数，该函数接受一个重采样的DataFrame作为参数，计算企鹅体重作为鳍长度的最小二乘拟合，并返回截距。 -
使用许多重采样的
adelie调用该函数并收集截距。 -
使用
plot_kde绘制截距的抽样分布。 -
计算标准误差和 90%置信区间。
-
检查从重采样得到的标准误差是否与
RegressionResult对象中的intercept_stderr属性一致——它可能略小一些。
练习 10.3
人的体质指数（BMI）是他们的体重（千克）除以他们的身高（米）的平方。在 BRFSS 数据集中，我们可以这样计算 BMI，在将身高从厘米转换为米之后。
heights_m = heights / 100
bmis = weights / heights_m**2
在这个定义中，高度是平方的，而不是提高到某个其他指数，这是因为观察——在统计学历史早期——平均体重大致与高度平方成正比。
为了验证这一点，我们可以使用 BRFSS 数据、最小二乘拟合和一些数学知识。假设体重与高度的不知指数(a)成正比。在这种情况下，我们可以写出：
[w = b h^a]
其中 (w) 是体重，(h) 是身高，(b) 是一个未知的比例常数。对等式两边取对数：
[\log w = \log b + a \log h]
因此，如果我们计算以对数变换的高度为函数的对数变换重量的最小二乘拟合，拟合线的斜率估计了未知的指数 (a)。
计算身高和体重的对数。你可以使用任何底数进行对数变换，只要两个变换使用相同的底数即可。计算最小二乘拟合。斜率接近 2 吗？
《Python 数据探索：Think Stats》第 3 版
版权所有 2024 艾伦·B·唐尼
代码许可：MIT 许可
文本许可：Creative Commons 署名-非商业性使用-相同方式共享 4.0 国际
多元回归
上一章中的线性最小二乘拟合是回归的一个例子，它是建模一组变量（称为响应变量或因变量）与另一组变量（称为解释变量或自变量）之间关系的一个更一般的问题。
在上一章的示例中，只有一个响应变量和一个解释变量，这被称为简单回归。在本章中，我们将继续探讨多元回归，其中包含多个解释变量，但仍然只有一个响应变量。如果有多个响应变量，那就是多元回归，我们在这本书中不会涉及。
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("Downloaded " + local)
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")
try:
import empiricaldist
except ImportError:
%pip install empiricaldist
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from thinkstats import decorate
StatsModels
在上一章中，我们使用了 SciPy 函数 linregress 来计算最小二乘拟合。此函数执行简单回归，但不能执行多元回归。为此，我们将使用 StatsModels，这是一个提供多种回归和其他分析的包。
作为第一个例子，我们将继续探索企鹅数据。
以下单元格从 Allison Horst 创建的存储库中下载数据。
Horst AM, Hill AP, Gorman KB (2020)。palmerpenguins：帕默群岛（南极）企鹅数据。R 包版本 0.1.0。allisonhorst.github.io/palmerpenguins/。doi: 10.5281/zenodo.3960218。
数据是作为导致这篇论文的研究的一部分收集的：Gorman KB, Williams TD, Fraser WR (2014)。南极企鹅（属 Pygoscelis）群落中的生态性别二态性和环境可变性。PLoS ONE 9(3):e90081。doi.org/10.1371/journal.pone.0090081
download(
"https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)
当我们加载数据时，我们将使用以下字典来给出不包含空格的列名，这将使它们更容易与 StatsModels 一起使用。
columns = {
"Body Mass (g)": "mass",
"Flipper Length (mm)": "flipper_length",
"Culmen Length (mm)": "culmen_length",
"Culmen Depth (mm)": "culmen_depth",
}
现在我们可以加载数据，删除缺失体重的行，并重命名列。
penguins = (
pd.read_csv("penguins_raw.csv")
.dropna(subset=["Body Mass (g)"])
.rename(columns=columns)
)
penguins.shape
(342, 17)
数据集包含三种企鹅物种。我们将仅使用阿德利企鹅。
adelie = penguins.query('Species.str.startswith("Adelie")').copy()
len(adelie)
151
在上一章中，我们计算了企鹅的鳍长和体重之间的最小二乘拟合。
flipper_length = adelie["flipper_length"]
body_mass = adelie["mass"]
作为提醒，以下是使用 linregress 实现的方法。
from scipy.stats import linregress
result_linregress = linregress(flipper_length, body_mass)
result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002514), np.float64(32.83168975115009))
StatsModels 提供了两个接口（APIs）——我们将使用“公式”API，它使用 Patsy 公式语言来指定响应变量和解释变量。以下公式字符串指定响应变量 mass 是一个解释变量 flipper_length 的线性函数。
formula = "mass ~ flipper_length"
我们可以将此公式传递给 StatsModels 函数 ols，以及数据。
import statsmodels.formula.api as smf
model = smf.ols(formula, data=adelie)
type(model)
statsmodels.regression.linear_model.OLS
名称 ols 代表“普通最小二乘法”，这表明此函数在最常见的或“普通”的假设下计算最小二乘拟合。
结果是一个表示模型的 OLS 对象。一般来说，模型是变量之间关系的简化描述。在这个例子中，它是一个线性模型，这意味着它假设响应变量是解释变量的线性函数。
fit 方法将模型拟合到数据，并返回一个包含结果的 RegressionResults 对象。
result_ols = model.fit()
# Technically it's a RegressionResultsWrapper
type(result_ols)
statsmodels.regression.linear_model.RegressionResultsWrapper
RegressionResults 对象包含大量信息，因此 thinkstats 提供了一个仅显示我们现在所需信息的函数。
from thinkstats import display_summary
display_summary(result_ols)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -2535.8368 | 964.798 | -2.628 | 0.009 | -4442.291 | -629.382 |
| flipper_length | 32.8317 | 5.076 | 6.468 | 0.000 | 22.801 | 42.862 |
| R-squared: | 0.2192 |
第一列包含截距和斜率，即模型的 系数。我们可以确认它们与我们从 linregress 得到的系数相同。
result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002514), np.float64(32.83168975115009))
第二列包含系数的标准误差 – 同样，它们与我们从 linregress 得到的值相同。
result_linregress.intercept_stderr, result_linregress.stderr,
(np.float64(964.7984274994059), np.float64(5.076138407990821))
下一个列报告 (t) 统计量，用于计算 p 值 – 我们可以忽略它们，因为 p 值在下一列，标记为 P>|t|。flipper_length 的 p 值四舍五入到 0，但我们可以这样显示。
result_ols.pvalues["flipper_length"]
np.float64(1.3432645947789321e-09)
并确认 linregress 计算了相同的结果。
result_linregress.pvalue
np.float64(1.3432645947790051e-09)
p 值非常小，这意味着如果实际上体重和鳍长之间没有关系，我们几乎不可能偶然看到与估计值一样大的斜率。
最后两列，标记为 [0.025 和 0.975]，报告了截距和斜率的 95% 置信区间。因此，斜率的 95% 置信区间是 [22.8, 42.9]。
最后一行报告了模型的 (R²) 值，大约为 0.22 – 这意味着如果我们使用鳍长来预测体重，与只使用平均体重相比，我们可以将均方误差降低约 22%。
从简单相关得到的 (R²) 值是相关系数 (r) 的平方。因此，我们可以将 ols 计算出的 rsquared 与 linregress 计算出的 rvalue 的平方进行比较。
result_ols.rsquared, result_linregress.rvalue**2
(np.float64(0.21921282646854878), np.float64(0.21921282646854875))
它们除了由于浮点近似产生的小差异外，都是相同的。
在我们进行多元回归之前，让我们计算一个更简单的回归，以 culmen_length 作为解释变量（culmen 是喙的顶部脊）。
formula = "mass ~ culmen_length"
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| 截距 | 34.8830 | 458.439 | 0.076 | 0.939 | -870.998 | 940.764 |
| culmen_length | 94.4998 | 11.790 | 8.015 | 0.000 | 71.202 | 117.798 |
| R-squared: | 0.3013 |
再次，斜率的 p 值非常小，这意味着如果实际上质量和喙长之间没有关系，我们不太可能偶然看到这么大的斜率。你可能注意到与截距相关的 p 值较大，但这不是问题，因为我们不关心截距是否可能为零。在这个模型中，截距接近零，但这只是一个巧合——它并不表明模型存在问题。
这个模型的 (R²) 值约为 0.30，所以如果我们使用喙长而不是鳍长作为解释变量，MSE 的减少会略高（鳍长的 (R²) 值为 0.22）。现在，让我们看看将它们结合起来会发生什么。
接下来是多元回归
这里是关于质量是鳍长和喙长线性函数的多元回归模型的 Patsy 公式。
formula = "mass ~ flipper_length + culmen_length"
这里是将此模型拟合到数据的结果。
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -3573.0817 | 866.739 | -4.122 | 0.000 | -5285.864 | -1860.299 |
| flipper_length | 22.7024 | 4.742 | 4.787 | 0.000 | 13.331 | 32.074 |
| culmen_length | 76.3402 | 11.644 | 6.556 | 0.000 | 53.331 | 99.350 |
| R-squared: | 0.3949 |
这个模型有三个系数：一个截距和两个斜率。与鳍长相关的斜率是 22.7，这意味着我们预计鳍长增加一毫米的企鹅，重量会增加 22.7 克——假设喙长相同。同样，我们预计喙长增加一毫米的企鹅，重量会增加 76.3 克——假设鳍长相同。
与两个斜率相关的 p 值都很小，这意味着两个解释变量的贡献不太可能偶然发生。
而 (R²) 值为 0.39，高于只有喙长（0.30）的模型和只有鳍长（0.22）的模型。因此，基于两个解释变量的预测比基于任何一个单独的解释变量都要好。
但它们并不像我们希望的那样好。如果鳍长减少了 MSE 的 22%，喙长减少了 30%，为什么两者加在一起不能总共减少 52% 呢？原因是解释变量彼此相关。
from thinkstats import corrcoef
corrcoef(adelie, "flipper_length", "culmen_length")
np.float64(0.32578471516515944)
一只鳍更长的企鹅平均来说喙也更长。解释变量包含一些彼此之间的信息，这意味着它们包含关于响应变量的部分相同信息。当我们向模型添加一个解释变量时，(R²) 的提高仅反映了新变量提供的新信息。
如果我们添加喙深作为第三个解释变量，我们会看到相同的模式。
formula = "mass ~ flipper_length + culmen_length + culmen_depth"
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -4341.3019 | 795.117 | -5.460 | 0.000 | -5912.639 | -2769.964 |
| flipper_length | 17.4215 | 4.385 | 3.973 | 0.000 | 8.756 | 26.087 |
| culmen_length | 55.3676 | 11.133 | 4.973 | 0.000 | 33.366 | 77.369 |
| culmen_depth | 140.8946 | 24.216 | 5.818 | 0.000 | 93.037 | 188.752 |
| R-squared: | 0.5082 |
这个模型有四个系数。所有的 p 值都很小，这意味着每个解释变量的贡献不太可能偶然发生。并且 (R²) 值约为 0.51，略好于具有两个解释变量的先前模型（0.39），并且优于具有单个变量的任何模型（0.22 和 0.30）。
但再次，增量改进小于我们可能希望的，因为 culmen 深度与其他两个测量值相关。
[
corrcoef(adelie, "culmen_depth", "flipper_length"),
corrcoef(adelie, "culmen_depth", "culmen_length"),
]
[np.float64(0.30762017939668534), np.float64(0.39149169183587634)]
这个例子演示了多元回归的常见用途，通过结合多个解释变量来做出更好的预测。另一个常见用途是在控制另一组变量的贡献的同时，量化一组变量的贡献。
控制变量
在第 第四章 中，我们看到了第一个孩子平均比其他孩子轻。在第 第九章 中，我们看到了出生体重与母亲年龄相关——年龄较大的母亲平均有更重的孩子。
这些结果可能相关。如果第一个孩子的母亲比其他孩子的母亲年轻——这似乎很可能是这样——那么这或许可以解释为什么他们的孩子体重较轻。我们可以通过估计第一个孩子和其他孩子在出生体重上的差异，同时控制母亲的年龄，来使用多元回归测试这个假设。
下载 NSFG 数据的说明在本书的笔记本中。
以下单元格下载数据文件并安装 statadict，这是我们读取数据所需的。
download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz")
try:
import statadict
except ImportError:
%pip install statadict
我们可以使用 get_nsfg_groups 来读取数据，选择活产，并将活产分为第一个孩子和其他孩子。
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
我们将使用 dropna 来选择具有有效出生体重、出生顺序和母亲年龄的行。
valid = live.dropna(subset=["agepreg", "birthord", "totalwgt_lb"]).copy()
现在，我们可以使用 StatsModels 来确认出生体重与年龄相关，并估计斜率——假设这是一种线性关系。
formula = "totalwgt_lb ~ agepreg"
result_age = smf.ols(formula, data=valid).fit()
display_summary(result_age)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 6.8304 | 0.068 | 100.470 | 0.000 | 6.697 | 6.964 |
| agepreg | 0.0175 | 0.003 | 6.559 | 0.000 | 0.012 | 0.023 |
| R-squared: | 0.004738 |
斜率很小，只有每年 0.0175 磅。所以如果两个母亲的年龄相差十年，我们预计他们的孩子的体重将相差 0.175 磅。但是 p 值很小，所以即使实际上没有关系，这个斜率——尽管很小——也很可能不会出现。
(R²) 值也较小，这意味着母亲的年龄作为预测变量并不是非常有用。如果我们知道母亲的年龄，我们预测婴儿体重的能力几乎不会有所提高。
这种小 p 值和小(R²)值的组合是常见的混淆来源，因为它似乎自相矛盾——如果关系在统计上显著，它似乎应该具有预测性。但这个例子表明，这并不矛盾——一个关系可以是统计上显著的，但并不非常适合预测。如果我们可视化结果，我们会看到原因。首先，让我们选择相关的列。
totalwgt = valid["totalwgt_lb"]
agepreg = valid["agepreg"]
要计算拟合线，我们可以从result_age中提取截距和斜率，但不必这样做。RegressionResults对象提供了一个我们可以使用的predict方法。首先，我们将计算agepreg的值范围。
agepreg_range = np.linspace(agepreg.min(), agepreg.max())
要使用predict，我们必须将解释变量的值放入一个DataFrame中。
df = pd.DataFrame({"agepreg": agepreg_range})
DataFrame中的列必须与解释变量具有相同的名称。然后我们才能将其传递给predict。
fit_ys = result_age.predict(df)
结果是一个包含预测值的Series。以下是它们的形状，以及数据的散点图。
plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)
plt.plot(agepreg_range, fit_ys, color="C1", label="linear model")
decorate(xlabel="Maternal age", ylabel="Birth weight (pounds)")
因为拟合线的斜率很小，所以我们几乎看不到年轻母亲和年长母亲预期出生体重的差异。在每一个母亲年龄，出生体重的变化要大得多。
接下来，我们将使用StatsModels来确认第一个出生的孩子比其他孩子轻。为了实现这一点，我们将创建一个布尔Series，对于第一个出生的孩子为True，对于其他人则为False，并将其添加为名为is_first的新列。
valid["is_first"] = valid["birthord"] == 1
from thinkstats import value_counts
# check the results
value_counts(valid["is_first"])
is_first
False 4675
True 4363
Name: count, dtype: int64
这是将出生体重作为响应变量，is_first作为解释变量的模型的公式。在 Patsy 公式语言中，C和变量名周围的括号表示它是分类的——也就是说，它代表像“第一个孩子”这样的类别，而不是像出生体重这样的测量值。
formula = "totalwgt_lb ~ C(is_first)"
现在我们可以拟合模型并显示结果，就像往常一样。
result_first = smf.ols(formula, data=valid).fit()
display_summary(result_first)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 7.3259 | 0.021 | 356.007 | 0.000 | 7.286 | 7.366 |
| C(is_first)[T.True] | -0.1248 | 0.030 | -4.212 | 0.000 | -0.183 | -0.067 |
| R-squared: | 0.00196 |
在结果中，标签C(is_first)[T.True]表示is_first是一个分类变量，系数与值True相关。True前面的T代表“治疗”——在控制实验的语言中，第一个出生的孩子被视为治疗组，其他孩子被视为参照组。这些指定是任意的——我们也可以考虑第一个出生的孩子为参照组，其他人为治疗组。但我们需要知道哪个是哪个，以便解释结果。
截距约为 7.3，这意味着参考组的平均体重为 7.3 磅。is_first 的系数为 -0.12，这意味着治疗组（第一个孩子）的平均体重比参考组轻 0.12 磅。我们可以通过直接计算来检查这两个结果。
others["totalwgt_lb"].mean()
np.float64(7.325855614973262)
diff_weight = firsts["totalwgt_lb"].mean() - others["totalwgt_lb"].mean()
diff_weight
np.float64(-0.12476118453549034)
除了这些系数外，StatsModels 还计算了 p 值、置信区间和 (R²)。与第一个孩子相关的 p 值很小，这意味着两组之间的差异在统计上具有显著性。而 (R²) 值很小，这意味着如果我们试图猜测一个婴儿的体重，知道它是否是第一个孩子并不会有多大帮助。
现在我们来看一下，出生体重的差异是否可能是由于母亲年龄的差异。平均而言，第一个孩子的母亲比其他母亲年轻约 3.6 岁。
diff_age = firsts["agepreg"].mean() - others["agepreg"].mean()
diff_age
np.float64(-3.5864347661500275)
并且出生体重作为年龄函数的斜率是每年 0.0175 磅。
slope = result_age.params["agepreg"]
slope
np.float64(0.017453851471802638)
如果我们将斜率乘以年龄差异，我们就能得到由于母亲年龄差异而导致的第一个孩子和其他孩子预期出生体重的差异。
slope * diff_age
np.float64(-0.0625970997216918)
结果是 0.063 磅，大约是观察到的差异的一半。所以看起来观察到的出生体重差异可以部分地用母亲年龄的差异来解释。
使用多元回归，我们可以同时估计母亲年龄和第一个孩子的系数。
formula = "totalwgt_lb ~ agepreg + C(is_first)"
result = smf.ols(formula, data=valid).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 6.9142 | 0.078 | 89.073 | 0.000 | 6.762 | 7.066 |
| C(is_first)[T.True] | -0.0698 | 0.031 | -2.236 | 0.025 | -0.131 | -0.009 |
| agepreg | 0.0154 | 0.003 | 5.499 | 0.000 | 0.010 | 0.021 |
| R-squared: | 0.005289 |
is_first 的系数为 -0.0698，这意味着在考虑了母亲年龄差异后，第一个孩子平均比其他孩子轻 0.0698 磅。这大约是未考虑母亲年龄差异时差异的一半。
并且 p 值为 0.025，这仍然被认为是统计上显著的，但它处于一个边缘范围，我们不能排除这种大小的差异可能偶然发生的可能性。
由于这个模型考虑了由于母亲年龄导致的体重差异，我们可以说它控制了母亲年龄。但它假设体重和母亲年龄之间的关系是线性的。所以让我们看看这是否正确。
非线性关系
为了检查 agepreg 的贡献是否可能是非线性的，我们可以在数据集中添加一个新列，其中包含 agepreg 的平方值。
valid["agepreg2"] = valid["agepreg"] ** 2
现在我们可以定义一个包含线性关系和二次关系的模型。
formula = "totalwgt_lb ~ agepreg + agepreg2"
我们可以用通常的方式拟合这个模型。
result_age2 = smf.ols(formula, data=valid).fit()
display_summary(result_age2)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 5.5720 | 0.275 | 20.226 | 0.000 | 5.032 | 6.112 |
| agepreg | 0.1186 | 0.022 | 5.485 | 0.000 | 0.076 | 0.161 |
| agepreg2 | -0.0019 | 0.000 | -4.714 | 0.000 | -0.003 | -0.001 |
| R-squared: | 0.00718 |
与二次项agepreg2相关的 p 值非常小，这表明它比我们预期的偶然性提供了更多关于出生体重的信息。并且该模型的(R²)值为 0.0072，高于线性模型的(R²)值（0.0047）。
通过估计agepreg和agepreg2的系数，我们实际上是在将抛物线拟合到数据中。为了验证这一点，我们可以使用RegressionResults对象为一系列母体年龄生成预测。
首先，我们将创建一个临时的DataFrame，其中包含名为agepreg和agepreg2的列，基于agepreg_range中的年龄范围。
df = pd.DataFrame({"agepreg": agepreg_range})
df["agepreg2"] = df["agepreg"] ** 2
现在，我们可以使用predict方法，将DataFrame作为参数传递，并返回一个预测的Series。
fit_ys = result_age2.predict(df)
下面是拟合的抛物线以及数据的散点图。
plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)
plt.plot(agepreg_range, fit_ys, color="C1", label="quadratic model")
decorate(xlabel="Maternal age", ylabel="Birth weight (pounds)")
曲率微妙，但暗示着最年轻和最年长的母亲的出生体重较低，而中间的较高。
二次模型比线性模型更好地捕捉这些变量之间的关系，这意味着它可以更有效地解释由于母体年龄差异导致的出生体重差异。因此，让我们看看当我们将is_first添加到二次模型中会发生什么。
formula = "totalwgt_lb ~ agepreg + agepreg2 + C(is_first)"
result = smf.ols(formula, data=valid).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 5.6923 | 0.286 | 19.937 | 0.000 | 5.133 | 6.252 |
| C(is_first)[T.True] | -0.0504 | 0.031 | -1.602 | 0.109 | -0.112 | 0.011 |
| agepreg | 0.1124 | 0.022 | 5.113 | 0.000 | 0.069 | 0.155 |
| agepreg2 | -0.0018 | 0.000 | -4.447 | 0.000 | -0.003 | -0.001 |
| R-squared: | 0.007462 |
通过更有效地控制母体年龄，我们估计第一个孩子与其他孩子之间的体重差异为 0.0504 磅，小于仅使用线性模型时的估计（0.0698 磅）。并且与is_first相关的 p 值为 0.109，这意味着这些组之间剩余的差异可能是由于偶然性。
我们可以得出结论，出生体重的差异至少部分地，可能全部地，是由母亲年龄的差异所解释的。
逻辑回归
线性回归基于一个模型，其中响应变量的期望值是解释变量的加权和以及一个截距。当响应变量是连续量，如出生体重或企鹅体重时，此模型是合适的，但当响应变量是离散量，如计数或类别时，则不合适。
对于这些类型的响应变量，我们可以使用广义线性模型或 GLMs。例如：
-
如果响应变量是计数，我们可以使用泊松回归。
-
如果它是只有两个类别的分类变量，我们可以使用逻辑回归。
-
如果它是具有两个以上类别的分类变量，我们可以使用多项逻辑回归。
-
如果它是有序的，并且类别可以按顺序排列，我们可以使用有序逻辑回归。
我们不会在这本书中涵盖所有这些——只涵盖最常用的逻辑回归。作为一个例子，我们将再次使用企鹅数据集，看看我们是否可以根据其体重和其他测量值来判断企鹅是雄性还是雌性。
StatsModels 提供了一个进行逻辑回归的函数——它被称为 logit，因为这是逻辑回归定义中出现的数学函数的名称。在我们能够使用 logit 函数之前，我们必须将响应变量转换，使其值为 0 和 1。
adelie["y"] = (adelie["Sex"] == "MALE").astype(int)
adelie["y"].value_counts()
y
0 78
1 73
Name: count, dtype: int64
我们将从以y作为响应变量和质量作为解释变量的简单模型开始。以下是创建和拟合模型的方法——参数disp=False抑制了关于拟合过程的输出。
model = smf.logit("y ~ mass", data=adelie)
result = model.fit(disp=False)
下面是结果。
display_summary(result)
| 系数 | 标准误 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| 截距 | -25.9871 | 4.221 | -6.156 | 0.000 | -34.261 | -17.713 |
| 质量 | 0.0070 | 0.001 | 6.138 | 0.000 | 0.005 | 0.009 |
| 伪 R 平方 | 0.5264 |
决定系数 (R²) 不适用于逻辑回归，但有一些替代方案被用作“伪 (R²) 值”。该模型的伪 (R²) 值约为 0.526，这本身并没有什么意义，但我们将用它来比较模型。
质量的系数为正，这意味着体重较重的企鹅更有可能是雄性。除此之外，系数不易解释——我们可以通过绘制预测图来更好地理解模型。我们将创建一个包含质量值范围的DataFrame，并使用predict来计算一系列预测值。
mass = adelie["mass"]
mass_range = np.linspace(mass.min(), mass.max())
df = pd.DataFrame({"mass": mass_range})
fit_ys = result.predict(df)
每个预测值都是企鹅为雄性的概率，这是其体重的函数。以下是预测值的外观。
plt.plot(mass_range, fit_ys)
decorate(xlabel="Mass (g)", ylabel="Prob(male)")
最轻的企鹅几乎肯定是雌性，而最重的企鹅很可能是雄性——在中间，体重为 3750 克的企鹅是雄性或雌性的可能性大致相等。
现在我们来看看如果我们添加其他测量值作为解释变量会发生什么。
formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"
model = smf.logit(formula, data=adelie)
result = model.fit(disp=False)
display_summary(result)
| 系数 | 标准误 | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| 截距 | -60.6075 | 13.793 | -4.394 | 0.000 | -87.642 | -33.573 |
| 质量 | 0.0059 | 0.001 | 4.153 | 0.000 | 0.003 | 0.009 |
| 翼展长度 | -0.0209 | 0.052 | -0.403 | 0.687 | -0.123 | 0.081 |
| 顶点长度 | 0.6208 | 0.176 | 3.536 | 0.000 | 0.277 | 0.965 |
| 顶点深度 | 1.0111 | 0.349 | 2.896 | 0.004 | 0.327 | 1.695 |
| 伪 R 平方 | 0.6622 |
该模型的伪 (R²) 值为 0.662，高于之前的模型（0.526）——因此，额外的测量包含区分雄性和雌性企鹅的信息。
culmen_length 和 depth 的 p 值很小，这表明它们比我们预期的偶然性贡献了更多信息。flipper_length 的 p 值很大，这表明如果你知道企鹅的体重和喙尺寸，flipper_length 不会提供额外的信息。
为了理解这个模型，让我们看看它的某些预测。我们将使用以下函数，它接受一系列质量和特定的culmen_length值。它将其他测量值设置为它们的平均值，计算作为质量函数的预测概率，并绘制结果。
def plot_predictions(mass_range, culmen_length, **options):
"""Plot predicted probabilities as a function of mass."""
df = pd.DataFrame({"mass": mass_range})
df["flipper_length"] = adelie["flipper_length"].mean()
df["culmen_length"] = culmen_length
df["culmen_depth"] = adelie["culmen_depth"].mean()
fit_ys = result.predict(df)
plt.plot(mass_range, fit_ys, **options)
这里是culmen_length三个值的结果：平均值的 1 个标准差以上，平均值，和平均值的 1 个标准差以下。
culmen_length = adelie["culmen_length"]
m, s = culmen_length.mean(), culmen_length.std()
plot_predictions(mass_range, m + s, ls="--", label="Above average culmen length")
plot_predictions(mass_range, m, alpha=0.5, label="Average culmen length")
plot_predictions(mass_range, m - s, ls=":", label="Below average culmen length")
decorate(xlabel="Mass (g)", ylabel="Prob(male)")
正如我们在更简单的模型中看到的，体重较重的企鹅更有可能是雄性。此外，在任何体重下，喙较长的企鹅更有可能是雄性。
这个模型比它看起来更有用——事实上，它与用于收集此数据集的原始研究论文中使用的模型相似。该研究的主要主题是性二形性，即雄性和雌性身体差异的程度。量化二形性的方法之一是使用测量值来分类雄性和雌性。在一个二形性较高的物种中，我们预计这些分类将更加准确。
原始论文是 Gorman KB, Williams TD, Fraser WR (2014)。南极企鹅（属 Pygoscelis）群落中的生态性二形性和环境可变性。PLoS ONE 9(3):e90081。 doi.org/10.1371/journal.pone.0090081
为了测试这个方法，让我们尝试在另一种物种上使用相同的模型。除了我们之前工作的阿德利企鹅外，数据集还包含 123 只金图企鹅的测量值。我们将使用以下函数来选择它们。
def get_species(penguins, species):
df = penguins.query(f'Species.str.startswith("{species}")').copy()
df["y"] = (df["Sex"] == "MALE").astype(int)
return df
gentoo = get_species(penguins, "Gentoo")
len(gentoo)
123
这里是逻辑回归模型的预测结果。
formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"
model = smf.logit(formula, data=gentoo)
result = model.fit(disp=False)
display_summary(result)
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -173.9123 | 62.326 | -2.790 | 0.005 | -296.069 | -51.756 |
| mass | 0.0105 | 0.004 | 2.948 | 0.003 | 0.004 | 0.017 |
| flipper_length | 0.2839 | 0.183 | 1.549 | 0.121 | -0.075 | 0.643 |
| culmen_length | 0.2734 | 0.285 | 0.958 | 0.338 | -0.286 | 0.833 |
| culmen_depth | 3.0843 | 1.291 | 2.389 | 0.017 | 0.554 | 5.614 |
| 伪(R²)： | 0.848 |
伪(R²)值为 0.848，高于我们用阿德利企鹅得到的值（0.662）。这意味着与阿德利企鹅相比，可以使用物理测量值更准确地分类金图企鹅，这表明金图企鹅的性二形性更高。
术语表
-
回归：一种估计系数的方法，使模型与数据拟合。
-
响应变量：回归模型试图预测的变量，也称为因变量。
-
解释变量：模型用来预测响应变量的变量，也称为独立变量。
-
简单回归：一个包含一个响应变量和一个解释变量的回归。
-
多元回归：一个包含多个解释变量但只有一个响应变量的回归。
-
系数：在回归模型中，系数是解释变量的截距和估计斜率。
-
分类变量：可以有一个离散值集的变量，通常不是数值。
-
控制变量：在回归中包含的变量，用于将解释变量的直接效应与间接效应分开。
-
广义线性模型：基于解释变量和响应变量之间不同数学关系的一系列回归模型。
-
逻辑回归：当响应变量只有两个可能值时使用的广义线性模型。
练习
练习 11.1
婴儿男孩比婴儿女孩重吗？为了回答这个问题，我们将再次使用 NSFG 数据。
拟合一个线性回归模型，以totalwgt_lb作为响应变量，以babysex作为分类解释变量——该变量的值为1表示男孩，2表示女孩。估计的体重差异是多少？这是否具有统计学意义？如果你控制了母亲的年龄——母亲的年龄是否解释了部分或全部的明显差异？
练习 11.2
特里弗斯-威尔拉德假说认为，对于许多哺乳动物，性别比例取决于“母体状况”——也就是说，像母亲的年龄、体型、健康状况和社会地位这样的因素。一些研究表明人类中存在这种效应，但结果不一。
查看 en.wikipedia.org/wiki/Trivers-Willard_hypothesis。
让我们看看母亲的年龄和生男孩的概率之间是否存在关系。拟合一个逻辑回归模型，以婴儿的性别作为响应变量，母亲的年龄作为解释变量。年龄较大的母亲更有可能还是不太可能生男孩？如果你使用母亲的年龄的二次模型会怎样？
为了在逻辑回归中使用babysex作为响应变量，我们将它重新编码为男孩的值为1，女孩的值为0。
valid["y"] = (valid["babysex"] == 1).astype(int)
练习 11.3
对于阿德利企鹅，拟合一个线性回归模型，预测企鹅的体重作为flipper_length、culmen_depth和作为分类变量的Sex的函数。如果我们控制了鳍长和喙深，雄性企鹅会重多少？为一系列鳍长、雄性和雌性企鹅生成并绘制预测图，其中culmen_depth设置为平均值。
练习 11.4
让我们看看帝企鹅在数据集中是否比其他企鹅物种的异质性更多或更少，这是通过模型的伪 (R²) 值来衡量的。使用 get_species 函数来选择帝企鹅，然后使用逻辑回归来拟合一个以性别作为响应变量，所有四个测量值作为解释变量的逻辑回归模型。伪 (R²) 值与其他模型相比如何？
chinstrap = get_species(penguins, "Chinstrap")
len(chinstrap)
68
《Think Stats：Python 中的探索性数据分析，第 3 版》
版权所有 2024 艾伦·B·唐尼
代码许可：MIT 许可证
浙公网安备 33010602011771号