2.1.8 Python案例：线性回归调优实战
1.残差图
为了观察残差的分布，我们先定义一个绘制“残差-预测值”及“残差-自变量”的函数：
def plot_resid(resid, x, y_hat, **kw):
plt.figure(figsize=(8, 8))
# "残差-X"散点图
ax1 = plt.subplot(211)
ax1.set_ylabel('resid', fontsize=15)
ax1.set_xlabel('x', fontsize=15)
ax1.scatter(x, resid, **kw)
ax1.set_title('Scatter of resid - x', fontsize=15)
# "残差-Y估计"散点图
ax2 = plt.subplot(212)
ax2.set_ylabel('resid', fontsize=15)
ax2.set_xlabel('y_hat', fontsize=15)
ax2.scatter(y_hat, resid, **kw)
ax2.set_title('Scatter of resid - y_hat', fontsize=15)
# 输出图形
plt.tight_layout()
plt.show()
建立三个模型，除了对信用卡支出和收入做回归，还分别对二者取对数后做回归，其中信用卡支出的对数包含在原数据表中，因此单独对收入取对数放入原表。这里需要注意在训练集上做的任何变换，在测试集和验证集中也需要做同样的变换。示例如下：
train['Income_ln'] = np.log(train['Income']) # 样本外数据做同样变换 oos_data['Income_ln'] = np.log(oos_data['Income']) # 对x和y各自取对数，比较结果 linear_model1 = ols(formula='avg_exp ~ Income', data=train).fit() linear_model2 = ols(formula='avg_exp_ln ~ Income', data=train).fit() linear_model3 = ols(formula='avg_exp_ln ~ Income_ln', data=train).fit()
对支出和收入不取对数时的模型残差分布：
plot_resid(linear_model1.resid,
train['Income'],
linear_model1.predict(train))
输出结果如图2-17所示。
图2-17 残差图1
图2-17 （续）
可以看到在简单线性回归中，“残差-自变量”和“残差-预测值”的分布十分类似（当然横轴坐标是不一样的）。根据2.1.6节中的论述可知，使用原始数据时，残差出现了异方差情况。在这种情况下，可使用加权最小二乘法，但是实际上最常用的是对被解释变量取对数。从业务角度来看，对信用卡支出和收入做对数转换也更加符合业务，因此我们再比较一下对X、Y都做对数转换时的残差图：
plot_resid(linear_model3.resid,
train['Income_ln'],
linear_model3.predict(train))
结果如图2-18所示。
图2-18 残差图2
可以看到，取对数后残差的分布更符合独立同分布的假设。
不妨再看一下残差的分布是否符合正态分布：
fig, ax = plt.subplots(figsize=(6, 6)) ax.set_xlabel(..., fontsize=15) ax.set_ylabel(..., fontsize=15) sm.qqplot(linear_model3.resid, fit=True, line='45', ax=ax, figure=fig) plt.show()
输出的QQ图如图2-19所示。
图2-19 QQ图
可以看到，取对数后残差的分布基本符合正态分布。
我们再比较一下三个模型的拟合优度，并且使用pprint让输出结果更加美观：
import pprint
r_squares = {'exp ~ Income': linear_model1.rsquared,
'ln(exp) ~ Income': linear_model2.rsquared,
'ln(exp) ~ ln(Income)': linear_model3.rsquared}
pprint.pprint(r_squares)
输出结果如下：
{'exp ~ Income' : 0.45429062315565294,
'ln(exp)~ Income': 0.403085555532965,
'ln(exp)~ ln(Income)': 0.4803927993893108}
可以发现，自变量和因变量都取对数后的拟合优度最高。
2.强影响点处理
从残差图和QQ图中可以看出，有两个数据点的残差偏离较大，我们使用学生化残差进行判断：
train['resid'] = linear_model3.resid student_resid = (train['resid'] - train['resid'].mean() ) / train['resid'].std() train[np.abs(student_resid) > 2]
因为样本量较小，我们将学生化残差绝对值大于2的样本输出，结果如下：
可以看到，73号和98号样本是离群点，对模型影响较大。此外，statemodels包提供了更多强影响点判断指标：
from statsmodels.stats.outliers_influence import OLSInfluence OLSInfluence(linear_model3).summary_frame().head()
每个样本都会获得一系列用于判断其是否是离群点的指标，具体如下：
我们选择任意一个指标去判断强影响点即可。这里选择删除学生化残差绝对值大于2的样本，重新建立模型，并观察残差：
train = train[np.abs(student_resid) <= 2].copy()
linear_model4 = ols('avg_exp_ln ~ Income_ln', train).fit()
plot_resid(linear_model4.resid,
train['Income_ln'],
linear_model4.predict(train))
此时的残差如图2-20所示。
图2-20 绝对值大于2的残差
图2-20 （续）
3.多重共线性的识别和处理
经过单变量线性回归的处理，我们基本对模型的性质有了一定的了解，接下来可以放入更多的连续型解释变量。在加入变量之前，要注意变量的函数形式的转变，比如当地房屋均价、当地平均收入的性质和个人收入一样，都需要取对数，具体如下：
# 对解释变量取对数
train['dist_home_val_ln'] = np.log(train['dist_home_val'])
train['dist_avg_income_ln'] = np.log(train['dist_avg_income'])
# 对本期数据也做同样变换
oos_data['dist_home_val_ln'] = np.log(oos_data['dist_home_val'])
oos_data['dist_avg_income_ln'] = np.log(oos_data['dist_avg_income'])
linear_model5 = ols('''avg_exp_ln ~ Age + Income_ln +
dist_home_val_ln + dist_avg_income_ln''', train).fit()
linear_model5.summary()
生成的模型概述如下：
从输出结果中可以获得许多重要信息，笔者已在图中用数字编号标示出来：
①回归方程的F检验中P值很小，说明方程整体显著。
②假如设定的显著性水平为0.05，则Age、Income_ln和dist_avg_income_ln的回归参数检验中P值大于显著性水平，说明这几个回归参数不显著，可以考虑剔除。但要注意，一旦剔除一个变量，其他变量的显著性水平可能会发生较大变化，尤其当变量间存在多重共线性时。
③当回归系数不显著时，其标准误差（相对于参数估计本身）比较大，对应的参数置信区间也比较大，也说明了回归系数并不很可靠。
④DW统计量接近2，说明残差的自相关较低。
为了识别多重共线性的自变量，我们根据前文的描述定义一个计算方差膨胀因子的函数，并对案例中的4个连续型自变量进行分析。
def vif(df, y_col):
x_cols = list(df.columns)
x_cols.remove(y_col)
formula = y_col + '~' + '+'.join(x_cols)
rsquared = ols(formula, df).fit().rsquared
return 1. / (1. - rsquared)
exog = train[['Age', 'Income_ln', 'dist_home_val_ln', 'dist_avg_income_ln']]
for y_col in exog.columns:
print('VIF of "%s" is %.2f' %(y_col, vif(df=exog, y_col=y_col)))
输出结果如下：
VIF of "Age"is 1.17 VIF of "Income_ln" is 36.98 VIF of "dist_home_val_ln" is 1.05 VIF of "dist_avg_income_ln"is 36.92
根据判断标准，Income_ln和dist_avg_income_ln的VIF大于10，意味着存在多重共线性自变量，此时我们可以删除其中一个。考虑到这两个变量一个为收入（对数），另一个为地区平均收入（对数），使用高出平均收入的比率代替其中一个变量也许是一个更好的选择，这样新自变量既保留了信息，又与其他自变量的多重共线性较低。
其中，高出地区平均收入已经在原数据集中给出，我们只需要计算比率即可。同样，对新自变量计算方差膨胀因子：
train['high_avg_ratio'] = train['high_avg'] / train['dist_avg_income']
oos_data['high_avg_ratio'] =oos_data['high_avg']/oos_data['dist_avg_income']
exog = train[['Age', 'high_avg_ratio', 'dist_home_val_ln',
'dist_avg_income_ln']]
for y_col in exog.columns:
print('VIF of "%s" is %.2f' %(y_col, vif(df=exog, y_col=y_col)))
输出结果如下：
VIF of "Age" is 1.17 VIF of "high_avg_ratio" is 1.13 VIF of "dist_home_val_ln"is 1.05 VIF of "dist_avg_income_ln"is 1.31
我们发现新自变量的引入使得多重共线性问题基本解决了。不妨利用之前定义的函数，再次对变量进行逐步筛选：
var_select_data = train[['avg_exp_ln', 'Age', 'high_avg_ratio',
'dist_home_val_ln', 'dist_avg_income_ln']]
linear_model6 = forward_select(data=var_select_data,
response='avg_exp_ln')
print('r_squared is ', linear_model6.rsquared)
输出结果如下：
aic is 23.81679370073742,continuing! aic is 20.83095227956072,continuing! forward selection over! final formula is avg_exp_ln ~ dist_avg_income_ln + dist_home_val_1n r_squared is 0.5520397736845986
上述算法中去除了两个变量，包括新生成的变量。这是因为笔者生成的变量仍然不够好，也许可以尝试一下“高出地区平均收入的对数”或者“收入与房价比”等变量，而这取决于你对业务的理解。在实际中，建模之前就需要对数据执行特征工程，其中很重要的一点就是尽量多地生成既有业务意义，又与Y有关的候选特征（变量），随后再进行变量筛选。生成新变量到筛选变量的过程是必要的。
此外，要注意在使用逐步法进行变量筛选时，有较大的可能会去除那些具有多重共线性的变量，但并不总是有效，因此识别多重共线性变量仍然是需要做的。
4.离散变量的处理
线性回归适用于自变量和因变量都是连续型的情况。如果自变量包含离散型，则需对其做哑变量变换（也称为one-hot编码）。这个过程在许多框架中要手动执行，而在stats-models框架下可以自动完成，我们只需要在formula中指定变量是离散型即可，示例如下：
formula = ''' avg_exp_ln ~ dist_avg_income_ln + dist_home_val_ln + C(gender) + C(Ownrent) + C(Selfempl) + C(edu_class) ''' linear_model7 = ols(formula, train).fit() linear_model7.summary()
使用形如“C（离散变量）”的方式在formula中指定离散型变量，模型能自动对该变量进行哑变量变换。仅包含模型参数及检验部分的截图如下：
可以看到，gender、Ownrent、Selfempl都是仅包含两个分类水平的离散型变量，其哑变量变换就是其本身。以Selfempl为例，其回归系数说明：当控制其他变量的时候，自谋职业者（对应1）相对于无职业者（对应0）的信用卡支出的对数增加了-0.3805，相当于信用卡支出约增长了-31.6%（e-0.3805-1）。
由于教育程度edu_class包含4个分类水平，因此生成了3个哑变量编码，分别代表教育程度为1、2、3的3个水平。其中，教育程度为0的样本作为对照组。
Ownrent变量的参数不显著，可以考虑删除。此外，我们还可以对离散型变量的交叉项进行建模，只需在formula中对两个离散型变量使用星号“*”连接，Statsmodels就会将这两个离散型变量及其交叉项放入模型。如果只想放入交叉项，则使用冒号“:”连接离散型变量即可：
formula = ''' avg_exp_ln ~ dist_avg_income_ln + dist_home_val_ln + C(Selfempl) + C(gender) * C(edu_class) ''' linear_model8 = ols(formula, train).fit() linear_model8.summary()
模型概述中的参数估计及检验部分截图如下：
可以看到，交叉项的系数显著，但是gender本身的系数不显著，可以考虑仅使用交叉项。要注意一点，如果一个离散型变量（或者交叉项）对应的哑变量中的部分系数显著，部分系数不显著，不能只删除不显著的部分，即一个原始离散型变量对应的哑变量要么都保留，要么都不保留。例如我们不能只删除教育程度为高中对应的哑变量，而只保留教育程度为小学、初中、大学对应的哑变量。不过，可以对系数不显著的变量进行重新编码，比如将初中和高中合并为中学水平。
关于模型概述中未展示的部分，读者可以自行生成并解读。
5.保存结果
模型经过检验后，我们通过样本外数据进行预测，并将结果保存了下来。其中，注意两点：
1）由于样本量较小，训练集中的edu_class的取值包括1、2、3，而测试集中的edu_class则只有0、1。由于小学及以下的样本极少，且和中学教育程度比较接近（都是低学历），因此将二者合并入小学教育程度。
2）建立的模型预测的是信用卡支出的对数值，需要进行转换获得信用卡支出的预测值，代码如下：
output_data_dir = './data_processed'
if not os.path.exists(output_data_dir):
os.mkdir(output_data_dir)
output_data_path = os.path.join(output_data_dir,
'creditcard_exp_linear_regression.csv')
# 将edu_class=0 概化为1
oos_data['edu_class'].replace({0: 1}, inplace=True)
oos_data['Pred_avg_exp'] = np.exp(linear_model8.predict(oos_data))
oos_data.to_csv(output_data_path, index=False)
这样，我们就将预测结果保存在了creditcard_exp_linear_regression.csv文件中。