第2章 直线回归分析
2.1 问题与数据
【例2-1】 某研究组测定了20名儿童体内血红蛋白与血铁的位置互换含量,数据见表2-1,试采用常规的SAS编程进行直线回归分析,以得到由血红蛋白含量预测血铁含量的回归方程。
表2-1 20名儿童体内的血红蛋白与血铁含量
【例2-2】 内容同例2-1资料,试采用高级SAS编程技术实现智能化的直线回归分析,以得到由血红蛋白含量预测血铁含量的回归方程。
【例2-3】 在突变实验中,用不同剂量的射线照射植物的种子,发现苗期高度与成活株之间有一定的联系。用X射线照射大麦的种子,记处理株第一叶平均高度与对照株高度的比值为x,存活百分率为y,结果如表2-2所示,试采用直线回归分析,检验两者之间是否存在一定的联系。
表2-2 大麦处理株第一叶平均高度与对照株高度的比值与存活百分率
【例2-4】 内容同例2-3资料,试采用高级SAS编程技术实现智能化的直线回归分析,以检验两者之间是否存在一定的联系。
2.2.1 对例2-1资料的分析
2.2 分析与解答
本资料中,研究者欲得到由血红蛋白含量预测血铁含量的回归方程。一般而言,研究一个定量变量随另一个定量变量的依存变化关系,若两变量基本呈直线关系,可采用直线回归分析,否则应进行相应的变量变换或采用合适的曲线回归分析。为演示如何借助SAS软件进行直线回归分析,此处直接采用此分析方法。SAS程序如下,程序名为nr2 1.sas。
datanr2_1; /*1*/ input x y@ @ ; cards; 13.50 518.70 13.00 467.30 11.00 469.80 14.30 456.60 12.50 448.70 12.50 424.10 11.80 405.60 11.50 446.00 11.00 416.70 10.70 430.80 10.20 409.80 10.00 384.10 9.50 356.30 9.40 388.60 8.80 325.90 6.30 292.80 7.30 332.80 7.80 283.00 7.30 312.50 7.00 294.70 ; run; ods html; proc gplot; /*2*/ plot y* x; run; proc reg; /*3*/ model y=x/r p; run; ods html close;
【程序说明】 程序第1步为建立数据集nr2 1,以x表示原因变量,以y表示结果变量。第2步为调用GPLOT过程,以绘制该资料的散点图(见图2-1)。第3步调用REG过程,进行含有截距项的直线回归分析,其中“r”选项和“p”选项分别用来进行残差分析和输出预测值。
图2-1 例2-1资料的散点图
【SAS输出结果及其解释】
由图2-1可以看出,散点基本上都在一条既不平行于x轴,也不平行于y轴的不太宽的条带中随机分布着,故可以采用直线回归分析。以下是直线回归分析的结果。
The REG Procedure
Model:MODEL1
Dependent Variable:y
Number of Observations Read 20
Number of Observations Used 20
这是模型的基本信息,共有观测点20个。
Analysis of Variance
这是对模型整体有无统计学意义进行假设检验的结果。F=98.06,对应的P<0.0001,说明拟合的直线回归方程在整体上是有统计学意义的。
Root MSE 27.60259 R-Square 0.8449
Dependent Mean 393.24000 Adj R-Sq 0.8363
Coeff Var 7.01927
这是评价模型拟合效果的一些统计量。“Root MSE”表示误差均方根,其值越小,反映模型拟合效果就越好。“R-Square”表示决定系数,其值越大,反映模型拟合效果越好,一般其值大于0.5时,模型方有实际的使用价值。“Adj R-Sq”表示校正的R2,其值越大反映模型拟合效果越好。
这是参数估计的结果。从前至后依次是变量名、自由度、参数估计值、标准误、t值及P值。可以看出,截距项和原因变量x的回归系数对应的P值分别为0.0007和<0.0001,即二者均有统计学意义。故可构造回归方程如下:
y = 116.58837 +26.93784x
这是模型对观测的预测值及残差分析的结果。由左至右,各列的含义依次为:观测号、结果变量的实际值、结果变量的预测值、结果变量预测值的标准误、实际值与预测值之间的残差、残差的标准误、学生化残差、学生化残差图和Cook's D统计量。其中,学生化残差大于2或Cook's D统计量大于0.5时,对应的观测可能是异常点,此时需判断有无过失误差,并结合专业知识综合判断该观测是否是异常点,如有可能,最好在此观测点上重新进行一次试验。本资料中,3号观测点可能是异常点,如无专业依据认为该观测点需要保留时,可删除该观测点,重新进行直线回归分析。学有余力的读者,可尝试进行分析。
Sum of Residuals 0
Sum of Squared Residuals 13714
Predicted Residual SS(PRESS)17289
这是回归方程预测效果的统计量。“Sum of Residuals”表示残差之和,“Sum of Squared Residu-als”表示残差平方和,“Predicted Residual SS(PRESS)”表示预测残差平方和。
2.2.2 对例2-2资料的分析
对例2-2资料采用高级SAS编程技术,实现智能化的直线回归分析,SAS程序如下,程序名为nr2 2.sas。
datanr2_2; /*1*/ input x y@ @ ; cards; 13.50 518.70 13.00 467.30 11.00 469.80 14.30 456.60 12.50 448.70 12.50 424.10 11.80 405.60 11.50 446.00 11.00 416.70 % linear_reg(nr2_2); /*2*/ 10.70 430.80 10.20 409.80 10.00 384.10 9.50 356.30 9.40 388.60 8.80 325.90 6.30 292.80 7.30 332.80 7.80 283.00 7.30 312.50 7.00 294.70 ; run;
【程序说明】 程序第1步为建立数据集nr2 2,以x表示原因变量,以y表示结果变量,这样可节省第2步中所调用宏程序中宏参数的数量。第2步为调用已编制好的宏程序linear reg,此宏程序可自动识别截距项是否有统计学意义,从而直接给出最终的分析结果,无须使用者调试修改。此宏包含一个宏参数,即数据集的名称,调用时需由使用者给定。
【SAS输出结果及其解释】
例2-2资料的散点图如图2-2所示。
图2-2 例2-2资料的散点图
由图2-2可以看出,散点基本上都在一条既不平行于x轴,也不平行于y轴的不太宽的条带中随机分布着,故可以采用直线回归分析。该程序之后输出的统计分析结果与2.2.1小节一致,不再单独列出,读者可参阅之前的内容。
图2-3是所得直线回归方程对例2-2资料的拟合效果图。可见,拟合效果较好。
图2-3 例2-2资料进行直线回归分析的拟合效果图
需要说明的是,由于例2-1资料拟合的直线回归方程中截距项有统计学意义,所以SAS程序nr2 2并未体现出其智能化的优点。为比较两种编程方式在智能化程度上的差异,下面再以例2-3资料进行说明。
2.2.3 对例2-3资料的分析
本资料中,研究者欲得到大麦处理株第一叶平均高度和对照株高度的比值与存活百分率之间的关系。一般而言,研究一个定量变量随另一个定量变量的依存变化关系,若两变量基本呈直线关系,可采用直线回归分析,否则应进行相应的变量变换或采用合适的曲线回归分析。为演示如何借助SAS软件进行直线回归分析,此处直接采用此分析方法。SAS程序如下,程序名为nr2 3.sas。
datanr2_3; /*1*/ input x y@ @ ; cards; 40 60 41 61 41 67 43 68 47 70 51 85 52 88 55 89 61 90 62 100 ; run; ods html; proc gplot; /*2*/ plot y* x; run; proc reg; /*3*/ model y=x; run; proc reg; /*4*/ model y=x/noint r p; run; ods html close;
【程序说明】 程序第1步为建立数据集nr2 3,以x表示原因变量,以y表示结果变量。第2步为调用GPLOT过程,以绘制该资料的散点图(见图2-4)。第3步调用REG过程,进行含有截距项的直线回归分析,查看下面的SAS输出结果可知,截距项与0之间的差异并无统计学意义,故仍需进行不含截距项的直线回归分析。第4步调用REG过程,进行不含截距项的直线回归分析,其中“noint”选项表示模型中不包含截距项,“r”选项和“p”选项分别用来进行残差分析和输出预测值。
图2-4 例2-3资料的散点图
【SAS输出结果及其解释】
由图2-4可以看出,散点基本上都在一条既不平行于x轴,也不平行于y轴的不太宽的条带中随机分布着,故可以采用直线回归分析。以下是直线回归分析的结果。
首先输出的是含有截距项时,直线回归分析的结果。
The REG Procedure
Model:MODEL1
Dependent Variable:y
Number of Observations Read 10
Number of Observations Used 10
这是模型的基本信息,共有观测点10个。
这是对模型整体有无统计学意义进行假设检验的结果。F=89.07,对应的P<0.0001,说明拟合的直线回归方程在整体上是有统计学意义的。
这是评价模型拟合效果的一些统计量。“Root MSE”表示误差均方根,其值越小,反映模型拟合效果就越好。“R-Square”表示决定系数,其值越大,反映模型拟合效果越好。“Adj R-Sq”表示校正的R2,其值越大反映模型拟合效果越好。
这是参数估计的结果。从前至后依次是变量名、自由度、参数估计值、标准误、t值及P值。可以看出,截距项对应的P值为0.7209>0.05,说明截距项与0的差异没有统计学意义。故应删除截距项,重新进行不含截距项的直线回归分析。
以下输出的是不含截距项时直线回归分析的结果。
Number of Observations Read 10
Number of Observations Used 10
Note:No intercept in model. R-Square is redefined.
这是模型的基本信息,共有观测点10个。所得的模型中不包含截距项,此时R2的计算方法有所调整。
这是对模型整体有无统计学意义进行假设检验的结果。F=3717.43,对应的P<0.0001,说明拟合的直线回归方程在整体上是有统计学意义的。
这是模型拟合效果的一些统计量。需要注意的是,此模型中不包含截距项,其R2的计算方法不同于含有截距项的回归模型,故两类模型不宜直接依据R2的大小比较它们的拟合效果。
这是参数估计的结果。从前至后依次是变量名、自由度、参数估计值、标准误、t值及P值。可以看出,原因变量x的回归系数对应的P<0.0001,有统计学意义。故可构造回归方程如下:
y = 1.57969x
这是模型对观测的预测值及残差分析的结果。由左至右,各列的含义依次为:观测号、结果变量的实际值、结果变量的预测值、结果变量预测值的标准误、实际值与预测值之间的残差、残差的标准误、学生化残差、学生化残差图和Cook's D统计量。其中,学生化残差大于2或Cook's D统计量大于0.5时,对应的观测可能是异常点,此时需判断有无过失误差,并结合专业知识综合判断该观测点是否是异常点,如有可能,最好在此观测点上重新进行一次试验。
Sum of Residuals -0.78764
Sum of Squared Residuals 150.52362
Predicted Residual SS(PRESS)190.70070
这是回归方程预测效果的统计量。“Sum of Residuals”表示残差之和,“Sum of Squared Residu-als”表示残差平方和,“Predicted Residual SS(PRESS)”表示预测残差平方和。
2.2.4 对例2-4资料的分析
对例2-4资料采用高级SAS编程技术,实现智能化的直线回归分析,SAS程序如下,程序名为nr2 4.sas。
datanr2_4; /*1*/ input x y@ @ ; cards; 40 60 41 61 41 67 43 68 47 70 51 85 52 88 55 89 61 90 62 100 ; run; % linear_reg(nr2_4); /*2*/
【程序说明】 程序第1步为建立数据集nr2 4,以x表示原因变量,以y表示结果变量,这样可节省第2步中所调用宏程序中宏参数的数量。第2步为调用已编制好的宏程序linear reg,此宏程序可自动识别截距项是否有统计学意义,从而直接给出最终的分析结果,无须使用者调试修改。此宏包含一个宏参数,即数据集的名称,调用时需由使用者给定。
【SAS输出结果及其解释】
例2-4资料的散点图如图2-5所示。
图2-5 例2-4资料的散点图
由图2-5可以看出,散点基本上都在一条既不平行于x轴,也不平行于y轴的不太宽的条带中随机分布着,故可以采用直线回归分析。该程序之后输出的统计分析结果与2.2.3小节中不含截距项时直线回归分析的结果一致,不再单独列出,读者可参阅之前的内容。
图2-6是所得直线回归方程对例2-4资料的拟合效果图,可见,拟合效果较好。
图2-6 例2-4资料的拟合效果图
2.3 计算原理
直线回归分析是用直线回归方程表示两个定量变量间依存关系的统计分析方法。此分析方法的基本思路是:首先由实际测量得到的数据计算反映两定量变量依赖关系的直线回归方程,即计算直线回归方程的截距a、斜率b;然后对回归方程进行假设检验,也就是检验总体回归系数是否为零;当然,还要检验总体截距是否为零,以决定回归直线是否通过坐标原点;最后结合专业知识,评价此直线回归方程是否有实用价值。
如果用X表示原因变量,Y表示结果变量,直线回归分析要求资料符合下列条件:
(1)X和Y之间的关系为线性关系;
(2)各个观察值之间相互独立;
(3)每个X值对应的Y服从正态分布;
(4)不同X值对应的Y之分布具有相同的方差,换句话说Y的方差与X无关。
直线回归分析的基本作用包括以下三条。
(1)描述两变量之间的依赖关系:若认为两变量间存在着直线回归关系,则可用直线回归方程来描述;
(2)利用回归方程进行预测:所谓预测就是把预报因子(原因变量X)代入回归方程对预报量(结果变量Y)进行估计;
(3)利用回归方程进行统计控制:统计控制是利用回归方程进行逆估计,如要求Y在一定范围内波动,可以通过控制X的取值来实现。
2.3.1 参数估计
进行直线回归分析的两个变量一般有原因变量和结果变量之分,当在专业上无法区分时,常把容易测量的变量看作原因变量,另一个较难测量的变量看作结果变量。结果变量Y与原因变量X之间的回归模型为:
式中,α表示回归模型的截距,β为回归模型的斜率,常称之为回归系数,ε为随机误差。如果用Y表示Y的估计值,用a与b分别表示α和β的样本估计值,则样本回归方程可以表示为:
式中, 就是由回归方程求得的结果变量Y的回归值,也称为Y的预测值,它是当X固定时,Y的总体均数μY|X的估计值。a,b可以根据最小二乘法原理求得,被称为总体参数α,β的最小二乘估计量。用最小二乘法求参数估计量的基本原理是使回归直线上各估计值 与观察值Y之差的平方和达到最小,估计值 与观察值Y之差通常称为残差,α,β的最小二乘估计量a,b的具体计算公式见式(2-3)和式(2-4):
如果根据专业知识需求过定点(X0,Y0)的直线回归方程,则只需将式(2-3)和式(2-4)中的X —换成X0、Y —换成Y0即可。其中,
,表示X的离均差平方和;
,表示Y的离均差平方和;
,表示X与Y的离均差积和。
2.3.2 假设检验
1.对总体回归系数的假设检验
在估计出样本回归系数b以后,仍然需要对总体回归系数β是否为0进行假设检验,检验的原假设和备择假设分别为H0:β=0;H1:β≠0,a=0.05。检验方法有方差分析和t检验两种。
在进行方差分析时,首先要对结果变量Y的离均差平方和进行分解。总的离均差平方和为观察值Y与均数Y —之差的平方和∑(Y - Y —)2,用SST表示。它被分解为两部分,一部分是估计值 与均数Y —之差的平方和∑( - Y —)2,称为回归平方和,用SSR表示,它反映了Y的总变异中可以用X与Y的直线关系解释的那部分变异,回归平方和越大,说明回归的效果越好;另一部分是观察值Y与估计值Y 之差的平方和∑(Y - )2,称为残差平方和或剩余平方和,用SSE表示,它反映了X对Y的线性影响之外的一切因素对Y的变异的作用,也就是总变异中无法用X解释的部分,残差平方和越小,说明直线回归的估计误差越小,回归的作用越明显。离均差平方和的分解可以用下式表示:
相应的自由度也可以分解为:
式中,νT=n-1,νR=1,νE=n-2。
在对离均差平方和进行分解之后,可以按下式计算检验统计量F:
式中,MSR与MSE分别称为回归均方和残差均方。求出F值后查F界值表,就可以得到P值,然后,根据具体情况作出统计和专业结论。
通过上述统计量,还可以引入决定系数(R2),它是回归平方和与总的离均差平方和之比:
决定系数也就是相关系数的平方,它是相关与回归分析中一个十分有用的统计量,反映了回归的贡献大小,说明了在Y的总变异中原因变量X所能解释的百分比,其取值越接近于1,说明回归的效果越好,即原因变量X对结果变量Y的贡献就越大。
对回归系数的检验也可以采用t检验,检验统计量为:
式(2-9)中,Sb是样本回归系数的标准误,SY·X为剩余标准差,是排除了X的影响后,单独Y方面的变异大小,常用它作为预报精确度的标志。因为它的单位与Y一致,容易在实际中进行比较和检验,所以,一个回归能否对解决实际问题有所帮助,只要比较SY·X与允许的偏差即可,故它是检验一个回归是否有效的极其重要的标志。
在对回归系数的假设检验中,方差分析和t检验是等价的,并且有t2 =F。
2.对总体截距的假设检验
在估计出样本的截距a 以后,仍然需要对总体截距 α 是否为零进行假设检验,检验的原假设和备择假设分别为H0:α=0;H1:α≠0;α(显著性水平,不是总体截距)=0.05。检验的方法为t 检验,计算公式如下:
式中,Sa为a的标准误,SY·X,lxx与式(2-10)、式(2-3)中相应符号的含义相同。
2.3.3 区间估计
1. 总体回归系数β的置信区间
总体回归系数β的100(1-α)%(此处的α为显著性水平,不是总体截距)置信区间可以按下式计算:
式中,Sb为样本回归系数的标准误,ν=n-2。
2. 总体截距α的置信区间
总体截距α的100(1-α)%(后一个α为显著性水平,不是总体截距)置信区间可以按下式计算:
式中,Sa为样本截距的标准误,ν=n-2。
3. 给定X=X0条件下Y的总体均数的置信区间
若记μY|X=X0 为给定X=X0条件下Y的总体均数,则它的100(1-α)%置信区间为:
式中, 为给定X=X0时Y的估计值, 为 的标准误,按照下式计算:
4. 给定X=X0条件下Y的个体值的预测区间
给定X=X0时,Y的个体值的100(1-α)%预测区间为:
预测区间反映了个体Y值的波动范围,解决了对结果变量进行预报的问题。式中SY0 为给定X=X0时Y值的标准差,其计算公式为:
2.4 分析步骤
2.4.1 分析要领
1.受试对象应具有同质性
在分析前,应明确欲分析的受试对象是否同质。例如,受试对象是否是相同的病情、病程、接受相同的处理等。由医学专业知识和常识可知,人的胸围与肺活量之间有一定的联系,是否一定可以用直线回归分析研究他们之间的关系呢?那要看测量谁的胸围与肺活量了,若测量同一批成年男性的胸围与肺活量,则可以加以研究;但是,若测量一组成年男性的胸围与肺活量,同时测另一组成年女性的胸围与肺活量,将这样的两组数据放在一起进行直线回归分析就要慎重了,因为男女胸围与肺活量的关系不一定相同。
另外,如果将两个或两个以上样本合并在一起进行相关和回归分析,有可能得出实际上并不存在相关或回归关系,也可能掩盖确实存在的相关与回归关系,如图2-7所示。因此,将不同的样本放在一起进行相关与回归分析时应慎重。
图2-7 两个不同样本资料对相关的影响
关于图2-7的说明:
(1)左图中圆圈和叉号分别代表性质不同的两类受试对象,每一类受试对象的散点呈一个“圆盘状”,表明这两类的内部两变量之间都没有线性关系,而盲目地将它们拼接在一起,碰巧出现了“正相关趋势”,也就很容易错误地得出两变量之间呈现直线相关或回归关系的结论;
(2)与左图情况正好相反,右图中将两组分别有直线相关或回归关系的资料合并成了一组无直线相关或回归关系的资料,即由于错误地将性质不同的两类受试对象放在一起分析,将每一类样本中实际存在的直线相关或回归关系掩盖掉了。
2.应以专业知识为依据
统计分析方法只能帮助人们揭示数据之间内在的规律,而不能创造规律。变量之间是否存在本质联系要靠专业知识来解释。例如,研究“随时间的推移,用小孩的身高去推测小树苗高度的直线回归分析”和“根据人的年龄去推测人被犬咬伤的发生率的直线回归分析”都是毫无专业根据的。
3.应绘制并正确分析散布图
在进行直线回归分析之前,应该绘制散布图。专业上有联系的两项指标之间的关系并非都是直线关系,事实上,如果两项指标之间呈一条弯曲度不大的“S”型或反“S”型曲线趋势,错误地用一条直线回归方程来描述,在统计学上往往会认为很有意义,即该直线回归方程是成立的,但生物学上是解释不通的(当结果变量是某种率时最易发生这种现象)。正确的做法是:在直角坐标系内绘制X与Y的散布图,如果n个点形成的散布图呈一条明显的曲线趋势时,宜拟合一条曲线回归方程;如果n个点在一条不太宽的长带内随机地分布着,不存在明显的曲线趋势且此带不平行于X轴也不垂直于X轴,可考虑进行直线回归分析;如果n个点形成的散布图近似于一个圆盘,说明X与Y之间无确切的变化趋势,几乎是相互独立的,不必硬把它们捏合在一起进行分析。
4.应善于发现可疑的异常点
通过细心观察散布图,可以发现可疑的异常点。异常点出现的原因有很多,如过失误差所致(仪器出故障、试剂失效、读写或录入错误)、某些个体与其他绝大多数个体不同质、突变点(这或许意味着产生“奇迹”的机会到了)。
5.应进行正确的计算
正确求解直线回归方程中参数的点估计值和区间估计值,并对参数进行假设检验。
6.应考察拟合优度
应考察拟合误差是否在专业上允许的范围内;应考虑用某种曲线回归方程取代直线回归方程,是否会大幅度提高拟合的效果。
7.应考察直线回归方程是否与实际情况相违背
直线回归方程中斜率的正负号在专业上是否能得到合理的专业解释;结果变量Y的预测值是否超出了专业上允许的范围。
8.应考察直线回归分析是否具有实际意义
设反映两定量变量之间呈相互关系密切程度的系数,即相关系数为r,当r2 > 0.5时才表明原因变量X对结果变量Y的贡献率大于50%,即结果变量Y的总变异中,有50%以上的信息可由原因变量X来控制。
2.4.2 基本步骤
第一步,看变量之间在专业上是否有联系,有联系方有必要进行分析。
第二步,看分析的对象是否具有同质性,若满足同质性要求则可进行下一步。
第三步,绘制散布图,一是要敏锐地发现是否存在可疑的“异常点”,若存在这种点,应考虑在此类点存在和不存在两类情形下分别做分析,并加以讨论;二是看散布图的变化趋势,若呈直线变化趋势则可进行直线回归分析,若呈某种曲线变化趋势,则应考虑进行曲线回归分析。
第四步,计算出截距和斜率,并对截距和斜率进行假设检验。
第五步,考察所得直线方程有无实际意义。结合专业知识及决定系数r2 是否大于0.5(r2 越接近于1越好,但至少应大于0.5)对此方程有无实用价值做出科学评价。
2.4.3 分析策略
直线回归分析较为简单,具体操作时,只有含有截距项和不含截距项两种情形。所以,可在两种情形下分别进行分析,然后从所得的两个直线回归模型中选出对资料拟合最好的那个模型。该策略的框图见图2-8。
图2-8 直线回归分析时最优模型的拟合策略框图
2.4.4 实施步骤
具体的实施步骤如下:
第一步,建立SAS数据集。为减少使用者需要填写的宏参数的数量,这里规定,以y表示结果变量,以x表示原因变量。
第二步,分别从含截距项和不含截距项两方面入手,进行直线回归分析。
第三步,根据对模型的假设检验结果,选择输出最终模型。若参数与0之间的差异全部有统计学意义的模型只有一个,则该模型即为最优模型;若两个模型中参数与0之间的差异全部都有统计学意义,则选择含有截距项的那个模型作为最优模型;若两个模型中参数与0的差异均不是全部有统计学意义,则给出有截距项的那个模型,作为参考模型,供使用者参考。
参考文献
[1] SAS Institute Inc. SAS/STAT 9.2 User's Guide.Cary,NC:SAS Institute Inc.,2008:3253-3474,5427-5640.
[2] 孙振球. 医学统计学. 北京:人民卫生出版社,2002:140-168,241-271.
[3] 柳青.中国医学统计百科全书(多元统计分册).北京:人民卫生出版社,2004:1-9.
[4] 胡良平. 心血管病科研设计与统计分析. 北京:人民军医出版社,2010:145-172.
[5] 胡良平. 面向问题的统计学——(1)科研设计与统计基础. 北京:人民卫生出版社,2012:564-576.