加性与广义加性模型回归分析

2018-09-11胡良平

四川精神卫生 2018年4期

关键词：加性计分因变量

胡良平

(1.军事科学院研究生院，北京 100850；2.世界中医药学会联合会临床科研统计学专业委员会，北京 100029

1 概述[1-2]

1.1 加性模型

将多重线性回归模型进行推广，使其表达式成为下面的式(1)形式：

Y=s0+s1(X1)+s2(X2)+…+sp(Xp)+ε

(1)

在式(1)中，sj(Xj),j=1,2,…,p,是P个“光滑函数”；误差“ε”满足如下条件：它的期望为0[E(ε)=0]、方差为σ2[Var(ε)=σ2]。为了使式(1)成为可以估计的，要求光滑函数si(Xi)必须满足如下的标准化条件：即期望为0，E[sj(Xj)]=0。式(1)中的P个光滑函数不以参数形式呈现，而以非参数形式呈现。

1.2 广义加性模型

在加性模型式(1)中，假定因变量y服从正态分布。然而，在很多场合下，因变量不服从正态分布，而可能服从其他某种分布。现假定式(1)中的因变量y具有下面的指数族分布密度，见式(2)：

(2)

在式(2)中，θ被称为“自然参数”，φ被称为“尺度参数”；以因变量y的平均值μ为自变量构造出的函数g(μ)被称为“连接函数”，它与协变量“X1、X2、…、Xp”之间建立了关系。下面的数量定义了“加性分量(或成分)”，见式(3)：

(3)

式(3)中，S1()、S2()、…、Sp()都是“光滑函数”。μ与η之间的关系由下式来定义：

g(μ)=η

(4)

最常用的连接函数为“典型连接函数”，即η=θ。

1.3 加性与广义加性模型回归分析应用的场合

广义线性模型强调对模型中参数的估计和推断，而广义加性模型则聚焦于如何用非参数法探测数据。换言之，广义加性模型更适合于探查数据并可视化因变量与自变量之间的关系。

1.4 加性与广义加性模型回归分析的计算原理

1.4.1 加性模型回归分析的计算原理

基于加性模型式(1)，可以构造如下的残差，见式(5)：

(5)

式(5)被称为“第k项光滑参数与因变量y”之间的“残差”，即Rk≈sk(Xk)。严格地说，应该有下式成立，见式(6)：

E(Rk|Xk)=sk(Xk)

(6)

1.4.2 后退拟合算法

1.4.2.1 未加权的后退拟合算法

未加权的后退拟合算法步骤如下：

第1步：初始化。

第2步：迭代。令m=m+1；让j从1到p循环，循环体内的计算公式为：

第3步：终止。

直到上式的计算结果不再下降，或满足收敛的临界值，就停止迭代。

值得注意的是：对于除正态分布之外的其他分布而言，具有权重的数值不稳定可能会引起收敛问题。即使当算法收敛时，各个个体函数并不需要彼此完全不同，即便对于同一个拟合曲面来说，由于协变量之间的依赖性会导致多于一个表达式出现。

1.4.2.2 加权的后退拟合算法

除了要对光滑器进行加权之外，加权的后退拟合算法具有与未加权的后退拟合算法相同的形式。在SAS的GAM过程步中，具体是在采用“局部计分过程”中且资料为非正态分布时，使用了“权重”。

GAM过程使用下面的“条件”作为后退拟合算法的“收敛临界值”：

此处，ε=10-8是缺省的界值。用户可以通过修改模型语句中的选项“EPSILON=”来改变此界值。

1.4.3 广义加性模型回归分析的计算原理

广义加性模型回归分析的计算原理比前面介绍的加性模型回归分析的计算原理稍复杂，其中，最关键的内容为“局部计分算法”。该算法的重要内容取决于与每个特定分布对应的“连接函数”。它们之间的关系见表1。

表1 局部计分算法涉及到的重要内容

注：表1摘录自SAS软件GAM过程帮助信息

由表1可知：一旦分布被指定，相应的“那些量”也就被定义了。于是，可按下面的步骤实施“局部计分算法”。

1.4.4 一般的局部计分算法

第1步：初始化。

第2步：迭代。令m=m+1；从前一次迭代中获得各变量的相应数值，这些变量分别是：预测量η、均值μ、权重w和校正后的因变量z：

第3步：终止。直到达到收敛临界值或离差不再减少时，停止迭代。这里所说的“离差”，实际上是广义线性模型中“RSS”的一个扩展或推广。

GAM过程使用下面的“条件”作为局部计分算法的“收敛临界值”：

此处，ε=10-8是缺省的界值。用户可以通过修改模型语句中的选项“EPSILON=”来改变此界值。

算法小结：广义加性模型的估计过程由两个循环构成。在局部计分算法(外循环)每一步内部，使用加权后退拟合算法(内循环)，直到收敛或RSS不再减少。然后，基于来自这个加权后退拟合算法得到的估计量，计算出一组新的权重，开始计分算法的下轮迭代。当达到收敛临界值或估计量的离差停止下降时，计分算法也就停止了。

2 基于加性模型回归分析解决实际问题[2]

2.1 问题与数据结构

【例1】下面是一个假设的来自化学试验的例子：每次试验，研究者将某种催化剂加入到某种化学溶液中，从而合成一种新化合物。其数据是测量溶液的温度(temperature)、加入的催化剂量(catalyst)和化学反应量(yield)的结果。试验数据的结构很简单，两个计量原因变量和一个计量结果变量及其取值，前6次试验数据见表2。

表2 不同“溶液温度”和“催化剂量”条件下化学反应量的测定结果

注：表中仅列出了少量数据，详细试验数据见下面的SAS程序，此处从略

【对数据结构的分析】对于每次试验而言，可以观测到3个计量的数据，即溶液的温度(temperature)、加入的催化剂量(catalyst)和化学反应的量(yield)。

【统计分析方法的选择】若希望考察化学反应的量(yield)是如何随溶液的温度(temperature)和催化剂量(catalyst)变化而变化的依赖关系，可选择多重线性回归分析。因本例中的因变量为计量变量，故可以考虑选用“加性模型回归分析”。若因变量是定性变量或计数变量，就可能需要选用“广义加性模型回归分析”。因篇幅所限，本文只介绍如何用SAS实现“加性模型回归分析”。

2.2 基于常规方法构建多重线性回归模型[3-4]

2.2.1 创建SAS数据集

创建一个名为“ExperimentA”的临时SAS数据集所需的数据步程序：

data ExperimentA;

format Temperature f4.0 Catalyst f6.3 Yield f8.3;

input Temperature Catalyst Yield @@;

x1=temperature;x2=Catalyst;y=Yield;

datalines;

80 0.005 6.039 80 0.010 4.719 80 0.015 6.301

80 0.020 4.558 80 0.025 5.917 80 0.030 4.365

80 0.035 6.540 80 0.040 5.063 80 0.045 4.668

80 0.050 7.641 80 0.055 6.736 80 0.060 7.255

80 0.065 5.515 80 0.070 5.260 80 0.075 4.813

80 0.080 4.465 90 0.005 4.540 90 0.010 3.553

90 0.015 5.611 90 0.020 4.586 90 0.025 6.503

90 0.030 4.671 90 0.035 4.919 90 0.040 6.536

90 0.045 4.799 90 0.050 6.002 90 0.055 6.988

90 0.060 6.206 90 0.065 5.193 90 0.070 5.783

90 0.075 6.482 90 0.080 5.222 100 0.005 5.042

100 0.010 5.551 100 0.015 4.804 100 0.020 5.313

100 0.025 4.957 100 0.030 6.177 100 0.035 5.433

100 0.040 6.139 100 0.045 6.217 100 0.050 6.498

100 0.055 7.037 100 0.060 5.589 100 0.065 5.593

100 0.070 7.438 100 0.075 4.794 100 0.080 3.692

110 0.005 6.005 110 0.010 5.493 110 0.015 5.107

110 0.020 5.511 110 0.025 5.692 110 0.030 5.969

110 0.035 6.244 110 0.040 7.364 110 0.045 6.412

110 0.050 6.928 110 0.055 6.814 110 0.060 8.071

110 0.065 6.038 110 0.070 6.295 110 0.075 4.308

110 0.080 7.020 120 0.005 5.409 120 0.010 7.009

120 0.015 6.160 120 0.020 7.408 120 0.025 7.123

120 0.030 7.009 120 0.035 7.708 120 0.040 5.278

120 0.045 8.111 120 0.050 8.547 120 0.055 8.279

120 0.060 8.736 120 0.065 6.988 120 0.070 6.283

120 0.075 7.367 120 0.080 6.579 130 0.005 7.629

130 0.010 7.171 130 0.015 5.997 130 0.020 6.587

130 0.025 7.335 130 0.030 7.209 130 0.035 8.259

130 0.040 6.530 130 0.045 8.400 130 0.050 7.218

130 0.055 9.167 130 0.060 9.082 130 0.065 7.680

130 0.070 7.139 130 0.075 7.275 130 0.080 7.544

140 0.005 4.860 140 0.010 5.932 140 0.015 3.685

140 0.020 5.581 140 0.025 4.935 140 0.030 5.197

140 0.035 5.559 140 0.040 4.836 140 0.045 5.795

140 0.050 5.524 140 0.055 7.736 140 0.060 5.628

140 0.065 6.644 140 0.070 3.785 140 0.075 4.853

140 0.080 6.006

;

run;

【SAS程序说明】数据中每行上有3次试验的结果，每次试验结果都有3个数据，即温度数值(temperature)、催化剂量(catalyst)与产量(yield)。

创建一个名为“ExperimentB”的临时SAS数据集的SAS数据步程序：

data ExperimentB;

set ExperimentA;

x3=x1*x1;x4=x2*x2;x5=x1*x2;x6=x3*x1;

x7=x4*x2;x8=x3*x2;x9=x4*x1;

run;

2.2.2 基于常规方法构建多重线性回归模型

利用下面的两个SAS过程步程序可以创建两个二重线性回归模型：

proc reg data=ExperimentA;

model y=x1 x2/r;

run;

记以上SAS程序创建的二重线性回归模型为模型(1)。

proc reg data=ExperimentA;

model y=x1 x2/noint r;

run;

记以上SAS程序创建的二重线性回归模型为模型(2)。

经比较，模型(1)优于模型(2)。具体方法详见下文，此处从略。

引入自变量的“二次项”，利用下面的两个SAS过程步程序可以创建两个多重线性回归模型：

proc reg data=ExperimentB;

model y=x1-x5/selection=backward sle=0.05 r;

run;

记以上SAS程序创建的多重线性回归模型为模型(3)。

proc reg data=ExperimentB;

model y=x1-x5/noint selection=backward sle=0.05 r;

run;

记以上SAS程序创建的多重线性回归模型为模型(4)。

经比较，模型(4)优于模型(3)。具体方法详见下文，此处从略。

引入自变量的“三次项”，利用下面的两个SAS过程步程序可以创建两个多重线性回归模型：

proc reg data=ExperimentB;

model y=x1-x9/selection=backward sle=0.05 r;

run;

记以上SAS程序创建的多重线性回归模型为模型(5)。

proc reg data=ExperimentB;

model y=x1-x9/noint selection=backward sle=0.05 r;

run;

记以上SAS程序创建的多重线性回归模型为模型(6)。

经比较，模型(5)优于模型(6)。具体方法详见下文，此处从略。

将模型(4)与模型(1)比较，得出模型(4)优于模型(1)。最后，需要将模型(5)与模型(4)作比较，具体方法如下：

模型(4)的有关信息为：SSε=128.48055(模型误差的离均差平方和)、dfε=108(误差的自由度)；

模型(5)的有关信息为：SSε=78.07028(模型误差的离均差平方和)、dfε=106(误差的自由度)。

利用下面的F检验对上述回归模型(5)与模型(4)进行拟合优度比较：

对应的F临界值F((2,106)(0.01))<4.82，因F=34.222>4.82，说明P<0.01，故需要选择参数多的回归模型(5)。模型(5)的输出结果如下：

方差分析源自由度平方和均方F值Pr > F模型588.6661117.7332224.08<.0001误差10678.070280.73651校正合计111166.73639

均方根误差0.85820R方0.5318因变量均值6.14077调整 R方0.5097变异系数13.97550

参数估计值变量自由度参数估计值标准误差t 值Pr > |t|Intercept1122.1032018.516226.59<.0001x11-3.459610.52226-6.62<.0001x310.033220.004826.89<.0001x411437.83096246.932355.82<.0001x61-0.000103550.00001460-7.09<.0001x71-17556.000003130.36119-5.61<.0001

根据最后的“参数估计值”，请读者写出相应的“五重线性回归模型”的表达式，此处从略。