多元线性回归

---

原理解析

多元线性回归通俗理解就是一个因变量与多个自变量之间的相关关系

---

scikit-learn中的线性回归使用

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn import datasets

# 这两行代码在画图时添加中文必须用
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 加载数据（数据为490行十四列,最后一列为y值（二分类时为0/1））
boston = datasets.load_boston()

X = boston.data
y = boston.target  
# 筛选数据集
X = X[y<50.0]
y = y[y<50.0]

X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=666)

# 多元线性回归
reg = LinearRegression()
reg.fit(X_train,y_train)

print(reg.coef_) # 回归系数
print(reg.intercept_)  # 公式中的误差项 W0（偏置）
# 判定系数R2的意义是由 X 引起的影响占总影响的比例来判断拟合程度
# reg.score效果和r2_score(y_text, y_pred))一样,即R方
print(reg.score(X_test,y_test)) # 模型准确度（拟合程度）

# 模型预测
y_pred = reg.predict(X_test)
print(y_pred)

# 估计标准误差评价
# (1) 评价测度
# 对于分类问题，评价测度是准确率，但这种方法不适用于回归问题。
# 我们使用针对连续数值的评价测度(evaluation metrics)。
# 这里介绍3种常用的针对线性回归的测度。
# (1)平均绝对误差(Mean Absolute Error, MAE)
# (2)均方误差(Mean Squared Error, MSE)
# (3)均方根误差(Root Mean Squared Error, RMSE)
# 这里我使用RMES。
rmse = sqrt(mean_squared_error(y_test, y_pred))  # 计算均方根误差
print(rmse)

# 做ROC曲线
plt.figure()
plt.plot(range(len(y_pred)), y_pred, 'b', label="predict")
plt.plot(range(len(y_pred)), y_test, 'r', label="test")
plt.legend(loc="upper right")  # 显示图中的标签
plt.xlabel("the number of sales")
plt.ylabel('value of sales')
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

---

用 statsmodels 包（不需要划分数据集）

import sys
import io  # 解决输出乱码问题
sys.stdout = io.TextIOWrapper(sys.stdout.buffer,encoding='gb18030')
import statsmodels.api as sm  # 可以解决单变量 / 多变量线性回归
# 模型可以做单变量/多变量
# y为因变量,X为多个自变量但需将X加入sm模型中
x=sm.add_constant(X)  # 若模型中有截距，必须有这一步
model = sm.OLS(y, x).fit()  # 生成模型并拟合模型
print(model.summary())     # 统计结果
predicts = model.predict() # 模型的预测值

plt.figure()
plt.plot(range(len(y)), predicts, 'b', label="predict")
plt.plot(range(len(y)), y, 'r', label="test")
plt.legend(loc="upper right")  # 显示图中的标签
plt.show()

datas = pd.read_excel(r'D:\datas\linear_regression.xlsx') 
model = sm.ols('不良贷款~各项贷款余额', datas).fit()
# 构建最小二乘模型并拟合，此时不用单独输入 x，y了，
# 而是将自变量与因变量用统计语言公式表示，将全部数据导入
print(model.summary()) # 输出回归结果
predicts = model.predict() # 模型的预测值
y = datas.iloc[:, 1] # 因变量为第 2 列数据
x = datas.iloc[:, 2] # 自变量为第 3 列数据
plt.scatter(x, y, label='实际值')
plt.plot(x, predicts, color = 'red', label='预测值')
plt.legend() 
plt.show() 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

---

用 scipy 包（scipy 目前只能做一元线性回归，也不能用来做预测）

import scipy.stats as st
import pandas as pd

datas = pd.read_excel(r'D:\datas\linear_regression.xlsx') 
y = datas.iloc[:, 1] # 因变量为第 2 列数据
x = datas.iloc[:, 2] # 自变量为第 3 列数据

# 线性拟合，可以返回斜率，截距，r 值，p 值，标准误差
slope, intercept, r_value, p_value, std_err = st.linregress(x, y)

print(slope)# 输出斜率
print(intercept) # 输出截距
print(r_value**2) # 输出 r^2

1
2
3
4
5
6
7
8
9
10
11
12
13
14

---

线性回归的显著性检验

1. 线性关系检验（ X 影响占比越高，预测效果就越好 ）

     结论：在多元线性回归中，只要有一个自变量系数不为零（即至少一个自变量系数与因变量有线性关系），我们就说这个线性关系是显著的。如果不显著，说明所有自变量系数均为零。

2. 回归系数检验（检验统计量使用t分布）

3. Python代码实现

    因为使用statsmodels包里的ols模块进行上述两种显著性检验很方便，所以直接对其进行分析。 （通过 model.summary()查看）

通过上面结果我们清楚看到：

• F统计量的p值非常小，拒绝原假设，说明线性关系显著
• 两个回归系数的t统计量p值均为0，拒绝原假设，说明回归系数也都显著

---

线性回归的诊断

（1）残差分析                （2）线性相关性检验
（3）多重共线性分析     （4）强影响点分析

<1>残差分析

1. 正态性检验

    干扰项（即残差），服从正态分布的本质是要求因变量服从变量分布。 因此，验证残差是否服从正态分布就等于验证因变量的正态分布特性。关于正态分布的检验通常有以下几种方法。

（1）直方图法：
    直方图法就是根据数据分布的直方图与标准正态分布对比进行检验，主要是通过目测。

import sys
import io  # 解决输出乱码问题
sys.stdout = io.TextIOWrapper(sys.stdout.buffer,encoding='gb18030')
from scipy.stats import *
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
import statsmodels.api as sm  # 可以解决单变量 / 多变量线性回归

x=sm.add_constant(X)  # 若模型中有截距，必须有这一步
model = sm.OLS(y, x).fit()  # 生成模型并拟合模型

residual = model.resid  # 模块下的残差函数
sns.distplot(residual,
             bins = 13,
             kde = False,
             color = 'blue',
             fit = norm)  # 拟合标准正态分布
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

（2）PP图和QQ图：

    PP图是对比正态分布的累积概率值和实际分布的累积概率值。

    QQ图是通过把测试样本数据的分位数与已知分布相比较，从而来检验数据的分布情况。对应于正态分布的QQ图，就是由标准正态分布的分位数为横坐标，样本值为纵坐标的散点图。

# pp图法
pq = sm.ProbPlot(residual)
pq.ppplot(line='45')

# qq图法
pq = sm.ProbPlot(residual)
pq.qqplot(line='q')
1
2
3
4
5
6
7

    pp图和qq图判断标准是：如果观察点都比较均匀的分布在直线附近，就可以说明变量近似的服从正态分布，否则不服从正态分布。（若样本点并不是十分均匀地落在直线上，有的地方有一些较大的偏差，则判断不是正态分布）

（3）Shapiro检验：

    这种检验方法均属于非参数方法，先假设变量是服从正态分布的，然后对假设进行检验。一般的数据量低于5000则可以使用Shapiro检验，大于5000的数据量可以使用K-S检验，这种方法在scipy库中可以直接调用：

# shapiro检验
import scipy.stats as stats
stats.shapiro(residual)
out:
(0.9539670944213867, 4.640808128e-06)
1
2
3
4
5

    上面结果两个参数：第一个是Shaprio检验统计量值，第二个是相对应的p值。可以看到，统计量W接近于1，但是 p值非常小，远远小于0.05，因此拒绝原假设，说明残差不服从正态分布。

2. 独立性检验 （通过 model.summary()查看）

3. 方差齐性检验

（1）BP检验法

# BP检验 
# residual: residual = model.resid  # 残差函数获得残差
# model = sm.OLS(y, x).fit()  # 生成模型并拟合模型
# model.model.exog 
bptest=sm.stats.diagnostic.het_breuschpagan(residual,model.model.exog)
out:
(0.16586685109032384, 0.6838114989412791,
 0.1643444790856123, 0.6856254489662914)
"""
上述参数:
第一个为:LM统计量值

第二个为: 响应的p值,0.68远大于显著性水平0.05,
因此接受原假设,即残差方差是一个常数

第三个为: F统计量值,用来检验残差平方项与自变量之间是否独立,
如果独立则表明残差方差齐性

第四个为: F统计量对应的p值,也是远大于0.05的,
因此进一步验证了残差方差的齐性。
"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

<3>多重共线性检验

1. 多重共线性产生的问题

    当回归模型中两个或两个以上的自变量彼此相关时，则称回归模型中存在多重共线性，也就是说共线性的自变量提供了重复的信息。

    那么这种多重共线性会有什么不好的影响吗？答案是会的，而且影响非常不好。总结一下就是：会造成回归系数，截距系数的估计非常不稳定，即整个模型是不稳定。

    这种不稳定的具体表现是：很可能回归系数原来正，但因为共线性而变为负。这对于一些自变量的可解释性来讲可能是致命的，因为得到错误系数无法解释正常发生的现象。

2. 多重共线性的检测

如果出现以下情况，可能存在多重共线性：

（1）模型中各对自变量之间显著性相关。

（2）当模型线性关系（F检验）显著时，几乎所有回归系数的 t 检验不显著。

（3）回归系数的正负号与预期的相反。

（4）方差膨胀因子（VIF）检测，一般认为VIF大于10，则存在严重的多重共线性。

这里主要说明一下（1）和（4），因为（2）和（3）一般通过观察即可。

(1) 相关系数检验

    主要是用到了DataFrame的corr()方法，默认皮尔逊相关，然后通过seaborn的heatmap可以可视化展示出来。

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 加载数据
df = pd.read_excel('../data/boston.xlsx')
print(df.corr()) # 计算所有的变量的两两相关性
# 只计算选择的两个变量的相关性
print(data['化妆品费'].corr(data['置装费'])) 
sns.heatmap(df,cmap='YlGnBu')
sns.pairplot(df)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14

  计算相关性用到的方法有pearson、spearman、kendall，具体区别如下表所示：

(4) 方差膨胀因子检验

方差膨胀因子的公式如下：VIF = 1 / ( 1 - R^2)

    VIF的公式是基于拟合优度R2的，其中VIF代表自变量X的方差膨胀系数，R代表把自变量X最为因变量，与其他自变量做回归时的R2。

# 自定义VIF方差膨胀因子计算
import statsmodels.formula.api as smf

def vif(df, col_i):
    cols = list(df.columns)
    cols.remove(col_i)
    cols_noti = cols
    formula = col_i + '~' + '+'.join(cols_noti)
    r2 = smf.ols(formula, df).fit().rsquared
    return 1./(1.-r2)

for i in df.columns:
    print(i, '\t', vif(df,col_i=i))

1
2
3
4
5
6
7
8
9
10
11
12
13
14

---

多重共线性的处理方法

多重共线性对于线性回归是种灾难，并且我们不可能完全消除，而只能利用一些方法来减轻它的影响。对于多重共线性的处理方式，有以下几种思路：

（1）提前筛选变量： 可以利用相关检验来或变量聚类的方法。注意：决策树和随机森林也可以作为提前筛选变量的方法，但是它们对于多重共线性帮助不大，因为如果按照特征重要性排序，共线性的变量很可能都排在前面。

（2）子集选择： 包括逐步回归和最优子集法。因为该方法是贪婪算法，理论上大部分情况有效，实际中需要结合第一种方法。

（3）收缩方法： 正则化方法，包括岭回归和LASSO回归。LASSO回归可以实现筛选变量的功能。

（4）维数缩减： 包括主成分回归（PCR）和偏最小二乘回归（PLS）方法。

LASSO回归

from sklearn.linear_model import Lasso,LassoCV,LassoLarsCV
# 加载数据
df = pd.read_excel('../data/boston.xlsx')
# 相关系数查看
print(df.corr())

# 通过设置不同的alpha值建立lasso实例
lasso = Lasso().fit(X_train,y_train)
lasso001 =Lasso(alpha=0.01).fit(X_train,y_train)
print('**********************************')
print("Lasso alpha=1")
print ("training set score:{:.2f}".format(lasso.score(X_train,y_train)))
print ("test set score:{:.2f}".format(lasso.score(X_test,y_test)))
print ("Number of features used:{}".format(np.sum(lasso.coef_!=0)))
1
2
3
4
5
6
7
8
9
10
11
12
13
14

---

强影响点分析

针对于强影响点分析，一般的有以下几种方法：
（1）学生化残差（SR） （2）Cook’s D统计量
（3）DFFITS统计量 （4）DFBETAS统计量

    对于这些指标我们可以通过statsmodels直接查找到，对于我们建立的模型model自动检测每个样本的指标值是多少，我们只需要设置相应的临界点来判断就可以完成检测了。以下是代码实现部分：

# 强离散点各个指标
from statsmodels.stats.outliers_influence import OLSInfluence
import statsmodels.api as sm

model = sm.OLS(yArr,xArr).fit()
OLSInfluence(model).summary_frame().head()

1
2
3
4
5
6
7

    当然，如果我们想单独获取某个指标，我们也可以这样操作：

# 单独获取各个指标
ol = model.get_influence()

leverage = ol.hat_diag_factor
dffits = ol.dffits
resid_stu = ol.resid_studentized_external
cook = ol.cooks_distance

1
2
3
4
5
6
7
8

---