某化学药品公司的研究人员研究了药片表面积和体积对药品在释放剂量控制中释放速度的作用(Drug Development and Industrial pharmacy, Vol. 28, 2002),准备了6个形状相似具有不同重量和厚度的药片,测量了每个药片的表面积与体积的比,利用溶解设备,每个药片放在900ml去离子水中,确定渗滤药品释放速度(药品释放啊的百分比除以时间平方根),试验数据在下表中列出:
从以上数据中,您可以得出什么有用的结论呢?当表面积/体积=1.2时,药片释放速度大概是多少?
一、图形化探索——散点图
上述案例中,两个变量都是数值型变量,同时,我们希望寻找两个变量之间的相互关系,因此,通过散点图来观察两个变量的关系是再恰当不过的选择了。
通过可视化工具我们得到下图,从图中我们可以明显的看出,药品释放速度与药片表面积/体积呈现明显的线性正相关关系,也就是说,随着表面积/体积的增加,药片释放速度呈增加趋势。
可视化在数据分析中的重要性不言而喻,尤其是在演示中扮演着举足轻重的角色,但是图形的缺点也非常明显:
(1)无法量化变量之间的相关关系;
(2)无法进行预测,而回归的主要目的就是为了进行预测。
所以,以上的分析只是线性回归分析最初始的步骤,它可以帮助我们理解数据之间的关系。接下来,我们将主要讨论如何量化变量之间的相关关系,如何通过最小二乘法获取回归方程的最优参数,如何检验模型的显著性,以及非常重要的残差诊断,和如何进行预测,包括预测的陷阱等。
二、它们相关吗?相关关系未必是因果关系
对于两个变量相关关系的量化,著名统计学家卡尔·皮尔逊曾给出了以下关于相关系数的数学定义,根据他的定义:
(1)当r的值越趋近于1或-1时,两个变量的线性关系越强;
(2)当r的值越趋近于0,两个变量的线性关系越弱;
(3)而当r值正好为1或-1时,则表明所有的点都精确地落在一条直线上。
然而,通过相关系数的数学公式来理解相关关系非常地不直观,所以,这里通过引用向量概念来更直观地解释两个变量的相关关系,如下图,以向量空间为例,向量
和向量
之间的夹角的余弦值可以通过以下公式获得(同样,在
空间,以下公式也适用):
由此可见,当两个变量的相关系数越趋近于1时,也就意味着向量之间的夹角越接近于0°,而当其越趋近于-1时,则意味着向量之间的夹角越接近于180°,而当变量之间的相关系数越趋近于0时,则意味着向量之间的夹角越接近于90°,也就是说两个向量之间线性无关。
依据以上的相关系数的公式,我们可以求得本案例中药片释放速度和药片表面积/体积的相关系数为:
由此可见,两个变量线性相关性很强。
但是,这里面有一点需要提醒大家,两个变量之间虽然存在很强的相关关系,但是,这并不意味着一个变量会影响另一个变量,更不能说明两者之间存在因果关系,这仅仅揭示了两个变量之间存在某种数学关系,而至于两个变量之间是否存在因果关系,则需要研究者根据自身的经验来进行判定。
三、为何最小二乘法?——高斯的谷神星轨道计算
在通过散点图和相关系数明确了两个变量之间存在很强的线性关系后,我们还需要进行线性回归,才可以更加明确地揭露变量之间的对应关系,这样做的目的是为了后面进行预测服务,为此我们假设变量和
之间存在以下线性关系:
接下来,需要通过已收集到的数据来估测出最佳的和
,以使得误差平方和取最小值:
以上就是数学史上闻名遐迩的最小二乘法,谈起最小二乘法,需要追溯到1801年,当时,意大利天文学家G. Piazzi发现了谷神星,他在6个星期中跟踪这颗小行星,但由于太阳的干扰,这颗小行星突然不见了,很多著名的天文学家都发表了文章,预测谷神星的轨道,高斯也发表了一个预测,但是他预测的轨道和其他人有相当大的差异,此后,谷神星在1801年12月7日被一个观测者再度发现,并在1802年1月1日又被另一个观测者发现,这两个情况均和高斯预测的位置十分接近,高斯立刻在天文学界赢得了威望,并在一段时间内,被公认为是知名的天文学家而不是数学家,它成功的关键就在于使用了最小二乘法。自此以后,最小二乘法在各个行业领域都得到了广泛的应用。
对和
的估测,学术界存在多种方法,但是都殊途同归,最多见的是通过将残差平方和分别对
和
求偏导数,并使之等于0,通过解方程以求得
和
,此外,也有通过求解超定方程组的方法获得
和
的估测,并且可以给出非常直观的几何解释,此外,也有通过最大似然法来推断
和
的做法。这里我们仅仅给出通过求偏导数来估测
和
的方法:
通过求解上述方程组可以得到:
依据上述公式,我们求得本案例的回归方程为:
四、回归方程显著性检验——决定系数与相关系数
到这里,我们已经知道了本案例两个变量之间具有很强的线性关系,同时,我们也通过最小二乘法估测了回归方程的参数,但是,我们还需要对回归方程的显著性进行检验,而回归方程显著性的检验实际上是一个方差分析的过程,其基本原理是,观测值的变异性(称为总离差平方和,记为
)可以分解为两个部分:
(1)由于和
之间存在线性关系,
的变化会引起
的变化,这种变化可以用回归线上拟合点的波动加以解释,称为回归平方和,记为
;
(2)的观测值与回归线上拟合点的差异是由于随机误差导致的,称为残差平方和,记为
。
因此,如果远大于
,则说明回归方程有效,否则,就认为回归方程无效。
根据上述原理,我们从以下恒等式出发,给出以下推导过程:
最后给出回归显著性检验的方差分析如下表:
对于本案例,根据查F分布的分位数表可得,当α=0.01、分子自由度
、分母自由度时
,
,由于
,因此,可以认为该回归方程有效。
此外,为了度量回归效果的好坏,还有一个重要的统计指标:决定系数,它的含义是回归方程解释观测数据变异的能力,其计算公式如下:
对于本案例,我们可以得出:
但是,从的定义可以看出,当多一个自变量加入模型中,不管这个变量影响是否显著,回归平方和都会增加,因此
也会增大,因此,从
看不出新增加的自变量是否有意义,所以统计学家提出了
,它把残差平方和
和总离差平方和
的自由度考虑在公式里,这样就使得系数的评估更加公正客观。
需要注意的是,大家可能已经发现,本案例中决定系数和相关系数正好相等(),但并不能就此将决定系数和相关系数当作同一个概念,这个相等只发生在一元线性回归中,因此,在解读时需要加以注意,决定系数是度量回归方程解释观测数据变异的能力,而相关系数是度量两个变量之间相关关系的统计量。
五、残差诊断——安斯库姆的忠告
1973年,统计学家F.J.Anscombe曾构造了以下四组奇特的数据集(如表4),我们可以通过上述的线性回归的方法对这个数据集进行回归,我们惊奇地发现了所有四组数据的基本统计量和回归参数等都惊奇地相似(如表5),那么是不是这四组变量中的和
之间都具有线性相关关系呢?答案是否定的,实际上,F.J.Anscombe构造这组数据集的目的就是提醒人们,单纯地通过统计量来进行判断是不可靠的,确认变量之间的关系需要借助可视化工具来辅助判断,否则很可能会得出荒唐的结论。
因此,为了更好地展示这四组数据的关系,我们给出了它们的散点图(如图2),从图形中可以看出,除了第一组变量之间存在线性关系外,其余几组变量的线性关系都不明显,尤其是第三组和第四组变量都存在较为明显的异常值现象,需要特别注意。
F.J.Anscombe数据集给我们的启示是,判断两个或者多个变量之间的关系需要借助可视化来进行确认,而残差诊断则是线性回归中非常必要的一个可视化判断的步骤,如果残差分布具有某种特殊的形状,不够随机的话,可能模型就不是很成功,需要对回归方程做一些调整,比如残差对于拟合值有类似图4-(3)的喇叭口形状,这时建议对取平方根、取倒数或者取对数来尝试解决。如果残差对于自变量的散点图出现类似图4-(2)的弯曲情况,这时建议在模型中添加
的平方项等等。
通常而言,残差诊断包括以下四个部分的诊断:
(1)残差对于观测顺序的散点图;
(2)残差对于拟合值的散点图;
(3)残差的正态概率图;
(4)残差对于各自变量的散点图。
由此,可以看出,本案例的残差无任何异常,并且残差分布也是正态的,模型可以用于进行下一步的预测。
六、预测的陷阱——没有调查没有发言权
建立回归模型最主要的目的是为了进行预测,但是仅仅给出一个确定并不是一个严谨的预测,很多时候,我们希望得到结果是这个值有多大概率处于某个区间内,因此,统计学上通过严谨的数学推导给出了以下两个区间估计,其中第一个式子是对于响应变量均值的100×(1-α)%的置信区间估计,而第二个式子是新观测值的100×(1-α)%预测区间,很显然地可以看出,均值的区间估计要比单值的区间估计要窄一些。
在预测阶段,经常性会落入一个陷阱,比如我们的回归模型是建立在的区间内的,但是现在需要预测
时
的值,可是我们知道
,那么此时是否可以用上述方法来进行预测呢?这是一个非常严重的错误,因为我们前面得出的回归模型并没有在这个区域内有观测点,所以,根本无从知晓现有模型是否在这个区域依然适用,因此,强烈不建议做超出观测范围的预测推断。
对于本案例,假设α=0.05,当表面积/体积=1.2时,我们可以得出药片释放速度的区间估计如下:
七、线性回归在工业现场的应用场景介绍
以上通过医药研发当中一个案例对一元线性回归进行了简单的介绍,而在实际工作中,大部分响应变量很多时候会受多个因子的影响,对于多元线性回归,其基本思路和分析步骤与一元线性回归类似。同时,也有很多时候,通过残差分析我们发现残差与因子之间存在弯曲现象,这时候需要加入因子的高阶项,甚至会出现因子之间存在强烈的交互作用,这时候就需要在回归模型中加入类似项,甚至,在一些化工领域三传一反的模型中需要加入指数函数,使得模型方程不具有线性特征,如:
此时,我们只要通过简单的变换就可以将问题转化为线性回归的问题,如对于第一个回归方程,我们可以将诸如的高阶项当作一个新的自变量即可,对于第二个式子,我们可以对两边取对数,对
进行变换就可以简化为线性回归的问题。
线性回归在工业现场应用案例非常之多,无论是在生产、运营或者研发等场景中,都有非常多的应用,如化工工艺参数设定对收率的影响,配方对产品性能指标的影响,广告费用对销售收入的影响,不同产线产量与总生产能耗的关系等等,此外,线性回归模型在化学分析、化工装置设计等各个场景也有非常多的应用。
参考文献:
[1].《统计学》,[美]William M. Mendenhall等著
[2].《线性回归分析导论》,[美]Douglas C. Montgomery等著
[3].《深入浅出统计学》,[美]Dawn Griffiths著
[4].《线性代数》,[美]Steven J. Leon著
[5].《六西格玛管理统计指南》,马逢时等著