一、数据集
本次实战用到的是高炉煤气循环发电（CombinedCyclePowerPlant，CCPP）数据集，该数据集来自于UCI主页（http://archive.ics.uci.edu/ml/datasets.html）。
自变量 | AT，表示高炉的温度 |
V，表示炉内的压力； | |
AP，表示高炉的相对湿度 | |
RH，表示高炉的排气量 | |
因变量 | 连续性，PE，表示高炉的发电量 |
二、加载数据集
1、读取数据并展示数据的统计量。
data = pd.read_excel('./CCPP.xlsx')
print(data.describe())
2、绘制各变量之间的散点图，观察各变量的相关关系。
sns.pairplot(ccpp)
plt.show()
从上面的散点图来看，似乎AP（相对湿度）和RH（排气量）与PE（发电量）之间并不存在明显的线性关系。
三、建模
本次实战采用OLS（ordinary least squares，普通最小二乘）回归模型。
1、拆分训练集（用于建模）和测试集（用于模型的评估）。
train, test = train_test_split(ccpp, test_size=0.2, random_state=14)2、建模
fit_train = sm.formula.ols('PE~AT+V+AP', data=ccpp).fit()
print(fit_train.summary())
使用 fit_train.summary()的参数解释：
part1：
1、Dep.Variable：因变量，此回归分析中为PE（发电量）。
2、model：方法，此回归分析中为普通最小二乘法（OLS）。
3、No.Observations：观察次数，为样本的大小。
4、Df Residuals：残差的自由度（Df），Df = N- K（N = 样本数量（观察数），k = 变量数 + 1）。
5、Df Models：模型的自由度（Df），Df = K - 1（K = 变量数 + 1）。
6、R-squared：R-平方值，是确定系数，告诉我们自变量可以解释多少百分比的自变量。
7、Adj.R-suquared：它在R-squared的基础上，加入了一个“惩罚项”，当向现有模型加入一个“无关自变量”时，Adj.R-suquared会给这个“无关自变量”一个惩罚，从而使得Adj.R-suquared的值不一定增加，防止了“虚假提升信息的产生”。
8、F-statistic：用于完成模型的显著性检验，模型的显著性检验是指构成因变量的线性组合是否有效。
9、AIC：赤池信息准则，是衡量统计模型拟合优良性(Goodness of fit)的一种标准，赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。即用最少的自变量达到最好的拟合效果。
part2：
1、Intercept：常数项，是回归线的截距。在回归中，我们忽略了一些对因变量没有太大影响的自变量，截距表明了这些遗漏变量的平均值和模型中存在的噪声。
2、coef（Coefficient term）：系数项。如果熟悉导数，可以将其视为Y相对于X的变化率。
3、std err：标准误差也称为标准偏差。
4、t： t 统计值,衡量系数统计显著程度的指标。
5、P > |t|：P值是一种概率，指的是在H0假设为真的前提下，样本结果出现的概率。当P值小于0.05时，说明模型是通过显著性检验的，需要拒绝原假设。
6、[0.025,0.975]：置信区间，表示我们的系数可能下降的范围（可能性为95%）。
三、回归模型的假设性检验
回归模型的检验目的是为了检验模型是否有效，包括了模型总体是否有效以及其中单个变量是否具有统计学意义。
1、模型的显著性检验--F检验，检测自变量是否真正影响到因变量的波动。
1）提出问题的原假设和备择假设：
原假设：模型的所有偏回归系数为0。
备择假设：模型的所有偏回归系数不全为0，即至少存在一个自变量可以构成因变量的线性组合。
2）构造F-statistic：模型拟合的越好，F-statistic值越大。
python计算F统计量的函数：回归模型拟合后，调用model.fvalue即可输出 f 的实际值，同时，可通过引入scipy.stats模块计算 f 的理论值。dfn为自由度，n为样本数量。
通常来说，F的实际值大于理论值，则会拒绝原假设，即模型是显著的，模型的所有偏回归系数不全为0。f的实际值也可通过model.summary得出。
from scipy.stats import f
print("f的实际值：", fit_train_1.fvalue)
p = fit_train_1.df_model
n = train.shape[0]
f_theroy = f.ppf(q=0.95, dfn=p, dfd=n-p-1)
print("f的理论值：", f_theroy)
在本次实战中，通过模型反馈的F统计值<F的理论值，说明模型是通过显著性检验的，说明需要拒绝原假设（即认为模型的所有回归系数都不全为0），但模型的显著性通过检验并不代表每个变量都对因变量是重要的，所以还需要进行偏回归系数的显著性检验。
2、回归系数的显著性检验--t检验，单个自变量在模型中是否有效。
t检验同样有原假设，即：变量不具有统计学意义。
t检验在python中的操作很简单，只需要通过调用模型摘要：model.summary()即可看到每个变量的p值。p小于0.05的话则说明在95%的置信水平下，该变量具有统计学意义。
通过上图的检验结果显示，P>|t|这一栏，结果都小于0.05，均通过显著性检验。说明温度AT，压力A，相对湿度AP都影响到排气量PE的变动，故不需要将其中的变量从模型中剔除。
四、变量筛选
多元线性回归能够按照一些方法筛选建立回归的自变量，这些方法包括：向前法，向后法，逐步法。这三种方法进入或者剔除变量的一个标准是AIC准则，即最小信息准则。AIC值越小说明模型效果越好。
x.join(y)方法，y是可迭代对象，x插入可迭代对象中间，形成字符串。
举例：'+'.join(['AE', 'V', 'AP'])，输出：'AT+V+AP'。（用+号将多个字符串连接在一起）。
def forward_select(data, response):
remaining = set(data.columns)
remaining.remove(response)
selected = []
current_score, best_new_score = float('inf'), float('inf')
while remaining:
aic_with_candidates = []
for candidate in remaining:
formula = '{} ~ {}'.format(response, ' + '.join(selected + [candidate]))
aic = sm.formula.ols(formula=formula, data=data).fit().aic
aic_with_candidates.append((aic, candidate))
aic_with_candidates.sort(reverse=True)
best_new_score, best_candidate = aic_with_candidates.pop()
if current_score > best_new_score:
remaining.remove(best_candidate)
selected.append(best_candidate)
current_score = best_new_score
print('aic is {}, continuing!'.format(current_score))
else:
print('forward selection over!')
break
formula = '{} ~ {}'.format(response,' + '.join(selected))
print('final formula is {}'.format(formula))
model = sm.formula.ols(formula=formula, data=data).fit()
return(model)
data_for_select = train[['PE', 'AT', 'V', 'AP']]
fit_train_1 = forward_select(data=data_for_select, response='PE')
print(fit_train_1.rsquared)
模型预测结果
pred = fit_train_1.predict(test[['AT', 'V', 'AP']])
print(test['PE'])
print(pred)
五、可视化
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = 'Microsoft YaHei'
plt.scatter(test['PE'], pred, label='观测点', color='b')
plt.plot([test['PE'].min(), test['PE'].max()], [pred.min(), pred.max()], 'r--', lw=2, label='拟合线')
plt.title('真实值VS.预测值')
plt.xlabel('预测值')
plt.ylabel('真实值')
plt.legend(loc='upper left')
plt.show()
六、回归模型的诊断
由于本线性回归模型的偏回归系数是通过最小二乘法（OLS）实现的，关于最小二乘法的使用是有一些假设前提的，具体是：
自变量与因变量之间存在线性关系；
自变量之间不存在多重共线性；
回归模型的残差服从正态分布；
回归模型的残差满足方差齐性（即方差为某个固定值）；
回归模型的残差之间互相独立。
1、线性相关检验
关于线性关系的判断，我们可以通过图形或Pearson相关系数来识别，具体Python代码如下：
# 绘制各变量之间的散点图
sns.pairplot(ccpp)
plt.show()
# 发电量与自变量之间的相关系数
print(ccpp.corrwith(ccpp['PE']))
从返回的结果来看，PE（发电量）与AT（温度）和V（压力）之间的相关系数较高，而PE与AP（相对湿度）和RH（排气量）之间的相关系数较小。
一般情况下，当Pearson相关系数低于0.4时，表明变量之间存在弱相关关系；当Pearson相关系数在0.4~0.6之间时，则说明变量之间存在中度相关关系；当相关系数在0.6以上时，则反映变量之间存在强相关关系。
经对比发现，PE与RH之间的为弱相关关系，故不考虑将该变量纳入模型。当然，变量之间不存在任何关系，可能是二次函数关系、对数关系等，所以一般还需要进行检验和变量转换。
2、多重共线性检验
自变量之间不能有强共线性，如果多元线性回归中出现这一问题，那么会造成对回归系数、截距系数的估计非常不稳定。关于多重共线性的检验可以使用方差膨胀因子（VIF）来鉴定，如果VIF>10，则说明变量存在多重共线性。一旦发现变量之间存在多重共线性的话，可以考虑删除变量和重新选择模型（岭回归法）。
方差膨胀系数(variance inflation factor，VIF)是衡量多元线性回归模型中复 (多重)共线性严重程度的一种度量。它表示回归系数估计量的方差与假设自变量间不线性相关时方差相比的比值。
使用pasty的模型函数来得到子列，再将它实例化（np.ndarray），类似于dataframe选择子列以后再使用.values属性，不过pasty函数的操作流更加流程化，且高效。而且pasty的结果可以作为参数传入最小二乘法函数，数据分析中用的比较多。
patsy.dmatirces的return_type默认参数是Design_Matrix。这里使用pasty的dmatrices函数将因变量PE，自变量AT，V，AP和截距项（值为1的1为数组）以数据框的形式组合起来。这里导入了statsmodels模块中提供计算VIF的函数。
import patsy
from statsmodels.stats.outliers_influence import variance_inflation_factor
y, X = patsy.dmatrices('PE~AT+V+AP', data=ccpp, return_type='dataframe')
# formula +0去掉截距项
# y, X = patsy.dmatrices('PE~AT+V+AP+0', data=ccpp, return_type='dataframe')
print(y)
print(X)
vif = pd.DataFrame()
vif['VIF Factor'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif['features'] = X.columns
print(vif)
结果显示，以上所有的VIF值都小于10，所以各自变量之间并不存在多重共线性。
3、异常点检测
在异常点监测之前，我们需要对现有的数据，进行线性回归模型的构造。具体代码如下：
# 构造PE与AT、V和AP之间的线性模型
fit = sm.formula.ols('PE~AT+V+AP', data=ccpp).fit()
print(fit.summary())
# 计算模型的RMSE值
pred = fit.predict()
print(np.sqrt(mean_squared_error(ccpp.PE, pred)))计算异常点监测之前的RMSE值：
开始异常点监测，关于异常点的监测方法，一般通过高杠杆值点（帽子矩阵）或DFFITS值、学生话残差、cook距离和covratio值来判断。这些值的具体计算脚本如下：
（1）当高杠杆值点（或帽子矩阵）> 2 ( p + 1) / n 时，则认为该样本可能存在异常（其中P为自变量的个数，n为观测的个数）；
（2）当DFFITS统计值 > 2sqrt ( ( p + 1) / n) 时, 则认为该样本可能存在异常；
（3）当学生化残差的绝对值大于2，则认为该样本点可能存在异常；
（4）对于cook距离来说，则没有明确的判断标准，一般来说，值越大则认为该样本存在异常的可能性越高；
（5）对于covratio值来说，如果一个样本的covratio值离数值1越远，则认为该样本存在异常的可能性越高。
outliers = fit.get_influence()
leverage = outliers.hat_matrix_diag
dffits = outliers.dffits[0]
resid_stu = outliers.resid_studentized_external
cook = outliers.cooks_distance[0]
covratio = outliers.cov_ratio
contat1 = pd.concat([pd.Series(leverage, name='leverage'),
pd.Series(dffits, name='dffits'),
pd.Series(resid_stu, name='resid_stu'),
pd.Series(cook, name='cook'),
pd.Series(covratio, name='covratio'), ], axis=1)
ccpp_outliers = pd.concat([ccpp, contat1], axis=1)
print(ccpp_outliers.head(5))
这里我们以学生化残差作为评判标准，因为其包含了帽子矩阵和DFFITS的信息。并且我们将3.7%的异常值作删除处理。
outliers_ratio = sum(np.where(np.abs(ccpp_outliers['resid_stu']) > 2, 1, 0)) / ccpp_outliers.shape[0]
print(outliers_ratio)
ccpp_outliers = ccpp_outliers.loc[np.abs(ccpp_outliers.resid_stu) <= 2, ]
结果显示，样本确实存在异常值点，且异常值的数量占了3.7%。对于异常值的处理，我们来考虑下面三种方法：
（1）当异常比例极低时（如5%以内），可以考虑直接删除；
（2）当异常比例较高时，可以考虑将异常值衍生为哑变量，即异常值对应到1，非异常值则对应到0；
（3）单独把异常值提取出来，另行建模。
删除异常值后，我们重新建模，对比删除异常值前的模型，删除后的模型效果更好，具体表现为：信息准则（AIC和BIC）均变小，同时RMSE（误差均方根）也由原来的4.89降低到4.26。
fit2 = sm.formula.ols('PE~AT+V+AP', data=ccpp_outliers).fit()
print(fit2.summary())
pred2 = fit2.predict()
print(np.sqrt(mean_squared_error(ccpp_outliers['PE'], pred2)))
4、正态性检验
当模型的残差服从正态性假设时，才能保证模型偏回归系数对应的t值和模型的F值是有效的。古模型建好后，要对模型的残差进行正态性检验。关于正态性检验由两类方法，即定性的图形法（直方图、PP图和QQ图）和定量的非参数法（Shapiro检验和K-S检验）。
resid函数会返回一个向量,其中包含了每个观测值的残差。
stats.norm.pdf(X，mu，sigma)，一个计算正态分布函数的概率密度函数的函数。
from scipy.stats import norm
from matplotlib import mlab
# 残差的正态性检验（直方图法）
resid = fit2.resid
plt.rcParams['font.sans-serif'] =['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
plt.hist(resid,
density=True,
bins=100,
color='steelblue',
edgecolor='k',
stacked=True)
plt.title('残差直方图')
plt.ylabel('密度值')
# 生成正态曲线的数据
x1 = np.linspace(resid.min(),resid.max(), 1000)
normed = norm.pdf(x1, resid.mean(), resid.std())
# 绘制正态分布曲线
plt.plot(x1, normed, 'r-', linewidth=2, label='正态分布曲线')
# 生成核密度曲线的数据
x2 = np.linspace(resid.min(), resid.max(), 1000)
kde = mlab.GaussianKDE(resid)
# 绘制核密度曲线
plt.plot(x2, kde(x2), 'k-', linewidth=2, label='核密度曲线')
plt.legend(loc='best')
plt.show()
由上图可以看到，核密度曲线与理论的正态分布曲线基本吻合，残差基本服从正态分布。
# 残差的正态性检验（PP图和QQ图法）
pp_qq_plot = sm.ProbPlot(resid)
pp_qq_plot.ppplot(line='45')
plt.title('P-P图')
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q图')
plt.show()
从PP图和QQ图来看，有一部分样本点并没有落在参考线上，但绝大多数样本点还是与参考线保持一致的步调。
from scipy import stats
standard_resid = (resid-np.mean(resid))/np.std(resid)
print(stats.kstest(standard_resid, 'norm'))
由于shapiro正态性检验对样本量的要求是5000以内；而本次数据集的样本量有
lamd = stats.boxcox_normmax(ccpp_outliers['PE'], method='mle')
ccpp_outliers['PE'] = stats.boxcox(ccpp_outliers['PE'], lamd)
fit3 = sm.formula.ols('PE~AT+V+AP', data=ccpp_outliers).fit()
print(fit3.summary())5、残差方差齐性检验
在线性回归模型中，如果模型表现的非常好的话，那么残差与拟合值之间不应该存在某些明显的关系或趋势。如果模型的残差确实存在一定的异方差的话，会导致估计出来的偏回归系数不具备有效性，甚至导致模型的预测也不准确。所以，建模后需要验证残差方差是否具有齐性，检验的方法有两种，一种是图示法，一种是统计验证法。
所谓方差齐性，也就是方差相等，在t检验和方差分析中，都需要满足这一前提条件。
在两组和多组比较中，方差齐性的意思很容易理解，无非就是比较各组的方差大小，看看各组的方差是不是差不多大小，如果差别太大，就认为是方差不齐，或方差不等。如果差别不大，就认为方差齐性或方差相等。当然，这种所谓的差别大或小，需要统计学的检验，所以就有了方差齐性检验。
判断：散点随机散布、无规律则表明满足基本假设，有明显规律或者呈现一定趋势，则有异方差性。
plt.style.use('ggplot')
# 标准化残差与预测值之间的散点图
plt.scatter(fit2.predict(), (fit2.resid - fit2.resid.mean()) / fit2.resid.std(),color='b', s=20)
plt.xlabel('预测值')
plt.ylabel('标准化残差')
# 添加水平参考线
plt.axhline(y=0, color='r', linewidth=2)
plt.show()
从图中看，并没有发现明显的规律或趋势（判断标准：如果残差在参考线两侧均匀分布，则意味着异方差性较弱；而如果呈现出明显的不均匀分布，则意味着存在明显的异方差性。），故可以认为没有显著的异方差性特征。
# ===========统计法完成方差齐性的判断===============
# White's Test
print(sm.stats.diagnostic.het_white(fit2.resid, exog=fit2.model.exog))
# Breusch-Pagan
print(sm.stats.diagnostic.het_breuschpagan(fit2.resid, exog_het=fit2.model.exog))
从检验结果来看，不论是White检验还是Breush-Pagan检验，P值都远远小于0.05这个判别界限，即拒绝原假设（残差方差为常数的原假设），认为残差并不满足齐性这个假设。如果模型的残差确实不服从齐性的话，可以考虑两类方法来结果，一种是模型变换法，另一种是加权最小二乘法。这里演示一下加权最小二乘法。
# 三种权重
w1 = 1 / np.abs(fit2.resid)
w2 = 1 / fit2.resid ** 2
ccpp_outliers['loge2'] = np.log(fit2.resid ** 2)
# 第三种权重
model = sm.formula.ols('loge2~AT+V+AP', data=ccpp_outliers).fit()
w3 = 1 / (np.exp(model.predict()))
# 建模
fit3 = sm.formula.wls('PE~AT+V+AP', data=ccpp_outliers, weights=w1).fit()
# 异方差检验
het3 = sm.stats.diagnostic.het_breuschpagan(fit3.resid, exog_het=fit3.model.exog)
# AIC
fit3.aic
fit4 = sm.formula.wls('PE~AT+V+AP', data=ccpp_outliers, weights=w2).fit()
het4 = sm.stats.diagnostic.het_breuschpagan(fit4.resid, exog_het=fit4.model.exog)
fit4.aic
fit5 = sm.formula.wls('PE~AT+V+AP', data=ccpp_outliers, weights=w3).fit()
het5 = sm.stats.diagnostic.het_breuschpagan(fit5.resid, exog_het=fit5.model.exog)
fit5.aic
# fit2模型
het2 = sm.stats.diagnostic.het_breuschpagan(fit2.resid, exog_het=fit2.model.exog)
fit2.aic
print('fit2模型异方差检验统计量：%.2f，P值为%.4f：' % (het2[0], het2[1]))
print('fit3模型异方差检验统计量：%.2f，P值为%.4f：' % (het3[0], het3[1]))
print('fit4模型异方差检验统计量：%.2f，P值为%.4f：' % (het4[0], het4[1]))
print('fit5模型异方差检验统计量：%.2f，P值为%.4f：\n' % (het5[0], het5[1]))
print('fit2模型的AIC：%.2f' % fit2.aic)
print('fit3模型的AIC：%.2f' % fit3.aic)
print('fit4模型的AIC：%.2f' % fit4.aic)
print('fit5模型的AIC：%.2f' % fit5.aic)
通过对比发现，尽管我们采用了三种不同的权重，但都没能通过残差方差齐性的显著性检验，但似乎fit4模型更加理想，相比于fit2来说，AIC信息更小（当然也可能产生过拟合问题）。
6、残差独立性检验
之所以要求残差是独立的，是因为要求因变量y是独立的，因为在模型中只有y和残差项是变量，而自变量x是已知的。如果再配上正态分布的假设，那就是独立同布于正态分布，关于残差的独立性检验我么可以通过Durbin-Watson统计量来测试。其实，在模型的summary信息中就包含了残差的Durbin-Watson统计量值，如果该值越接近于2，说明残差是独立的。一般而言，在实际的数据集中，时间序列的样本之间可能会存在相关性，而其他数据集样本之间基本还是独立的。
下面对5个模型产生的预测值和实际值作散点图，如果散点图与预测线特别紧密，则认为模型拟合的非常棒。
model_list = {'fit': fit,
'fit2': fit2,
'fit3': fit3,
'fit4': fit4,
'fit5': fit5}
for i, model in enumerate(model_list.items()):
plt.subplot(2, 3, i+1)
if i==0:
plt.scatter(model[1].predict(), ccpp['PE'], color='b', edgecolors='white', s=20)
plt.plot([model[1].predict().min(), model[1].predict().max()],
[ccpp.PE.min(), ccpp.PE.max()],
'r-', linewidth=3)
else:
plt.scatter(model[1].predict(), ccpp_outliers['PE'], color='b', edgecolors='white', s=20)
plt.plot([model[1].predict().min(), model[1].predict().max()],
[ccpp_outliers.PE.min(), ccpp_outliers.PE.max()],
'r-', linewidth=3)
plt.title('模型：%s' % model[0])
plt.xlabel('预测值')
plt.ylabel('实际值')
plt.show()
