python统计分析——线性回归模型的结果_python线性回归模型报告

参考资料：python统计分析【托马斯】

        首先，我们从一个小数据集查看英国不同地区的烟草和酒精购买之间的相关性。该数据集有趣的特点是北爱尔兰报告为异常值。尽管如此，我们将使用次数据集描述计算线性回归的两个工具。我们将使用statsmodels和sklearn模块来计算线性回归，同时使用pandas进行数据管理并使用matplotlib进行绘图。代码如下：

# 导入库
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression
from scipy import stats
 
# 录入数据
region=['North','Yorkshire','Northeast','East_Midlands','West_Midlands',
        'East_Anglia','Southeast','Southwest','Wales','Scotland','Northern_Ireland']
Alcohol=[6.47,6.13,6.19,4.89,5.63,4.52,5.89,4.79,5.27,6.08,4.02]
Tobacco=[4.03,3.76,3.77,3.34,3.47,2.92,3.20,2.71,3.53,4.51,4.56]
 
# 生成dataFrame数据集
df=pd.DataFrame({
    'Region':region,
    'Alcohol':Alcohol,
    'Tobacco':Tobacco
})
 
# 数据可视化展示
df.plot('Tobacco','Alcohol',style='o')
plt.xlabel('Tobacco')
plt.ylabel('Alcohol')
plt.title('Sales in Several UK Regions')

# 暂时忽略掉异常值（上图的最后一个数据点），进行模型的拟合
result=smf.ols("Alcohol~Tobacco",df[:-1]).fit()
print(result.summary())

1、Df Model：指模型的自由度，即预测变量或解释变量的个数。

2、Df Residuals：指残差的自由度，是观测数减去模型自由度，再减去1（考虑到偏移量）。

3、R-squared：指决定系数，可以通过result.rsquared得到。表示y变量的方差中由于x变量的方差所引起的比例。对于简单线性回归来说，R-squared是样本相关系数的平方。对于带截距项的多从线性回归（也包括简单线性回归）来说：

        Adj. R-squared：调整后决定系数。为了评估模型的质量，许多研究人员更偏好使用调整后的决定系数，该参数会惩罚模型中过多的参数。表示为：

其中，k为估计的参数个数（含截距）。

4、F-statistic：F统计量。

        原假设：

        备择假设：，至少有一个j成立。

        F统计量：，服从自由度为的F分布。

        我们可以在python中直接检验

N=result.nobs # nobs:No. Observations
k=result.df_model+1
dfm,dfe=k-1,N-k
F=result.mse_model/result.mse_resid
p=1-stats.f.cdf(F,dfm,dfe)
print('F-statistic:{:.3f},p-value:{:.5f}'.format(F,p))

5、Log-Likelihood：对数似然函数

        统计学中一个非常常见的方法就是：最大似然估计的思想。基本思想与OLS（最小二乘）方法截然不同：在最小二乘法中，模型是常数，响应误差是可变的；相比之下，在最大似然法中，数据响应值被认为是常数，并且将模型的可能性最大化。python中具体计算如下：

N=result.nobs
SSR=result.ssr
s2=SSR/N
L=(1/np.sqrt(2*np.pi*s2))**N*np.exp(-SSR/(s2*2))
print('ln(L)=',np.log(L))

6、AIC和BIC

        为了判断模型的质量，应首先进行数据的可视化检验。此外，还可以用一些数据标准来评估统计模型的质量。这些标准代表了平衡模型准确性与简洁性的各种方法。

        我们已经见过校正后决定系数，和决定系数相比，它在模型有许多回归变量的时候会降低。

        其他常见的标准是Akaike信息准则（AIC）和Schwartz/贝叶斯信息准则（BIC）。两种方法引入了模型复杂度的惩罚，不过AIC对复杂性的惩罚没有BIC那么厉害。AIC可以通过下式获得：

其中，k为模型参数的个数，为对数似然。

        Schwartz/贝叶斯信息准则（BIC）的计算如下：

其中N是观测次数。相当于result.nobs。

        AIC是不同模型之间信息损失的相对估计。BIC最初是使用贝叶斯论证剔除的，并不涉及信息的想法。这两种措施仅在尝试在不同模型之间做决定时才使用。我们应该选择具有较低AIC或BIC值得模型。

7、系数

        线性回归的系数/权重包含在result.params中，并作为pandas的Series对象返回，因为我们使用的是pandas数据框作为输入。效果如下：

        当然我们也可以通过计算值得获得系数：

计算过程如下：

# 在X矩阵中增加一个常数列，用于获得截距
df['Ones']=np.ones(len(df))
 
# 取值-1的索引，是因为最后一个异常点被排除在外
Y=df.Alcohol[:-1]
X=df[['Tobacco','Ones']][:-1]
 
coef_=np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),Y)
print(coef_)

8、标准误

        为了获得系数的标准误，我们将计算协方差-方差矩阵，也叫作协方差矩阵，通过下式来估计预测变量的系数β：

这里，是方差或残差的均方误差。标准误是这个协方差矩阵对角线上元素的平方根。

X=df.Tobacco[:-1]
 
# 加上一列1作为常数截距项
X=np.vstack((np.ones(X.size),X))
 
# 将numpy数组转化为矩阵
X=np.matrix(X)
 
# 进行矩阵乘法再求逆
C=np.linalg.inv(X*X.T)
 
# 乘以残差的均方误差
C=result.mse_resid*C
 
# 计算平方根
SE=np.sqrt(C)
 
print(SE)

此结果与result.summary()输出结果的表2中的标准误相一致。

9、t统计量

        我们使用t检验来检查给定预测变量的系数为0的无效假设，也就是说给定的预测变量对响应变量有没有影响。

        假定模型的残差是大致正态分布在0的时候，t检验能够大体上让我们评估不同预测变量的重要性。如果残差并非正态分布在0附近，这说明变量之间存在一些非线性关系，那么每个预测变量的重要程度就不能用t检验来评估了。更进一步说，也许最好试着修正我们的模型，使得残差在0周围正态性聚集。

        t统计量由感兴趣的预测变量的系数（或因子）与其响应的标准误之比给出。如果β是我们预测变量或因素的系数向量，SE是我们的标准误差，那么t统计量是：

接续上面代码和结果：

i=1
beta=result.params[i]
se=SE[i,i]
t=beta/se
print("t=",t)

        一旦我们得到了t统计量，根据我们使用这些代码得到的误差的正态性假设，我们可以计算观察到一个统计量至少和我们已经观察到的同样极端的概率。

N=result.nobs
k=result.df_model+1
dof=N-k
p_onesided=1-stats.t(dof).cdf(t)
p=p_onesided*2
print('p={0:3f}'.format(p))

10、置信区间

        执行区间使用标准误，t检验的p值和具有N-k自由度的t检验的临界值构建。较小的置信区间表面，我们对估计的系数或常数项的值有信心；较大的置信区间表明估计的项存在更多的不确定性或方差。

        置信区间可以通过下式得到：

其中，β_i是估计的系数之一，z是一个t统计量需要获得小于α置信水平的概率的临界值，SE时标准误。接续上面的代码如下：

i=0
 
# 估计的系数和它的方差
beta,c=result.params[i],SE[i,i]
 
# t统计量临界值
N=result.nobs
P=result.df_model
dof=N-P-1
z=stats.t(dof).ppf(0.975)
 
# 置信区间
print(beta-z*c,beta+z*c)

11、偏度和峰度

        偏度和峰度是指一种分布的形状。偏度是一种分布不对称的度量，峰度是其曲率的度量。python代码实现如下：

S=stats.skew(result.resid,bias=True)
 
K=stats.kurtosis(result.resid,fisher=False,bias=True)
 
print("Skewness:{:.3f}, Kurtosis:{:.3f}".format(S,K))

12、Omnibus检验

        Omnibus检验使用偏度和峰度检验一个分布是否为正态分布。如果我们得到一个非常小的P(Omnibus)值，那么说明分布不符合正态分布。使用statsmodels.stats.normaltest()函数。的代码如下：

(k2,p)=stats.normaltest(result.resid)
print("Omnibus: {0:.3f},p={1:.3f}".format(k2,p))

13、Durbin-Watson检验

        Durbin-Watson检验用来检测残差中的自相关（彼此检测一定时间的值之间的关系），在这里间隔是1。

DW=np.sum(np.diff(result.resid.values)**2)/result.ssr
print('Durbin-Watson:{:.5f}'.format(DW))

14、Jarque-Bera检验

        Jarque-Bera检验是另一个偏度（S）和峰度（K）的检验。无效假设是：分布是正态的，即偏度和超值峰度都是0，或着说偏度是0且峰度为3。但是小样本情况下，即使分布为正态分布Jarque-Bera也倾向于拒绝无效假设。

        用自由度为2的卡方分布计算JB统计量：

JB=(N/6)*(S**2+0.25*(K-3)**2)
p=1-stats.chi2(2).cdf(JB)
print("JB-statistic:{:.5f},p-value:{:.5f}".format(JB,p))

15、条件数：

        条件数测量时一个函数的输出对其输入的灵敏度。当两个预测变量高度相关时（即多重共线性），这些预测变量的系数或因子可能会在数据获或模型的小变化中不稳定地波动。理想情况下，类似的模型应该是类似的，即具有近似相等的系数。多重共线性可能导致数值矩阵求逆失败或产生不准确的结果。解决这个问题的方法是岭回归技术。

        我们首先得到预测变量（包括常数1向量）的乘积的特征值，然后通过最大特征值和最小特征值得比值的平方根来计算条件数。如果条件数大于30，那么回归可能存在多重共线性。

X=np.matrix(X)
EV=np.linalg.eig(X*X.T)
print(EV)
 
CN=np.sqrt(EV[0].max()/EV[0].min())
print('Condition No.: {:.5f}'.format(CN))

16、关于异常值：在实践中，应该在丢弃异常值之前弄明白它们，因为它们可能变得非常重要。它们可能意味着新的趋势，或者一些可能的灾难性事件。