1.引言

在统计学与机器学习的广阔领域中，线性回归作为一种基础而强大的预测技术，其核心在于通过输入变量（或称预测器、自变量）来估计输出变量（响应变量、因变量）的连续值。本章聚焦于线性回归的一个基本但意义深远的实例——单变量线性回归，即仅使用一个连续型自变量x来预测另一个连续型因变量y。

具体而言，我们构建了一个线性模型，其形式为$y_i=a+b{x}_i+\epsilon$，其中$y_i$代表第i个观测值的因变量，$\xi$是自变量的相应观测值，$a$和$b$是模型的参数（分别称为截距和斜率），而$\xi$代表误差项，用于捕捉模型无法完全解释的变异性。我们的目标是基于给定的数据集$(x_i,y_i)$（其中$i$从1到$n$，$n$为样本数量），找到最佳的$a‘$和$b$值，使得模型能够最准确地预测或拟合数据。

为了深入理解这一过程，我们将通过一个实际案例来逐步展示：从数据准备、模型拟合、结果可视化（包括数据点和拟合线的展示），到模型解释（解读斜率和截距的含义）。此外，为了验证模型的稳健性和拟合效果，我们还将采用模拟数据的方法，通过人为生成的数据集来重复模型拟合过程，以此检查拟合流程的有效性和准确性。

尤为值得一提的是，单变量线性回归不仅是复杂统计模型的基础，而且它还隐含了简单比较的逻辑。在某些特定情境下，当斜率$b$接近或等于某一固定值时（如0或1），线性回归模型就退化为一种更为直观的比较分析，从而展示了其作为数据分析工具箱中通用而灵活的工具之一的重要地位。

2.单变量线性回归的示例

本节以经济增长与美国总体选举的关系展示单变量回归模型的基本特性。

在探讨总统选举中选票份额与经济表现之间关系的统计模型时，我们聚焦于由著名政治学家道格拉斯·希布斯（Douglas Hibbs）提出的“面包与和平”理论框架。这一理论简约而有力，其核心在于通过经济增长这一单一变量来预测选举结果，并巧妙地考虑了战争时期的特殊情况，如1952年史蒂文森（Stevenson）因战争而遭遇的选举挫折，以及1968年汉弗莱（Humphrey）在越南战争阴影下的失利。尽管通过整合更多信息（如现任总统的执政状态、民意调查结果等）能够进一步提升预测精度，但这一基础模型已展现出其非凡的预测能力。

为了具体实现这一模型，我们采用了线性回归的方法，并将结果可视化于图2.2中。此图直观地描绘了经济增长（自变量x）与现任政党选票份额（因变量y）之间的关系，通过一条拟合线展现了两者之间的统计关联。数据来源于精心整理的$hibbs.dat$文件，其中包含了历次选举中现任政党在两党制投票中的得票百分比（$vote$）以及前几年的平均个人收入增长（$growth$），作为分析的关键指标。

在R环境中，我们首先读取并绘制了这些数据，以便直观感受经济增长与选举结果之间的潜在联系。随后，利用$sta{n}_{g}lm$函数（来自$rstanarm$包，一个结合了贝叶斯统计与广义线性模型的强大工具）拟合了线性模型$M1$，其中$vote$被建模为$growth$的函数。通过执行$print(M1)$，我们在R控制台中获得了模型的详细输出，这些输出包括参数估计值（如斜率和截距，它们分别量化了经济增长对选票份额的直接影响）、拟合优度指标（评估模型对数据解释的充分程度）、以及残差分析（揭示模型未能解释的变异性）。

深入理解这些输出结果，对于评估模型的有效性、识别潜在的数据模式以及探索经济增长与选举结果之间复杂关系的本质至关重要。尽管第8章将更深入地探讨模型的拟合过程与细节，但当前阶段的分析已为我们提供了关于这一重要社会现象的初步见解，揭示了经济增长在塑造选举结果方面不可忽视的力量。

图2.1 道格拉斯·希布斯的“面包与和平”投票与经济模型。自1952年以来的总统选举按前一届*任期末经济表现（通过通货膨胀调整后的平均个人收入增长来衡量）的顺序列出。对决以现任政党的候选人对阵其他政党的候选人列出。经济表现越好，现任政党的候选人表现得越好，最大的例外是1952年（朝鲜战争）和1968年（越南战争）。

图2.2 （a）根据经济预测选举：图7.1中的数据以散点图的形式表示，每年有一个数据点，（b）带有线性拟合的数据，$y=46.3+3.0x$。

3.利用经济预测分析总统选举投票

在探讨如何通过经济增长率来预测总统选举投票份额的过程中，我们深入分析了投票结果与经济状况之间的关联。本例中，投票份额作为结果变量，而经济增长率则作为预测变量。值得注意的是，我们的模型默认包含了一个截距项，这意味着回归线并非简单地从原点出发，而是形式为 $y=a+bx$，其中 $a$ 是截距，$b$ 是斜率，而 $x$ 代表经济增长率。

3.1.截距项的意义与调整

在某些特定情境下，我们可能希望拟合一条截距为零的直线，即 $y=bx$。这可以通过在R中使用 $sta{n}_{g}lm$ 函数时添加 $-1$ 参数来实现，如 $sta{n}_{g}lm(vote-1+growth)$。然而，在本例中，由于我们希望模型能够反映经济增长为零时现任政党可能遭受的惩罚性影响，因此保留截距项是有意义的。

3.2.模型系数解读

模型拟合后，我们得到了两个关键系数：截距和增长率的斜率，以及它们的不确定性度量（通过中位数和MAD_SD表示）和残差标准差（sigma）。具体来说，拟合的直线方程为 $y=46.3+3.0x$，这意味着：

• 当经济增长率为零时，现任政党预计将获得46.3%的选票，这通常预示着选举的失利。
• 经济增长率每增加一个百分点，现任政党的预期选票份额将增加3.0个百分点。

斜率的标准误差较小（0.7），表明这一估计具有较高的可靠性，且95%置信区间为[1.6, 4.4]，与零显著分离，进一步支持了经济增长对选举结果的重要影响。

3.3.残差标准差与模型预测能力

残差标准差（σ=3.9）反映了模型预测与实际结果之间的平均差异，约68%的观测值将落在拟合线±3.9个百分点之内。这表明模型虽然提供了有价值的预测信息，但并非完全准确，选举结果仍受多种因素影响。

3.4案例分析

回顾2008年选举，当时的经济增长率约为0.1%，根据我们的模型，预测现任政党（共和党）将获得约46.6%的选票，而奥巴马（*党）则有望获得53.4%的选票。这一预测与实际情况相符，展示了模型在特定条件下的有效性。

在2016年选举前，我们基于约2%的经济增长率进行了预测。根据模型，希拉里·克林顿的预期投票份额为52.3%。然而，这个预测值本身并不直接告诉我们希拉里获胜的概率。为了评估这一预测的不确定性，我们构建了图2.3，展示了以52.3为中心、标准差为3.9的正态分布图，从而量化了预测的不确定性。

图2.3 基于2%的经济增长率，对希拉里·克林顿在2016年两党投票中所占百分比的预测分布。曲线显示了一个以预测点 46.3 + 3.0 × 2.0 = 52.3 为中心的正态分布，标准差为3.9，这是基于模型拟合的估计不确定性。阴影区域表示基于模型的克林顿获胜的概率，即 1 - pnorm(50, 52.3, 3.9)，或者0.72。

为了更直观地解释如何利用投票份额预测来估算希拉里·克林顿在2016年普选中的获胜概率，我们可以进一步探讨这一过程的逻辑和实用性。在图7.3中，通过将预测的投票份额分布（以52.3%为中心，标准差为3.9%）与50%的获胜阈值进行比较，我们得到了一个概率值，即希拉里赢得选举的概率为72%。这一计算不仅展示了统计模型的应用，还揭示了概率分布在预测不确定性中的关键作用。

为什么预测投票份额而非直接预测获胜者？

预测投票份额而非直接预测获胜者的原因，根植于选举复杂性和结果不确定性的本质。以下是通过分析三类选举情况来阐述这一点的深度解析：

1. 势均力敌的选举：
   在势均力敌的选举中，如1960年的肯尼迪与尼克松对决，两位候选人的支持率极为接近，使得直接预测获胜者变得极为困难且不可靠。这种情况下，预测投票份额（即每位候选人可能获得的选票百分比）比单纯判断谁将获胜更加合理和科学。它提供了关于选举结果可能性的更细致、更全面的信息，有助于理解选举的紧张程度和潜在的不确定性。

2. 具有竞争性的选举：
   对于像2008年奥巴马与麦凯恩这样的选举，虽然一方被普遍认为是更有可能的获胜者，但另一方仍保留着相当程度的竞争力。在这种情境下，预测投票差距（即两位候选人之间的选票百分比差异）以及任何一方获胜的概率，比简单地宣布某位候选人将获胜更为准确和有用。这不仅有助于评估选举的紧张程度，还能为选民、媒体和政策制定者提供关于选举结果可能性的重要参考。

3. 概率视角的优越性：
   通过预测投票份额并计算获胜概率，我们能够以更加量化和精确的方式理解选举结果的不确定性。这种方法不仅有助于减少主观臆断和偏见的影响，还能为决策制定提供更加科学、客观的依据。此外，概率分布还允许我们考虑极端情况的可能性（尽管它们可能较为罕见），从而进一步完善我们对选举结果的理解。

预测投票份额而非直接预测获胜者，是出于对选举复杂性和结果不确定性的充分考虑。通过预测投票份额并计算获胜概率，我们能够以更加全面、细致和量化的方式理解选举结果的可能性，为选民、媒体和政策制定者提供更加准确、有用的信息。这种方法不仅提高了预测的准确性和可靠性，还促进了对选举结果更深入、更科学的理解和分析。

4.利用模拟数据检验模型拟合过程

虽然前一个例子简单到可以通过绘图直观判断直线是否穿过数据点，但更普遍地说，在已知真实情况的受控条件下进行拟合并检查拟合效果是一种良好的实践。本文将通过选举模型来演示这一过程。

步骤1：构建虚拟世界
我们首先为模型中的所有参数设定真实的值。鉴于我们已经对数据进行了模型拟合，我们将假设这些特定的参数值是真实的。即，我们假设公式 $y=a+bx+\text{error}$ 正确无误，其中误差来自均值为0、标准差为 $\sigma$ 的正态分布。然后，我们使用数据集中已有的预测变量 $x$ 来检查这些预测变量是否能够复现与我们观察到的 $y$ 一致的分布。

a <- 46.3
b <- 3.0
sigma <- 3.9
x <- hibbs$growth
n <- length(x)

步骤2：生成模拟数据
接下来，我们生成一个假数据向量 $y$ 并将其放入数据框中：

y <- a + b*x + rnorm(n, 0, sigma)
fake <- data.frame(x, y)

步骤3：拟合模型并比较拟合值与预设值
然后，我们对这些数据执行回归分析。在拟合过程中，不使用预设的真实值 $\alpha$、$\beta$ 和 $\sigma$。

it <- stan_glm(y ~ x, data=fake)
print(fit)

回归分析的输出结果如下：

Median MAD_SD
(Intercept) 44.4 1.7
x 3.2 0.7
Auxiliary parameter(s):
Median MAD_SD
sigma 4.0 0.8

将估计得到的系数与预设的真实值46.3和3.0进行比较，可以看出拟合结果是合理的：虽然估计值并不完全精确，但它们落在了可接受的误差范围内。

为了更正式地进行比较，我们可以从回归对象中提取系数的估计值和标准误差。为了编程简便，这里我们仅提取斜率 $b$ 的值：

b_hat <- coef(fit)["x"]
b_se <- se(fit)["x"]

接着，我们检查 $b$ 的真实值是否落在通过加减一个或两个标准误差得到的68%和95%置信区间内：

cover_68 <- abs(b - b_hat) < b_se
cover_95 <- abs(b - b_hat) < 2*b_se
cat(paste("68% coverage: ", cover_68, "\n"))
cat(paste("95% coverage: ", cover_95, "\n"))

步骤4：将模拟过程嵌入循环中
我们已经确认了置信区间在单次模拟中的有效性，但它们是否具有正确的覆盖概率——也就是说，这些区间是否如预期那样按照宣称的百分比包含了真实值？为了验证这一点，我们将模拟过程、模型拟合和覆盖率检查嵌入一个循环中，并重复执行1000次。这样的循环也可以使用R语言中的 replicate 函数来隐式执行，如文档中第72页和74页所示。这里我们选择直接编写循环代码：

n_fake <- 1000
cover_68 <- rep(NA, n_fake)
cover_95 <- rep(NA, n_fake)
for (s in 1:n_fake){
y <- a + b*x + rnorm(n, 0, sigma)
fake <- data.frame(x, y)
fit <- stan_glm(y ~ x, data=fake, refresh=0) # suppress output on console
b_hat <- coef(fit)["x"]
b_se <- se(fit)["x"]
cover_68[s] <- abs(b - b_hat) < b_se
cover_95[s] <- abs(b - b_hat) < 2*b_se
}
cat(paste("68% coverage: ", mean(cover_68), "\n"))
cat(paste("95% coverage: ", mean(cover_95), "\n"))

这个过程需要几分钟来运行，完成后在控制台上将显示如下内容：

68% coverage: 0.628
95% coverage: 0.928

这表明平均而言，63%的68%置信区间和93%的95%置信区间包含了真实参数值，这与理论值0.68和0.95相差不远。覆盖率稍低，部分是因为我们使用的是正态分布的标准误差界限，而实际上，鉴于样本量只有16，我们应该使用具有14个*度的t分布来进行更准确的推断。为了谨慎起见，我们可以重新检查适当t区间的覆盖率：

n_fake <- 1000
cover_68 <- rep(NA, n_fake)
cover_95 <- rep(NA, n_fake)
t_68 <- qt(0.84, n - 2)
t_95 <- qt(0.975, n - 2)
for (s in 1:n_fake){
y <- a + b*x + rnorm(n, 0, sigma)
fake <- data.frame(x, y)
fit <- stan_glm(y ~ x, data=fake, refresh=0)
b_hat <- coef(fit)["x"]
b_se <- se(fit)["x"]
cover_68[s] <- abs(b - b_hat) < t_68 * b_se
cover_95[s] <- abs(b - b_hat) < t_95 * b_se
}
cat(paste("68% coverage: ", mean(cover_68), "\n"))
cat(paste("95% coverage: ", mean(cover_95), "\n"))

这次模拟给出了预期的结果，即大约68%的68%置信区间和大约95%的95%置信区间包含了真实的参数值。

5. 将比较转化为回归模型

为了与基本统计方法相衔接，本文将展示如何将简单的平均值和比较操作视为线性回归的特殊情况。这种更为通用的公式化方法不仅有助于在概念上统一这些看起来不同的概念，而且在实践中，它将成为构建更复杂模型的一个关键部分，这将在后续章节中展示。

要将比较转化为回归模型，我们需要引入指标变量的概念。指标变量是一种特殊的预测变量，它可以取值为1或0，用于指示数据点是否属于某个特定类别。例如，我们可以为性别创建一个指标变量，其中男性为1，女性为0；或者为种族创建多个指标变量，如“白人”、“黑人”和“西班牙裔”，每个指标变量对于属于特定种族的个体取值为1，否则为0。

为了简化说明，我们将使用模拟数据来展示比较操作与指标回归之间的联系。

5.1.估计均值等同于对常数项进行回归

示例：假设我们从一个均值为2.0、标准差为5.0的总体中模拟20个观测值：

n_0 <- 20
y_0 <- rnorm(n_0, mean = 2.0, sd = 5.0)
fake_0 <- data.frame(y_0)
print(y_0)

在一次特定的模拟中，生成的数字可能如下：

-0.3, 4.1, -4.9, 3.3, 6.4, 7.2, 10.7, -4.6, 4.7, 6.0, 1.1, -6.7, 10.2,
9.7, 5.6, 1.7, 1.3, 6.2, -2.1, 6.5

将这些数值视为一个随机样本，我们可以通过 mean(y_0) 来估计总体均值，其标准误差为 sd(y_0)/sqrt(n_0)。例如，我们得到的估计值可能是3.3，标准误差为1.1。

使用最小二乘回归对常数项进行拟合，我们可以得到相同的结果：

fit_0 <- stan_glm(y_0 ~ 1, data = fake_0, 
                 prior_intercept = NULL, prior = NULL, prior_aux = NULL)
print(fit_0)

这将产生以下输出：

Median MAD_SD
(Intercept) 3.3 1.1
Auxiliary parameter(s):
Median MAD_SD
sigma 5.3 0.9

上述设置分配了一个平坦的先验，以便 stan_glm 能够产生传统的最小二乘回归估计。在后续的讨论中，我们将更详细地讨论先验分布；在这里，我们使用简单的最小二乘估计来证明简单平均值和只有截距的回归之间的等价性。

5.2.估计差异等同于对指标变量进行回归

接下来，我们添加一个新的组别：来自均值为8.0、标准差为5.0的总体的30个观测值：

n_1 <- 30
y_1 <- rnorm(n_1, mean = 8.0, sd = 5.0)

我们可以直接比较每组的平均值，并计算相应的标准误差：

diff <- mean(y_1) - mean(y_0)
se_0 <- sd(y_0) / sqrt(n_0)
se_1 <- sd(y_1) / sqrt(n_1)
se <- sqrt(se_0^2 + se_1^2)

在特定的模拟中，这些计算返回的差异值为4.2，其标准误差为1.3，这与模拟中构建的总体差异6.0一致。

或者，我们可以将问题构建为回归问题，通过将数据合并到一个单独的向量 $y‘$中，并创建一个指标变量 $x$：

x_i <- ifelse(observation_i 在组 0, 0, 1)

在R中，我们可以这样编写：

n <- n_0 + n_1
y <- c(y_0, y_1)
x <- c(rep(0, n_0), rep(1, n_1))
fake <- data.frame(x, y)
fit <- stan_glm(y ~ x
上一篇：漏洞挖掘基本概念及步骤介绍

下一篇：基于上一篇博客，用阻塞队列实现异步下单-引言