线性回归模型

看书的时候发现书上代码的一些错误，故自己整理一下，但还是有不明白的地方（PP图的错误）。

1. 一元线性回归模型

数学表达式：y = a + bx + ε

  1. 图形绘制
# 工作年限与收入之间的散点图
# 导入第三方模块
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 导入数据集
income = pd.read_csv('Salary_Data.csv')
# 绘制散点图
sns.lmplot(x = 'YearsExperience', y = 'Salary', data = income, ci = None)
# 显示图形
plt.show()
  1. 拟合线求解
# 简单线性回归模型的参数求解
# 样本量
n = income.shape[0]
# 计算自变量、因变量、自变量平方、自变量与因变量乘积的和
sum_x = income.YearsExperience.sum()
sum_y = income.Salary.sum()
sum_x2 = income.YearsExperience.pow(2).sum()
xy = income.YearsExperience * income.Salary
sum_xy = xy.sum()
# 根据公式计算回归模型的参数
b = (sum_xy-sum_x*sum_y/n)/(sum_x2-sum_x**2/n)
a = income.Salary.mean()-b*income.YearsExperience.mean()
# 打印出计算结果
print('回归参数a的值：',a)
print('回归参数b的值：',b)
  1. 使用第三方模块求解
# 导入第三方模块
import statsmodels.api as sm
# 利用收入数据集，构建回归模型
fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit()
# 返回模型的参数值
fit.params

2. 多元线性回归模型

数学表达式： y = Xβ + ε

β： p × 1 的一维向量，代表多元线性回归模型的偏回归系数
ε：n × 1 的一维向量，代表每一个样本的误差项，期望为0，正态分布

  1. 回归模型的参数求解
import statsmodels.api as sm 
from sklearn import model_selection
# 导入数据
Profit = pd.read_excel('Predict to Profit.xlsx')
# 将数据集拆分为训练集和测试集
train, test = model_selection.train_test_split(Profit, test_size = 0.2, random_state=1234)
# 根据train数据集建模（数据state为字符型变量，建模时进行特殊处理。）
model = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + C(State)', data = train).fit()
print('模型的偏回归系数分别为：\n', model.params)
# 删除test数据集中的Profit变量，用剩下的自变量进行预测
test_X = test.drop(labels = 'Profit', axis = 1)
pred = model.predict(exog = test_X)
print('对比预测值和实际值的差异：\n',pd.DataFrame({'Prediction':pred,'Real':test.Profit}))
  1. 自己定义哑变量

字符型变量通过 C() 被转化为了哑变量，既衍生出两个变量，另一个做对照组。但这是随机的，如需将New York作为参照：

# 生成由State变量衍生的哑变量 (1列变3列 0，1数据)
dummies = pd.get_dummies(Profit.State)
# 将哑变量与原始数据集水平合并
Profit_New = pd.concat([Profit,dummies], axis = 1)
# 删除State变量和California变量（因为State变量已被分解为哑变量，New York变量需要作为参照组）
Profit_New.drop(labels = ['State','New York'], axis = 1, inplace = True)

# 拆分数据集Profit_New
train, test = model_selection.train_test_split(Profit_New, test_size = 0.2, random_state=1234)
# 建模
model2 = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + Florida + California', data = train).fit()
print('模型的偏回归系数分别为：\n', model2.params)

获得各变量的参数即可写出回归方程。

3. 回归模型的假设检验

一、 F检验

  1. 提出假设
    原假设 H0：所有偏回归系数都为0
    备择假设 H1：不都为0

就F检验而言，研究者更倾向于通过数据来推翻H0，从而接受备择假设。

  1. 构造统计量
    ESS：误差平方和（真实值与预测值）
    RSS：回归离差平方和（预测值与总体平均值）
    TSS：总离差平方和（真实值与总体平均值）
    TSS = RSS + ESS


    image.png

可以构造F统计量

image.png
  1. 计算统计量

手动计算法

# 导入第三方模块
import numpy as np
# 计算建模数据中，因变量的均值
ybar = train.Profit.mean()
# 统计变量个数和观测个数
p = model2.df_model
n = train.shape[0]
# 计算回归离差平方和
RSS = np.sum((model2.fittedvalues-ybar) ** 2)
# 计算误差平方和
ESS = np.sum(model2.resid ** 2)
# 计算F统计量的值
F = (RSS/p)/(ESS/(n-p-1))
print('F统计量的值：',F)

直接利用 fvalue 获得F统计量值

# 返回模型中的F值
model2.fvalue
  1. 对比F统计量与理论F分布的值
# 导入模块
from scipy.stats import f
# 计算F分布的理论值
F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1)
print('F分布的理论值为：',F_Theroy)

# 模型的概览信息
model2.summary()

二、回归系数的显著性检验 -- t 检验

  1. 提出假设
    H0：第 j 个变量的偏回归系数为0，既该变量不是因变量的影响因素。
    H1：第 j 个变量的偏回归系数不为0

  2. 构造统计量


    image.png
image.png

  1. model.summary() 中查看t，及p>[t]
    每个t统计量值都对应了概率值p，当p<0.05时表示拒绝原假设，即该项对因变量有显著影响。如图所示，只有 Intercept 和RD_Spend有显著影响。


    model2.summary().png

4. 回归模型诊断

统计学家提出，只有模型满足一些前提的情况下，模型才是合理的。

  1. 误差项 ε 服从正态分布
  2. 无多重共线性
  3. 线性相关性
  4. 误差项 ε 的独立性
  5. 方差齐性

除此之外，需要对异常点进行识别和处理

1. 正态性检验

要求误差项服从正态分布，就是要求因变量服从正态分布

图形法1：直方图法

# 正态性检验
# 直方图法
# 导入第三方模块
import scipy.stats as stats
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘制直方图
sns.distplot(a = Profit_New.Profit, bins = 10, fit = stats.norm, norm_hist = True,
             hist_kws = {'color':'steelblue', 'edgecolor':'black'}, 
             kde_kws = {'color':'black', 'linestyle':'--', 'label':'核密度曲线'}, 
             fit_kws = {'color':'red', 'linestyle':':', 'label':'正态密度曲线'})
# 显示图例
plt.legend()
# 显示图形
plt.show()
image.png

直观上，可以认为利润变量服从正态分布。

图形法2：PP图与QQ图

PP图：比对正态分布的累计概率值与实际分布的累计概率值。
QQ图：比对正态分布的分位数与实际分布的分位数
判断标准：散点比较均匀地散落在直线上，则说明变量服从正态分布，否则就认为不服从正态分布。

# 残差的正态性检验（PP图和QQ图法）
pp_qq_plot = sm.ProbPlot(Profit_New.Profit)
# 绘制PP图
pp_qq_plot.ppplot(line = "45")
plt.title('P-P图')
# 绘制QQ图
pp_qq_plot.qqplot(line = 'q')
plt.title('Q-Q图')
# 显示图形
plt.show()

这里的PP图有些问题，但目前还不知道问题出在哪里，如果有大佬看到这篇笔记的话希望可以解释一下。


image.png

QQ图是正常的


image.png

非参数法：shapiro检验和K-S检验

原假设 H0：变量服从正态分布
若数据量低于5000，使用shapiro检验法（stats.shapiro）较合理
否则，推荐使用K-S检验法（stats.kstest）

数据量小于5000，适应shapiro检验法

# 导入模块
import scipy.stats as stats
stats.shapiro(Profit_New.Profit)

输出结果1：统计量值
输出结果2：p值，大于0.05则接受原假设。

K-S检验法 代码示例

stats.kstest(Profit_New.Profit, cdf = "norm")

构造一组数据，应用K-S检验

# 生成正态分布和均匀分布随机数
rnorm = np.random.normal(loc = 5, scale=2, size = 10000)
runif = np.random.uniform(low = 1, high = 100, size = 10000)
# 正态性检验
KS_Test1 = stats.kstest(rvs = rnorm, args = (rnorm.mean(), rnorm.std()), cdf = 'norm')
KS_Test2 = stats.kstest(rvs = runif, args = (runif.mean(), runif.std()), cdf = 'norm')
print(KS_Test1)
print(KS_Test2)
2. 多重共线性检验

多重共线性是指自变量之间存在较高的线性相关关系。
可以使用方差膨胀因子VIF来鉴定，若 VIF > 10 ，则说明变量间存在多重共线性，若 VIF > 100，则说明变量间存在严重的多重共线性。

# 导入statsmodels模块中的函数
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 自变量X(包含RD_Spend、Marketing_Spend和常数列1)
X = sm.add_constant(Profit_New.iloc[:, [0, 2]]) #pandas 取消了.ix，可以使用iloc代替。

# 构造空的数据框，用于存储VIF值
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
3. 线性相关性检验

确保自变量和因变量之间存在线性关系，可以使用Pearson相关系数进行识别。

image.png

Pearson相关系数可以使用“corrwith”方法，可以计算任意变量间的相关系数。

# 计算数据集Profit_New中每个自变量与因变量利润之间的相关系数
Profit_New.drop('Profit', axis = 1).corrwith(Profit_New.Profit)

相关系数:
0.8-1.0 极强相关
0.6-0.8 强相关
0.4-0.6 中等程度相关
0.2-0.4 弱相关
0.0-0.2 极弱相关或无相关

使用 pairplot 画图

# 散点图矩阵
# 导入模块
import matplotlib.pyplot as plt
import seaborn
# 绘制散点图矩阵
seaborn.pairplot(Profit_New.iloc[:,0:4])
# 显示图形
plt.show()
image.png

模型修正

只保留 RD_Spend 和 Marketing_Spend 两个自变量

# 模型修正
model3 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = train).fit()
# 模型回归系数的估计值
model3.params

最终结果：
Profit = 51902.11 + 0.79RD_Spend + 0.02Marketing_Spend

4. 异常值检测
  1. 帽子矩阵

考察第i个样本对预测值（y hat）的影响大小。

  1. DFFITS 准则
  2. Cook 距离
# 异常值检验
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]

# 合并各种异常值检验的统计量值
contat1 = 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数据的行索引
train.index = range(train.shape[0])
# 将上面的统计量与train数据集合并
profit_outliers = pd.concat([train,contat1], axis = 1)
profit_outliers.head()

设定，标准化残差大于2时，即认为对应的数据点为异常值。计算异常值的比例。

# 计算异常值数量的比例
outliers_ratio = profit_outliers[abs(profit_outliers["resid_stu"]) > 2].shape[0]/profit_outliers.shape[0]
outliers_ratio

删除异常值或将其转化为哑变量

# 挑选出非异常的观测点
none_outliers = profit_outliers[abs(profit_outliers["resid_stu"]) <=2]

应用无异常值的数据集重新建模

model4 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = none_outliers).fit()
model4.param

Profit = 51827.42 + 0.80RD_Spend + 0.02Marketing_Spend

5. 独立性检验

因变量 y 或残差项 ε 的独立性检验，通常使用 Durbin-Watson 统计量来测试。
DW值接近2，则说明满足残差的独立性假设。

在 summary 中可以直接查看 Durbin-Watson 的值
model4.summary()
6. 方差齐性检验

模型的残差项不随自变量的变动而呈现某种趋势

  1. 图形法

绘制残差和自变量之间的散点图，就可以发现两者间是否存在某种趋势。

# 方差齐性检验
# 设置第一张子图的位置
ax1 = plt.subplot2grid(shape = (2,1), loc = (0,0))
# 绘制散点图
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', linestyles = '--')
# 添加x轴和y轴标签
ax1.set_xlabel('RD_Spend')
ax1.set_ylabel('Std_Residual')

# 设置第二张子图的位置
ax2 = plt.subplot2grid(shape = (2,1), loc = (1,0))
# 绘制散点图
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 = 'red', linestyles = '--')
# 添加x轴和y轴标签
ax2.set_xlabel('Marketing_Spend')
ax2.set_ylabel('Std_Residual')

# 调整子图之间的水平间距和高度间距
plt.subplots_adjust(hspace=0.6, wspace=0.3)
# 显示图形
plt.show()
image.png

所有点都均匀分布在参考线 y = 0 附近，所以满足方差齐性。

  1. BP 检验

假设残差的方差是一个常数，通过构建拉格朗日乘子LM统计量，实现方差齐性检验。
可以借助 statsmodels 中的 het_breuschpagan 函数实现。

#BP检验
sm.stats.diagnostic.het_breuschpagan(model4.resid, exog_het = model4.model.exog)

5. 回归模型预测

# 模型预测
# model4对测试集的预测
pred4 = model4.predict(exog = test.iloc[:,[0, 2]])
# 绘制预测值与实际值的散点图
plt.scatter(x = test.Profit, y = pred4)
# 添加斜率为1，截距项为0的参考线
plt.plot([test.Profit.min(),test.Profit.max()],[test.Profit.min(),test.Profit.max()],
        color = 'red', linestyle = '--')
# 添加轴标签
plt.xlabel('实际值')
plt.ylabel('预测值')
# 显示图形
plt.show()
image.png
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 216,324评论 6赞 498
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 92,356评论 3赞 392
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 162,328评论 0赞 353
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,147评论 1赞 292
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,160评论 6赞 388
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,115评论 1赞 296
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,025评论 3赞 417
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,867评论 0赞 274
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,307评论 1赞 310
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,528评论 2赞 332
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,688评论 1赞 348
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,409评论 5赞 343
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,001评论 3赞 325
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,657评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,811评论 1赞 268
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,685评论 2赞 368
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,573评论 2赞 353