Python——线性回归模型的应用

前言：
线性回归模型属于经典的统计学模型，该模型的应用场景是根据已知的变量（即自变量）来预测某个连续的数值变量（即因变量）。例如餐厅根据媒体的营业数据（包括菜谱价格、就餐人数、预订人数、特价菜折扣等）预测就餐规模或营业额；网站根据访问的历史数据（包括新用户的注册量、老用户的活跃度、网站内容的更新频率等）预测用户的支付转化率；医院根据患者的病历数据（如体检指标、药物复用情况、平时的饮食习惯等）预测某种疾病发生的概率。

由于针对具体操作相关文档太多，所以本文内容涉及具体操作较少，主要是讲方法。

本文内所用到的包：

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import scipy.stats as stats
import statsmodels.api as sm
from scipy.stats import chi2_contingency

1.简单线性回归模型

• 也被称为一元线性回归模型，是指模型中只含有一个自变量和一个因变量
• 一般可以通过散点图刻画两个变量之间的关系，并基于散点图绘制简单线性拟合线，进而使变量之间的关系体现地更加直观

以经典的刹车距离为例：

ccpp=pd.read_csv('cars.csv')
sns.set(font=getChineseFont(8).get_name())
sns.lmplot(x='speed',y='dist',data=ccpp,
               legend_out=False,#将图例呈现在图框内
               truncate=True#根据实际的数据范围，对拟合线做截断操作
               )
plt.show()

从散点图来看，自变量speed与因变量dist之间存在明显的正相关，即刹车速度越大，刹车距离越长。图内的阴影部分为拟合线95%的置信区间，每个散点都是尽可能地围绕在拟合线附近。

通过ols函数求得回归模型的参数解’y~s’表示简单线性回归模型

fit=sm.formula.ols('dist~speed',data=ccpp).fit()
print(fit.params)#Intercept：-17.579095,speed：3.932409

因此，关于刹车距离的简单线性模型为：
dist=-17.579095+3.932409speed
也就是说，刹车速度每提升一个单位，将促使刹车距离增加3.93个单位。

2.多元线性回归模型

• 实际应用中，简单线性回归模型并不多见，因为影响因变量的自变量往往不止一个
• 用于构建多元线性回归模型的数据实际上由两部分组成：一个是一维的因变量y；另一个是二维矩阵的自变量x

以利润表Profit为例，研究影响利润的因素
表结构如下：

profit=pd.read_csv(‘Profit.csv’,sep=",")

fit=sm.formula.ols('Profit~RD_Spend+Administration+Marketing_Spend',data=profit).fit()
print(fit.params)#Intercept:50122.192990,RD_Spend:0.805715,Administration:-0.026816,Marketing_Spend:0.027228

不考虑模型的显著性和回归系数的显著性，得到的回归模型可以表示为：
Profit=50122.192990+0.805715RD_Spend-0.026816Administration+0.027228Marketing_Spend

但是，在实际应用中，单纯的利用ols函数将偏回归系数计算出来，并构造一个多元线性回归模型，得到的结果往往不是理想的，这时候需要借助于统计学中的F检验法和t检验法，完成模型的显著性检验和回归系数的显著性检验。

2.1 应用F检验法完成模型的显著性检验

步骤如下：

1. 提出问题的原假设和备择假设；
2. 在原假设的条件下，构造统计量F；
3. 根据样本信息，计算统计量F的值；
4. 对比统计量的值和理论值，如果计算的统计量值超过理论的值，则拒绝原假设，否则需要接受原假设

假设认为模型的所有偏回归系数全为0（即认为没有一个自变量可以构成因变量的线性组合）
通常在实际的应用中将概率值p与0.05做比较，小于0.05则拒绝原假设，否则接受原假设

2.2应用t检验法完成回归系数的显著性检验

模型通过了显著性检验，只能说明模型关于因变量的线性组合是合理的，但不能说明每个自变量对因变量都具有显著意义，所以还要对模型的回归系数做显著性检验。

• 只有当回归系数通过了t检验，才可以认为模型的系数是显著的
• t检验的出发点就是验证每一个自变量是否能够成为影响因变量的重要因素
• t检验的原假设是假定第j变量的回归系数为0，即认为该变量不是因变量的影响因素
• t统计量大于理论的t分布值，则拒绝原假设，否则接受原假设；也可以根据概率值P判断是否需要拒绝原假设

#返回模型概览
print(fit.summary())

由图可知：F-statistic:296.0，Prob (F-statistic):4.53e-30，F统计量值为296.0，对应的概率值P远远小于0.05，说明应该拒绝原假设，认为模型是显著的。

在各自变量的t统计中，Administration和Marketing_Spend变量所对应的概率值p大于0.05，说明不能拒绝原假设，认为该变量是不显著的，无法认定其实影响Profit的重要因素

• 对于F检验来说，如果无法拒绝原假设，则认为模型是无效的，通常解决办法是增加数据量、改变自变量或选择其他模型；
• 对于t检验，如果无法拒绝原假设，则认为对应的自变量与因变量之间不存在线性关系，通常的解决办法是剔除该变量或修正该变量（如选择对于的数学转换函数对其修正处理）

根据返回的fit模型的概览信息，由于Administration和Marketing_Spend变量的t检验结果是不显著的，故可以探索其余因变量Profit之间的散点关系，如果确实没有线性关系，可将其从模型中剔除。

sns.lmplot(x='Administration',y='Profit',data=profit,
               legend_out=False,#将图例呈现在图框内
               fit_reg=False#不显示拟合曲线
               )
sns.lmplot(x='Marketing_Spend',y='Profit',data=profit,
               legend_out=False,#将图例呈现在图框内
               fit_reg=False#不显示拟合曲线
               )
plt.show()

图中自变量Administration、Marketing_Spend与因变量Profit没有呈现明显的线性关系，故可以认为两者不存在互相依赖关系

#将Administration、Marketing_Spend变量从模型中剔除
fit2 = sm.formula.ols('Profit~RD_Spend',data=profit).fit()
print(fit2.params)#Intercept：49032.899141，RD_Spend ：0.854291
print(fit2.summary())#Prob (F-statistic):3.50e-32；P>|t|：0.000

对模型fit重新调整后，得到的新模型fit2仍然通过了显著性检验，而且每个自变量所对应的系数也是通过显著性检验的。

最终得到的模为：

Profit=49032.899141+0.854291RD_Spend

该回归模型中的系数解释为：在其他条件不变的情况下RD_Spend每增加一个单位，将使Profit增加0.854291个单位。

3.基于回归模型识别异常点

• 回归模型计算过程会依赖于自变量的均值，均值的最大弊端是其容易受到异常点（或极端值）的影响
• 建模数据中存在异常点，一定程度上会影响到建模的有效性
• 对于现行回归模型来说，通常利用帽子矩阵、DFFITS准则、学生化残差或Cook距离进行异常点检测
• 使用以上4中方法判别数据集的第i个样本是否为异常点，前提是已构建好一个线性回归模型，然后基于由get_influence方法获得4种统计量的值

继续使用上面的数据。

#异常值检验
outliers=fit2.get_influence()
#高杠杆值点（帽子矩阵）
leverage=outliers.hat_matrix_diag
#diffits值
dffits=outliers.dffits[0]
#学生化残差
resid_stu=outliers.resid_studentized_external
#cook距离
cook=outliers.cooks_distance[0]

#合并以上4种异常值检验的统计量值
concat_result=pd.concat([pd.Series(leverage,name='leverage'),pd.Series(dffits,name='diffits'),
                             pd.Series(resid_stu,name='resid_stu'),pd.Series(cook,name='cook')],axis=1)
#将上面的concat_result结果与profit数据集合并
raw_outliers=pd.concat([profit,concat_result],axis=1)

前5行数据集打印结果如下：

简单起见，这里使用学生化残差，当学生化残差大于2时，即认为对应的数据点为异常值

#计算异常值数量的比例
outliers_ratio=sum(np.where(np.abs(raw_outliers.resid_stu)>2,1,0))/raw_outliers.shape[0]
print(outliers_ratio)#0.04

结果显示，通过学生化残差识别出了异常值，并且异常值比例为4%。由于异常值比例非常小，故可以考虑将其直接从数据集中删除，由此继续建模将会得到更加稳健且合理的模型

#通过筛选的方法，将异常点排除
none_outliers=raw_outliers.loc[np.abs(raw_outliers.resid_stu)<=2,]
#应用无异常值的数据集重新建模
fit3=sm.formula.ols('Profit~RD_Spend',data=profit).fit()
#返回模型的概览信息
print(fit3.params)#Intercept:49032.899141,RD_Spend:0.854291
print(fit3.summary())

排除异常点之后得到的模型fit3，不管是模型的显著性检验还是系数的显著性检验，各自的概率P值均小于0.05，说明他们均通过显著性检验。

模型fit3为：Profit=49032.899141+0.854291RD_Spend

从fit3模型来看RD_Spend（研发成本）每增加一个单位的成本，所带来的Profit（收益）提升要明显比其他的高，所以将更多的成本花费在研发上是个不错的选择。

以原始数据profit为例，根据fit3模型重新预测各成本下的收益预测值：

pred=fit3.predict(profit[['RD_Spend']])
#对于实际值与预测值的比较
df=pd.concat([pd.Series(profit.Profit/100,name='real'),pd.Series(pred/100,name='prediction')],axis=1)
df['误差绝对值']=np.abs((df['real']-df['prediction'])/100)
print(df.head(10))

从结果上看有的预测值比较接近实际值，有的预测测偏离实际值较远，但从总体上来说，预测值与实际值之间的差异并不是特别大。

4.含有离散变量的回归模型

• 虚拟变量也称为哑变量，专门用来解决离散型变量无法量化的问题
• 解决思路为：根据离散变量的值，衍生出多个"0-1"值的新变量 如：0表示不属于当前状态，1表示属于当前状态
• 假设有未婚、已婚、离婚三个哑变量，违背了数据的非多重共线性假设（即哑变量之间存在非常高的相关性），非已婚非离婚为未婚状态，所以在构建哑变量处理离散型自变量时，哑变量的个数应该为n-1，其中n表示离散型自变量的不同值个数
• 对于线性回归模型来说，从所有哑变量中删除某个哑变量时，被删除的哑变量便成为参照变量（因为可以将参照变量用于对比其他变量对因变量的影响）

以经典的泰坦尼克号数据为例：

titanic=pd.read_csv('Titanic_all.csv',sep=",")
print(titanic.dtypes)

12个变量中涉及5个离散型变量和7个数值型变量。根据知己情况可知，船舱等级Pclass应为离散型变量（3等、2等、1等，并非数值）

查看各变量的缺失比例

print(titanic.isnull().sum(axis=0)/titanic.shape[0])

Cabin缺失比例超过77%，Embarked缺失比例仅为0.22%，二者数据删除。

在这里我们需要利用年龄的非缺失数据构造多元线性回归模型，进而预测缺失比例为19.87%的乘客年龄

如需基于其余变量来预测年龄变量Age,至少有5个变量与年龄无关（乘客编号、姓名、票号信息、座位号信息和登船地点），删除。

1.删除无意义的变量

titanic.drop(labels=['PassengerId','Name','Ticket','Cabin','Embarked'],axis=1,inplace=True)

2.哑变量转换

titanic.Pclass=titanic.Pclass.astype(str)#Pclass变量转换为字符型变量
#将Pclass和Sex变量作哑变量处理
dummies=pd.get_dummies(data=titanic[['Pclass','Sex']])
#将titanic数据集与dummies数据集进行合并
titanic_new=pd.concat([titanic,dummies],axis=1)
#根据哑变量原则，需要从新生成的0和1变量中提出两个变量作为参照变量（以性别中的男性为参照变量，以船舱等几种的Pclass_3作为参照变量）
titanic_new.drop(labels=['Pclass','Sex','Pclass_3','Sex_male'],axis=1,inplace=True)

3.将数据拆分为两部分

• 在清洗后的数据集titanic_new中，仅有Age变量存在缺失值
• 为预测该变量的缺失值，需要将数据集按照年龄是否缺失拆分为两个数据子集，分别用于线性回归模型的构造和基于构造好的数据集对其进行预测

missing=titanic_new.loc[titanic.Age.isnull(),]
no_missing=titanic_new.loc[~titanic.Age.isnull(),]

4.构建多元线性回归模型

 #提取出所有的自变量名称
predictors=no_missing.columns[~no_missing.columns.isin(['Age'])]
#构造多元线性回归模型的"类"
lm_Class=sm.OLS(endog=no_missing.Age,#指定模型中的因变量
                            exog=no_missing[predictors]#指定模型中的自变量
                            )
#基于"类"，对模型进行拟合
lm_model=lm_Class.fit()
print(lm_model.summary())

根据F检验的结果可知模型是显著的，但是从t检验的结果来看，仅有船舱等级Pclass和性别Sex是通过显著性检验的。

绘制其他变量与年龄之间的散点图

sns.pairplot(no_missing[['Survived','SibSp','Parch','Fare','Age']])

唯有SibSp与Age之间存在一定趋势性，将出SibSp之外其他变量从模型中剔除。

基于散点图结果重新构造多元线性回归模型

predictors2=['SibSp','Pclass_1','Pclass_2','Sex_female']
lm_Class2=sm.OLS(endog=no_missing.Age,#指定模型中的因变量
                            exog=no_missing[predictors2]#指定模型中的自变量
                            )
lm_model2=lm_Class2.fit()
print(lm_model2.summary())

在新的模型中，F检验的统计量和t检验的统计量所对应的概率P值均小于0.05，说明他们通过了显著性检验。

最后再利用可视化的QQ图，验证模型的残差是否服从正态分布的假设

sns.set(font=getChineseFont(8).get_name())
pp_qq_plot=sm.ProbPlot(lm_model2.resid)
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q图')
plt.show()

图形中的散点基本上都围绕在斜线附近，可以判断模型lm_model2的残差服从正态分布，最终证明了lm_model2模型是合理的。

5.未知年龄的预测

pred_Age=lm_model2.predict(exog=missing[predictors2])
#将年龄的预测值替换missing数据集中的Age变量
missing.loc[:,'Age']=pred_Age
#print(missing.head(5))
#将结果插补到原始数据中
titanic.loc[:,'Age']=pd.concat([no_missing.Age,pred_Age],axis=0)
print(titanic.head(10))