第七章线性回归预测模型 |
您所在的位置:网站首页 › 百分比数据分析模型 › 第七章线性回归预测模型 |
前言
线性回归模型属于经典的统计学模型,该模型的应用场景是根据已知的变量(自变量)来预测某个连续的数值变量(因变量)。例如,餐厅根据每天的营业数据(包括菜谱价格、就餐人数、预定人数、特价菜折扣等)预测就餐规模或营业额;网站根据访问的历史数据(包括新用户的注册量、老用户的活跃度、网页内容的更新频率等)预测用户的支付转化率;医院根据患者的病历数据(如体检指标、药物服用情况、平时的饮食习惯等)预测某种疾病发生的概率。 站在数据挖掘的角度看待线性回归模型,它属于一种有监督的学习算法,即在建模过程中必须同时具备自变量x和因变量y。本章的重点就是介绍有关线性回归模型背后的数学原理和应用实战,通过本章内容的学习,读者将会掌握如下内容: 一元线性回归模型的实战;多元线性回归模型的系数推导;线性回归模型的假设检验;线性回归模型的诊断;线性回归模型的预测。 7.1 一元线性回归模型一元线性回归模型也被称为简单线性回归模型,是指模型中只含有一个自变量和一个因变量,用来建模的数据集可以表示成{(x1,y1),(x2,y2),…,(xn,yn)}。其中,xi表示自变量x的第i个值,yi表示因变量y的第i个值,n表示数据集的样本量。当模型构建好之后,就可以根据其他自变量x的值,预测因变量y的值,该模型的数学公式可以表示成: y=a+bx+ε如上公式所示,该模型特别像初中所学的一次函数。其中,a为模型的截距项,b为模型的斜率项,ε为模型的误差项。模型中的a和b统称为回归系数,误差项ε的存在主要是为了平衡等号两边的值,通常被称为模型无法解释的部分。 为了使读者理解简单线性回归模型的数学公式,这里不妨以收入数据集为例,探究工作年限与收入之间的关系。在第6章的数据可视化部分已经介绍了有关散点图的绘制,下面将绘制工作年限与收入的散点图,并根据散点图添加一条拟合线: # 工作年限与收入之间的散点图 # 导入第三方模块 import pandas as pd import matplotlib.pyplot as plt import seaborn as sns # 导入数据集 income = pd.read_csv(r'D:\projects\notebook\data\Salary_Data.csv') # 绘制散点图 sns.lmplot(x = 'YearsExperience', y = 'Salary', data = income, ci = None) # 显示图形 plt.show()结果: 拟合线的求解 本节的内容就是关于简单线性回归模型的求解,即如何根据自变量x和因变量y,求解回归系数a和b。前面已经提到,误差项ε是为了平衡等号两边的值,如果拟合线能够精确地捕捉到每一个点(所有的散点全部落在拟合线上),那么对应的误差项ε应该为0。按照这个思路来看,要想得到理想的拟合线,就必须使误差项ε达到最小。由于误差项是y与a+bx的差,结果可能为正值或负值,因此误差项ε达到最小的问题需转换为误差平方和最小的问题(最小二乘法的思路)。误差平方和的公式可以表示为: 第一步:展开平方项
结果: 回归参数a的值: 25792.200198668666 回归参数b的值: 9449.962321455081如上所示,利用Python的计算功能,最终得到模型的回归参数值。你可能会觉得麻烦,为了计算回归模型的参数还得人工写代码,是否有现成的第三方模块可以直接调用呢?答案是肯定的,这个模块就是statsmodels,它是专门用于统计建模的第三方模块,如需实现线性回归模型的参数求解,可以调用子模块中的ols函数。有关该函数的语法及参数含义可见下方: ols(formula, data, subset=None, drop_cols=None) formula:以字符串的形式指定线性回归模型的公式,如’y~x’就表示简单线性回归模型。data:指定建模的数据集。subset:通过bool类型的数组对象,获取data的子集用于建模。drop_cols:指定需要从data中删除的变量。这是一个语法非常简单的函数,而且参数也通俗易懂,但该函数的功能却很强大,不仅可以计算模型的参数,还可以对模型的参数和模型本身做显著性检验、计算模型的决定系数等。接下来,利用该函数计算模型的参数值,进而验证手工方式计算的参数是否正确: # 导入第三方模块 import statsmodels.api as sm # 利用收入数据集,构建回归模型 fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit() # 返回模型的参数值 fit.params结果: Intercept 25792.200199 YearsExperience 9449.962321 dtype: float64如上结果所示,Intercept表示截距项对应的参数值,YearsExperience表示自变量工作年限对应的参数值。对比发现,函数计算出来的参数值与手工计算的结果完全一致,所以,关于收入的简单线性回归模型可以表示成: Salary=25792.20+9449.96YearsExperience 7.2 多元线性回归模型读者会不会觉得一元线性回归模型比较简单呢?它反映的是单个自变量对因变量的影响,然而实际情况中,影响因变量的自变量往往不止一个,从而需要将一元线性回归模型扩展到多元线性回归模型。 如果构建多元线性回归模型的数据集包含n个观测、p+1个变量(其中p个自变量和1个因变量),则这些数据可以写成下方的矩阵形式: 根据线性代数的知识,可以将上式表示成y=Xβ+ε。其中,β为p×1的一维向量,代表了多元线性回归模型的偏回归系数;ε为n×1的一维向量,代表了模型拟合后每一个样本的误差项。 7.2.1 回归模型的参数求解 在多元线性回归模型所涉及的数据中,因变量y是一维向量,而自变量X为二维矩阵,所以对于参数的求解不像一元线性回归模型那样简单,但求解的思路是完全一致的。为了使读者掌握多元线性回归模型参数的求解过程,这里把详细的推导步骤罗列到下方: 第一步:构建目标函数 第二步:展开平方项 J(β)=(y-Xβ)'(y-Xβ) =(y'-β'X')(y-Xβ) =(y'y-y'Xβ-β'X'y+β'X'Xβ)由于上式中的y’Xβ和β’X’y均为常数,并且常数的转置就是其本身,因此y’Xβ和β’X’y是相等的。接下来,对目标函数求参数β的偏导数。 第三步:求偏导 第四步:计算偏回归系数的值 X'Xβ=X'y β=(X'X)-1X'y经过如上四步的推导,最终可以得到偏回归系数β与自变量X、因变量y的数学关系。这个求解过程也被成为“最小二乘法”。基于已知的偏回归系数β就可以构造多元线性回归模型。前文也提到,构建模型的最终目的是为了预测,即根据其他已知的自变量X的值预测未知的因变量y的值。 7.2.2 回归模型的预测 如果已经得知某个多元线性回归模型y=β0+β1x1+β2x2+…+βpxn,当有其他新的自变量值时,就可以将这些值带入如上的公式中,最终得到未知的y值。在Python中,实现线性回归模型的预测可以使用predict“方法”,关于该“方法”的参数含义如下: predict(exog=None, transform=True) exog:指定用于预测的其他自变量的值。transform:bool类型参数,预测时是否将原始数据按照模型表达式进行转换,默认为True。接下来将基于statsmodels模块对多元线性回归模型的参数进行求解,进而依据其他新的自变量值实现模型的预测功能。这里不妨以某产品的利润数据集为例,该数据集包含5个变量,分别是产品的研发成本、管理成本、市场营销成本、销售市场和销售利润,数据集的部分截图如表7-1所示。 结果: 模型的偏回归系数分别为: Intercept 58581.516503 C(State)[T.Florida] 927.394424 C(State)[T.New York] -513.468310 RD_Spend 0.803487 Administration -0.057792 Marketing_Spend 0.013779 dtype: float64 对比预测值和实际值的差异: Prediction Real 8 150621.345801 152211.77 48 55513.218079 35673.41 14 150369.022458 132602.65 42 74057.015562 71498.49 29 103413.378282 101004.64 44 67844.850378 65200.33 4 173454.059691 166187.94 31 99580.888894 97483.56 13 128147.138396 134307.35 18 130693.433835 124266.90如上结果所示,得到多元线性回归模型的回归系数及测试集上的预测值,为了比较,将预测值和测试集中的真实Profit值罗列在一起。针对如上代码需要说明三点: 为了建模和预测,将数据集拆分为两部分,分别是训练集(占80%)和测试集(占20%),训练集用于建模,测试集用于模型的预测。由于数据集中的State变量为非数值的离散变量,故建模时必须将其设置为哑变量的效果,实现方式很简单,将该变量套在C()中,表示将其当作分类(Category)变量处理。对于predict“方法”来说,输入的自变量X与建模时的自变量X必须保持结构一致,即变量名和变量类型必须都相同,这就是为什么代码中需要将test数据集的Profit变量删除的原因。对于输出的回归系数结果,读者可能会感到疑惑,为什么字符型变量State对应两个回归系数,而且标注了Florida和New York。那是因为字符型变量State含有三种不同的值,分别是California、Florida和New York,在建模时将该变量当作哑变量处理,所以三种不同的值就会衍生出两个变量,分别是State[Florida]和State[New York],而另一个变量State[California]就成了对照组。 正如建模中的代码所示,将State变量套在C()中,就表示State变量需要进行哑变量处理。但是这样做会存在一个缺陷,那就是无法指定变量中的某个值作为对照组,正如模型结果中默认将State变量的California值作为对照组(因为该值在三个值中的字母顺序是第一个)。如需解决这个缺陷,就要通过pandas模块中的get_dummies函数生成哑变量,然后将所需的对照组对应的哑变量删除即可。为了使读者明白该解决方案,这里不妨重新建模,并以State变量中的New York值作为对照组,代码如下: # 生成由State变量衍生的哑变量 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)结果: 模型的偏回归系数分别为: Intercept 58068.048193 RD_Spend 0.803487 Administration -0.057792 Marketing_Spend 0.013779 Florida 1440.862734 California 513.468310 dtype: float64如上结果所示,从离散变量State中衍生出来的哑变量在回归系数的结果里只保留了Florida和California,而New York变量则作为了参照组。以该模型结果为例,得到的模型公式可以表达为: Profit=58068.05+0.80RD_Spend-0.06Administation+0.01Marketing_Spend+1440.86Florida+513.47California虽然模型的回归系数求解出来了,但从统计学的角度该如何解释模型中的每个回归系数呢?下面分别以研发成本RD_Spend变量和哑变量Florida为例,解释这两个变量对模型的作用:在其他变量不变的情况下,研发成本每增加1美元,利润会增加0.80美元;在其他变量不变的情况下,以New York为基准线,如果在Florida销售产品,利润会增加1440.86美元。 关于产品利润的多元线性回归模型已经构建完成,但是该模型的好与坏并没有相应的结论,还需要进行模型的显著性检验和回归系数的显著性检验。在下一节,将重点介绍有关线性回归模型中的几点重要假设检验。 7.3 回归模型的假设检验模型的显著性检验是指构成因变量的线性组合是否有效,即整个模型中是否至少存在一个自变量能够真正影响到因变量的波动。该检验是用来衡量模型的整体效应。回归系数的显著性检验是为了说明单个自变量在模型中是否有效,即自变量对因变量是否具有重要意义。这种检验则是出于对单个变量的肯定与否。 模型的显著性检验和回归系数的显著性检验分别使用统计学中的F检验法和t检验法,接下来将介绍有关F检验和t检验的理论知识和实践操作。 7.3.1 模型的显著性检验——F检验 在统计学中,有关假设检验的问题,都有一套成熟的步骤。首先来看一下如何应用F检验法完成模型的显著性检验,具体的检验步骤如下: (1)提出问题的原假设和备择假设。 (2)在原假设的条件下,构造统计量F。 (3)根据样本信息,计算统计量的值。 (4)对比统计量的值和理论F分布的值,如果计算的统计量值超过理论的值,则拒绝原假设,否则需接受原假设。 下面将按照上述四个步骤对构造的多元线性回归模型进行F检验,进一步确定该模型是否可用,详细操作步骤如下: 步骤一:提出假设 步骤二:构造统计量 为了使读者理解F统计量的构造过程,可以先观看图7-2,然后掌握总的离差平方和、回归离差平方和与误差平方和的概念与差异。 结果: F统计量的值: 174.63721716844674为了验证手工计算的结果是否正确,可以通过fvalue“方法”直接获得模型的F统计量值,如下结果所示,经过对比发现,手工计算的结果与模型自带的F统计量值完全一致: model2.fvalue结果: 174.63721715703542步骤四:对比结果下结论 最后一步所要做的是对比F统计量的值与理论F分布的值,如果读者手中有F分布表,可以根据置信水平(0.05)和自由度(5,34)查看对应的分布值。为了简单起见,这里直接调用Python函数计算理论分布值: # 导入模块 from scipy.stats import f # 计算F分布的理论值 F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1) print('F分布的理论值为:',F_Theroy)结果: F分布的理论值为: 2.502635007415366如上结果所示,在原假设的前提下,计算出来的F统计量值174.64远远大于F分布的理论值2.50,所以应当拒绝原假设,即认为多元线性回归模型是显著的,也就是说回归模型的偏回归系数都不全为0。 7.3.2 回归系数的显著性检验——t检验 模型通过了显著性检验,只能说明关于因变量的线性组合是合理的,但并不能说明每个自变量对因变量都具有显著意义,所以还需要对模型的回归系数做显著性检验。关于系数的显著性检验,需要使用t检验法,构造t统计量。接下来按照模型显著性检验的四个步骤,对偏回归系数进行显著性检验。 步骤一:提出假设 步骤二:构造统计量
步骤三:计算统计量 如果读者对t统计量值的计算比较感兴趣,可以使用如上公式完成统计量的计算,这里就不手工计算了。为了方便起见,可以直接调用summary“方法”,输出线性回归模型的各项指标值: # 模型的概览信息 model2.summary()结果: 步骤四:对比结果下结论 在第二部分的内容中,含有每个偏回归系数的t统计量值,它的计算就是由估计值coef和标准误std err的商所得的。同时,每个t统计量值都对应了概率值p,用来判别统计量是否显著的直接办法,通常概率值p小于0.05时表示拒绝原假设。从返回的结果可知,只有截距项Intercept和研发成本RD_Spend对应的p值小于0.05,才说明其余变量都没有通过系数的显著性检验,即在模型中这些变量不是影响利润的重要因素。 7.4 回归模型的诊断当回归模型构建好之后,并不意味着建模过程的结束,还需要进一步对模型进行诊断,目的就是使诊断后的模型更加健壮。统计学家在发明线性回归模型的时候就提出了一些假设前提,只有在满足这些假设前提的情况下,所得的模型才是合理的。本节的主要内容就是针对如下几点假设,完成模型的诊断工作: 误差项ε服从正态分布。无多重共线性。线性相关性。误差项ε的独立性。方差齐性。除了上面提到的五点假设之外,还需要注意的是,线性回归模型对异常值是非常敏感的,即模型的构建过程非常容易受到异常值的影响,所以诊断过程中还需要对原始数据的观测进行异常点识别和处理。接下来,结合理论知识与Python代码逐一展开模型的诊断过程。 7.4.1 正态性检验 虽然模型的前提假设是对残差项要求服从正态分布,但是其实质就是要求因变量服从正态分布。对于多元线性回归模型y=Xβ+ε来说,等式右边的自变量属于已知变量,而等式左边的因变量为未知变量(故需要通过建模进行预测)。所以,要求误差项服从正态分布,就是要求因变量服从正态分布,关于正态性检验通常运用两类方法,分别是定性的图形法(直方图、PP图或QQ图)和定量的非参数法(Shapiro检验和K-S检验),接下来通过具体的代码对原数据集中的利润变量进行正态性检验。 1.直方图法 # 正态性检验 # 直方图法 # 导入第三方模块 import scipy.stats as stats # 中文和负号的正常显示 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 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()结果: 2.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()结果: 3.Shapiro检验和K-S检验 这两种检验方法均属于非参数方法,它们的原假设被设定为变量服从正态分布,两者的最大区别在于适用的数据量不一样,若数据量低于5000,则使用shapiro检验法比较合理,否则使用K-S检验法。scipy的子模块stats提供了专门的检验函数,分别是shapiro函数和kstest函数,由于利润数据集的样本量小于5000,故下面运用shapiro函数对利润做定量的正态性检验: # 导入模块 import scipy.stats as stats stats.shapiro(Profit_New.Profit)结果: (0.9793398380279541, 0.537902295589447)如上结果所示,元组中的第一个元素是shapiro检验的统计量值,第二个元素是对应的概率值p。由于p值大于置信水平0.05,故接受利润变量服从正态分布的原假设。 为了应用K-S检验的函数kstest,这里随机生成正态分布变量x1和均匀分布变量x2,具体操作代码如下: # 生成正态分布和均匀分布随机数 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)结果: KstestResult(statistic=0.0054919289086649, pvalue=0.9236250840080669) KstestResult(statistic=0.05906629005008068, pvalue=9.941845404073026e-31)如上结果所示,正态分布随机数的检验p值大于置信水平0.05,则需接受原假设;均匀分布随机数的检验p值远远小于0.05,则需拒绝原假设。需要说明的是,如果使用kstest函数对变量进行正态性检验,必须指定args参数,它用于传递被检验变量的均值和标准差。 如果因变量的检验结果不满足正态分布时,需要对因变量做某种数学转换,使用比较多的转换方法有: 关于多重共线性的检验可以使用方差膨胀因子VIF来鉴定,如果VIF大于10,则说明变量间存在多重共线性;如果VIF大于100,则表名变量间存在严重的多重共线性。方差膨胀因子VIF的计算步骤如下: (1)构造每一个自变量与其余自变量的线性回归模型,例如,数据集中含有p个自变量,则第一个自变量与其余自变量的线性组合可以表示为: x1=c0+α2x1+…+αpxp+ε(2)根据如上线性回归模型得到相应的判决系数R2,进而计算第一个自变量的方差膨胀因子VIF: Python中的statsmodels模块提供了计算方差膨胀因子VIF的函数,下面利用该函数计算两个自变量的方差膨胀因子: # 导入statsmodels模块中的函数 from statsmodels.stats.outliers_influence import variance_inflation_factor # 自变量X(包含RD_Spend、Marketing_Spend和常数列1) X = sm.add_constant(Profit_New.ix[:,['RD_Spend','Marketing_Spend']]) # 构造空的数据框,用于存储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结果: 7.4.3 线性相关性检验 线性相关性检验,顾名思义,就是确保用于建模的自变量和因变量之间存在线性关系。关于线性关系的判断,可以使用Pearson相关系数和可视化方法进行识别,有关Pearson相关系数的计算公式如下: 结果: RD_Spend 0.978437 Administration 0.205841 Marketing_Spend 0.739307 California -0.083258 Florida 0.088008 dtype: float64如上结果所示,自变量中只有研发成本和市场营销成本与利润之间存在较高的相关系数,相关性分别达到0.978和0.739,而其他变量与利润之间几乎没有线性相关性可言。通常情况下,可以参考表7-3判断相关系数对应的相关程度: 结果: 结果: Intercept 51902.112471 RD_Spend 0.785116 Marketing_Spend 0.019402 dtype: float64如上结果所示,返回的是模型两个自变量的系数估计值,可以将多元线性回归模型表示成: Profit=51902.11+0.79RD_Spend+0.02Marketing_Spend7.4.4 异常值检验 由于多元线性回归模型容易受到极端值的影响,故需要利用统计方法对观测样本进行异常点检测。如果在建模过程中发现异常数据,需要对数据集进行整改,如删除异常值或衍生出是否为异常值的哑变量。对于线性回归模型来说,通常利用帽子矩阵、DFFITS准则、学生化残差或Cook距离进行异常点检测。接下来,分别对这四种检测方法做简单介绍。 1.帽子矩阵 帽子矩阵的设计思路就是考察第i个样本对预测值: 2.DFFITS准则 DFFITS准则借助于帽子矩阵,构造了另一个判断异常点的统计量,该统计量可以表示成如下公式: 结果: 结果: 0.02564102564102564如上结果所示,通过标准化残差监控到了异常值,并且异常比例为2.5%。对于异常值的处理办法,可以使用两种策略,如果异常样本的比例不高(如小于等于5%),可以考虑将异常点删除;如果异常样本的比例比较高,选择删除会丢失一些重要信息,所以需要衍生哑变量,即对于异常点,设置哑变量的值为1,否则为0。如上可知,建模数据的异常比例只有2.5%,故考虑将其删除。 # 挑选出非异常的观测点 none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu) |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |