谷恒明，胡良平,2*

(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心，北京 100850；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專(zhuān)業(yè)委員會(huì)，北京 100029*通信作者：胡良平，E-mail：lphu812@sina.com)

復(fù)雜曲線回歸分析及其應(yīng)用

谷恒明1，胡良平1,2*

(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心，北京 100850；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專(zhuān)業(yè)委員會(huì)，北京 100029*通信作者：胡良平，E-mail：lphu812@sina.com)

本文目的是介紹曲線方程的兩種擬合方法(即直線化法與直接擬合法)之間的區(qū)別。通過(guò)對(duì)一個(gè)實(shí)例中散點(diǎn)圖為“S型曲線”的資料采用直線化法與直接擬合法擬合Logistic曲線回歸方程，對(duì)所得到的擬合結(jié)果進(jìn)行比較，發(fā)現(xiàn)基于粗估值再采用SAS中的NLIN過(guò)程直接擬合Logistic曲線回歸方程，可以得到更精確的擬合效果。

非線性回歸分析；Logistic曲線回歸方程；迭代計(jì)算；相關(guān)指數(shù)

*Correspondingauthor:HuLiangping,E-mail:lphu812@sina.com)

1 概 述

1.1 問(wèn)題的提出

本刊本期的前一篇文章《簡(jiǎn)單曲線回歸分析及其應(yīng)用》介紹了用曲線直線化的方法實(shí)現(xiàn)幾個(gè)常用初等函數(shù)曲線的曲線回歸或曲線擬合的方法。然而，人們可能會(huì)提出下面兩個(gè)問(wèn)題。問(wèn)題一：曲線直線化方法與直接進(jìn)行曲線擬合的效果一樣嗎？問(wèn)題二：是否所有的函數(shù)曲線都可以通過(guò)曲線直線化構(gòu)建曲線回歸方程？

對(duì)第二個(gè)問(wèn)題的回答也是否定的。其理由如下：事實(shí)上，能夠通過(guò)直線化來(lái)實(shí)現(xiàn)曲線擬合的曲線類(lèi)型是很少的，絕大多數(shù)函數(shù)曲線是無(wú)法通過(guò)簡(jiǎn)單直線化法來(lái)實(shí)現(xiàn)曲線擬合的。例如下面的兩個(gè)函數(shù)曲線就無(wú)法通過(guò)直線化法進(jìn)行曲線擬合。

二項(xiàng)型指數(shù)函數(shù)是由兩個(gè)指數(shù)函數(shù)項(xiàng)相加而構(gòu)成的函數(shù)表達(dá)式。此函數(shù)表達(dá)式所描繪出的曲線稱(chēng)為二項(xiàng)型指數(shù)曲線，見(jiàn)式(1)。

y=A×eαt+B×e-βt

(1)

又例如：三項(xiàng)型指數(shù)函數(shù)是由三個(gè)指數(shù)函數(shù)項(xiàng)相加而構(gòu)成的函數(shù)表達(dá)式。此函數(shù)表達(dá)式所描繪出的曲線稱(chēng)為三項(xiàng)型指數(shù)曲線。此曲線常用于研究三室模型藥物靜脈注射或二室模型藥物血管外給藥后血藥濃度與時(shí)間的關(guān)系，見(jiàn)式(2)。

y=Ne-kx+Le-αx+Me-βx

(2)

1.2 解決方案

SAS軟件提供了NLIN和NLMIXED過(guò)程，可用于非線性曲線的擬合。但在擬合時(shí)需要讀者給出未知參數(shù)的初始值。若初始值偏離真實(shí)值較遠(yuǎn)，則可能擬合過(guò)程無(wú)法收斂或擬合曲線并非最優(yōu)。所以，較為合適的方法是通過(guò)曲線直線化求出參數(shù)的粗估值，然后將其作為初始值引入NLIN或NLMIXED過(guò)程進(jìn)行迭代計(jì)算，以尋找擬合效果更好的曲線回歸方程。然而，對(duì)于類(lèi)似上面的式(1)和(2)那樣無(wú)法實(shí)施直線化的函數(shù)方程式，估計(jì)參數(shù)的粗估值可能就帶有很大的盲目性了。具體如何實(shí)現(xiàn)，請(qǐng)參見(jiàn)相關(guān)參考文獻(xiàn)[1-2]。

由于直接進(jìn)行曲線擬合問(wèn)題難度較大，再考慮到文章篇幅有限，故本文僅介紹一個(gè)最常用的曲線擬合的實(shí)例。

2 實(shí)例解析

2.1 問(wèn)題與數(shù)據(jù)

【例1】某研究者欲分析某縣瘧疾發(fā)病的季節(jié)性特點(diǎn)，觀測(cè)了某縣1961年-1996年瘧疾的月累計(jì)發(fā)病率(1/10萬(wàn))，結(jié)果見(jiàn)表1。試對(duì)該資料進(jìn)行曲線擬合。

表1 某縣瘧疾月累計(jì)發(fā)病率

2.2 如何用SAS實(shí)現(xiàn)計(jì)算

【分析與解答】圖1散點(diǎn)圖顯示，資料的趨勢(shì)近似S型曲線，故可采用Logistic曲線方程擬合此資料。曲線上限為130，下限為0，上下漸近線之間的垂直間隔為130，SAS程序如下：

%let ul=130; /*1*/

dataabc; /*2*/

do x=1 to 12;

input y@@;

z=log((&ul-y)/y);

output;

end;

cards;

0.555 1.295 2.751 11.116

24.879 43.476 69.297 97.037

114.631 121.645 124.412 125.619

;

run;

ods html;

proc gplot; /*3*/

plot y*x/haxis=0 to 12 vaxis=0 to &ul by 10;

symbol value=dot;

run;

proc reg data=abc outest=est(keep=Intercept x);

model z=x; /*4*/

run;quit;

data set1; /*5*/

set est;

call symput('a1',exp(Intercept));

call symput('b1',x);

run;

proc nlin data=abc; /*6*/

parms K=&ul a=&a1 b=&b1;

model y=K/(1+a*exp(b*x));

output out=set2 p=yp r=resid;

run;

proc sql; /*7*/

create table set3 (sum num,css

num);

insert into set3 select

1.人力資源配置優(yōu)化。成立區(qū)域檢修公司，通過(guò)統(tǒng)一調(diào)度檢修人員，可以解決區(qū)域公司內(nèi)部檢修人員工作量不均衡，同時(shí)解決單個(gè)企業(yè)檢修維護(hù)部組建人員和技術(shù)力量不足的問(wèn)題。

sum(resid**2),css(y) from set2;

quit;

data set4; /*8*/

set set3;

r2=1-sum/css;

run;

proc print data=set2 round; /*9*/

run;

proc print data=set4; /*10*/

run;

ods html close;

【程序說(shuō)明】程序中第1步是設(shè)立宏變量ul，表示曲線的上、下漸近線的垂直間距，其值為130，以便于后面調(diào)用。第2步是創(chuàng)建數(shù)據(jù)集，x表示月份，y表示瘧疾月累計(jì)發(fā)病率，z是對(duì)y進(jìn)行l(wèi)ogit變換所得的變量，以便后面進(jìn)行曲線直線化分析。第3步是繪制散點(diǎn)圖，以發(fā)現(xiàn)散點(diǎn)分布趨勢(shì)，本資料散點(diǎn)趨勢(shì)近似S型曲線。第4步對(duì)變量z和x進(jìn)行直線回歸分析，將所得直線的截距和斜率保存在數(shù)據(jù)集est中。第5步對(duì)數(shù)據(jù)集est進(jìn)行操作，將以自然對(duì)數(shù)為底、以截距為指數(shù)的數(shù)值賦給宏變量a1，將斜率賦給宏變量b1。當(dāng)然，也可直接將截距賦給宏變量a1，將斜率賦給宏變量b1，但第6步中model語(yǔ)句給出的曲線方程需改為“y=K/(1+exp(a+b*x))”。第6步是調(diào)用NLIN過(guò)程來(lái)擬合Logistic曲線，宏變量ul、a1、b1的值分別賦給變量K、a、b作為初值，并將所得曲線對(duì)各觀測(cè)擬合的預(yù)測(cè)值及殘差輸出到數(shù)據(jù)集set2中，預(yù)測(cè)值以變量yp表示，殘差以resid表示。第7步將數(shù)據(jù)集set2的內(nèi)容輸出到output窗口。第7步是調(diào)用sql過(guò)程，建立一個(gè)表set3，包括sum和css兩個(gè)變量，然后將數(shù)據(jù)集set2中所有觀測(cè)resid的平方和賦給set3中的sum變量，將set2中所有觀測(cè)y變量的離均差平方和賦給css。第8步是對(duì)數(shù)據(jù)集set3進(jìn)行操作，計(jì)算新的變量r2，即相關(guān)指數(shù)。第9步是將數(shù)據(jù)集set2的內(nèi)容輸出到output窗口，round選項(xiàng)用來(lái)指定數(shù)據(jù)最多保留兩位小數(shù)。第10步是將數(shù)據(jù)集set4的內(nèi)容輸出到output窗口。

2.3 SAS輸出結(jié)果及其解釋

【主要輸出結(jié)果及解釋】

圖1 y隨x變化的散點(diǎn)圖

圖1是SAS程序繪制的散點(diǎn)圖?？梢钥闯觯Ⅻc(diǎn)近似S型曲線趨勢(shì)，上限為130，下限為0。

AnalysisofVarianceDependentVariable：zSourceDFSumofSquaresMeanSquareFValuePr>FModel1102．44535102．44535658．11<．0001Error101．556660．15567CorrectedTotal11104．00201

RootMSE0．39455R-Square0．9850DependentMean0．50213AdjR-Sq0．9835CoeffVar78．57356

ParameterEstimatesVariableDFParameterEstimateStandardErrortValuePr>|t|Intercept16．003770．2428324．72<．0001x1-0．846400．03299-25．65 <．0001

這是曲線直線化分析的結(jié)果，即變量z與x之間直線回歸分析的結(jié)果。截距為6.00377，斜率為-0.8464，它們與0之間的差異均有統(tǒng)計(jì)學(xué)意義，檢驗(yàn)統(tǒng)計(jì)量t值分別為24.72、-25.65，P均<0.0001。直線回歸方程為z=6.00377-0.8464x，轉(zhuǎn)換成y與x的關(guān)系，即為：y=130/(1+e6.00377-0.8464x)。此時(shí)，殘差平方和為194.02219，相關(guān)指數(shù)R2=0.99740。

The NLIN Procedure Dependent Variable y Method: Gauss-Newton

NOTE: Convergence criterion met

Note： An intercept was not specified for this model

這是NLIN過(guò)程擬合曲線模型的有關(guān)信息。迭代方法為Gauss-Newton法，以K=130，a=405，b=-0.8464為初始值，迭代了5次，得到了收斂的結(jié)果。

SourceDFSumofSquaresMeanSquareFValueApproxPr>FModel376037．425345．811770．1<．0001Error919．38072．1534UncorrectedTotal1276056．8

ParameterEstimateApproxStdErrorApproximate95%ConfidenceLimitsK127．81．0737125．4130．2a394．665．1746247．1542．0b-0．88790．0263-0．9474-0．8284

ApproximateCorrelationMatrixKabK1．0000000-0．50866920．6148707a-0．50866921．0000000-0．9798514b0．6148707-0．97985141．0000000

這是NLIN過(guò)程擬合模型的結(jié)果。對(duì)模型的檢驗(yàn)結(jié)果為F=11770.1，P<0.0001，說(shuō)明所擬合的模型是有統(tǒng)計(jì)學(xué)意義的。模型的三個(gè)參數(shù)取值分別為：K=127.8，a=394.6，b=-0.8879。所以，本資料所擬合的Logistic曲線回歸方程為：

Obsxyzypresid110．565．450．78-0．23221．304．601．88-0．59332．753．834．49-1．734411．122．3710．380．745524．881．4422．602．286643．480．6943．83-0．357769．30-0．1371．46-2．178897．04-1．0896．500．5399114．63-2．01112．761．871010121．65-2．68121．160．491111124．41-3．10124．99-0．581212140．62-3．36126．64-1．02

Obssumcssr2119．380730827．910．99937

這是曲線擬合效果的結(jié)果。首先給出了曲線回歸方程對(duì)各觀測(cè)的擬合結(jié)果，yp列為曲線對(duì)各觀測(cè)的預(yù)測(cè)值，resid為預(yù)測(cè)值與觀測(cè)值之間的殘差。然后給出了對(duì)擬合整體效果的評(píng)價(jià)情況，殘差平方和為19.3807，相關(guān)指數(shù)R2=0.99937。

由此可知，在粗估值基礎(chǔ)上直接進(jìn)行Logistic曲線擬合，可獲得更好的擬合效果。

[1] 胡良平, 高輝. 非線性回歸分析與SAS智能化實(shí)現(xiàn)[M].北京: 電子工業(yè)出版社, 2013: 68-119.

[2] 胡良平. 科研設(shè)計(jì)與統(tǒng)計(jì)分析[M].北京: 軍事醫(yī)學(xué)科學(xué)出版社, 2012: 399-426.

Complexcurveregressionanalysisanditsapplication

GuHengming1,HuLiangping1,2*

(1.ConsultingCenterofBiomedicalStatistics,AcademyofMilitaryMedicalSciences,Beijing100850,China;2.SpecialtyCommitteeofClinicalScientificResearchStatisticsofWorldFederationofChineseMedicineSocieties,Beijing100029,China

The purpose of this paper is to introduce the difference between the two fitting methods (linearization and direct fitting) of curve equation. By fitting the Logistic curve regression equation with the linearized method and the direct fitting method for an example in which the scatter plot is "S-shaped curve". Compared the fitting results producing from the two methods mentioned above, it was found that the more accurate fitting result can be gotten by means of providing the crude estimators of the parameters in the regression equation and then directly fitting Logistic curve regression equation through the NLIN procedure in SAS/STAT.

Nonlinear regression analysis; Logistic curve regression equation; Iteration calculation; Correlation index

國(guó)家高技術(shù)研究發(fā)展計(jì)劃課題資助(2015AA020102)

R195.1

A

10.11886/j.issn.1007-3256.2017.06.004

2017-12-03)

陳 霞)