多元线性回归模型检验和预测

一、概述

（F检验）显著性检验：检测自变量是否真正影响到因变量的波动。

（t检验）回归系数检验：单个自变量在模型中是否有效。

二、回归模型检验

检验回归模型的好坏常用的是F检验和t检验。F检验验证的是偏回归系数是否不全为0（或全为0），t检验验证的是单个自变量是否对因变量的影响是显著的（或不显著）。

F检验和t检验步骤：

• 提出问题的原假设和备择假设
• 在原假设的条件下，构造统计量
• 根据样本信息，计算统计量的值
• 对比统计量的值和理论F分布的值，计算统计量的值超过理论值，则拒绝原假设，否则接受原假设

待检验数据集先构造成向量的模式：

可表示成如下形式：

　　                   

其中β 为n×1的一维向量。

1) 假设

F检验假设：

 

 t检验假设：

H0为原假设，H1为备择假设。F检验拒绝原假设的条件为计算的F检验的值大于查到的理论F值。t检验可以通过P值和拟合优度判断变量的显著性及模型组合的优劣。

2）计算过程

F检验计算过程：

 上图为假设其中一个点所在的平面，由以上点计算出ESS(误差平方和)，RSS（回归离差平方和）,TSS（总的离差平方和）。

计算公式为：

其中ESS和RSS都会随着模型的变化而发生变化（估计值变动）。而TSS衡量的是因变量和均值之间的离差平方和，不会随着模型的变化而变化。

由以上公式构造F统计量：

 

 （n为数据集向量的行数，p为列数（自变量的个数），n为离差平方和RSS的自由度，n-p-1为误差平方和ESS的自由度），模型拟合越好，ESS越小，RSS越大，F值越大。

3）python中计算F值得方法：可以直接通过model.fvalue得到f值或者通过计算得到。

from sklearn import model_selection
import statsmodels.api as sm
import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt
from sklearn import preprocessing

Profit = pd.read_excel(r'Predict to Profit.xlsx')
Profit.head()

dummies = pd.get_dummies(Profit.State,prefix = 'State')
Profit_New = pd.concat([Profit,dummies],axis=1)
Profit_New.drop(labels = ['State','State_New York'],axis =1,inplace = True)

train , test = model_selection.train_test_split(Profit_New,test_size = 0.2,random_state=1234)
model = sm.formula.ols('Profit~RD_Spend+Administration+Marketing_Spend+State_California+State_Florida',data = train).fit()

ybar = train.Profit.mean()
p = model.df_model  #自变量个数
n= train.shape[0]    #行数，观测个数
RSS = np.sum((model.fittedvalues - ybar)**2)  #计算离差平方和 估计值model.fittedvalues
ESS = np.sum(model.resid**2)  #计算误差平方和，误差表示model.resid
F = RSS/p/(ESS/(n-p-1))
print( F, model.fvalue)

　　结果如下：

 求理论F值，使用Scipy库子模块stats中f.ppf

from scipy.stats import f
F_Theroy = f.ppf(q=0.95,dfn = p,dfd = n-p-1)
F_Theroy

　　

 由于计算的F值远远大于理论F值，所以拒绝原假设，证明多元线性回归方程是显著的，偏回归系数不全为0，即所有自变量联合起来的组合确实对利润有显著性影响。

4) t检验手工计算比较复杂，套用python的验计算方式，model.summary()

 参数：R-squared为拟合优度（反映模型对样本的拟合程度）,Adj.R-squared为修正的可决系数。

 通过t检验得出，只有RD_Spend的P值小于0.05，其余变量均大于0.05，说明其余变量不是影响利润Profit的重要因素。

三、回归模型诊断

实际应用中，因变量若为数值型变量，可以考虑线性回归模型。模型需满足如下假设：

　　因变量（残差）服从正态分布，自变量之间不存在多重共线性，自变量和因变量之间存在线性关系，用于建模的数据不存在异常值，残差项满足方差异性和独立性。

• 正态性检验（直方图，pp图或qq图，shapiro检验和K-S检验）
• 多重共线性检验（VIF>10，存在多重共线性，>100则变量间存在严重的多重共线性）
• 线性相关性（数据框的corrwith方法或者seaborn模块的pairplot函数作图观察）|p|~>=0.8 高度相关  ，0.5《 |p|《0.8 中度相关，0.3《 |p|《0.5 弱相关，|p|《0.3 几乎不相关
• 异常值检验（帽子矩阵、DFFITS准则、学生化残差、Cook距离）
• 残差独立性检验（summary()方法观察Durbin-Watson值），DW值在2左右，则残差项之间不相关，偏离2较远，则不满足残差独立性假设
• 方差齐性检验（图形法和BP检验）

具体方法过程：

1）正态性检测（模型的前提假设是残差服从（0，1）正态分布，实质上由方程y=Xβ+ε，因变量也应该服从正态分布），检测的是因变量是否符合正态分布

• 直方图  scipy.stats seaborn.distplot

import seaborn as sns
import scipy.stats as stats
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

sns.distplot(a = Profit_New.Profit,
            bins =10,
            fit = stats.norm,
            norm_hist = True,
            hist_kws = {'color':'green','edgecolor':'black'},
            kde_kws = {'color':'black','linestyle':'--','label':'核密度曲线'},
            fit_kws = {'color':'red','linestyle':':','label':'正态密度曲线'}
           )
plt.legend()
plt.show()

　　

 通过观察核密度曲线和理论正态密度曲线的拟合程度，2条曲线近似吻合，可认为因变量Profit服从正态分布。

• PP图和QQ图  ppplot qqplot

pq_plot = sm.ProbPlot(Profit_New.Profit)
pq_plot.ppplot(line='45')
plt.title('PP图')
pq_plot.qqplot(line='q')
plt.title('QQ图')
plt.show()

　　

 

 观察到散点均均匀的分布在直线上，说明因变量近似服从正态分布。

• Shapiro检验和K-S检验（因变量数据量>5000 use K-S检验法，<5000 use Shapiro检验）

stats.shapiro(Profit_New.Profit)

　　输出为

(0.9793398380279541, 0.537902295589447)，第一个值为shapiro检验的统计量值，第二个值为对应概率p值,p值大于0.05置信水平，利润变量服从正态分布假设。

K-S检验 KSValue = stats.kstest(rvs = datay,args = (datay.mean(),datay.std(),cdf = 'norm'))
需要注意的是KS检验，必须指定args参数，来传递检验变量的均值和标准差。

2）多重共线性检验 （自变量之间的检测）

检测模型中自变量之间存在较高的线性相关关系，若存在会给模型带来严重的后果。其检验可以通过方差膨胀因子VIF来鉴定。（VIF>10，存在多重共线性，>100则变量间存在严重的多重共线性）

消除多重共线性，可以通过不同自变量之间的组合观察vif值，确定最佳组合。R语言中消除多重共线性，确定最优组合使用模型的step方法（逐步回归）。

计算方式：

X1为第一个自变量。构造x1与其余自变量的线性回归模型，得到回归模型的判决系数R2，得到第一个自变量的膨胀因子，依次类推计算剩余变量的VIF。

Python 的statsmodels模块可直接计算VIF值

from statsmodels.stats.outliers_influence import variance_inflation_factor
X = sm.add_constant(Profit_New.ix[:,['RD_Spend','Marketing_Spend']])
vif = pd.DataFrame()
vif['features'] = X.columns
vif["VIF Factor"] = [variance_inflation_factor(X.values,i) for i in range(X.shape[1])]

　　

得到vif的值：

 计算所有的变量:

X = sm.add_constant(Profit_New.drop('Profit',axis =1).ix[:])
vif = pd.DataFrame()
vif['features'] = X.columns
vif["VIF Factor"] = [variance_inflation_factor(X.values,i) for i in range(X.shape[1])]
vif

　　

vif值:

 第二种模型得到的const的值大于10，说明构建模型的变量之间存在多重共线性，需要删除变量重新选择模型。第一种模式选择2个自变量（'RD_Spend','Marketing_Spend'）是符合自变量之间无线性相关要求的。

3）线性相关性检验 （各变量之间线性相关性检测）

• 图形法 seaborn pariplot函数

sns.pairplot(Profit_New.ix[:,['RD_Spend','Administration','Marketing_Spend','Profit']]) #去除哑变量
plt.show()

　　散点图矩阵图示：

 观察散点图分布得出，RD_Spend，Marketing_Spend两个变量和因变量Profit之间存在线性关系，Administration与Profit之间散点图呈水平趋势，说明2者之间不相关。

• 数据框的corrwith方法(该方式优点是可以计算任意指定变量之间的相关系数)

判断条件：

# Profit 和各变量之间的相关系数
Profit_New.drop(['Profit','State_California','State_Florida'],axis =1).corrwith(Profit_New.Profit)

　　

 

 结果看出：自变量RD_Spend和Marketing_Spend与因变量之间存在较高的线性相关性，而其他变量与利润之间线性相关性比较弱。可忽略其线性相关关系。

需要注意的是通过corrwith检测变量之间不存在线性相关关系，结合seaborn.pariplot观察，可能存在其他关系如2次方，3次方，倒数关系等等（在变量通过t检验被确认为对因变量影响是显著的情况下，还需要对模型再次检查）。

通过如下方式重新确定模型为： Profit = 51902.112471 + 0.785116RD_Spend +  0.019402Marketing_Spend 

model3 = sm.formula.ols('Profit~RD_Spend + Marketing_Spend ',data = train).fit()
model3.params

　　

 4）异常值检测  （帽子矩阵、DFFITS准则、学生化残差、Cook距离）

确认方式：|学生化残差| > 2，则对应的数据点为异常值。

对异常值得处理：异常值检测 model.get_influence() ， 学生化残差值  model.get_influence().resid_studentized_external

确认异常样本的比例，比例<5%，考虑可以直接删除异常值，>5%，可以生成哑变量，对于异常点，设置哑变量为1，否则为0。

# 异常值检测
outliers = model3.get_influence()   

# 帽子矩阵
leverage = outliers.hat_matrix_diag
# dffits值
dffits= outliers.dffits[0]
# 学生化残差
resid_stu = outliers.resid_studentized_external
# cook距离
cook = outliers.cooks_distance[0]


# 合并各种异常值检验的统计量值
contatl = pd.concat([pd.Series(leverage, name = 'leverage'),
                     pd.Series(dffits, name = 'dffits'),
                     pd.Series(resid_stu, name = 'resid_stu'),
                     pd.Series(cook, name = 'cook')
                    ],axis =1 )

train.index = range(train.shape[0])
profit_outliers = pd.concat([train,contatl ],axis =1)
profit_outliers.head()

　

 使用非异常值的数据点进行建模：

# 计算异常值比例
outliers_ratio = np.sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/profit_outliers.shape[0]
outliers_ratio  # 0.02564

none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu)<=2,]
model4 = sm.formula.ols('Profit~RD_Spend+Marketing_Spend',data = none_outliers).fit()
model4.params

再次更新模型为：Profit = 51827.416821 + 0.797038RD_Spend + 0.01774Marketing_Spend

5）残差独立性检验

通过mdoel.summary函数观察Durbin-Watson的值，在2左右说明残差之间是不相关的，偏离2较远，说明残差之间是不独立的。

model4.summary()，DW的值为2.065，说明模型残差满足独立性假设。

　

6）残差方差齐性检验 （图形法和BP法）

 目的： 残差项的方差不随自变量的变动而呈现某种趋势，即残差项的方差不受自变量变化而影响

 图形法：绘制残差和自变量之间的散点图

ax1 = plt.subplot2grid(shape = (2,1) , loc = (0,0)) # 设置第一张子图位置
# 散点图绘制
# ax1.scatter(none_outliers.RD_Spend, none_outliers.resid_stu )  #学生化残差与自变量散点图
ax1.scatter(none_outliers.RD_Spend, (model4.resid - model4.resid.mean())/model4.resid.std()) #标准化残差和自变量散点图
# 添加水平参考线
ax1.hlines(y = 0 ,
          xmin = none_outliers.RD_Spend.min(),
           xmax = none_outliers.RD_Spend.max(),
           color = 'red',
           linestyle = '--'
          )
ax1.set_xlabel('RD_Spend')
ax1.set_ylabel('Std_Residual')
ax2 = plt.subplot2grid(shape = (2,1) , loc = (1,0))
# ax1.scatter(none_outliers.Marketing_Spend, none_outliers.resid_stu )
ax2.scatter(none_outliers.Marketing_Spend, (model4.resid - model4.resid.mean())/model4.resid.std())
ax2.hlines(y = 0 ,
          xmin = none_outliers.Marketing_Spend.min(),
           xmax = none_outliers.Marketing_Spend.max(),
           color = 'magenta',
           linestyle = '--'
          )
ax2.set_xlabel('Marketing_Spend')
ax2.set_ylabel('Std_Residual')

# 调整2子图之间距离
plt.subplots_adjust(hspace = 0.6,wspace = 0.3)
plt.show()

　注意的是采用学生化残差或标准化残差绘制和自变量的散点图均可以，2者图形类似。

原因是学生化残差计算公式为：

 图形结果如下：

 根据结果可知标准化残差并没有随着自变量变化而呈现某种趋势，图形中所有点均匀的分布在y=0两侧，可以认为残差项满足方差齐性假设。

方差齐性检测完整代码：

from sklearn import model_selection
import statsmodels.api as sm
import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt

Profit = pd.read_excel(r'Predict to Profit.xlsx')

dummies = pd.get_dummies(Profit.State,prefix = 'State')
Profit_New = pd.concat([Profit,dummies],axis=1)
Profit_New.drop(labels = ['State','State_New York'],axis =1,inplace = True)
train , test = model_selection.train_test_split(Profit_New,test_size = 0.2,random_state=1234)
model3 = sm.formula.ols('Profit~RD_Spend + Marketing_Spend ',data = train).fit()
# 异常值检测
outliers = model3.get_influence()   

# 帽子矩阵
leverage = outliers.hat_matrix_diag
# dffits值
dffits= outliers.dffits[0]
# 学生化残差
resid_stu = outliers.resid_studentized_external
# cook距离
cook = outliers.cooks_distance[0]
contatl = pd.concat([pd.Series(leverage, name = 'leverage'),
                     pd.Series(dffits, name = 'dffits'),
                     pd.Series(resid_stu, name = 'resid_stu'),
                     pd.Series(cook, name = 'cook')
                    ],axis =1 )

train.index = range(train.shape[0])
profit_outliers = pd.concat([train,contatl ],axis =1)

outliers_ratio = np.sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/profit_outliers.shape[0]

none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu)<=2,]
model4 = sm.formula.ols('Profit~RD_Spend+Marketing_Spend',data = none_outliers).fit()

ax1 = plt.subplot2grid(shape = (2,1) , loc = (0,0)) # 设置第一张子图位置
# 散点图绘制
# ax1.scatter(none_outliers.RD_Spend, none_outliers.resid_stu )
ax1.scatter(none_outliers.RD_Spend, (model4.resid - model4.resid.mean())/model4.resid.std())
# 添加水平参考线
ax1.hlines(y = 0 ,
          xmin = none_outliers.RD_Spend.min(),
           xmax = none_outliers.RD_Spend.max(),
           color = 'red',
           linestyle = '--'
          )
ax1.set_xlabel('RD_Spend')
ax1.set_ylabel('Std_Residual')
ax2 = plt.subplot2grid(shape = (2,1) , loc = (1,0))
# ax1.scatter(none_outliers.Marketing_Spend, none_outliers.resid_stu )
ax2.scatter(none_outliers.Marketing_Spend, (model4.resid - model4.resid.mean())/model4.resid.std())
ax2.hlines(y = 0 ,
          xmin = none_outliers.Marketing_Spend.min(),
           xmax = none_outliers.Marketing_Spend.max(),
           color = 'magenta',
           linestyle = '--'
          )
ax2.set_xlabel('Marketing_Spend')
ax2.set_ylabel('Std_Residual')

# 调整2子图之间距离
plt.subplots_adjust(hspace = 0.6,wspace = 0.3)
plt.show()

BP法（拉格朗日乘子检验）：statsmodels模块het_breuschpagan函数  https://programtalk.com/python-examples/statsmodels.stats/

sm.stats.diagnostic.het_breuschpagan(model4.resid, exog_het = model4.model.exog)　

结果得到4个值

(1.4675103668308342, 0.48010272699006384, 0.7029751237162462, 0.5019659740962872)

第一个值为LM统计量，第二个为p值，第三个为F统计量，第四个为F统计量的p值，通过观察p值均大于0.05，证明接受残差方差齐性原假设，即残差不受自变量的影响而变化。

若果残差不满足方差齐性假设，观察残差和自变量的关系，通过模型变换或加权最小二乘法对模型加以处理。

四、回归模型预测

pred4 = model4.predict(exog = test.ix[:,['RD_Spend','Marketing_Spend']])
plt.scatter(x = test.Profit ,y = pred4)
plt.plot([ test.Profit.min(),test.Profit.max()],
         [test.Profit.min(),test.Profit.max()],
        color = 'red',
         linestyle = '--'
)
plt.xlabel('实际值')
plt.ylabel('预测值')
plt.show()　　

　真实值和预测值之间差异：

print(pd.DataFrame({"prediction":pred4,'real':test.Profit}))

　　

图形展示效果：