简单线性回归建模背后的基本目标是从成对的X值和Y值(即X和Y测量值)组成的二维平面中找到最吻合的直线。一旦用最小方差法找到这条直线,就可以执行各种统计测试,以确定这条直线与观测到的Y值的偏离量吻合程度。
线性方程(y=mx+b)有两个参数必须根据所提供的X和Y数据估算出来,它们是斜率(m)和y轴截距(b)。一旦估算出这两个参数,就可以将观测值输入线性方程,并观察方程所生成的Y预测值。
要使用最小方差法估算出m和b参数,就要找到m和b的估计值,使它们对于所有的X值得到的Y值的观测值和预测值最小。观测值和预测值之差称为误差(yi-(mxi+b)),并且,如果对每个误差值都求平方,然后求这些残差的和,其结果是一个被称为预测平方差的数。使用最小方差法来确定最吻合的直线涉及寻找使预测方差最小的m和b的估计值。
可以用两种基本方法来找到满足最小方差法的估计值m和b。第一种方法,可以使用数值搜索过程设定不同的m和b值并对它们求值,最终决定产生最小方差的估计值。第二种方法是使用微积分找到用于估算m和b的方程。我不打算深入讨论推导出这些方程所涉及的微积分,但我确实在SimpleLinearRegression类中使用了这些分析方程,以找到m和b的最小平方估计值(请参阅SimpleLinearRegression类中的getSlope()和getYIntercept方法)。
即使拥有了可以用来找到m和b的最小平方估计值的方程,也并不意味着只要将这些参数代入线性方程,其结果就是一条与数据良好吻合的直线。这个简单线性回归过程中的下一步是确定其余的预测方差是否可以接受。
可以使用统计决策过程来否决“直线与数据吻合”这个备择假设。这个过程基于对T统计值的计算,使用概率函数求得随机大的观测值的概率。正如第1部分所提到的,SimpleLinearRegression类生成了为数众多的汇总值,其中一个重要的汇总值是T统计值,它可以用来衡量线性方程与数据的吻合程度。如果吻合良好,则T统计值往往是一个较大的值;如果T值很小,就应该用一个缺省模型代替您的线性方程,该模型假定Y值的平均值是最佳预测值(因为一组值的平均值通常可以是下一个观测值的有用的预测值)。
要测试T统计值是否大到可以不用Y值的平均值作为最佳预测值,需要计算随机获得T统计值的概率。如果概率很低,那就可以不采用平均值是最佳预测值这一无效假设,并且相应地可以确信简单线性模型是与数据良好吻合的。(有关计算T统计值概率的更多信息,请参阅第1部分。)
回过头讨论统计决策过程。它告诉您何时不采用无效假设,却没有告诉您是否接受备择假设。在研究环境中,需要通过理论参数和统计参数来建立线性模型备择假设。
您将构建的数据研究工具实现了用于线性模型(T测试)的统计决策过程,并提供了可以用来构造理论和统计参数的汇总数据,这些参数是建立线性模型所需要的。数据研究工具可以归类为决策支持工具,供知识工作者在中小规模的数据集中研究模式。
从学习的角度来看,简单线性回归建模值得研究,因为它是理解更高级形式的统计建模的必由之路。例如,简单线性回归中的许多核心概念为理解多次回归(MultipleRegression)、要素分析(FactorAnalysis)和时间序列(TimeSeries)等建立了良好的基础。
简单线性回归还是一种多用途的建模技术。通过转换原始数据(通常用对数或幂转换),可以用它来为曲线数据建模。这些转换可以使数据线性化,这样就可以使用简单线性回归来为数据建模。所生成的线性模型将被表示为与被转换值相关的线性公式。
概率函数
在前一篇文章中,我通过交由R来求得概率值,从而避开了用PHP实现概率函数的问题。我对这个解决方案并非完全满意,因此我开始研究这个问题:开发基于PHP的概率函数需要些什么。
我开始上网查找信息和代码。一个两者兼有的来源是书籍NumericalRecipesinC中的概率函数。我用PHP重新实现了一些概率函数代码(gammln.c和betai.c函数),但我对结果还是不满意。与其它一些实现相比,其代码似乎多了些。此外,我还需要反概率函数。
幸运的是,我偶然发现了JohnPezzullo的InteractiveStatisticalCalculation。John关于概率分布函数的网站上有我需要的所有函数,为便于学习,这些函数已用JavaScript实现。
我将StudentT和FisherF函数移植到了PHP。我对API作了一点改动,以便符合Java命名风格,并将所有函数嵌入到名为Distribution的类中。该实现的一个很棒的功能是doCommonMath方法,这个库中的所有函数都重用了它。我没有花费力气去实现的其它测试(正态测试和卡方测试)也都使用doCommonMath方法。
这次移植的另一个方面也值得注意。通过使用JavaScript,用户可以将动态确定的值赋给实例变量,譬如:
varPiD2=pi()/2
在PHP中不能这样做。只能把简单的常量值赋给实例变量。希望在PHP5中会解决这个缺陷。
请注意清单1中的代码并未定义实例变量—这是因为在JavaScript版本中,它们是动态赋予的值。
清单1.实现概率函数
<?php//Distribution.php//CopyrightJohnPezullo//ReleasedundersametermsasPHP.//PHPPortandOO'fyingbyPaulMeagherclassDistribution{functiondoCommonMath($q,$i,$j,$b){$zz=1;$z=$zz;$k=$i;while($k<=$j){$zz=$zz*$q*$k/($k-$b);$z=$z+$zz;$k=$k+2;}return$z;}functiongetStudentT($t,$df){$t=abs($t);$w=$t/sqrt($df);$th=atan($w);if($df==1){return1-$th/(pi()/2);}$sth=sin($th);$cth=cos($th);if(($df%2)==1){return1-($th+$sth*$cth*$this->doCommonMath($cth*$cth,2,$df-3,-1))/(pi()/2);}else{return1-$sth*$this->doCommonMath($cth*$cth,1,$df-3,-1);}}functiongetInverseStudentT($p,$df){$v=0.5;$dv=0.5;$t=0;while($dv>1e-6){$t=(1/$v)-1;$dv=$dv/2;if($this->getStudentT($t,$df)>$p){$v=$v-$dv;}else{$v=$v+$dv;}}return$t;}functiongetFisherF($f,$n1,$n2){//implementedbutnotshown}functiongetInverseFisherF($p,$n1,$n2){//implementedbutnotshown}}?>
既然您已经用PHP实现了概率函数,那么开发基于PHP的数据研究工具剩下的唯一难题就是设计用于显示分析结果的方法。
简单的解决方案是根据需要将所有实例变量的值都显示到屏幕上。在第一篇文章中,当显示燃耗研究(BurnoutStudy)的线性方程、T值和T概率时,我就是这么做的。能根据特定目的而访问特定值是很有帮助的,SimpleLinearRegression支持此类用法。
然而,另一种用于输出结果的方法是将输出的各部分系统化地进行分组。如果研究用于回归分析的主要统计软件包的输出,就会发现它们往往是用同样的方式对输出进行分组的。它们往往有摘要表(SummaryTable)、偏离值分析(AnalysisOfVariance)表、参数估计值(ParameterEstimate)表和R值(RValue)。类似地,我创建了一些输出方法,名称如下:
- showSummaryTable()
- showAnalysisOfVariance()
- showParameterEstimates()
- showRValues()
- Y轴截距没有有意义的解释,或者
- 输入值可能是经过转换的,而您可能需要取消对它们的转换以获取最终的解释。
所有这些方法都假定输出媒介是网页。考虑到您有可能希望用非网页的其它媒介输出这些汇总值,所以我决定将这些输出方法包装在一个继承了SimpleLinearRegression类的类中。清单2中的代码旨在演示输出类的通用逻辑。为了使通用逻辑更突出,所以除去了实现各种show方法的代码。
清单2.演示输出类的通用逻辑
<?php//HTML.php//Copyright2003,PaulMeagher//DistributedunderGPLinclude_once"slr/SimpleLinearRegression.php";classSimpleLinearRegressionHTMLextendsSimpleLinearRegression{functionSimpleLinearRegressionHTML($X,$Y,$conf_int){SimpleLinearRegression::SimpleLinearRegression($X,$Y,$conf_int);}functionshowTableSummary($x_name,$y_name){}functionshowAnalysisOfVariance(){}functionshowParameterEstimates(){}functionshowFormula($x_name,$y_name){}functionshowRValues(){}}?>
这个类的构造函数只是SimpleLinearRegression类构造函数的包装器。这意味着如果您想显示SimpleLinearRegression分析的HTML输出,则应该实例化SimpleLinearRegressionHTML类,而不是直接实例化SimpleLinearRegression类。其优点是不会有许多未使用的方法充斥SimpleLinearRegression类,并且可以更自由地定义用于其它输出媒介的类(也许会对不同媒介类型实现同一API)。
图形输出
迄今为止,您已经实现的输出方法都以HTML格式显示汇总值。它也适合于用GIF、JPEG或PNG格式显示这些数据的分布图(scatterplot)或线图(lineplot)。
与其亲自编写生成线图和分布图的代码,我认为最好使用名为JpGraph的基于PHP的图形库。JpGraph正由JohanPersson积极开发,其项目网站这样描述它:
无论是对于只有最少代码的“以快捷但不恰当方式获得的”图形,还是对于需要非常细粒度控制的复杂专业图形,JpGraph都可以使它们的绘制变得简单。JpGraph同样适用于科学和商业类型的图形。
JpGraph分发版中包含大量可以根据特定需求进行定制的示例脚本。将JpGraph用于数据研究工具非常简单,只需找到功能与我的需求类似的示例脚本,然后对该脚本进行改写以满足我的特定需求即可。
清单3中的脚本是从样本数据研究工具(explore.php)中抽取的,它演示了如何调用该库以及如何将来自于SimpleLinearRegression分析的数据填入Line和Scatter类。这段代码中的注释是JohanPersson编写的(JPGraph代码库的文档化工作做得很好)。
清单3.来自于样本数据研究工具explore.php的函数的详细内容
<?php//Snippetextractedfromexplore.phpscriptinclude("jpgraph/jpgraph.php");include("jpgraph/jpgraph_scatter.php");include("jpgraph/jpgraph_line.php");//Createthegraph$graph=newGraph(300,200,'auto');$graph->SetScale("linlin");//Setuptitle$graph->title->Set("$title");$graph->img->SetMargin(50,20,20,40);$graph->xaxis->SetTitle("$x_name","center");$graph->yaxis->SetTitleMargin(30);$graph->yaxis->title->Set("$y_name");$graph->title->SetFont(FF_FONT1,FS_BOLD);//makesurethattheX-axisisalwaysatthe//bottomattheplotandnotjustatY=0whichis//thedefaultposition$graph->xaxis->SetPos('min');//Createthescatterplotwithsomenicecolors$sp1=newScatterPlot($slr->Y,$slr->X);$sp1->mark->SetType(MARK_FILLEDCIRCLE);$sp1->mark->SetFillColor("red");$sp1->SetColor("blue");$sp1->SetWeight(3);$sp1->mark->SetWidth(4);//Createtheregressionline$lplot=newLinePlot($slr->PredictedY,$slr->X);$lplot->SetWeight(2);$lplot->SetColor('navy');//Addthepltostotheline$graph->Add($sp1);$graph->Add($lplot);//...andstroke$graph_name="temp/test.png";$graph->Stroke($graph_name);?><imgsrc='<?phpecho$graph_name?>'vspace='15'>?>
该数据研究工具由单个脚本(explore.php)构成,该脚本调用SimpleLinearRegressionHTML类和JpGraph库的方法。
该脚本使用了简单的处理逻辑。该脚本的第一部分对所提交的表单数据执行基本验证。如果这些表单数据通过验证,则执行该脚本的第二部分。
该脚本的第二部分所包含的代码用于分析数据,并以HTML和图形格式显示汇总结果。清单4中显示了explore.php脚本的基本结构:
清单4.explore.php的结构
<?php//explore.phpif(!empty($x_values)){$X=explode(",",$x_values);$numX=count($X);}if(!empty($y_values)){$Y=explode(",",$y_values);$numY=count($Y);}//displayentrydataentryformifvariablesnotsetif((empty($title))OR(empty($x_name))OR(empty($x_values))OR(empty($y_name))OR(empty($conf_int))OR(empty($y_values))OR($numX!=$numY)){//Omittedcodefordisplayingentryform}else{include_once"slr/SimpleLinearRegressionHTML.php";$slr=newSimpleLinearRegressionHTML($X,$Y,$conf_int);echo"<h2>$title</h2>";$slr->showTableSummary($x_name,$y_name);echo"<br><br>";$slr->showAnalysisOfVariance();echo"<br><br>";$slr->showParameterEstimates($x_name,$y_name);echo"<br>";$slr->showFormula($x_name,$y_name);echo"<br><br>";$slr->showRValues($x_name,$y_name);echo"<br>";include("jpgraph/jpgraph.php");include("jpgraph/jpgraph_scatter.php");include("jpgraph/jpgraph_line.php");//Thecodefordisplayingthegraphicsisinlineinthe//explore.phpscript.Thecodeforthesetwolineplots//finishesoffthescript://Omittedcodefordisplayingscatterpluslineplot//Omittedcodefordisplayingresidualsplot}?>
为了演示如何使用数据研究工具,我将使用来自假想的火灾损失研究的数据。这个研究将主要住宅区火灾损失的金额与它们到最近消防站的距离关联起来。例如,出于确定保险费的目的,保险公司会对这种关系的研究感兴趣。
该研究的数据如图1中的输入屏幕所示。
图1.显示研究数据的输入屏幕
数据被提交之后,会对它进行分析,并显示这些分析的结果。第一个显示的结果集是TableSummary,如图2所示。
图2.TableSummary是所显示的第一个结果集
TableSummary以表格形式显示了输入数据和其它列,这些列指出了对应于观测值X的预测值Y、Y值的预测值和观测值之间的差以及预测Y值置信区间的下限和上限。
图3显示了TableSummary之后的三个高级别数据汇总表。
图3.显示了TableSummary之后的三个高级别数据汇总表
AnalysisofVariance表显示了如何将Y值的偏离值归为两个主要的偏离值来源,由模型解释的方差(请看Model行)和模型不能解释的方差(请看Error行)。较大的F值意味着该线性模型捕获了Y测量值中的大多数偏离值。这个表在多次回归环境中更有用,在那里每个独立变量都在表中占有一行。
ParameterEstimates表显示了估算的Y轴截距(Intercept)和斜率(Slope)。每行都包括一个T值以及观测到极限T值的概率(请看Prob>T列)。斜率的Prob>T可用于否决线性模型。
如果T值的概率大于0.05(或者是类似的小概率),那么您可以否决该无效假设,因为随机观测到极限值的可能性很小。否则您就必须使用该无效假设。
在火灾损失研究中,随机获得大小为12.57的T值的概率小于0.00000。这意味着对于与该研究中观测到的X值区间相对应的Y值而言,线性模型是有用的预测器(比Y值的平均值更好)。
最终报告显示了相关性系数或R值。可以用它们来评估线性模型与数据的吻合程度。高的R值表明吻合良好。
每个汇总报告对有关线性模型和数据之间关系的各种分析问题提供了答案。请查阅Hamilton、Neter或Pedhauzeur编写的教科书,以了解更高级的回归分析处理。
要显示的最终报告元素是数据的分布图和线图,如图4所示。
图4.最终报告元素—分布图和线图
大多数人都熟悉线图(如本系列中的第一幅图)的说明,因此我将不对此进行注释,只想说JPGraph库可以产生用于Web的高质量科学图表。当您输入分布或直线数据时,它也做得很好。
第二幅图将残差(观测的Y、预测的Y)与您预测的Y值关联起来。这是研究性数据分析(ExploratoryDataAnalysis,EDA)的倡导者所使用的图形示例,用以帮助将分析人员对数据中的模式的检测和理解能力提到最高程度。行家可以使用这幅图回答关于下列方面的问题:
- 可能的非正常值或影响力过度的例子
- 可能的曲线关系(使用转换?)
- 非正态残差分布
- 非常量误差方差或异方差性
数学库体系结构
对数学的业余爱好使我在最近几个月中保持着对数学库的浓厚兴趣。此类研究推动我思考如何组织我的代码库以及使其预期在未来能不断增长。
我暂时采用清单5中的目录结构:
清单5.易于增长的目录结构
phpmath/burnout_study.phpexplore.phpfire_study.phpnavbar.phpdist/Distribution.phpfisher.phpstudent.phpsource.phpjpgraph/etc...slr/SimpleLinearRegression.phpSimpleLinearRegressionHTML.phptemp/
例如,未来有关多次回归的工作,将涉及扩展这个库以包括matrix目录,该目录用来容纳执行矩阵操作(这是对于更高级形式的回归分析的需求)的PHP代码。我还将创建一个mr目录,以容纳实现多次回归分析输入方法、逻辑和输出方法的PHP代码。
请注意这个目录结构包含一个temp目录。必须设置该目录的许可权,使explore.php脚本能够将输出图写到该目录。在尝试安装phpmath_002.tar.gz源代码时请牢记这一点。此外,请在JpGraph项目网站上阅读安装JpGraph的指示信息(请参阅参考资料)。
最后提一点,如果采取以下作法,可以将所有软件类移到Web根目录之外的文档根目录:
- 使某个全局PHP_MATH变量有权访问非Web根目录位置,并且
- 确保在所有需要或包括的文件路径前面加上这个已定义的常量作为前缀。
您学到了什么?
在本文中,您了解了如何使用SimpleLinearRegression类开发用于中小规模的数据集的数据研究工具。在此过程中,我还开发了一个供SimpleLinearRegression类使用的本机概率函数,并用HTML输出方法和基于JpGraph库的图形生成代码扩展该类。
从学习的角度来看,简单线性回归建模是值得进一步研究的,因为事实证明,它是理解更高级形式的统计建模的必由之路。在深入学习更高级的技术(如多次回归或多变量方差分析)之前,对于简单线性回归的透彻理解将使您受益匪浅。
即使简单线性回归只用一个变量来说明或预测另一个变量的偏离值,在所有的研究变量之间寻找简单线性关系仍然常常是研究性数据分析的第一步。仅因为数据是多元的并不意味着就必须使用多元工具研究它。实际上,在开始时使用简单线性回归这样的基本工具是着手探究数据模式的好方法。
(举报)