目录
回归分析
当变量之间存在相互依存的关系时，还可以进行回归分析(regression)。
回归分析的研究领域有：线性回归、非线形回归、定性自变量回归、离散因变量回归等。
线性回归
基本原理
相关分析中的两个变量之间的地位时对等的，即变量A与变量B相关等价于变量B与变量A相关，相关分析的两个变量均为随机变量；而回归分析中要确定自变量和因变量，通常情况下只有因变量时随机变量，可以利用回归分析来对研究对象进行预测与控制。
回归分析往往时通过一条拟合的曲线来表示模型的建立。以线性回归为例，设 \(y\) 表示因变量，\(x\) 表示自变量，则有如下线形回归模型
\[y = \alpha+\beta x+ \varepsilon
\]
其中，\(\alpha,\beta\) 是回归模型的参数，称为回归系数（\(\alpha\) 也可称为截距项），\(\varepsilon\) 是随机误差项或随机扰动项，反映了除 \(x,y\)之间线性关系之外的随机因素或不可观测的因素。
在回归分析只能够，对 \(\varepsilon\) 有如下最为常见的基本经典假设：
- \(\varepsilon\) 的期望值为0；
- \(\varepsilon\) 对于所有 \(x\) 而言具有同方差性；
- \(\varepsilon\) 是服从正态分布且相互独立的随机变量
如果存在多个自变量，回归模型可以写作
\[y = \alpha+\beta_1 x_1+\beta_2 x_2+\cdots+\beta_i x_i+ \varepsilon
\]
因为上述直线中含有随机误差项，所以回归模型反映的直线是不确定的。回归分析的主要目的就是要从这些不确定的直线中，找出一条最能够代表数据原始信息的直线，来描述因变量和自变量之间的关系，这条直线称之为回归方程。
如只有一个自变量的情况下，可对模型左右两边取 \(x\) 的条件期望并根据 \(\varepsilon\) 的经典假设，可得
\[E(y|x) = \alpha+\beta x+0
\]
然后通过一定的参数估计方法，可得到估计的直线方程如下
\[\hat{y} = \hat{\alpha}+\hat{\beta}x
\]
因为回归模型的参数往往是未知的，只能依靠样本数据去进行参数估计，所以可用 \(\hat{\alpha},\hat{\beta}\) 分别表示回归模型中 \(\alpha,\beta\) 的参数估计值。方程 \(\hat{y} = \hat{\alpha}+\hat{\beta}x\) 也可称为估计的回归方程。
\(\hat{\alpha}\) 是估计的回归直线在 \(y\) 轴上的截距，\(\hat{\beta}\) 是直线的斜率，表示自变量 \(x\) 每变动一个单位时，因变量 \(y\) 的平均变动值，\(\hat{y}\) 是 \(y\) 的估计值。
参数估计的普通最小二乘法
把某个因变量的观测值记为 \(y_i\)，其期望值为 \(\bar{y}\)，其估计值(即在回归直线上的值)记为 \(\hat{y}\)。因变量的离差为 \(y_i-\bar{y}\)，可以把离差分解为 \(y_i-\hat{y}_i\)（残差）和 \(\hat{y}_i-\bar{y}\)（回归离差），即 \(y_i-\bar{y}=(y_i-\hat{y}_i)+(\hat{y}_i-\bar{y})\)。对于因变量的观测值而言，为了避免正负号影响，对3个差值求平方后，发现（可证明的）
\[\sum{(y_i-\bar{y})^2}=\sum{(y_i-\hat{y}_i)^2}+ \sum{(\hat{y}_i-\bar{y})^2}
\]
即，总离差平方和（SST）等于残差平方和（SSE）与回归离差平方和（SSR）。
如果残差越小，\(y_i\) 就越往回归直线靠近。对于所有的因变量而言，残差平方和越小，观测值就越往回归直线靠近。当残差平方和（SSE）达到极小值是，估计处理啊的回归直线能在最大程度上代表原始数值的信息。
\[\sum{(y_i-\hat{y}_i)^2}=\sum{[y_i-(\hat{\alpha}+\hat{\beta}x)]^2} \Rightarrow min
\]
当如上时，可以利用微分求极值的方法，分别求 \(\hat{\alpha},\hat{\beta}\)的偏微粉，并使之同时为0，然后求解联立方程组便可计算出\(\hat{\alpha},\hat{\beta}\)的具体数值。
这种参数估计的方法和过程是在随机误差项 \(\varepsilon\) 满足基本经典假设的情况下进行的，通常被称之为普通最小二乘法。普通最小二乘法同样适用于多个自变量和因变量之间的模型参数估计，在大多数的回归分析模型中，普遍默认使用普通最小二乘法对线性回归模型进行估计。
回归方程的检验及模型预测
- 回归方程的拟合优度判定
回归方程的拟合优度主要用于判定 \(\hat{y}_i\)估计 \(y\) 的可靠性问题，可用来衡量模型的解释程度。拟合优度判定是建立在模型参数估计时对总离差平方和（SST）分解的基础上的，SST可分解为残差平方和（SSE）与回归离差平方和（SSR），通常使用SSR占SST的比重来判定模型的解释能力，即
\[R^2 =\frac{\sum{(\hat{y}_i-\bar{y})^2}}{\sum{(y_i-\bar{y})^2}}=1-\frac{\sum{(y_i-\hat{y}_i)^2}}{\sum{(y_i-\bar{y})^2}}
\]
其中，\(R^2\) 表示拟合优度的判定系数或决定系数，取值范围[0,1]，越接近1，说明变量之间的相互依存关系越密切，接近于函数关系，两变量之间的相关程度越高，回归方程的拟合程度越好。
但是，\(R^2\) 的数值与自变量的数目有关，即自变量的个数越多 \(R^2\) 越大，这在一定策划高难度上削弱了 \(R^2\) 的评价能力，因此可考虑剔除自变量数目影响的修正 \(R^2\)。
- 回归方程总体显著性检验
利用最小二乘法拟合出来的回归方程，都是由样本数据进行的，需要对回归方程整体的显著性进行检验，来确定回归方程对总体进行推断是否显著。
回归方程总体显著性检验主要是检验因变量和自变量之间的线性关系是否显著。其原假设与备择假设：
\[H_0:\beta_1=\beta_2=\cdots=\beta_i=0;H_1:\beta_i不全为0
\]
对于回归方程总体显著性检验可用 \(F\) 检验并计算出用于检验判定的 \(P\) 值来进行。如果 \(P\) 值小于理论显著水平 \(\alpha\) 值，可以认为在显著水平 \(\alpha\) 条件下，回归方程总体显著。
- 回归方程系数显著性检验
如果模型的线性关系显著，还应对模型参数估计的结果即回归方程的系数及性能显著性检验，用于考察单个自变量与因变量的线性关系是否成立。其原假设与备择假设：
\[H_0:\beta_i=0(i=1,2,...,k);H_1:\beta_i \neq 0(i=1,2,...,k)
\]
回归方程系数显著性检验要求所有估计出来的回归系数分别进行检验（截距项通常不进行显著性检验），可以利用 \(t\) 检验进行。\(t\) 检验可以计算出每个回归系数所对应的 \(t\) 统计量值及其检验 \(P\) 值。如果某个系数对应的 \(P\) 值小于理论显著性水平 \(\alpha\) 值，可认为在显著性水平 \(\alpha\) 条件下，该回归系数是显著的。
有些情况下，没有任何关联的变量之间进行回归分析也可能得到显著的检验结果，会对分析过程造成不良影响。因此，在进行回归分析之前，必须考虑好变量之间的关系及其所代表的经济含义。
- 回归方程的预测
回归预测是一种有条件的预测，根据估计出来的回归方程，在给定自变量数值的条件下，对因变量进行预测，其预测基本公式为
\[\hat{y}_f = \hat{\alpha}+\hat{\beta}x_f
\]
其中，\(x_f\) 是另外给定的自变量的值，\(\hat{y}_f\) 为根据回归方程计算出来的预测值。
一元线性回归
模型
\[y = \alpha + \beta x + \varepsilon
\]
基本步骤
1.根据变量之间的关系，判断其是否有线性关系
2.若是线性关系，利用OLS方法或其他方法进行回归模型的参数估计
3.根据参数估计的结果进行检验
检验过程中：
1.拟合优度判定：若是拟合优度的判定系数较小，说明建立的回归方程较差，可能还有其他重要因素没有加入到模型中，考虑增加有重要影响的自变量；
2.总体显著性检验：如果不显著，说明变量之间的线性关系不明显，不适合做线性回归；
3.系数显著性检验：在拟合优度判定系数比较高、方差总体显著的清咖滚下，对回归系数进行己拿烟那，通过显著性检验的回归系数才对因变量有解释能力
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 研究文盲率与犯罪率之间关系，收集不同区域的信息
murder = pd.read_csv(\'./data/murder.csv\')
print(murder.head(5))
# region illiteracy murder
# 0 4 0.5 11.5
# 1 8 0.5 2.3
# 2 8 0.5 1.7
# 3 4 0.6 5.3
# 4 4 0.6 5.0
# 数据预处理：值标签
murder[\'region\'] = murder[\'region\'].astype(\'category\')
murder[\'region\'].cat.categories = [\'East North Central\', \'East South Central\', \'Middle Atlantic\', \'Mountain\',
\'New England\', \'Pacific\', \'South Atlantic\', \'West North Central\',
\'West South Central\']
murder[\'region\'].cat.set_categories = [\'East North Central\', \'East South Central\', \'Middle Atlantic\', \'Mountain\',
\'New England\', \'Pacific\', \'South Atlantic\', \'West North Central\',
\'West South Central\']
# 1.线性关系分析
# 散点图
plt.scatter(murder[\'illiteracy\'], murder[\'murder\'], c=\'blue\', alpha=0.75)
plt.xlabel(\'illiteracy\')
plt.ylabel(\'murder\')
plt.show()
# 从图上可知一定程度上呈线性关系
# 2.一元回归分析
# 方法一：stats.linregress
murder_model1 = stats.linregress(murder[\'illiteracy\'], murder[\'murder\']) # 自变量在前，因变量在后
print(murder_model1) # LinregressResult(slope=4.257456742653117,
# intercept=2.396775611095853, rvalue=0.7029751986841699, pvalue=1.2579116392501817e-08, stderr=0.6217132752343032)
print(\'R^2\', murder_model1.rvalue ** 2) # 0.49417412996504817
# 拟合优度判定系数R^2=0.494，拟合度不高，但是本例中勉强可接受；
# 没有返回用于回归方程总体显著性检验的F统计量及其显著性水平；
# 截矩项的参数估计值为2.396，通常不对截距项做显著性检验;
# 自变量illiteracy对应的回归系数为4.257，其对应的p值<0.0001,在给定的显著水平条件下非常显著；
# 方法二：ols
from statsmodels.formula.api import ols
formula = \'murder~illiteracy\'
murder_model2 = ols(formula, data=murder).fit()
res = murder_model2.summary2()
print(res) # 与下面的结果一样
# 类的方式调用
x = murder[\'illiteracy\']
x = sm.add_constant(x) # 构建模型的自变量数组，包含截距项，若无此行表明无截距项
y = murder[\'murder\']
murder_model3 = sm.OLS(y, x).fit()
res = murder_model3.summary2()
print(res)
# Results: Ordinary least squares
# =================================================================
# Model: OLS Adj. R-squared: 0.484
# Dependent Variable: murder AIC: 241.4099
# Date: 2020-06-08 16:10 BIC: 245.2340
# No. Observations: 50 Log-Likelihood: -118.70
# Df Model: 1 F-statistic: 46.89
# Df Residuals: 48 Prob (F-statistic): 1.26e-08
# R-squared: 0.494 Scale: 7.0367
# -------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------
# Intercept 2.3968 0.8184 2.9285 0.0052 0.7512 4.0424
# illiteracy 4.2575 0.6217 6.8479 0.0000 3.0074 5.5075
# -----------------------------------------------------------------
# Omnibus: 0.452 Durbin-Watson: 2.000
# Prob(Omnibus): 0.798 Jarque-Bera (JB): 0.452
# Skew: 0.210 Prob(JB): 0.798
# Kurtosis: 2.797 Condition No.: 4
# =================================================================
# 给出了拟合优度系数R^2=0.484
# 回归方程总体显著性检验：F=46.89,P_v=1.26e-08，总体上非常显著
# 回归方程系数显著性检验：t=6.8479,P_v=0,系数非常显著
# 截距为2.3968,系数为4.2575
# statsmodels中提供了Rainbow检验方法可以处理回归方程的总体显著性检验
# 该检验方法原假设为：回归建模被正确的进行线性建模
res = sm.stats.diagnostic.linear_rainbow(murder_model2)
print(res) # (0.6152086826947939, 0.8811091720318103)
# P_v=0.88，不能决绝原假设，可以认为线性建模是正确的
# 3.残差分析
# 根据回归方程估计出来的残差项应当符合经典假设。
# 通常可用残差图来判定残差是否与变量有关，用P-P图或Q-Q图来判定残差是否符合正态分布
# a.判定残差是否符合正态分布
sm.qqplot(murder_model2.resid, fit=True, line=\'45\')
plt.show()
# 类ProbPlot的实例对象也可绘制P-P图或Q-Q图
sm.ProbPlot(murder_model2.resid, stats.t, fit=True).ppplot(line=\'45\')
sm.ProbPlot(murder_model2.resid, stats.t, fit=True).qqplot(line=\'45\')
plt.show()
# P-P图或Q-Q图在残差符合正态假定条件下，散点图看起来应该项一条截距为0、斜率为1的直线。
# 本例中数据大部分散落在45度线上，表明残差基本符合正态分布
# b.判定残差是否与变量有关
# 方法一：手动
fig = plt.gcf()
fig.set_size_inches(8, 4)
plt.subplot(121)
plt.plot(murder[\'illiteracy\'], murder_model2.resid, \'o\')
plt.xlabel(\'illiteracy\')
plt.ylabel(\'residual\')
plt.subplot(122)
plt.plot(murder_model2.fittedvalues, murder_model2.resid, \'o\')
plt.xlabel(\'predicted_value\')
plt.ylabel(\'residual\')
plt.show()
# 自变量与残差之间的关系不明显，基本无关，符合\varepsilon对于所有x而言具有同方差性的假定；
# 残差大体均匀地分布在[-6,6]之间，其均值与0接近，符合\varepsilon零均值的假定
# 方法二：函数
# 用于回归诊断的图形其他画法
fig = plt.figure(figsize=(12, 8)) # 因图形重叠会导致部分文字叠加在一起，故事先定义好图像大小
from statsmodels.graphics.regressionplots import plot_regress_exog
plot_regress_exog(murder_model2, 1, fig=fig)
# 参数1：数据，参数2：模型中的第几个自变量，参数3：表示在指定的对象上绘图
plt.show()
# 预测值与观察值对比图形
from statsmodels.graphics.regressionplots import plot_fit
plot_fit(murder_model2, 1) # 参数2：模型中第几个自变量
plt.show()
# 拟合曲线和因变量预测置信区间的图形
from statsmodels.sandbox.regression.predstd import wls_prediction_std
prstd, interval_l, interval_u = wls_prediction_std(murder_model2, alpha=0.5)
# 计算模型预测值的标准差和置信区间
fig = plt.subplots(figsize=(7, 4))
plt.plot(murder[\'illiteracy\'], murder[\'murder\'], \'o\', label=\'data\')
plt.plot(murder[\'illiteracy\'], murder_model2.fittedvalues, \'r--\', label=\'OLS\')
plt.plot(murder[\'illiteracy\'], interval_u, \'r--\')
plt.plot(murder[\'illiteracy\'], interval_l, \'r--\')
plt.legend(loc=\'best\')
plt.show()
# 至此，得回归方程
# y = 2.3968 + 4.2575*illiteracy
多元线性回归
模型
\[y = \beta_0+\beta_1 x_1+\beta_2 x_2+\cdots+\beta_n x_n+ \varepsilon
\]
\(\varepsilon\) 仍然服从零均值、相互独立且同方差服从正态分布等经典假定。
对上述模型左右两边同时取条件期望，得
\[E(y|x) = \beta_0+\beta_1 x_1+\beta_2 x_2+\cdots+\beta_n x_n+0
\]
仍然采用普通最小二乘法，可估计出参数 \(\beta_0,\beta_1,\cdots,\beta_n\) 的估计值 \(\hat{\beta}_0,\hat{\beta}_1,\cdots,\hat{\beta}_n\)，即可得到多元回归方程
\[\hat{y} = \hat{\beta}_0+\hat{\beta}_1 x_1+\hat{\beta}_2 x_2+\cdots+\hat{\beta}_n x_n
\]
对于多元回归方程的拟合度有限判定、方程总体显著性检验与回归系数显著性检验过程与一元回归方程的检验过程一样。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 对员工的当前薪酬及其影响因素进行回归分析，\alpha=0.1
salary = pd.read_csv(\'./data/salary_r.csv\')
print(salary.head(5))
# position ID Gender ... Begin_Salary Experience Age
# 0 1 1 1 ... 27000 144.0 55.0
# 1 1 34 1 ... 39990 175.0 58.0
# 2 1 18 1 ... 27510 70.0 51.0
# 3 1 200 1 ... 34980 9.0 44.0
# 4 1 199 1 ... 27480 69.0 49.0
# 数据中有较多缺失值，需做预处理。
# 使用statsmodels会自动不考虑有缺失值的行
salary = salary.dropna() # 删除有缺失值的行
# 数据预处理：值标签
salary[\'position_valuelabel\'] = salary[\'position\']
salary[\'Gender_valuelabel\'] = salary[\'Gender\']
salary[\'position_valuelabel\'] = salary[\'position_valuelabel\'].astype(\'category\')
salary[\'position_valuelabel\'].cat.categories = [\'经理\', \'主管\', \'普通员工\']
salary[\'position_valuelabel\'].cat.set_categories = [\'经理\', \'主管\', \'普通员工\']
salary[\'Gender_valuelabel\'] = salary[\'Gender_valuelabel\'].astype(\'category\')
salary[\'Gender_valuelabel\'].cat.categories = [\'女\', \'男\']
salary[\'Gender_valuelabel\'].cat.set_categories = [\'女\', \'男\']
# 性别和职位都是定性变量，对于定性自变量的回归有其特殊的处理方法，此处暂不考虑定性变量对因变量的影响
# 1.多元回归分析
from statsmodels.formula.api import ols
formula = \'Current_Salary~Education+Begin_Salary+Experience+Age\'
salary_model1 = ols(formula, data=salary).fit()
res = salary_model1.summary2()
print(res)
# Results: Ordinary least squares
# ===================================================================
# Model: OLS Adj. R-squared: 0.781
# Dependent Variable: Current_Salary AIC: 9243.8648
# Date: 2020-06-08 17:46 BIC: 9264.3664
# No. Observations: 446 Log-Likelihood: -4616.9
# Df Model: 4 F-statistic: 398.8
# Df Residuals: 441 Prob (F-statistic): 5.53e-145
# R-squared: 0.783 Scale: 5.8068e+07
# -------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------
# Intercept 703.3118 3147.8045 0.2234 0.8233 -5483.2503 6889.8740
# Education 490.5022 180.7222 2.7141 0.0069 135.3185 845.6860
# Begin_Salary 1.8936 0.0722 26.2203 0.0000 1.7517 2.0356
# Experience -11.6277 5.9764 -1.9456 0.0523 -23.3735 0.1182
# Age -74.9263 53.2311 -1.4076 0.1600 -179.5445 29.6918
# -------------------------------------------------------------------
# Omnibus: 194.312 Durbin-Watson: 1.887
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 1481.311
# Skew: 1.695 Prob(JB): 0.000
# Kurtosis: 11.259 Condition No.: 160242
# ===================================================================
# * The condition number is large (1e+05). This might indicate
# strong multicollinearity or other numerical problems.
# R^2=0.781，拟合优度系数较高，拟合程度高，解释性好；
# 总体显著性检验p \approx 0，回归方程总体显著；
# 给定\alpha=0.1，可以判定，Education,Begin_Salary,Experience的p值均<0.1，因而在模型中显著
# Age的p值>0.1, 说明年龄对当前薪酬的影响不显著
# 去掉对于回归系数不显著的变量
formula = \'Current_Salary~Education+Begin_Salary+Experience\'
salary_model2 = ols(formula, data=salary).fit()
res = salary_model2.summary2()
print(res)
# Results: Ordinary least squares
# =====================================================================
# Model: OLS Adj. R-squared: 0.780
# Dependent Variable: Current_Salary AIC: 9264.6475
# Date: 2020-06-08 17:52 BIC: 9281.0577
# No. Observations: 447 Log-Likelihood: -4628.3
# Df Model: 3 F-statistic: 529.5
# Df Residuals: 443 Prob (F-statistic): 4.72e-146
# R-squared: 0.782 Scale: 5.8206e+07
# ---------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------
# Intercept -2662.7078 2010.0497 -1.3247 0.1860 -6613.1256 1287.7100
# Education 507.1855 180.1850 2.8148 0.0051 153.0619 861.3092
# Begin_Salary 1.8942 0.0723 26.1975 0.0000 1.7521 2.0363
# Experience -18.1933 3.7141 -4.8985 0.0000 -25.4927 -10.8939
# ---------------------------------------------------------------------
# Omnibus: 192.374 Durbin-Watson: 1.873
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 1429.155
# Skew: 1.680 Prob(JB): 0.000
# Kurtosis: 11.090 Condition No.: 102573
# =====================================================================
# * The condition number is large (1e+05). This might indicate
# strong multicollinearity or other numerical problems.
# 拟合优度系数、总体显著性检验、回归系数显著性检验通过
# 但上述提示有多重共线问题或其他问题，经检查，不存在数据问题，
# 但是变量之间存在多重共线 问题，如受教育年限越高代表学历越高其起薪也会相应增高
# 解决方案参考计量经济学相关资料
# 2.对模型进行诊断
# a.判定残差是否符合正态分布
sm.qqplot(salary_model2.resid, fit=True, line=\'45\')
plt.show()
# 类ProbPlot的实例对象也可绘制P-P图或Q-Q图
sm.ProbPlot(salary_model2.resid, stats.t, fit=True).ppplot(line=\'45\')
sm.ProbPlot(salary_model2.resid, stats.t, fit=True).qqplot(line=\'45\')
plt.show()
# P-P图或Q-Q图在残差符合正态假定条件下，散点图看起来应该项一条截距为0、斜率为1的直线。
# 本例中数据大部分散落在45度线上，表明残差基本符合正态分布
# b.判定残差是否与变量有关
fig = plt.gcf()
fig.suptitle(\'Residual by Regressors & Predicted for Current_Salary\')
fig.set_size_inches(6, 6)
plt.subplot(221)
plt.plot(salary[\'Education\'], salary_model2.resid, \'o\')
plt.xlabel(\'Education\')
plt.ylabel(\'Residual\')
plt.subplot(222)
plt.plot(salary[\'Begin_Salary\'], salary_model2.resid, \'o\')
plt.xlabel(\'Begin_Salary\')
plt.subplot(223)
plt.plot(salary[\'Experience\'], salary_model2.resid, \'o\')
plt.xlabel(\'Experience\')
plt.ylabel(\'Residual\')
plt.subplot(224)
plt.plot(salary_model2.fittedvalues, salary_model2.resid, \'o\')
plt.xlabel(\'Predicted\')
plt.subplots_adjust(hspace=0.3, wspace=0.35)
plt.show()
# 自变量与残差之间的关系不明显，基本无关，符合\varepsilon对于所有x而言具有同方差性的假定；
# 残差大体均匀地分布在[-2000,2000]之间，其均值与0接近，符合\varepsilon零均值的假定
# 用于回归诊断的图形其他画法
fig = plt.figure(figsize=(12, 8)) # 因图形重叠会导致部分文字叠加在一起，故事先定义好图像大小
from statsmodels.graphics.regressionplots import plot_regress_exog
plot_regress_exog(salary_model2, 1, fig=fig)
# 参数1：数据，参数2：模型中的第几个自变量，参数3：表示在指定的对象上绘图
plt.show()
# 预测值与观察值对比图形
from statsmodels.graphics.regressionplots import plot_fit
plot_fit(salary_model2, 1) # 参数2：模型中第几个自变量
plt.show()
# 至此得回归方程
# y = -2662.7078 + 507.1855*Education + 1.8942*Begin_Salary - 18.1933*Experience
# 其他方法进行多元回归：sklearn
from sklearn import linear_model
salary_model3 = linear_model.LinearRegression()
salary_model3.fit(salary[[\'Education\', \'Begin_Salary\', \'Experience\']], salary[\'Current_Salary\'])
print(\'intercept:{0},coefficients:{1}\'.format(salary_model3.intercept_, salary_model3.coef_))
# intercept:-2708.6776308644767,coefficients:[512.4962845 1.89389266 -18.22300788]
含有定性自变量的线性回归
在影响因变量的诸多因素中，除了定量变量的影响之外，有时候还有一些定性因素的影响。定性因素对因变量的影响在进行回归分析的过程中，需要进行特殊的处理，即应当把定性变量转化为虚拟变量（或哑变量）之后再引入回归模型中进行分析。
虚拟变量的设定即是把对变量的定性描述转换成定量数据进行描述。如对性别设定为0、1代表男女。
设定虚拟变量时应当遵循如下原则：
- 对于有 \(k\) 个表现值的定性变量，只设定 \((k-1)\) 个虚拟变量；
- 虚拟变量的值通常用0、1来表示；
- 对于每个样本而言，同一个定性变量对应虚拟变量的值之和不超过1。
如性别变量，有2个表现值，则只需设定1个虚拟变量(gender)即可，0表示女，1表示男；而职位变量，有3个表现值，需设定2个虚拟变量（position1,position2）。虚拟变量中，其具体数值表示的含义，可以根据需求进行指定。
# python中虚拟变量处理
from patsy.contrasts import Treatment
contrast = Treatment(reference=3).code_without_intercept([1, 2, 3])
print(contrast)
# ContrastMatrix(array([[1., 0.], [0., 1.], [0., 0.]]), [\'[T.1]\', \'[T.2]\'])
# 将类Treatment实例化的reference参数，可以指定对照组即参考对象(所有虚拟变量均为0)的位置。
# 上式中的分类变量有3个谁能够，分别可用1，2，3来表示，参考属性设置为3
对于上式中的回归方程
\[Current\_Salary=\alpha+\beta_1\times Gender+\beta_2\times Education+\beta_1\times Postion1+\beta_4\times Postion2\\+\beta_5\times Begin\_Salary+\beta_6\times Experience +\beta_7\times Agge
\]
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 对员工的当前薪酬及其影响因素进行回归分析，\alpha=0.1
salary = pd.read_csv(\'./data/salary_r.csv\')
print(salary.head(5))
# position ID Gender ... Begin_Salary Experience Age
# 0 1 1 1 ... 27000 144.0 55.0
# 1 1 34 1 ... 39990 175.0 58.0
# 2 1 18 1 ... 27510 70.0 51.0
# 3 1 200 1 ... 34980 9.0 44.0
# 4 1 199 1 ... 27480 69.0 49.0
# 数据中有较多缺失值，需做预处理。
# 使用statsmodels会自动不考虑有缺失值的行
salary = salary.dropna() # 删除有缺失值的行
# 数据预处理：值标签
salary[\'position_valuelabel\'] = salary[\'position\']
salary[\'Gender_valuelabel\'] = salary[\'Gender\']
salary[\'position_valuelabel\'] = salary[\'position_valuelabel\'].astype(\'category\')
salary[\'position_valuelabel\'].cat.categories = [\'经理\', \'主管\', \'普通员工\']
salary[\'position_valuelabel\'].cat.set_categories = [\'经理\', \'主管\', \'普通员工\']
salary[\'Gender_valuelabel\'] = salary[\'Gender_valuelabel\'].astype(\'category\')
salary[\'Gender_valuelabel\'].cat.categories = [\'女\', \'男\']
salary[\'Gender_valuelabel\'].cat.set_categories = [\'女\', \'男\']
# 1.多元回归分析
from statsmodels.formula.api import ols
# 性别和职位都是定性变量，做虚拟变量处理
# 方法一：使用Treatment
from patsy.contrasts import Treatment
formula = \'Current_Salary~Education+Begin_Salary+Experience-Age+C(position, Treatment(reference=3))+C(Gender)\'
salary_model1 = ols(formula, data=salary).fit()
res = salary_model1.summary2()
print(res)
# Results: Ordinary least squares
# =================================================================================================
# Model: OLS Adj. R-squared: 0.813
# Dependent Variable: Current_Salary AIC: 9175.9545
# Date: 2020-06-08 21:04 BIC: 9204.6567
# No. Observations: 446 Log-Likelihood: -4581.0
# Df Model: 6 F-statistic: 323.8
# Df Residuals: 439 Prob (F-statistic): 1.02e-157
# R-squared: 0.816 Scale: 4.9647e+07
# -------------------------------------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------------------------------------
# Intercept 3908.9330 2155.2003 1.8137 0.0704 -326.8598 8144.7258
# C(position, Treatment(reference=3))[T.1] 11501.1261 1543.7087 7.4503 0.0000 8467.1481 14535.1041
# C(position, Treatment(reference=3))[T.2] 6691.6927 1692.3071 3.9542 0.0001 3365.6621 10017.7234
# C(Gender)[T.1] 1970.6471 813.6961 2.4218 0.0158 371.4232 3569.8711
# Education 506.7274 172.5026 2.9375 0.0035 167.6939 845.7609
# Begin_Salary 1.3221 0.0988 13.3870 0.0000 1.1280 1.5162
# Experience -22.5668 3.7782 -5.9728 0.0000 -29.9925 -15.1411
# -------------------------------------------------------------------------------------------------
# Omnibus: 206.434 Durbin-Watson: 1.962
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 1762.413
# Skew: 1.786 Prob(JB): 0.000
# Kurtosis: 12.060 Condition No.: 130698
# =================================================================================================
# * The condition number is large (1e+05). This might indicate strong multicollinearity
# or other numerical problems.
# 拟合优度系数、总体显著性检验、回归系数显著性检验通过
#
# 方法二：手动赋值，变量值多时不便处理
salary.loc[(salary.position_valuelabel == \'经理\'), \'position1\'] = 1
salary.loc[(salary.position_valuelabel == \'经理\'), \'position2\'] = 0
salary.loc[(salary.position_valuelabel == \'主管\'), \'position1\'] = 0
salary.loc[(salary.position_valuelabel == \'主管\'), \'position2\'] = 1
salary.loc[(salary.position_valuelabel == \'普通员工\'), \'position1\'] = 0
salary.loc[(salary.position_valuelabel == \'普通员工\'), \'position2\'] = 0
print(salary.head())
# position ID Gender ... Gender_valuelabel position1 position2
# 0 1 1 1 ... 男 1.0 0.0
# 1 1 34 1 ... 男 1.0 0.0
# 2 1 18 1 ... 男 1.0 0.0
# 3 1 200 1 ... 男 1.0 0.0
# 4 1 199 1 ... 男 1.0 0.0
# 方法三：pandas的get_dummies
dm = pd.get_dummies(salary[\'position_valuelabel\'], prefix=\'position_\') # 将指定变量处理为虚拟变量
salary = salary.join(dm) # 将虚拟变量加入原数据集中
res = salary[[\'position_valuelabel\', \'position1\', \'position2\', \'position__主管\', \'position__普通员工\']].iloc[[0, 330, 440]]
print(res)
# position_valuelabel position1 position2 position__主管 position__普通员工
# 0 经理 1.0 0.0 0 0
# 340 普通员工 0.0 0.0 0 1
# 465 主管 0.0 1.0 1 0
formula = \'Current_Salary~Education+Begin_Salary+Experience-Age+position1+position2+Gender\'
salary_model2 = ols(formula, data=salary).fit()
res = salary_model2.summary2()
print(res)
# Results: Ordinary least squares
# =====================================================================
# Model: OLS Adj. R-squared: 0.813
# Dependent Variable: Current_Salary AIC: 9175.9545
# Date: 2020-06-08 21:30 BIC: 9204.6567
# No. Observations: 446 Log-Likelihood: -4581.0
# Df Model: 6 F-statistic: 323.8
# Df Residuals: 439 Prob (F-statistic): 1.02e-157
# R-squared: 0.816 Scale: 4.9647e+07
# ---------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------
# Intercept 3908.9330 2155.2003 1.8137 0.0704 -326.8598 8144.7258
# Education 506.7274 172.5026 2.9375 0.0035 167.6939 845.7609
# Begin_Salary 1.3221 0.0988 13.3870 0.0000 1.1280 1.5162
# Experience -22.5668 3.7782 -5.9728 0.0000 -29.9925 -15.1411
# position1 11501.1261 1543.7087 7.4503 0.0000 8467.1481 14535.1041
# position2 6691.6927 1692.3071 3.9542 0.0001 3365.6621 10017.7234
# Gender 1970.6471 813.6961 2.4218 0.0158 371.4232 3569.8711
# ---------------------------------------------------------------------
# Omnibus: 206.434 Durbin-Watson: 1.962
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 1762.413
# Skew: 1.786 Prob(JB): 0.000
# Kurtosis: 12.060 Condition No.: 130698
# =====================================================================
# * The condition number is large (1e+05). This might indicate
# strong multicollinearity or other numerical problems.
# 2.对模型进行诊断
# a.判定残差是否符合正态分布
sm.qqplot(salary_model2.resid, fit=True, line=\'45\')
plt.show()
# 类ProbPlot的实例对象也可绘制P-P图或Q-Q图
sm.ProbPlot(salary_model2.resid, stats.t, fit=True).ppplot(line=\'45\')
sm.ProbPlot(salary_model2.resid, stats.t, fit=True).qqplot(line=\'45\')
plt.show()
# P-P图或Q-Q图在残差符合正态假定条件下，散点图看起来应该项一条截距为0、斜率为1的直线。
# 本例中数据大部分散落在45度线上，表明残差基本符合正态分布
# b.判定残差是否与变量有关
# 方式一：手动处理
fig = plt.gcf()
fig.suptitle(\'Residual by Regressors & Predicted for Current_Salary\')
fig.set_size_inches(6, 6)
plt.subplot(241)
plt.plot(salary[\'Education\'], salary_model2.resid, \'o\')
plt.xlabel(\'Education\')
plt.ylabel(\'Residual\')
plt.subplot(242)
plt.plot(salary[\'Begin_Salary\'], salary_model2.resid, \'o\')
plt.xlabel(\'Begin_Salary\')
plt.subplot(243)
plt.plot(salary[\'Experience\'], salary_model2.resid, \'o\')
plt.xlabel(\'Experience\')
plt.subplot(244)
plt.plot(salary[\'position1\'], salary_model2.resid, \'o\')
plt.xlabel(\'position1\')
plt.subplot(245)
plt.plot(salary[\'position2\'], salary_model2.resid, \'o\')
plt.xlabel(\'position2\')
plt.ylabel(\'Residual\')
plt.subplot(246)
plt.plot(salary[\'Gender\'], salary_model2.resid, \'o\')
plt.xlabel(\'Gender\')
plt.subplot(247)
plt.plot(salary_model2.fittedvalues, salary_model2.resid, \'o\')
plt.xlabel(\'Predicted\')
plt.ylabel(\'Residual\')
plt.subplots_adjust(hspace=0.3, wspace=0.35)
plt.show()
# 方式二：函数处理
fig = plt.figure(figsize=(12, 8)) # 因图形重叠会导致部分文字叠加在一起，故事先定义好图像大小
from statsmodels.graphics.regressionplots import plot_regress_exog
plot_regress_exog(salary_model2, 1, fig=fig)
# 参数1：数据，参数2：模型中的第几个自变量，参数3：表示在指定的对象上绘图
plt.show()
# 自变量与残差之间的关系不明显，基本无关，符合\varepsilon对于所有x而言具有同方差性的假定；
# 残差大体均匀地分布在[-2000,2000]之间，其均值与0接近，符合\varepsilon零均值的假定
# 至此得回归方程
# y = 3908.93+506.73*Education+1.32*Begin_Salary-22.57*Experience
# +11501.13*position1+6691.69*position2+1970.65*Gender
# 虚拟变量为0时表示女性且职位为普通员工，相应的性别和岗位影响通过虚拟变量取值来得出不同的回归方程
非线性回归
很多情况下，变量之间可能存在非线性关系。若是用。新型回归来拟合曲线，则会丢失数据之间的大量有用信息，甚至会得出错误的结论。此时，应使用非线性回归分析来对数据之间的统计关系进行考察。
可线性化的非线形分析
变量之间的非线性形式较多，常见的有：幂函数形式、对数形式、指数形式、逻辑斯蒂形式、抛物线形式等。
非线性形式的变量关系一般可以通过变量代换或转换的方式转化为线性关系。
对于幂函数形式
\[y=\alpha x^{\beta} + \varepsilon
\]
模型两边取对数，可得
\[\ln y = \ln \alpha + \beta \ln x + \varepsilon
\]
令 \(y^{\'}=\ln y,x^{\'}=\ln x,\alpha^{\'}=\ln \alpha\)，可得线性模型如下：
\[y^{\'}=\alpha^{\'}+\beta x^{\'} + \varepsilon
\]
对数模型
\[y=\alpha+\beta\ln x + \varepsilon \to y=\alpha+\beta x^{\'}+\varepsilon
\]
指数模型
\[y=\alpha e^{\beta x}+ \varepsilon \to \ln y=\ln \alpha+\beta x + \varepsilon \to y^{\'}=\alpha^{\'}+\beta x + \varepsilon
\]
逻辑斯蒂模型
\[y = \frac{1}{1+\alpha e^{\pm\beta x}} + \varepsilon \to \ln(\frac{y}{1-y})=-\ln \alpha \pm \beta x \to ...
\]
类似地，在存在多个自变量情形下的非线性回归，也可以按照上述变量转换和代换的方式把多元非线性模型转化为多元线性模型。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 对25个电商平台在某时期的注册用户数和销售收入进行回归分析
eb = pd.read_csv(\'./data/electronic_business.csv\')
print(eb.head(5))
# registration sales
# 0 16.9 67.5
# 1 15.3 88.3
# 2 4.5 77.4
# 3 14.7 84.6
# 4 4.0 77.2
# 1.散点图展示收入与注册数关系
plt.scatter(eb[\'registration\'], eb[\'sales\'], marker=\'o\')
plt.xlabel(\'registration\')
plt.ylabel(\'sales\')
plt.show()
# 图形显示接近抛物线。
# 2.建立抛物线模型，把抛物线模型转换为多元线性模型
from statsmodels.formula.api import ols
foramula = \'sales~registration+np.square(registration)\'
eb_model = ols(foramula, data=eb).fit()
res = eb_model.summary2()
print(res)
# Results: Ordinary least squares
# ========================================================================
# Model: OLS Adj. R-squared: 0.968
# Dependent Variable: sales AIC: 148.0375
# Date: 2020-06-08 22:27 BIC: 151.6942
# No. Observations: 25 Log-Likelihood: -71.019
# Df Model: 2 F-statistic: 360.7
# Df Residuals: 22 Prob (F-statistic): 1.53e-17
# R-squared: 0.970 Scale: 19.520
# ------------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# ------------------------------------------------------------------------
# Intercept 9.1189 3.4492 2.6437 0.0148 1.9656 16.2721
# registration 20.7981 0.8147 25.5296 0.0000 19.1086 22.4876
# np.square(registration) -1.0394 0.0390 -26.6602 0.0000 -1.1202 -0.9585
# ------------------------------------------------------------------------
# Omnibus: 0.921 Durbin-Watson: 1.490
# Prob(Omnibus): 0.631 Jarque-Bera (JB): 0.777
# Skew: 0.095 Prob(JB): 0.678
# Kurtosis: 2.157 Condition No.: 780
# ========================================================================
# 拟合优度判定R^2=0.968, 拟合很好，解释性高；
# 回归方程总体显著性检验：F= 360.7，P值近似0，总体显著；
# 回归方程回归系数显著性检验：p近似0，非常显著
# 3.对模型进行诊断
# a.判定残差是否符合正态分布
sm.qqplot(eb_model.resid, fit=True, line=\'45\')
plt.show()
# 类ProbPlot的实例对象也可绘制P-P图或Q-Q图
sm.ProbPlot(eb_model.resid, stats.t, fit=True).ppplot(line=\'45\')
sm.ProbPlot(eb_model.resid, stats.t, fit=True).qqplot(line=\'45\')
plt.show()
# P-P图或Q-Q图在残差符合正态假定条件下，散点图看起来应该项一条截距为0、斜率为1的直线。
# 本例中数据大部分散落在45度线上，表明残差基本符合正态分布
# b.判定残差是否与变量有关
fig = plt.figure(figsize=(12, 8)) # 因图形重叠会导致部分文字叠加在一起，故事先定义好图像大小
from statsmodels.graphics.regressionplots import plot_regress_exog
# plot_regress_exog(eb_model, 1, fig=fig)
# # 参数1：数据，参数2：模型中的第几个自变量，参数3：表示在指定的对象上绘图
# plt.show()
plot_regress_exog(eb_model, 2, fig=fig)
plt.show()
# 自变量与残差之间的关系不明显，基本无关，符合\varepsilon对于所有x而言具有同方差性的假定；
# 残差大体均匀地分布在[-6,6]之间，其均值与0接近，符合\varepsilon零均值的假定
# 得回归方程
# y = 9.12 + 20.7981*registration - 1.04*registration^2
非线性回归模型
对常见非线性模型进行线性转换后用线性回归的参数估计方法进行参数估计虽然较简单，但是有时估计效果不理想。当对因变量 \(y\) 作变换时，由于线性回归的最小二乘估计是对变换后的 \(y\) 而不是直接对 \(y\) 进行估计，再次基础上估计的曲线可能会造成拟合效果并不理想。此外，有些时候，变量间的曲线关系比较明显，关系式已知，但是难以用变量变换或代换的方式将其线性化，这时，可考虑直接使用非线性最小二乘估计方法来估计模型参数。
此外，非线性回归模型有一种情况：模型中至少有一个参数不是线性的，该模型也可称之为非线性模型。
非线性模型的参数一般可以用最小二乘及迭代算法进行估计，主要估计方法有：最速下降法（Steepest-Descent）或梯度法（Gradient Method）、牛顿法（Newton Method）、修正高斯-牛顿法（Modified Gauss-Newton Method）和麦夸特法（Marquardt Method）等。
一般而言，非线性曲线的拟合程度均较高，在考虑实际数据的拟合问题时，一般的分析过程往往不会给出模型检验结果。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 对电商平台在某时期的注册用户数和销售收入进行回归分析
eb_ex = pd.read_csv(\'./data/eb_extended.csv\')
print(eb_ex.head(5))
# registration sales
# 0 16.9 67.5
# 1 15.3 88.3
# 2 4.5 77.4
# 3 14.7 84.6
# 4 4.0 77.2
# 1.观察变量关系
plt.scatter(eb_ex[\'registration\'], eb_ex[\'sales\'], marker=\'o\')
plt.xlabel(\'registration\')
plt.ylabel(\'sales\')
plt.show()
# 由图形可知：
# 在注册用户<20是，变量之间关系可以用抛物线描述；在注册用户>20时，变量之间关系可以用指数形式
# 2.建立分段模型
from scipy.optimize import curve_fit
# 抛物线形式
def func1(x, alpha1, beta1, beta2):
return alpha1 + beta1 * x + beta2 * np.square(x)
# 指数形式
def func2(x, alpha2, beta3, gamma):
return alpha2 * (gamma ** (beta3 * x))
# 3.拟合参数
# 进行非线性最小二乘曲线拟合
m1 = curve_fit(func1, eb_ex[eb_ex[\'registration\'] <= 20][\'registration\'],
eb_ex[eb_ex[\'registration\'] <= 20][\'sales\'])
m2 = curve_fit(func2, eb_ex[eb_ex[\'registration\'] > 20][\'registration\'],
eb_ex[eb_ex[\'registration\'] > 20][\'sales\'])
print(m1, m2)
# 返回使离差f(x)-y平方和最小的参数估计值以及参数估计值的协方差
# (array([ 9.11887122, 20.79806299, -1.03935516]),
# array([[ 1.18971759e+01, -2.54291948e+00, 1.10578766e-01],
# [-2.54291948e+00, 6.63677073e-01, -3.11621801e-02],
# [ 1.10578766e-01, -3.11621801e-02, 1.51985540e-03]]))
# (array([3.3457962 , 1.27874618, 1.02971756]),
# array([[ 1.25324968e+00, -1.92397588e+05, 4.53702683e+03],
# [-1.92397590e+05, 3.21232066e+10, -7.57513956e+08],
# [ 4.53702686e+03, -7.57513956e+08, 1.78633286e+07]]))
# 得到分段非线性模型如下
# sales = 9.119 + 20.798*registration - -1.039*registration^2; if registration<=20;
# sales = 3.346 + 1.279*e^{1.030*registration}; if registration>20
# 4.预测
# 使用回归方程对给定用户注册数量对销售进行预测
eb_ex.loc[eb_ex[\'registration\'] <= 20, \'predicted\'] = func1(eb_ex[eb_ex[\'registration\'] <= 20][\'registration\'],
m1[0][0], m1[0][1], m1[0][2])
eb_ex.loc[eb_ex[\'registration\'] > 20, \'predicted\'] = func2(eb_ex[eb_ex[\'registration\'] > 20][\'registration\'],
m2[0][0], m2[0][1], m2[0][2])
print(eb_ex.head(5))
# registration sales predicted
# 0 16.9 67.5 63.755910
# 1 15.3 88.3 84.026587
# 2 4.5 77.4 81.663213
# 3 14.7 84.6 90.256142
# 4 4.0 77.2 75.681441
# 根据模型预测进行曲线拟合
eb_ex_sorted = eb_ex.sort_values(by=\'registration\')
plt.plot(eb_ex_sorted[\'registration\'], eb_ex_sorted[\'sales\'], \'o\',
eb_ex_sorted[\'registration\'], eb_ex_sorted[\'predicted\'], \'r-\', lw=3)
plt.show()
# 通过图形知，非线性分段方程较好地拟合了样本数据
多项式回归
如果有多个自变量的情况下，无法绘制散点图，同时很难对模型形式进行预估，这时需要使用多项式回归。根据数学的相关理论，任何曲线均可以使用多项式进行逼近，这种逼近的分析过程即多项式回归。
多项式回归类似于可线性化的非线性模型，可通过变量代换的方式使用普通最小二乘对参数进行估计。
设有因变量 \(y\) 和自变量 \(x\)，它们之间的关系为 \(n\) 次多项式的关系，则有如下模型
\[y = \alpha + \beta_1 x + \beta_2 x^2 +\cdots+\beta_n x^n + \varepsilon
\]
令 \(x_1=x, x_2=x^2,\cdots,x_n=x^n\)，则多项式模型就转化为如下的多元线性模型
\[y = \alpha + \beta_1 x_1 + \beta_2 x_2 +\cdots+\beta_n x_n + \varepsilon
\]
这样就可以按照多元线性回归模型进行分析了。对于多元的多项式模型
\[y = \alpha + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1^2 + \beta_4 x_1x_2 + \beta_5 x_2^2 +\cdots + \beta_m x_m^n + \varepsilon
\]
做变量代换，令 \(z_1=x_1,z_2=x_2,z_3=x_1^2,z_4=x_1x_2,z_5=x_2^2,\cdots,z_m=x_m^n\)，则有
\[y = \alpha + \beta_1 z_1 + \beta_2 z_2 + \beta_3 z_3 + \beta_4 z_4 + \beta_5 z_5 +\cdots + \beta_m z_m + \varepsilon
\]
转化之后的模型同样可以按照多元性回归模型进行分析。
多项式回归当阶数过高时，待估参数过多，在样本量不大的情况下会比较困难，这时多项式回归的一大缺陷。因此，一般的多项式回归模型很少应用到三阶以上。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 吉尼系数是基尼根据洛伦兹曲线所定义用于判断收入分配公平程度的指标，值在0和1之间，越往1靠近则表示收入分配公平性越低。
# 拟合洛伦兹曲线的方法主要有广义二次洛伦兹曲线法、指数曲线法等
# 利用3008年人均家庭收入数据，计算出收入累计百分比(cincome)和人口累计百分比(cpop)，对洛伦兹曲线进行拟合
lorenz = pd.read_csv(\'./data/lorenz.csv\')
print(lorenz.head(5))
# cpop cincome
# 0 0.018 0.000000e+00
# 1 0.019 9.701260e-07
# 2 0.020 3.366135e-06
# 3 0.020 5.785387e-06
# 4 0.021 8.521446e-06
# 由于事先不知道曲线的形状，因此无法像非线形回归模型那样为模型设置参数，故考虑使用多项式回归
# plt.scatter(lorenz[\'cpop\'], lorenz[\'cincome\'], marker=\'o\')
# plt.xlabel(\'cpop\')
# plt.ylabel(\'cincome\')
# plt.show()
# 1.拟合多项式函数
lorenz_est = np.polyfit(lorenz[\'cpop\'], lorenz[\'cincome\'], 2)
# 参数1：自变量，参数2：因变量，参数3：阶数
print(lorenz_est) # 多项式系数
# [ 1.02722918 -0.20824691 0.02550764]
lorenz_poly = np.poly1d(lorenz_est)
print(lorenz_poly) # 多项式回归方程
# 1.027 x^2 - 0.2082 x + 0.02551
# 2.模型检验
# 多项式回归模型转换后按照多元线性回归模型进行分析
from statsmodels.formula.api import ols
foraula = \'cincome~cpop+np.square(cpop)\'
lorenz_model1 = ols(foraula, data=lorenz).fit()
res = lorenz_model1.summary2()
print(res)
# Results: Ordinary least squares
# ===================================================================
# Model: OLS Adj. R-squared: 0.994
# Dependent Variable: cincome AIC: -6659.2479
# Date: 2020-06-09 09:20 BIC: -6643.6872
# No. Observations: 1322 Log-Likelihood: 3332.6
# Df Model: 2 F-statistic: 1.105e+05
# Df Residuals: 1319 Prob (F-statistic): 0.00
# R-squared: 0.994 Scale: 0.00037922
# -------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------
# Intercept 0.0255 0.0017 14.9508 0.0000 0.0222 0.0289
# cpop -0.2082 0.0077 -26.9092 0.0000 -0.2234 -0.1931
# np.square(cpop) 1.0272 0.0074 139.4493 0.0000 1.0128 1.0417
# -------------------------------------------------------------------
# Omnibus: 732.937 Durbin-Watson: 0.002
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 8902.124
# Skew: 2.320 Prob(JB): 0.000
# Kurtosis: 14.835 Condition No.: 24
# ===================================================================
# R^2比较高，模拟优度比较高；F统计量比较大，P值接近0，模型总体显著；回归系数p值接近0，非常显著。
# 回归方程为：y = 0.0255 - 0.2082*cpop + 1.0272*cpop^2
# 拟合曲线
# 方法一
from statsmodels.graphics.regressionplots import plot_fit
plot_fit(lorenz_model1, 1)
plt.show()
# 方法二
plt.plot(lorenz[\'cpop\'], lorenz[\'cincome\'], \'ro\',
lorenz[\'cpop\'], lorenz_model1.fittedvalues, linewidth=2)
plt.show()
# 3阶模拟
def cube(x):
return x**3
formula = \'cincome~cpop+np.square(cpop)+cube(cpop)\'
lorenz_model2 = ols(formula, data=lorenz).fit()
res = lorenz_model2.summary2()
print(res)
# Results: Ordinary least squares
# ===================================================================
# Model: OLS Adj. R-squared: 0.998
# Dependent Variable: cincome AIC: -7995.1774
# Date: 2020-06-09 09:32 BIC: -7974.4298
# No. Observations: 1322 Log-Likelihood: 4001.6
# Df Model: 3 F-statistic: 2.034e+05
# Df Residuals: 1318 Prob (F-statistic): 0.00
# R-squared: 0.998 Scale: 0.00013794
# -------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------
# Intercept -0.0238 0.0015 -16.4017 0.0000 -0.0267 -0.0210
# cpop 0.3350 0.0122 27.3831 0.0000 0.3110 0.3589
# np.square(cpop) -0.2877 0.0277 -10.3769 0.0000 -0.3421 -0.2333
# cube(cpop) 0.8611 0.0179 48.0430 0.0000 0.8259 0.8963
# -------------------------------------------------------------------
# Omnibus: 943.668 Durbin-Watson: 0.005
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 23537.954
# Skew: 2.997 Prob(JB): 0.000
# Kurtosis: 22.784 Condition No.: 134
# ===================================================================
# 拟合优度、总体显著性、回归系数显著性更好
# 回归方程：y = -0.0238 + 0.3350*cpop - 0.2877*cpop^2 + 0.8611*cpop^3
# 拟合曲线
plot_fit(lorenz_model2, 1)
plt.show()
分位数回归
前面的回归分析都是通过对模型求其条件期望的形式，然后利用最小误差平方（即最小二乘）进行参数估计。估计方程为
\[\hat{y}^{OLS} = \hat{\beta}_0^{OLS} + \hat{\beta}_1^{OLS} x_1 + \hat{\beta}_2^{OLS} x_2+ \cdots + \hat{\beta}_n^{OLS} x_n
\]
其模型的实质是随机变量的条件均值函数。因此，回归分析得到的常见结论是：\(x\) 变动一个单位，\(y\) 平均变动 \(\beta\) 个单位。但是当数据具有尖峰厚尾的分布特征或有离群点（即异常值）时，模型稳健型较差。
分位数回归(quantile regression)便是对以古典条件均值模型为基础的延伸，它用几个分位函数来估计整体模型。这种模型依据因变量的条件分位数对自变量进行回归，可得到所有分位数下的回归模型，能够得到更加稳健的估计结果。
分位数其实质是顺序数据划分为 \(k\) 个等分，每个等分点上的数值。因此，实值随机变量 \(y\) 的右连续分布函数为
\[F(y)=P(Y\leq y)
\]
则 \(y\) 的 \(\tau\) 分位数为 \(Q(\tau)=\inf(y:F(y)\geq \tau)\)，对回归模型取条件中位数，可得
\[Median(y|x) = \beta_0 + \beta_1 x_1 + \beta_2 x_2+\cdots + \beta_n x_n + Median(\varepsilon)
\]
然后通过最小化绝对误差损失函数
\[\min \sum|\mu| = \min|y-(\beta_0 + \beta_1 x_1 + \beta_2 x_2+\cdots + \beta_n x_n)|
\]
可得到最小绝对偏差（LAD，least abasolute deviation）的估计方程
\[\hat{y}^{LAD} = \hat{\beta}_0^{LAD} + \hat{\beta}_1^{LAD} x_1 + \hat{\beta}_2^{LAD} x_2+ \cdots + \hat{\beta}_n^{LAD} x_n
\]
上述方程的估计系数则反映了各个自变量对因变量中位数的影响。因为中位数也是分位数的一种特殊形式，更进一步，可以取回归模型的条件 \(\tau\) 分位数，即
\[Q_{y|x}(\tau)=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n + Q_{\varepsilon|x}(\tau)
\]
通过最小化不对称损失函数
\[\min\sum{\rho_{\tau}(y-(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n))}
\]
的方式，得到分位数回归方程
\[\hat{y}^{\tau} =\hat{\beta}_0^{\tau} + \hat{\beta}_1^{\tau} x_1 + \hat{\beta}_2^{\tau} x_2+ \cdots + \hat{\beta}_n^{\tau} x_n
\]
在分位数回归分析中，损失函数可进行如下定义
\[\rho_{\tau}(\varepsilon) = \varepsilon(\tau-I(\varepsilon))
\]
其中，\(I(\varepsilon)\) 为示性函数，即
\[I(\varepsilon)= \begin{cases}
0 ,& \varepsilon \geq 0 \\
1 ,& \varepsilon < 0
\end{cases}
\]
则有
\[\rho_{\tau}(\varepsilon) = \begin{cases}
\tau\varepsilon ,& \varepsilon \geq 0 \\
\varepsilon(\tau-1) ,& \varepsilon < 0
\end{cases},且\rho_{\tau}(\varepsilon)\geq 0
\]
因此，可将损失函数改写为
\[\rho_{\tau}(\varepsilon) = \varepsilon(\tau-I(\varepsilon))=\begin{cases}
\tau\varepsilon ,& \varepsilon \geq 0 \\
\varepsilon(\tau-1) ,& \varepsilon < 0
\end{cases}
= \tau\varepsilon I_{(\varepsilon \geq 0)} + (\tau-1)\varepsilon I_{\varepsilon<0}
\]
最小化不对成损失函数
\[\min\sum{\rho_{\tau}(y-\xi)} = \min \bigg[\sum_{y_i \geq \xi}{\tau|y-\xi|} + \sum{(\tau-1)|y-\xi|}\bigg]
\]
其中，\(\xi = (\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n)\)。于是，求 \(\tau\) 分位点的参数估计如下
\[\hat{\beta}^{\tau} = \arg\min_{\beta\in R^m}(\sum{\rho_{\tau}(y-\beta x)})
\]
其中，\(\beta, x\) 分别为回归系数和自变量的向量。
在 \(\tau \in (0,1)\)上，可得到 \(y\) 关于 \(x\) 的条件分布曲线簇。分位数回归可用的估计方法包括：单纯型法、内点法(预处理后内点法)、平滑法及 Majorize-Mininize法等。
与一般的回归分析过程一样，分位数回归模型进行参数估计后，也需要对模型及性能评价以及进行显著性检验。Koenker与Machado依据最小二乘中拟合优度 \(R^2\) 的计算思想，提出了一种拟合优度统计量可以用来判定模型的拟合程度，还给出了拟似然比检验对回归系数的显著性进行判定。除此之外，分位数回归中系数是随分位水平变化的，因此可对于不同分位点上的模型斜率相等性进行检验（即异方差性检验）,也可对不同分位点的模型关于中位数是否对称进行对称性检验。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 对公司人员的当前薪酬及其影响因素进行分位数回归分析
# 在现实状况中，员工的现有工资水平是有等级之分的，这些不同等级或层次是分位数的一种具体体现。
salary = pd.read_csv(\'./data/salary_r.csv\')
print(salary.head(5))
# position ID Gender ... Begin_Salary Experience Age
# 0 1 1 1 ... 27000 144.0 55.0
# 1 1 34 1 ... 39990 175.0 58.0
# 2 1 18 1 ... 27510 70.0 51.0
# 3 1 200 1 ... 34980 9.0 44.0
# 4 1 199 1 ... 27480 69.0 49.0
# 数据中有较多缺失值，需做预处理。
# 使用statsmodels会自动不考虑有缺失值的行
salary = salary.dropna() # 删除有缺失值的行
# 分位数回归分析
from statsmodels.formula.api import quantreg
formula = \'Current_Salary~Begin_Salary+Education+Experience\'
salary_qt = quantreg(formula, data=salary)
salary_qtmodel1 = salary_qt.fit(q=0.1) # q表示分位点，此处10%分位点上的模型和拟合结果
res = salary_qtmodel1.summary()
print(res)
# QuantReg Regression Results
# ==============================================================================
# Dep. Variable: Current_Salary Pseudo R-squared: 0.4469
# Model: QuantReg Bandwidth: 3071.
# Method: Least Squares Sparsity: 1.726e+04
# Date: Tue, 09 Jun 2020 No. Observations: 446
# Time: 11:15:52 Df Residuals: 442
# Df Model: 3
# ================================================================================
# coef std err t P>|t| [0.025 0.975]
# --------------------------------------------------------------------------------
# Intercept 20.1218 1435.641 0.014 0.989 -2801.408 2841.652
# Begin_Salary 1.4562 0.056 25.805 0.000 1.345 1.567
# Education 302.4687 134.098 2.256 0.025 38.919 566.018
# Experience -13.2062 2.268 -5.822 0.000 -17.664 -8.748
# ================================================================================
# 在给定的显著性水平\alpha=0.1下，在0.1分位点下模型的所有参数估计结果均显著
salary_qtmodel2 = salary_qt.fit(q=0.85, max_iter=20000)
res = salary_qtmodel2.summary()
print(res)
# QuantReg Regression Results
# ==============================================================================
# Dep. Variable: Current_Salary Pseudo R-squared: 0.6465
# Model: QuantReg Bandwidth: 3103.
# Method: Least Squares Sparsity: 2.767e+04
# Date: Tue, 09 Jun 2020 No. Observations: 446
# Time: 11:33:33 Df Residuals: 442
# Df Model: 3
# ================================================================================
# coef std err t P>|t| [0.025 0.975]
# --------------------------------------------------------------------------------
# Intercept -5684.6774 2716.874 -2.092 0.037 -1.1e+04 -345.082
# Begin_Salary 2.4948 0.082 30.530 0.000 2.334 2.655
# Education 440.0576 241.952 1.819 0.070 -35.461 915.576
# Experience -12.9493 5.476 -2.365 0.018 -23.711 -2.188
# ================================================================================
# 在给定的显著性水平\alpha=0.1下，在0.85分位点下模型的所有参数估计结果均显著
# 对比两个模型的结果，对于高薪和底薪的员工而言，起始薪酬、工作经验、受教育水平的影响是非常显著的。
# 对比各分位点上回归模型的参数估计值
import warnings
warnings.filterwarnings(\'ignore\') # 分位数回归执行过程中有些分位点估计参数时会有收敛循环提示，屏蔽这些信息
quanties = np.arange(0.05, 1, 0.05)
def fit_model(v, q):
res = salary_qt.fit(q=q)
return [q, res.params[\'Intercept\'], res.params[v]] + res.conf_int(alpha=0.05).loc[v].tolist()
# 自动生成各自变量在不同分位点上的参数估计值及其95%置信区间
names = locals()
for i in [\'Intercept\', \'Begin_Salary\', \'Education\', \'Experience\']:
names[\'params_%s\' % i] = [fit_model(i, x) for x in quanties]
names[\'params_%s\' % i] = pd.DataFrame(
names[\'params_%s\' % i], columns=[\'quantile\', \'Intercept\', \'Beta\', \'Beta_L\', \'Beta_U\'])
# 会自动生成内存变量：params_Intercept,params_Begin_Salary,params_Education,params_Experience
print(params_Intercept.head(5))
# quantile Intercept Beta Beta_L Beta_U
# 0 0.05 29.668429 29.668429 -2633.298978 2692.635836
# 1 0.10 20.121792 20.121792 -2801.408031 2841.651615
# 2 0.15 32.929636 32.929636 -2615.891367 2681.750640
# 3 0.20 689.525970 689.525970 -2000.781722 3379.833661
# 4 0.25 21.784563 21.784563 -2853.827528 2897.396655
print(params_Begin_Salary.head(5))
# quantile Intercept Beta Beta_L Beta_U
# 0 0.05 29.668429 1.290738 1.153871 1.427606
# 1 0.10 20.121792 1.456163 1.345258 1.567067
# 2 0.15 32.929636 1.505495 1.401500 1.609491
# 3 0.20 689.525970 1.551411 1.446781 1.656041
# 4 0.25 21.784563 1.630912 1.522986 1.738838
# 绘制图形
# 2*2共享x轴坐标的图形
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True)
fig.set_size_inches(6, 6)
fig.suptitle(\'Estimated Paramter by Quantile for Current_Salary with 95% confidence limits\')
ax1.plot(params_Intercept[\'quantile\'], params_Intercept[\'Beta\'])
ax1.fill_between(params_Intercept[\'quantile\'], params_Intercept[\'Beta_L\'], params_Intercept[\'Beta_U\'],
color=\'b\', alpha=0.1)
ax1.set_ylabel(\'Intercept\')
ax2.plot(params_Begin_Salary[\'quantile\'], params_Begin_Salary[\'Beta\'])
ax2.fill_between(params_Begin_Salary[\'quantile\'], params_Begin_Salary[\'Beta_L\'], params_Begin_Salary[\'Beta_U\'],
color=\'b\', alpha=0.1)
ax2.set_ylabel(\'Begin_Salary\')
ax3.plot(params_Education[\'quantile\'], params_Education[\'Beta\'])
ax3.fill_between(params_Education[\'quantile\'], params_Education[\'Beta_L\'], params_Education[\'Beta_U\'],
color=\'b\', alpha=0.1)
ax3.set_xlabel(\'quantiles\')
ax3.set_ylabel(\'Education\')
ax4.plot(params_Experience[\'quantile\'], params_Experience[\'Beta\'])
ax4.fill_between(params_Experience[\'quantile\'], params_Experience[\'Beta_L\'], params_Experience[\'Beta_U\'],
color=\'b\', alpha=0.1)
ax4.set_xlabel(\'quantiles\')
ax4.set_ylabel(\'Experience\')
fig.subplots_adjust(wspace=0.3) # 将横向图形之间宽度调宽以便美观
plt.show()
# 由图可知各个自变量第因变量条件分位数影响：
# 对于起始薪酬自变量而言，当前收入越高(分位点越高)的阶层，其收入受起始薪酬的影响越大，有呈线性关系的趋势
# 受教育年限在中低收入阶层影响较稳定，但对高收入阶层而言，该自变量体现出的学历影响则有显著提升
# 工作经历对中等薪酬阶层的影响较小，对于高薪阶层而言，工作经历的影响稍稍偏大