04_一元线性回归
第4章 一元线性回归
4.1 一元线性回归模型
从总体中抽取n个样本，一元线性回归模型如下：
\[y_i=\alpha + \beta x_i + \epsilon_i \quad(i=1，\dots，n)
\]
- \(n\)：样本容量
- \(\alpha + \beta x_i\)：总体回归线(population regression line) / 总体回归函数(population regression function，PRF)
模型本身也被称为数据生产过程(Data Generation Process，DGP)
4.2 OLS估计量的推导
任务
根据观测值 \(\{x_i,y_i \}_{i=1}^n\) 来估计总体回归线（PRL）。
核心思路
- 找到一条直线，让它离所有观测值最近。
- 给定一条直线：\(\hat y_i=\hat\alpha + \hat\beta x_i\)
- 计算观测值与直线的距离，：\(e_i \equiv y_i - \hat\alpha - \hat\beta x_i\)
- 观测值与估计量的差，就是残差（Residual）
- 当残差平方和最小时，直线离所有观测值最近。
实现方法
- OLS最小二乘法(ordinary least squares)
- 选择 \(\hat\alpha\) 和 \(\hat\beta\)，使残差平方和最小。
推导过程
- 目标函数：
\[\min_{\hat \alpha,\hat \beta}\sum_{i=1}^ne_i^2=\sum_{i=1}^n(y_i-\hat \alpha-\hat \beta x_i)^2
\]
- 最小的条件是一阶偏导均为0：
\[\begin{cases}
\cfrac{\partial }{\partial \hat \alpha}\sum_{i=1}^ne_i^2=-2\sum_{i=1}^n(y_i-\hat \alpha-\hat \beta x_i)=0\\
\cfrac{\partial }{\partial \hat \beta}\sum_{i=1}^ne_i^2=-2\sum_{i=1}^n(y_i-\hat \alpha-\hat \beta x_i)x_i=0
\end{cases}\]
- 得到正规方程组（normal equation）：
\[\begin{cases}
n\hat \alpha+\hat \beta\sum_{i=1}^n x_i =\sum_{i=1}^n y_i\\
\hat \alpha\sum_{i=1}^n x_i+\hat \beta\sum_{i=1}^n x_i^2=\sum_{i=1}^n x_i y_i
\end{cases}\]
- 求解上述二元一次方程组 （\(\overline x \equiv \frac{1}{n}\sum_{i=1}^n x_i\)，\(\overline y \equiv \frac{1}{n}\sum_{i=1}^n y_i\)）：
\[\begin{cases}
\hat \alpha = \overline y - \hat \beta \overline x \\
\hat \beta=\frac{\sum_{i=1}^n( x_i - \overline x )( y_i - \overline y)}{\sum_{i=1}^n ( x_i - \overline x )^2}
\end{cases}\]
- \(x_i\) 不能是常数项，否则分母为零。
- 残差平方和：sum of squared residuals（SSR）/ residuals squared sum（RSS）
- 普通最小二乘法：ordinary least squares（OLS）
- \(\hat y_i=\hat\alpha + \hat\beta x_i\) ：
- 样本回归线 sample regression line（SRL）
- 样本回归函数 sample regression function （SRF）
4.3 OLS的正交性
OLS残差（\(e\)）与解释变量（\(x\)、\(1\)）及拟合值（\(\hat y\)）是正交的。
- \(e_i \equiv y_i - \hat y_i\)
- 常数项可视为取值为1的解释变量，向量表示是 \(1\) 。
两个矩阵正交：\(x' y = 0\)
性质：
\[1'e = \sum_{i=1}^ne_i =0
\]
- 上式\(e_i\)展开后，求和并除于n，可以推导出，$$\overline y = \overline {\hat y}$$
\[x'e = \sum_{i=1}^nx_ie_i =0
\]
\[\hat y'e = \sum_{i=1}^n\hat y_ie_i =0
\]
OLS估计量是最佳无偏估计的前提就是OLS的正交性
4.4 平方和分解公式
若回归方程有常数项，则有：\(TSS = ESS + RSS\)
\[\sum_{i=1}^n(y_i-\overline y)^2=\sum_{i=1}^n(\hat y_i-\overline y)^2+\sum_{i=1}^ne_i^2
\]
若回归方程无常数项，上式不成立。
TSS：总平方和 \(\sum_{i=1}^n(y_i-\overline y)^2\)
ESS：可解释平方和 \(\sum_{i=1}^n(\hat y_i-\overline y)^2\)
RSS：残差平方和 \(\sum_{i=1}^ne_i^2\)
4.5 拟合优度
- 目的：
- 用来衡量样本回归线离所有样本点究竟有多近？
- 如果模型的可解释平方和（ESS）占总平方和（TSS）的比重越大，则拟合程度越高
- 定义：
- 拟合优度（goodness of fit）\(R^2\)：$$R^2=\frac{ESS}{TSS}=1-\frac{RSS}{TSS}$$ ^4adc1c
- 拟合优度也称可决系数（conefficient of determination）
title:注意
$R^2$只是反映拟合程度的好坏，除此无太多意义。
评估回归方程是否显著，仍应使用F检验。
4.6 无常数项的回归
偶尔会用到无常数项的回归。
无常数项回归也称为过原点回归（regression through the origin）
- 其他性质与有常数项OLS一致
- 但不适宜用\(R^2\)来度量拟合优度
- 因为没有常数项，所以使用非中心\(R^2\)来判断拟合优度。
\[R_{uc}^2=\frac{\sum_{i=1}^n\hat y_i^2}{\sum_{i=1}^n y_i^2}
\]
\(R^2\)和\(R_{uc}^2\)定义不同，无可比性，但在大多数统计软件中都被称为‘R-squared’。
4.7 一元线性回归的python命令及实例
使用statsmodels进行计量回归的基本步骤：
- 导入所需的库
- 读取数据->grilic.dta
- 动机与模型
- 定义变量
- 模型拟合
- 检验
import pandas as pd
import statsmodels.api as sm
# OLS模型使用 sm.OLS()构建
# 读取数据
grilic = pd.read_stata('../2_Data/Data-2e/grilic.dta')
grilic.head()
# 定义变量
y = grilic['lnw']
X = grilic[['s']]
X_const = sm.add_constant(X) # 增加常数项
# 含常数项的一元线性回归模型
model_const = sm.OLS(y,X_const)
grilic_ols_const = model_const.fit()
print(grilic_ols_const.summary())
# 不含常数项的一元线性回归模型
model = sm.OLS(y,X)
grilic_ols = model.fit()
print(grilic_ols.summary())
# 绘图可以更直观的看到样本回归线
# 使用statsmodels 的 plot_partregress 函数绘制样本回归线
sm.graphics.plot_partregress(endog='lnw',
exog_i='s',
exog_others=[ ],
data=grilic,
obs_labels=False # 不显示样本点的标签
)
# 或者使用seaborn画不同主题的图案
import seaborn as sns
sns.set_theme(color_codes=True)
sns.regplot(x=grilic['s'], y = grilic['lnw'])
statsmodel官方文档中提到了一种 Rainbow test的检验方法，用于检验线性模型的显著性，是对模型本身的建议。
- 原假设：使用全部样本拟合的模型和使用样本的中心子集拟合的模型是一样的
- 备择假设：是不一样的。
彩虹检验对许多不同形式的非线性具有检测能力。
[statsmodels.stats.diagnostic.linear_rainbow - statsmodels 0.14.1](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.linear_rainbow.html)
4.8 Stata命令运行结果的存储与调用
略
4.9 总体回归函数与样本回归函数：蒙特卡洛模拟
蒙特卡洛模拟(Monte Carlo Method)：通过计算机模拟，从总体抽取大量随机样本的计算方法。
对这样一个数据生成过程DGP或总体回归函数PRF：$$y_i = 1+2x_i+\epsilon_i \quad (i =1,\dots,20)$$
- 其中：
- 解释变量 \(x_i \sim N(3, 4^2)\)
- 扰动项 \(\epsilon \sim N(0, 9^2)\)
随机抽样60个样本，并进行回归，看看SRF和PRF的区别。
# 导入必要库
import numpy as np
import statsmodels.api as sm
import seaborn as sns
# 设置随机种子以获得可重现的结果（可选）
np.random.seed(62223)
# 生成60个x和e的样本
x = np.random.normal(3, 4, 60)
e = np.random.normal(0, 9, 60)
y =1 + 2*x + e
# 设置常数项
X_const = sm.add_constant(x)
model = sm.OLS(y, X_const)
results = model.fit()
print(results.summary())
# 绘图
sns.set_theme(color_codes=True)
sns.regplot(x = x, y = y, color = 'b', label='SRL')
sns.lineplot(x = x, y = 1 + 2*x, color = 'r',label='PRL')
习题
4.7 代码实现
import pandas as pd
import statsmodels.api as sm
import seaborn as sns
# 读取数据
galton = pd.read_stata('../2_Data/Data-2e/galton.dta')
galton.head(30)
# 计算child列和parent列的基本统计特征
child_stats = [galton['child'].mean(),
galton['child'].std(),
galton['child'].min(),
galton['child'].max()]
parent_stats = [galton['parent'].mean(),
galton['parent'].std(),
galton['parent'].min(),
galton['parent'].max()]
print('child_stats: 平均值为{}, 标准差为{}, 最小值为{}, 最大值为{}'.format(*child_stats))
print('parent_stats: 平均值为{}, 标准差为{}, 最小值为{}, 最大值为{}'.format(*parent_stats))
# 可能理解错了，应该是统计回归的特征吧
# 模型:child = \alpha + \beta * parent + \epsilon
# 定义变量
x = galton['parent']
y = galton['child']
x = sm.add_constant(x)
print('OLS统计回归结果:')
res = sm.OLS(y, x).fit().summary()
print(res)
# 绘制散点图与线性拟合图画在一起
sns.set_theme(color_codes=True)
sns.regplot(x=galton['parent'], y=galton['child'], ci=None)
（1）结果
child_stats:
平均值为68.0884780883789,
标准差为2.5179409980773926,
最小值为61.70000076293945,
最大值为73.69999694824219
parent_stats:
平均值为68.30818939208984,
标准差为1.787333369255066,
最小值为64.0,
最大值为73.0
![[习题4.7-1.png]]
（2）结果
![[习题4.7-2.png]]
(5)
# (5)定义父母身高与父母平均身高差为parent_dev变量，定义子女身高与父母身高差为gengap变量，将gengap对parent_dev进行回归分析，并绘制回归曲线。
galton_dev = pd.DataFrame()
galton_dev['gengap'] = galton['child'] - galton['parent']
galton_dev['parent_dev'] = galton['parent'] - galton['parent'].mean()
y_dev = galton_dev['gengap']
x_dev = galton_dev['parent_dev']
x_dev = sm.add_constant(x_dev)
res = sm.OLS(y_dev, x_dev).fit().summary()
print(res)
sns.set_theme(color_codes=True)
dev = sns.regplot(x=galton_dev['parent_dev'],
y=galton_dev['gengap'],
ci=None)
print(dev)
结果：\(\hat \alpha = -0.2197\) , \(\hat \beta = -0.3537\)，SRL是条斜率为负的直线。![[习题4.7-5.png]]
意味着：高个子父母的子女的平均代差更多的出现负值，也就是子女比父母矮。
