线性回归Statsmodels模型报告（重在Statsmodels.OLS.summary()）

最近看了一下Statsmodels.OLS，即用Statsmodels使用最小二乘法获得线性回归的系数、截距，主要有一个model.summary()，其中有一些参数想深入弄明白，将学习结果分享：如果用python，有很多种方法实现线性回归（带不带常数项截距都无所谓）：从计算原理上来分：一般经常使用正规方程（最小二乘法）和梯度下降。，其本质还是上述两种方法的不同库的实现；

TravelLight92

1573人浏览 · 2024-05-26 20:52:15

TravelLight92 · 2024-05-26 20:52:15 发布

写作原因：最近看了一下Statsmodels.OLS，即用Statsmodels使用最小二乘法获得线性回归的系数、截距，主要有一个model.summary()，其中有一些参数想深入弄明白，将学习结果分享：

如果用python，有很多种方法实现线性回归（带不带常数项截距都无所谓）：

从计算原理上来分：一般经常使用正规方程（最小二乘法）和梯度下降。

比如博客园这篇文章：机器学习算法的Python实现（一）：线性回归 - Geek-Thomas - 博客园 (cnblogs.com)

说了八种方法，其本质还是上述两种方法的不同库的实现；

本文内容：

这里只用最小二乘法，毕竟咱们目的是研究模型报告，但是过程还是要走一遍做做样子，写三种方法的计算实现，带上代码和对比思考：

1.sklearn.linear_model ，里面各种线性回归的库,具体计算有最小二乘、梯度下降、批量、随机梯度下降等；

2.numpy.linalg,也是最小二乘法，相当于借助numpy逐步实现（也算半个手动挡）

3.Statsmodels.OLS，也是最小二乘法，其中的summary()模型报告是本文的重点；

当然还有纯手动挡自己写方法计算的。

一、回顾最小二乘法

原理是通过矩阵形式，求MSE关于系数矩阵的偏导，令其为0，求解：

推导公式比较简单，但是敲latex太长嫌麻烦，此处从上链接那里截个图，该公式随便都能搜到。

梯度下降原理略，本文不讲这个。

二、构造数据与使用四种方法具体过程

我们考虑一个带常数项、二个系数、噪音的数据，因为三个维度可以画出来。

2.1生成数据及展示

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import time
import random
from collections import Counter
import itertools

np.set_printoptions(suppress=True)
plt.rcParams['font.sans-serif'] = ['SimHei']	# 显示中文
plt.rcParams['axes.unicode_minus'] = False		# 显示负号
np.random.seed(1999)
# 生成数据
x1 = np.linspace(1,30,num=100)
x2 = np.linspace(-10,10,num=100)
w1 = 2
w2 = 3
noise = np.random.normal(0,1,size=100)
y= x1*w1 + x2*w2 + noise +5
y_dream = x1*w1 + x2*w2+5
# print(y)
# 先画个3D图
ax = plt.subplot(111,projection='3d')
ax.scatter(x1,x2,y,color='navy',s=10,label='生成的数据点')
ax.plot(x1,x2,y_dream,color='red',lw=2,label='理想的拟合曲线（非最佳）')
ax.set_xlabel('x1',fontsize=20)
ax.set_ylabel('x2',fontsize=20)
ax.set_zlabel('y',fontsize=20)
plt.title('原数据3D展示图')
plt.legend(loc='upper left')
plt.show()

差不多数据就是这样，噪音稍微弄小点，一开始设置标准差设得太大，发现对初始的点影响极大；由于存在噪音，所以我们认为的这个红色的线，极大可能并不是最好的。

2.2sklearn.linear_model

from sklearn.linear_model import LinearRegression
modle_1 = LinearRegression(fit_intercept=True)
X=np.c_[x1,x2]
modle_1.fit(X,y)
coef_1 = modle_1.coef_
intercept_1 = modle_1.intercept_
print('系数是:',coef_1)
print('截距是:',intercept_1)

#系数是: [2.75724514 1.90154838]
#截距是: -6.760579051034362

2.3 手算

此时我们要手动添加一个截距项，系数全部为1，可知我们把系数放在了后面，注意要用np.linalg.pinv，如果用np.linalg.inv可能在矩阵奇异的情况下，解不出来或者与不太准确！

X_intercept = np.c_[x1,x2,np.ones(100)]
# 对应了最小二乘法的几个部分
theta = np.linalg.pinv(X_intercept.T.dot(X_intercept)).dot(X_intercept.T).dot(y)
print(theta)

# array([ 2.58969124,  2.89297533, 17.43027338])

2.4 statsmodels OLS方法

OLS means：（ordinary least square）

import statsmodels.api as sm
import statsmodels.formula.api as smf
# 我们使用自带的添加截距系数方法add_constant
X_new= sm.add_constant(X)
X_new[:10]

看来这个方法，把截距项的系数放在了前面

model_3 = sm.OLS(y,X_new)
result = model_3.fit()
print(result.params)
# 前面是截距，后面两个是系数
# [-0.08619059  2.32663944  2.52592665]

2.5 三个方法比较

我们构建的公式y=2*x_1 + 3*x_2 + 5 +noise

但由于噪音干扰，实际上最佳的曲线（MSE最小）已经不是原本的曲线了。

首先看系数和截距：

# 三组系数和截距
print('Sklearn前面系数,最后截距',coef_1,intercept_1)
print('Numpy.linalg前面系数,最后截距',theta)
print('statsmodels.OLS','切片调下位置同上',result.params[1:],result.params[0])

Sklearn前面系数,最后截距 [2.75724514 1.90154838] -6.760579051034362
Numpy.linalg前面系数,最后截距 [ 2.32663944  2.52592665 -0.08619059]
statsmodels.OLS 切片调下位置同上 [2.32663944 2.52592665] -0.08619059391555109

可以看到numpy.linalg和statsmodels.OLS方法结果一样，说明都是矩阵求逆再乘两下得出的结果，sklearn的默认纯线性回归(不带正则项的)，里面没有solver可选（岭回归Ridge就有），看来十有八九计算实现并不是用的正规方程，毕竟正规方程在系数多的情况下，很可能求不出来。

那么sklearn会不会结果比他们两个差呢？

residual_df = pd.DataFrame()
residual_df['y_real'] = y

# sklearn学习原数据
y_sklean = modle_1.predict(X)
residual_df['y_sklean'] =y_sklean

# 手算原数据
y_manual = x1*theta[0] + x2*theta[1] +  theta[2]
residual_df['y_manual'] =y_manual
# X_intercept.dot(theta) # 或者用矩阵乘法

# OLS模型计算结果
y_ols = X_new.dot(result.params)
residual_df['y_ols'] =y_ols
residual_df
# 或者用result.fittedvalues

乍一看，结果貌似完全一样，我们仔细看看：

residual_df.iloc[0,1],residual_df.iloc[0,2],residual_df.iloc[0,3]

# (-23.01881766462413, -23.018817664624127, -23.018817664624144)

结果的很多位后的不一样，还是稍微有一点点的差异；

from sklearn.metrics import mean_squared_error,mean_absolute_error
pd.set_option( 'display.precision',12) # 浮点数精度
error_df = pd.DataFrame()
# 计算绝对值误差
error_df.loc['Skleanr','MAE'] = mean_absolute_error(y,y_sklean)
error_df.loc['numpy','MAE'] = mean_absolute_error(y,y_manual)
error_df.loc['OLS','MAE'] = mean_absolute_error(y,y_ols)
# 计算均方误差
error_df.loc['Skleanr','MSE'] = mean_squared_error(y,y_sklean)
error_df.loc['numpy','MSE'] = mean_squared_error(y,y_manual)
error_df.loc['OLS','MSE'] = mean_squared_error(y,y_ols)
error_df

基本是一样的...，讲道理应该不是完全一样，不过差异肯定在小数点后面很多位。

三、model.summary()

终于要进入正题：

其实本文主要研究的是这个：

可以看到能分成三部分，上面是一些基础信息，中间是系数部分，下面是一些检验量。

其中一眼就能看出来的不讲：

(1)R-squared 即R方，越接近1证明拟合得越好

sum of squares of the regression,即预测数据与原始数据均值之差的平方和。

$SSR=\sum\nolimits_{i=1}\nolimits^{n}(\hat{y_i}-\bar{y})^2$

Total sum of squares，即原始数据和均值之差的平方和

$SST=\sum\nolimits_{i=1}\nolimits^{n}(y_i-\bar{y})^2$

The sum of squares due to error,拟合数据和原始数据对应点的误差的平方和

$SSE=\sum\nolimits_{i=1}\nolimits^{n}(y_i-\hat{y_i})^2$

$SST=SSR+SSE$

$R-squared=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SSE}{SST}$

# 验证R方
ssr = ((y_ols - y.mean())**2).sum() # 118359
sst = ((y - y.mean())**2).sum() # 118438
sse = ((y - y_ols)**2).sum() # 78.57
print('R方')
print(ssr/sst,result.rsquared) # 模型属性可查看

R方
0.9993365473073768 0.9993365473073765

（2）Adj. R-squared 调整R方，根据自由度和样本数修正了一下

注意这里面减去模型中的变量，我们此处DF model =1

# 验证调整R方
1-(100-1)*(1-result.rsquared)/(100-1-1),result.rsquared_adj

(0.9993297773819416, 0.9993297773819416)

（3）F-statistic，可参考

统计|如何理解多元线性回归的F检验的作用与目的-CSDN博客

f_stats = (ssr/1)/sse*(100-1-1) # 注意k=1
print(f_stats) 
print(result.fvalue)
# 完美匹配，最后一点误差可忽略
147614.11435208278
147614.11435208275

(4)Prob (F-statistic)

自然是查表得到的，要手算也行

from scipy.stats import f
# 懂点统计学都知道，曲线下面积即两个面积的和=1
p_value = f.cdf(f_stats,1,98) # 累计分布
print(p_value)

p_value_2= f.sf(f_stats,1,98) # 剩下的面积
print(p_value_2)             
print(result.f_pvalue)

0.9999999999999999
1.4929610240226694e-157
1.492961024022682e-157

(5）Log-Likelihood，似然函数的LOG

可参考https://www.le.ac.uk/users/dsgp1/COURSES/MATHSTAT/13mlreg.pdf

其中T为样本数，本文为100；sigma方即方差，是除以样本数100的，最右边公式可以抵消，sigma方即是总体的误差SSE除了一个样本数。

# T为样本数
L = -(100/2)*np.log(2*np.pi)-(100/2)*np.log(sse/100) - 100/2
print(L),print(result.llf)

-129.8399875471549
-129.8399875471549

(6) AIC,BIC

可参考：

模型中AIC和BIC以及loglikelihood的关系-腾讯云开发者社区-腾讯云

aic= -2*result.llf + 2*2
bic = -2*result.llf + 2*np.log(100)
print(aic),print(result.aic)
print(bic),print(result.bic)

263.6799750943098
263.6799750943098
268.890315466286
268.890315466286

至此第一部分的手动核算已经完毕；