第3章 可直线化的曲线回归分析
在生物医学试验中,两变量间往往并不呈直线关系,而是呈非线性关系或称曲线关系,这时拟合恰当的曲线比拟合直线更为合适。如细菌繁殖量与培养时间的关系,疾病发病率与时间的关系,毒物剂量与效应之间的关系等均非直线关系。本章介绍如何借助SAS软件实现两变量可直线化的曲线回归分析。
3.1 问题与数据
【例3-1】 上海医科大学微生物教研室以已知浓度的免疫球蛋白A(IgA,μg)做火箭电泳,测得火箭的高度(mm),结果如表3-1所示,试拟合其曲线回归方程。
表3-1 IgA浓度与火箭高度
【例3-2】 基于例3-1分析的结果,将所得的最终曲线与原散点图绘制在同一直角坐标系中,以直观展现其拟合效果。若存在多条拟合曲线,试进行回归曲线的拟合效果比较。
【例3-3】 某研究者收集了一只红铃虫在不同区域上的产卵数y和温度x之间的7组观测数据,结果如表3-2所示。试拟合其曲线回归方程。
表3-2 不同温度下红铃虫的产卵数
【例3-4】 基于例3-3分析的结果,将所得的最终曲线与原散点图绘制在同一直角坐标系中,以直观展现其拟合效果。若存在多条拟合曲线,试进行回归曲线的拟合效果比较。
【例3-5】 某研究组调查某县1980~1996年月累计麻疹的发病情况,观察指标为发病率(1/105),结果如表3-3所示,试拟合其曲线回归方程。
表3-3 某县1980~1996年月累计麻疹的发病情况
【例3-6】 基于例3-5分析的结果,将所得的最终曲线与原散点图绘制在同一直角坐标系中,以直观展现其拟合效果。若存在多条拟合曲线,试进行回归曲线的拟合效果比较。
3.2 分析与解答
3.2.1 对例3-1资料的分析
对于例3-1资料,绘制其散点图(见图3-1)可发现,散点的变化趋势近似于对数曲线、双曲线或幂函数曲线(可参见3.3.1小节),故可采用对数变换法或倒数法使其直线化。SAS程序如下,设程序名为nr3 1.sas。
图3-1 IgA浓度与火箭高度之间的散点图
data nr3_1; /*1*/ input x y@ @ ; x1=log10(x); x2=1/x; x3=log(x); y1=1/y; y2=log(y); cards; 0.2 7.6 0.4 12.3 0.6 15.7 0.8 18.2 1.0 18.7 1.2 21.4 1.4 22.6 1.6 23.8 ; run; ods html; proc gplot; /*2*/ plot y* x/haxis=0 to 1.6 by 0.4 vaxis=0 to 25 by 5; symbol value=dot; run; proc reg; /*3*/ model y=x1; model y1=x2; model y2=x3; run; ods html close;
【程序说明】 程序第1步是建立数据集,x代表IgA浓度,y代表火箭高度,并对两变量进行各种变换,以便进行曲线直线化分析。第2步是绘制散点图,以判断散点趋势。第3步是进行各种回归分析,包括对数曲线、双曲线和幂函数曲线的直线化分析。
【SAS输出结果及其解释】
图3-1是根据原始资料绘制的IgA浓度和火箭高度之间数据对的散点图,可发现其趋势与对数曲线、双曲线或幂函数曲线很相似,故对两变量进行相应变换。
以上是结果变量火箭高度(y)和IgA浓度(x)经对数变换后产生的新变量(x1)之间进行直线回归分析的结果,给出了拟合的直线回归方程中的截距项和斜率的估计值,以及对它们与零的差别是否具有统计学意义进行假设检验的结果。可以看出,截距和斜率与零的差别均有统计学意义。因此,直线回归方程为: = 19.74512 +17.90734x1。因x1 = log10(x),故 = 19.74512 +17.90734 log10(x)。
以上是结果变量火箭高度(y)经倒数变换后产生的新变量(y1)和IgA浓度(x)经倒数变换后产生的新变量(x2)之间进行直线回归分析的结果。可以看出,截距和斜率与0的差别均有统计学意义。因此,直线回归方程为:y 1 = 0.03029 +0.02029x2。因x2=1/x,y1=1/y,故y = x/(0.03029x +0.02029)。
以上是结果变量火箭高度(y)经对数变换后产生的新变量(y2)和IgA浓度(x)经对数变换后产生的新变量(x3)之间进行直线回归分析的结果。可以看出,截距和斜率与零的差别均有统计学意义。因此,直线回归方程为: = 2.96139 +0.53667x3。因x3 = lnx,y2 = lny,故 = exp(2.96139 +0.53667lnx)= 19.32481x0.53667。
3.2.2 对例3-2资料的分析
基于3.2.1小节的分析结果,得到了三条回归曲线,现比较一下三条回归曲线的拟合效果,并将三条曲线与原散点绘制在同一直角坐标系中,以直观展现其拟合效果。SAS程序如下,设程序名为nr3 2.sas。
% let n=8; /*1*/ data nr3_2; /*2*/ input x y@ @ ; y1=19.74512+17.90734* log10(x); y2=x/(0.03029* x+0.02029); y3=19.32481* x** 0.53667; residual1=y1-y; residual2=y2-y; residual3=y3-y; ssrs1=residual1** 2; proc gplot data=nr3_2; /*3*/ plot y* x y1* x y2* x y3* x/ haxis=0 to 1.6 by 0.4 vaxis=0 to 25 by 5 overlay legend=legend1; run; proc sql; /*4*/ create table test1 as select css(y)as total,sum(ssrs1) ssrs2=residual2** 2; ssrs3=residual3** 2; cards; 0.2 7.6 0.4 12.3 0.6 15.7 0.8 18.2 1.0 18.7 1.2 21.4 1.4 22.6 1.6 23.8 ; run; symbol1 color = green value =dot; symbol2 color=red value=point line=1 interpol=spline; symbol3 color=blue value=point line=2 interpol=spline; symbol4 color=green value=point line=3 interpol=spline; legend1 across=1 cborder=black label=none value=('实际值''对数曲线' '双曲线''幂曲线') position=(bottom right inside) mode=protect; as sum_resi1,sum(ssrs2) as sum_resi2,sum(ssrs3) as sum_resi3 from nr3_2; quit; data test; /*5*/ set test1; s1=sqrt(sum_resi1/(&n-2)); s2=sqrt(sum_resi2/(&n-2)); s3=sqrt(sum_resi3/(&n-2)); r_square1=1-sum_resi1/total; r_square2=1-sum_resi2/total; r_square3=1-sum_resi3/total; run; ods html; proc print data=nr3_2; /*6*/ run; proc print data=test noobs; /*7*/ run; ods html close;
【程序说明】 第1步是设立宏变量n,其含义为样本个数。第2步是创建数据集,并计算对数曲线y = 19.74512 +17.90734 log10(x)拟合资料的残差residual1及残差平方ssrs1,计算双曲
线y = x/(0.03029x+0.02029)拟合资料的残差residual2及残差平方ssrs2,计算幂函数曲线y =19.32481 × x0.53667 拟合资料的残差residual3及残差平方ssrs3。第3步是将三条曲线与原散点绘制在同一直角坐标系中,以直观展现它们对资料的拟合效果。第4步是计算变量y的离均差平方和,以及三条曲线拟合资料的残差平方和。第5步是计算各曲线拟合资料的剩余标准差(含义见3.3.4小节)和相关指数(含义见3.3.5部分)。第6步是将第2步所得到的数据集输出到output窗口。第7步输出数据集test至output窗口。
【SAS输出结果及其解释】
图3-2是对数曲线、双曲线和幂函数曲线对例3-1资料拟合效果的展示图。
图3-2 3种曲线对例3-1资料的拟合效果图
以上是3条曲线对结果变量的拟合情况,为节省篇幅,对后9列均人为保留两个小数位数。其中,x,y分别代表原始资料中的原因变量与结果变量,residual1,residual2,residual3分别表示对数曲线、双曲线和幂函数曲线拟合资料的残差,即估计值与真实值之差,ssrs1,ssrs2,ssrs3为对应的残差平方。
以上是衡量曲线拟合情况的几个统计量的值,为节省篇幅,均人为保留两个小数位数。其中,total为变量y的离均差平方和(计算相关指数用),sum resi1,sum resi2,sum resi3依次表示对数曲线、双曲线和幂函数曲线拟合资料的残差平方和,s1,s2,s3为对应的剩余标准差,r square1,r square2,r square3分别代表对应的相关指数。
若想比较3个曲线回归方程拟合此资料的效果是否有统计学差异,可进行曲线拟合优度的比较。三个曲线回归方程所含参数个数相同,其拟合资料的残差平方和分别为1.65,1.61和4.51,选出最大残差平方和与最小残差平方和进行比较。每个曲线回归方程中均有两个待估计参数,故其剩余的自由度均为8-2=6。
F = 4.51/1.61 = 2.80
查方差齐性检验用的F临界值表可知,F0.05/2,6,6 =5.82,因F=2.80<5.82,所以这两个曲线回归方程拟合此资料的效果没有统计学差异。当然,由于参比的是残差平方和最大和最小的两个回归曲线,故这两个曲线回归方程拟合此资料没有统计学差异,也就意味着三个曲线回归方程拟合此资料的效果均无统计学差异。
虽然三个曲线回归方程拟合资料的效果无统计学差异,但在实际应用时,如无专业上的特殊要求,一般仍倾向于选择拟合效果最好的曲线方程。对本资料而言,由于双曲线回归方程对应的残差平方和1.61,相比其他两个回归方程更小一些,故此资料采用双曲线来拟合效果更好。
3.2.3 对例3-3资料的分析
对于例3-3资料,绘制其散点图(见图3-3)可发现,散点的变化趋势近似于指数曲线,故采用指数变换法使其直线化。SAS程序如下,设程序名为nr3 3.sas。
图3-3 产卵数与温度之间的散点图
data nr3_3; /*1*/ input x y@ @ ; y1=log10(y); cards; 21 7 23 11 25 21 27 24 29 66 32 115 35325 ; run; ods html; proc gplot; /*2*/ plot y* x/haxis=21 to 35 by 2 vaxis=0 to 350 by 50; symbol value=dot; run; proc reg; /*3*/ model y1=x; run; ods html close;
【程序说明】 程序第1步是建立数据集,x代表温度,y代表产卵数,对y进行对数变换,以便进行曲线直线化分析。第2步是绘制散点图,以判断散点的趋势。第3步是对变换后的数据进行直线回归分析,本质上是指数曲线的曲线直线化。
【SAS输出结果及其解释】
图3-3是根据原始资料绘制的产卵数和温度之间数据对的散点图,可发现其趋势与指数曲线很相似,故对产卵数进行以10为底的对数变换。
以上是结果变量产卵数(y)经对数变换后产生的新变量(y1)和温度(x)之间进行直线回归分析的结果。可以看出,截距和斜率与0的差别均有统计学意义。因此,直线回归方程为: =-1.67168 +0.11814x。因y1 = log10(y),故
y = 10-1.67168+0.11814 x = 0.0213 × 100.11814 x
3.2.4 对例3-4资料的分析
基于3.2.3小节的分析结果,得到了1条指数回归曲线,现计算该回归曲线的拟合效果,并将这条曲线与原散点绘制在同一直角坐标系中,以直观展现其拟合效果。SAS程序如下,设程序名为nr3 4.sas。
% let n=7; /*1*/ data nr3_4; /*2*/ input x y@ @ ; y1=0.0213* 10** (0.11814* x); residual1=y1-y; ssrs1=residual1** 2; cards; 21 7 23 11 25 21 27 24 29 66 32 115 35325 ; run; symbol1 color=green value=dot; symbol2 color=red value=point line=1 interpol=spline; legend1 across=1 cborder=black label=none value=('实 际 值' '指数曲线') position=(bottom right inside) mode=protect; proc gplot data=nr3_4; /*3*/ plot y* x y1* x/haxis =21 to 35 by 2 vaxis=0 to 350 by 50 overlay legend=legend1; run; proc sql; /*4*/ create table test1 as select css(y) as total,sum(ssrs1) as sum_resi1 from nr3_4; quit; data test; /*5*/ set test1; s1=sqrt(sum_resi1/(&n-2)); r_square1=1-sum_resi1/total; run; ods html; proc print data=nr3_4; /*6*/ run; proc print data=test noobs; /*7*/ run; ods html close;
【程序说明】 第1步是设立宏变量n,其含义为样本个数。第2步是创建数据集,并计算曲线y = 0.0213 × 100.11814x拟合资料的残差residual1及残差平方ssrs1。第3步是将所得曲线与原散点绘制在同一直角坐标系中,以直观展现该曲线的拟合效果。第4步是计算变量y的离均差平方和,以及此曲线拟合资料的残差平方和。第5步是计算该曲线拟合资料的剩余标准差和相关指数。第6步是将第2步所得到的数据集输出到output窗口。第7步输出数据集test至output窗口。
【SAS输出结果及其解释】
图3-4是指数曲线y = 0.0213 100.11814 x对例3-3资料拟合效果的展示图。
图3-4 指数曲线对例3-3资料的拟合效果图
以上是所得指数曲线对结果变量的拟合情况。其中,x,y分别代表原始资料中的原因变量与结果变量,residual1表示指数曲线拟合资料的残差,即估计值与真实值之差,ssrs1为对应的残差平方。
从残差平方和、剩余标准差及相关指数这三个统计量的值可以看出,所得曲线对资料的拟合效果比较理想。
3.2.5 对例3-5资料的分析
对于例3-5资料,绘制其散点图(见图3-5)可发现,散点的变化趋势近似于“S”型,故采用logit变换使其直线化。SAS程序如下,设程序名为nr3 5.sas。
图3-5 麻疹累计发病率与月份之间的散点图
data nr3_5; /*1*/ input x y@ @ ; y1=log((19.34-y)/y); cards; 1 0.279 2 0.920 3 2.078 4 4.393 5 8.144 6 12.174 7 17.306 8 18.435 9 18.728 10 18.965 11 19.230 12 19.328 ; run; ods html; proc gplot; /*2*/ plot y* x/haxis=0 to 12 by 3 vaxis=0 to 20 by 2; symbol value=dot; run; proc reg; /*3*/ model y1=x; run; ods html close;
【程序说明】 程序第1步是建立数据集,x代表月份,y代表麻疹的累计发病率,对y进行logit变换,以便进行曲线直线化分析。第2步是绘制散点图,以判断散点趋势。第3步是对变换后的数据进行直线回归分析,本质上是Logistic曲线的曲线直线化。
【SAS输出结果及其解释】
图3-5是根据原始资料绘制的麻疹累计发病率和月份之间数据对的散点图,可发现其趋势近似“S”型,故对麻疹累计发病率进行logit变换。观察散点图,可发现下渐近线约为y=0,上渐近线约为y=19.34,两条渐近线之间的距离约为19.34,所以取L=0,K=19.34(L和K的含义参见3.3.1小节),对其进行变量变换。
以上是结果变量麻疹累计发病率(y)经logit变换后产生的新变量(y1)和月份(x)之间进行直线回归分析的结果。可以看出,截距和斜率与0的差别均有统计学差异。因此,直线回归方程为: = 5.09934 -0.97293x。因y1 = ln[(19.34-y)/y],故y =19.34/[1+exp(5.09934-0.97293x)] =19.34/[1 +163.91369exp(-0.97293x)]。
3.2.6 对例3-6资料的分析
基于3.2.5小节的分析结果,得到了1条Logistic曲线,现计算该回归曲线的拟合效果,并将这条曲线与原散点绘制在同一直角坐标系中,以直观展现其拟合效果。SAS程序如下,设程序名为nr3 6.sas。
% let n=12; /*1*/ data nr3_6; /*2*/ input x y@ @ ; y1=19.34/(1 +163.91369* exp(- 0.97293* x)); residual1=y1-y; ssrs1=residual1** 2; cards; 1 0.279 2 0.920 3 2.078 4 4.393 5 8.144 6 12.174 7 17.306 8 18.435 9 18.728 10 18.965 11 19.230 12 19.328 ; run; symbol1 color=green value=dot; symbol2 color = red value = point line=1 interpol=spline; legend1 across =1 cborder =black label=none value=('实 际 值' 'logistic曲线 ') position = (bottom right in- side) mode=protect; proc gplot data=nr3_6; /*3*/ plot y* x y1* x/haxis=0 to 12 by 3 vaxis=0 to 20 by 2 overlay leg- end=legend1; run; proc sql; /*4*/ create table test1 as select css(y) as total,sum(ssrs1) as sum_resi1 from nr3_6; quit; data test; /*5*/ set test1; s1=sqrt(sum_resi1/(&n-2)); r_square1=1-sum_resi1/total; run; ods html; proc print data=nr3_6; /*6*/ run; proc print data=test noobs; /*7*/ run; ods html close;
【程序说明】 第1步是设立宏变量n,其含义为样本个数。第2步是创建数据集,并计算曲线y =19.34/[1+163.91369exp(-0.97293x)]拟合资料的残差residual1及残差平方ssrs1。第3步是将所得曲线与原散点绘制在同一直角坐标系中,以直观展现该曲线的拟合效果。第4步是计算变量y的离均差平方和,以及此曲线拟合资料的残差平方和。第5步是计算该曲线拟合资料的剩余标准差和相关指数。第6步是将第2步所得到的数据集输出到output窗口。第7步输出数据集test至output窗口。
【SAS输出结果及其解释】
图3-6是Logistic曲线对例3-5资料拟合效果的展示图。
图3-6 logistic曲线对例3-5资料的拟合效果图
以上是所得Logistic曲线对结果变量的拟合情况。其中,x,y分别代表原始资料中的原因变量与结果变量,residual1表示指数曲线拟合资料的残差,即估计值与真实值之差,ssrs1为对应的残差平方。
从残差平方和、剩余标准差及相关指数这三个统计量的值可以看出,资料拟合比较理想。
3.3 计算原理
若散点图显示某问题中的两定量变量之间呈现某种简单的曲线关系,且这种曲线可以通过对其中一个或两个变量取某种简单的数学变换,对变换后的变量而言,其散布图就呈现直线变化关系了。我们把这样描述曲线关系的做法称之为曲线直线化法。此法仅适合某些简单的曲线拟合问题,如对数曲线、指数曲线、双曲线、幂函数曲线、S型或反S型曲线等。
3.3.1 常见曲线类型及其形状
可直线化的常见曲线类型有:
(1)对数曲线:y=a+blg(x)或y=a+blg(x+k)
(2)指数曲线:y=aebx或y=a10bx
(3)幂函数曲线:y=axb(a>0)
(4)双曲线: = a +(a>0)
(5)Logistic曲线:y = L +
曲线呈“S”型或反“S”型。其中L是下渐近线的纵坐标,K是上下两条平行的渐近线之间的距离,a,b是曲线拟合的待定常数,e是自然对数的底数。
3.3.2 变量变换方法
表3-4列出了常见的曲线直线化的变量转换方法。
表3-4 常见曲线的回归方程及其直线化时的变量转换方法
3.3.3 拟合优度检验
对某一资料采用曲线直线化法进行分析,散点图形态不太明显时,需要对这个资料拟合多种可能的曲线类型,这就牵扯到多个曲线回归方程拟合效果的比较了。当曲线回归方程中所含待定常数(常称为参数)个数相同时,欲比较多个曲线回归方程拟合同一资料的效果是否存在统计学差异,须做曲线拟合优度的相互比较,其方法为F检验,F统计量的计算公式为:
各曲线方程的残差平方和的计算公式为:
式中,y为结果变量的实际观测值,y 为由曲线回归方程算得的结果变量的估计值。
各曲线回归方程对应的剩余自由度为n-k,其中n为试验数据组数,即受试对象的个数;k为曲线拟合中待估计参数的个数。在本章讨论的曲线直线化法分析中,k一般取值为2,因为它包含截距与斜率两个参数。当然,若曲线直线化后拟合的直线中不含截距项,则k=1。
然后,查方差齐性检验用的F临界值表即可得到相应的P值。
当曲线回归方程中所含待定常数(常称为参数)个数不相同时,欲比较直线回归方程与曲线回归方程之间或多个曲线回归方程之间拟合同一资料的效果是否存在统计学差异,须做曲线拟合优度的逐步比较,其方法也是F检验,F统计量的计算公式为:
式中,SS多是参数较多的曲线模型拟合资料的残差平方和,SS少是参数较少的曲线模型拟合资料的残差平方和,残差平方和的计算公式见式(3-2)。v多是参数较多的曲线模型中剩余的自由度,v少是参数较少的曲线模型中剩余的自由度,其计算方法也是样本个数减去曲线模型中参数的个数。
然后,查方差分析用的F临界值表即可得到相应的P值。
曲线拟合优度的逐步比较可检验参数个数增加后,曲线的拟合优度是否有提高。在比较多条参数个数不同的曲线之间的拟合优度是否存在差异时,应按照参数个数从少到多的顺序依次选取两个参数个数相邻的曲线模型进行比较。如果假设检验的结果认为不同参数个数的曲线对资料的拟合效果无统计学差异,则应以参数个数少的曲线为准。反之,如果参数个数多的曲线对资料的拟合效果比参数个数少的曲线好,则可进一步比较参数个数多的这条曲线与参数个数更多(且参数个数与其相邻)的曲线谁的拟合效果好,以此类推。
3.3.4 残差分析
用曲线直线化法拟合的曲线方程有无实用价值,不能仅看直线化后的直线回归方程是否有统计学意义,更要看结果变量的观测值与由曲线回归方程算得的估计值之间的残差以及总的残差平方和或剩余标准差是否较小,同时还应考虑曲线回归方程中参数的个数是否较少。因为参数过多时,曲线回归方程显得很复杂,不便于应用,更重要的是减少了误差项的自由度。有时回归方程中的参数个数增多,残差平方和虽略有减少,但残差的均方却并没有减少甚至可能增大。所以,使用剩余标准差较残差平方和更为妥当。
残差平方和的计算公式同式(3-2)。
剩余标准差的计算公式为:
SS残和S剩都可以用来衡量曲线拟合的好坏,它们越小,说明拟合效果越好。n和k意义同前。
3.3.5 相关指数
在两变量可直线化的曲线回归分析中,可以利用相关指数(correlation index)来衡量曲线拟合的好坏,记为R2。其计算公式为:
如果SS残占SS总的比例很小,说明估计值与实际观察值很接近,曲线对资料的拟合效果较好。即R2越接近于1,曲线对资料的拟合效果越好。SS残的计算公式同式(3-2),SS总为总的离均差平方和,其公式为:
式中,n为受试对象个数。
3.4 分析步骤
可直线化曲线回归分析一般是先通过变量变换的方法将两变量之间的曲线关系转换为直线关系,然后再用最小二乘法求出参数的估计值,得到直线回归方程,最后对变量进行相应的逆变换,得到两个原变量之间的曲线关系。此分析方法的关键是选定合适的曲线模型和寻找适当的变量变换方法。其基本步骤可概括如下:
(1)绘制两变量间的散点图;
(2)根据各点在图中的散布趋势,并结合常见曲线图形的形状和专业知识,选定一种或几种最可能的曲线模型;
(3)根据所选定的曲线方程的特点,进行相应的变量变换,使变换后的两新变量呈直线关系;
(4)建立两新变量之间关系的直线回归方程,并对其假设检验;
(5)将变量还原,得到用原变量表达的曲线回归方程;
(6)若对同一资料拟合了多个可能的模型,需进行曲线拟合优度检验,看它们对资料的拟合效果是否存在统计学差异;
(7)选择最合适的曲线回归方程,并进行残差分析,考察所拟合的曲线回归方程在专业上是否成立,是否值得应用。
参考文献
[1] 胡良平. Windows SAS 6.12 &8.0实用统计分析教程. 北京:军事医学科学出版社,2001:387-399.
[2] 胡良平. 现代统计学与SAS应用. 北京:军事医学科学出版社,1996:235-239.
[3] 金丕焕. 医用统计方法(第2版).上海:复旦大学出版社,2003:317-324.
[4] 张意坚. 未免疫婴儿锡克氏试验阴性率的logistic曲线拟合. 中国卫生统计,1993,10(1):19.
[5] 娄冬华. 曲线直线化与非线性回归. 疾病控制杂志,1999,3(3):224-225.
[6] 郭祖超. 医用数理统计方法(第3版). 北京:人民卫生出版社,1988:573-631.
[7] 杨树勤. 中国医学百科全书(医学统计学). 上海:上海科学技术出版社,1985:170-173.
[8] 胡良平. 面向问题的统计学——(1)科研设计与统计基础. 北京:人民卫生出版社,2012:598-611.