非线性回归分析与SAS智能化实现
上QQ阅读APP看本书,新人免费读10天
设备和账号都新为新人

第6章 二项型指数曲线模型

多项型指数,是由两个或多个以e为底的指数项相加而构成的函数表达式。此函数表达式所描绘出的曲线称为多项型指数曲线。此曲线在药代动力学中应用较多,常用于研究多室模型中血药浓度与时间的关系,是药物动力学的重要研究内容,应用价值非常高。

常见的多项型指数曲线主要有二项型指数曲线和三项型指数曲线。本章主要介绍二项型指数曲线模型,它是两个以e为底的指数项相加而成的。在药代动力学中,单室模型药物血管外给药或二室模型药物静脉注射后,其血药浓度 时间曲线(简称“药 时”曲线)是一条二项型指数曲线。

6.1 问题与数据

【例6-1】 根据有关专业知识,已知某药物为双室模型药物,静脉注射100mg后,测得各时间点的血药浓度结果见表6-1。试拟合该药物的“药 时”曲线。

表6-1 某药物静脉注射后各时间点的血药浓度

6.2 分析与解答

已知药物为双室模型药物,且采用静脉注射,故其“药 时”曲线可由二项型指数曲线模型来拟合。SAS程序如下,设程序名为nr6 1.sas。

    data nr6_1;
      input x y@ @ ;
      cards;
       0.165  65.03
       0.500  28.69
       1.000  10.04
       1.500  4.93
       3.000  2.29
       5.000  1.36
       7.500  0.71
      10.000  0.38
      ;
    run;
    % binomial(nr6_1);

【程序说明】 程序中调用了自编的二项型指数曲线拟合的专用宏程序binomial,此宏程序含有一个宏参数,即数据集的名称。使用时给出宏参数的取值并替换程序左侧部分cards语句之后的数据即可,操作起来简单方便。

【SAS输出结果及其解释

本程序运行后,输出结果如下。

图6-1为例6-1资料的散点图。当无预实验、参考文献或专业知识可辅助确定应选用何种曲线模型拟合资料时,可先绘制散点图,从中判断可能的曲线形态,从而选择合适的曲线模型。当然,就例6-1而言,由专业知识即可确定应选用二项型指数曲线模型来拟合资料。

图6-1 例6-1资料的散点图

该资料拟合二项型指数曲线所得的回归方程为:

y = 94.3687*exp(-2.701*x) +4.7954*exp(-0.2525*x)

残差平方和为:0.0009

这是程序输出的所得最终模型的结果。可见,宏程序拟合该资料效果最好的二项型指数曲线方程(原理参见6.3节)为:

y = 94.3687e-2.7010x +4.7954e-0.2525x

拟合资料的残差平方和为0.0009,拟合效果非常好。

为展示该宏程序对资料的拟合优劣程度,下面给出不同散点组合时采用残数法(参阅6.3.3小节)所拟合的曲线方程及其拟合效果,见表6-2。表中“散点组合”列中的两位数字含义为两个计算阶段所选用的散点个数,如散点组合“34”,表示第一阶段由后至前选择3个散点,第二阶段在剩余散点中由后至前选择4个散点。

表6-2 不同散点组合时的残数法与非线性最小二乘法拟合的回归方程

可见,非线性最小二乘法所得的曲线方程对资料的拟合效果是最好的。图6-2可形象地展示非线性最小二乘法所得的曲线方程对资料的拟合效果。

图6-2 非线性最小二乘法所得二项型指数曲线方程对例6-1资料的拟合效果图

6.3 计算原理

6.3.1 曲线方程

二项型指数曲线模型的方程如下:

上述方程中,α>0,β>0,另外有一个潜在的条件,即αβ。否则,上式可转换为:

此为一般的指数曲线方程。所以,αβ之间必然存在着大小之分。

6.3.2 常见曲线形态

二项型指数曲线常见形态见图6-3,其中图6-3(a)和图6-3(c)较为常用。

图6-3 二项型指数曲线的常见形态

6.3.3 计算方法

现假设α>β,则当x足够大时,Ae-αx将趋向于0。式(6-1)可简化为:

对式(6-3)进行以e为底的对数变换,则有:

对结果变量y进行自然对数变换,得到新变量lny,然后以其为纵轴变量,以时间变量x为横轴变量,在平面直角坐标系内绘制散点图,记为lny-x图。若在此散点图中,由后向前选取连续的几个(3个或3个以上)基本呈直线趋势的散点,即可拟合一条回归直线。此回归直线的斜率为-β,截距为lnB。然后,可求出β值(等于斜率的相反数)和B值(等于exp(截距))。

β值和B值代入式(6-3)可得到各散点上结果变量的预测值。同样,可计算其他散点(即建立回归直线时所用散点之前的那些点,称为外推点)上结果变量的预测值,称为外推值。

对式(6-1)进行变换,可得:

以外推点的坐标值代入式(6-5)或式(6-6)。其中,y为外推点的实际值,Be-βx 为外推点的外推值,二者的差值是残数值,记为y。对式(6-5)和式(6-6)进行自然对数变换,可得:

两式中均存在一个默认假设,式(6-7)默认y - Be-βx > 0,式(6-8)默认Be-βx - y > 0。

作lny -x图,自后向前选取连续的几个(3个或3个以上)基本呈直线趋势的散点,若多数外推点(>总数据点数的1/2)实际值与外推值之差大于零,则可舍弃另一部分外推点,仅以二者之差大于零的这些外推点按照式(6-7)进行分析;反之,若多数外推点(>总数据点数的1/2)的实际值与外推值之差小于或等于零,则可仅以二者之差小于零的这些外推点按照式(6-8)进行分析。

拟合回归直线后,所得残数线的截距和斜率应为-α和lnA(对应式(6-7))或-α和ln(-A)(对应式(6-8))。然后,可求出 α 值(等于斜率的相反数)和 A 值(等于exp(截距)或 -exp(截距))。

至此,可计算出全部参数α,β,A,B的值。

6.4 分析步骤

6.4.1 分析策略

以上求解模型的方法称为残数法,它是把二项型指数曲线模型分解成两个指数成分,然后对这两个指数成分通过分段曲线直线化的方式得到相应指数成分的参数估计值,也就可以得到拟合的二项型指数曲线模型了。

但是,此法在计算过程中需要进行曲线直线化,而后再采用最小二乘法使变量变换后所得新变量的离均差平方和最小,这并不一定能使原结果变量的离均差平方和最小,所以其模型的拟合精度仍存在一定的提高空间。

以残数法计算所得的参数估计值为初始值,再借助于SAS软件STAT模块中的NLIN过程,采用非线性最小二乘法继续进行拟合,可得到拟合效果更好的曲线模型。当然,NLIN过程求解时使用的迭代算法,决定了其所得的曲线模型可能仅仅是局部最优的。为在一定程度上破解局部最优的困境,尽可能地得到全局最优解,并解决残数法求解过程中的散点选择与组合问题,可借助SAS软件的强大编程功能,对所有可能的散点组合均予以考虑,从而得到不同散点组合下的局部最优解,然后从这些局部最优解中选择出最优的,这一策略简称“优中选优”策略。

上述策略的框图如图64所示。

图6-4 二项型指数曲线模型分析策略框图

6.4.2 实施步骤

首先,以专业知识为依据,当判定某资料可能属于二项型指数曲线范畴时,再采用以下步骤实施二项型指数曲线拟合。

第一步,针对实际数据,建立相应的SAS数据集。以y表示结果变量,以x表示原因变量。

第二步,借助SAS宏,实现散点的全面组合。同时,在此宏内,需要计算不同散点组合情形下模型中参数的初始值和NLIN过程得到的最优解。

第三步,将各种局部最优解集合在一起,从中选择残差平方和最小或相关指数最大的那个模型,作为最终的模型。

具体编程时,主要的难点在第二步,具体问题和解决方法如下:

一是散点的全面组合问题,可通过在SAS宏内使用双do循环的方式解决。

二是每种散点组合情形下得到的参数初始值自动赋给NLIN过程的问题,可在拟合回归直线后,借助call symput语句,将所得截距和斜率经相应计算后,赋给相应的宏变量,供NLIN过程调用。

三是计算α值和A值时,式(6-7)和式(6-8)的自动判断与选择问题。借助数据步,用sign函数来判断残数值(设定为实际值减外推值)与0的关系,并对sign函数进行累积求和。若累积和大于0,则残数值大于0的外推点要多于小于0的外推点,此时应删去残数值小于或等于0的外推点,然后采用式(6-7)继续进行分析;反之,则残数值小于0的外推点多于残数值大于0的外推点,此时应删去残数值大于或等于0的外推点,然后采用式(6-8)继续进行分析。

参考文献

[1] 薛仲三. 医学统计方法和原理(内部资料). 北京:军事医学科学,1984:276-287.

[2] 梁文权. 生物药剂学与药物动力学(第2版). 北京:人民卫生出版社,2006:164-240.

[3] SAS Institute Inc. SAS/STAT 9.2 User's Guide. Cary,NC:SAS Institute Inc.,2008:4261-4336.