半变系数再生散度混合效应模型的影响分析 Influence Analysis of Semivarying Coefficient Reproductive Dispersion Mixed Models

设为首页加入收藏期刊导航网站地图

期刊菜单

文章导航

Statistical and Application 统计学与应用 , 2012, 1, 37-43

http://dx.doi.org/10.12677/sa.2012.12008 Published Online December 2012 (http://www.abtbus.com/journal/sa.html)

Influence Analysis of Semivarying Coefficient

Reproductive Dispersion Mixed Models

Rong Jiang, Xiaohan Yang, Weimin Qian

Department of Mathematics, Tongji University, Shanghai

Email: jrtrying@126.com, xiaohyang@tongji.edu.cn, wmqian2003@yahoo.com.cn

Received: Nov. 19 th , 2012; revised: Nov. 26th, 2012; accepted: Dec. 13th, 2012

Abstract: This paper proposes several case-deletion as well aslocal influence measures for assessing the in-

fluence of an observation for semivarying coefficient reproductive dispersion mixed models. The essential

idea is to treat the latent random effects in the m odel as missing data and estimate unknown parameters by

acceleration of Monte Carlo EM algorithm. On the basis of the Q-func tion which is associated with the con-

ditional expectation of the complete-data log-likelihood, we generate generalized Cook Distance. Moreover,

three different perturbation schemes are discussed. Finally ，one real illustrative example is presented to prove

the methodology.

Keywords: Semivarying Coefficient Reproductive Dispersion Mixed Models; Local Influences; Penalized

Spline; Cook Distance; Accelerati on of Monte Carlo EM Algorithm

半变系数再生散度混合效应模型的影响分析

姜荣，杨筱菡，钱伟民

同济大学应用数学系，上海

Email: jrtrying@126.com, xiaohyang@tongji.edu.cn, wmqian2003@yahoo.com.cn

收稿日期： 2012 年11月19日；修回日期： 2012 年11月26日；录用日期： 2012 年12月13日

摘要：本文把随机效应看作缺失数据并利用P-样条拟合非参数部分，应用Monte Carlo EM加速算

法得到半变系数再生散度混合效应模型的未知参数的估计，同时利用 Q 函数，得到了模型的广义Cook

距离。此外，本文还研究了三种不同扰动情形的局部影响分析，得到了相应的影响矩阵。最后，通过一

个实际例子验证了所提出的诊断统计量的有效性。

关键词：半变系数再生散度混合效应模型；局部影响；P-样条；广义Cook距离；Monte Carlo EM加

速算法

1. 引言

统计诊断从 20世纪 70 年代中期受到统计学家的

广泛关注，经过近 40 年的发展，异常点识别、残差

分析、影响分析和数据变换等内容现已成为统计诊断

的主要课题。特别地，基于数据删除模型和局部影响

的诊断分析方法现已成为统计诊断的通用方法，它们

可广泛地应用于各种统计模型的影响分析。例如，线

性模型 (Cook and Weisberg[1])，非线性回归模型(Seber

and W ild [2] )，半参数非线性模型(姜荣，邵明江，钱伟

民 [3])，线性混合效应模型(Beckman et al.[4])，半参数广

义线性混合效应模型 (张浩，朱仲义 [5] ) 。Jorgensen[6]

首次提出了再生散度模型 (RDM)，并指出广义线性模

型的理论可以推广到以 RDM 为随机误差的模型。唐

年胜和韦博成 [7] 研究了非线性再生散度随机效应

半变系数再生散度混合效应模型的影响分析

模型，讨论了该模型的几何结构、渐近性质和统计诊

断等问题。

变系数模型在生物医学、公共卫生、经济、农业、

制造业、道路安全等众多领域的数据分析中有广泛的

应用。本文研究的半变系数再生散度混合效应模型是

变系数模型的推广，对于此类模型的参数和非参数的

估计关键在于条件期望的计算， Lin and Zhang[8]提出

将随机效应当成参数从而用条件众数代替条件期望，

但是这种方法对于非正态的模型估计效果很差。本文

根据 McCulloch[9] 将随机效应看作缺失数据，进而引

入 EM 算法，并在 E 步中使用MCMH方法来计算条

件期望，再利用P-样条对非参数部分进行逼近。EM算

法和 Monte Carlo EM算法，其收敛速度都是线性的，

被缺损信息的倒数所控制，当缺损数据的比例很高时，

收敛速度就非常缓慢。 Monte Carlo EM加速算法(罗季

[10] )在后验众数附近具有二次收敛速度。本文应用

Monte Carlo EM 加速算法估计全部未知参数。并通过

一个实际例子验证了所提出的诊断统计量的有效性。

2. 主要结果

2.1. 模型介绍

假设第 i个接受试验单元第j次的观察值关于

随机效应的条件密度为：



 

,,,; exp;

~0,

ij iijijij

ij i

TTT

ijijijijij i

pyb aydyu

uxwv zb

  



















其中



,,,,, 1,,;1,,

ijij ijijiji

wvyzimjn 



表示 n个

独立的观察数据点，



 ij

u 为未知函数，是位置参数，



是散度参数，为已知函数，为已知单

位偏差度函数，



; a



; d









为联系函数。根据唐年胜和韦博

成 [7] 不妨设假设









 为典则联系，即ij。uu

2.2. 非参数函数的P- 样条估计

对于未知单变量函数









，本文采用 P- 样条估

计。根据 Yu, et al.[11]，假设：

 

01 1

K l

iililri

wwww

 

 r





 





为 K 个样条节点，且为整数。Yu, et al.1l

[11] 详细研究了节点的选择方法，对于光滑函数(取

2 l



并固定选取 5~10 个节点)。通常情况下，取预测

等分位点为节点。如果函数有不连续点，则在

其附近要有一个节点；如果函数有限多极大值和极小

值，则需要取 10 个以上的节点。设样条系数为



,,, T

 



 ，样条基为：

 

变量的



其中



r r





 





,, , T

iiiK

w ww





 1, , ,

Bw w l

l l

则函数









的样条

 

wBw





 。函数为

阵的形式，

将上述向量结合写成矩



,,, T

xx x ,



,,, T

www w 

X 

Vv 



,,, T

TT T

Yyyy 



,,, T

v v ,





ag, ,,di m

zz z，其中





,,,

i in

y y 

 

Bw Bw 

, i

z 类似的定







,,, T

Bw 

,,,

iii

xwv



和

yy 



义。









wBw





 ，则模型可写成如下矩阵形式：



; exp; ay dyu

 





,, ,2

~0,

T T

pyb

uBXV Zb

 















由 Yu, et al.[11]，上述模型的惩罚对数似然函数为：









,, ,,, ,,

PL YL YnK



 















1 1

,, ,

log, ,,

exp d

m i

ij i

i j

iii

pyb

bbb

 

 





 





















其中：表示观测到的数据集，

Y 0



 是光滑参数，K

为对

选

取光

是与节点 t 有关的矩阵，这里取K角矩阵，且只

取最后 K 个对角线元素的值为1，其它为0。

[ 注 ]：参照Yu, et al.[11]，我们可通过GCV方法

滑参数



， GCV的具体算法可以参见Yu, et al.[11]

的 (21)式。具体计算时，可用格子点方法获得最优的



。

2.3. 模型的估计

可根据文献 McCulloch[9]和Zhu, et al.[12]，用EM

算法对模型的参数和非参数进行估计。具体的做法是

半变系数再生散度混合效应模型的影响分析

将随机效应 i

b 看作缺失数据m

Y ，并用



com

YYY 表示完全数据，则完全数据的惩罚对数似

然函数为：







 

, ,

 

log ;;

11 1

22 2

m i

ij ij ij

m TT

aydy u

bbnK



 





















 







 

利用 EM算法求解。标准的 EM 算法包含E步和 M

步，给定初值







；







 



step:Max











其中， m表示EM算法中迭的

次数。

step: mm

EPLYY

 





,,, T

TT T





E 步是对分布

代



pby



求条件期望，而 M步

是求解





使得



m 达到最大。在相对较弱的

条件下，通过反复的迭





代，





能收敛到



的极大似

然估计



。根据以上的算法，可以得到每次迭代

计算



的极大似然估计方程为



：



















条件分布

PL Y



 

















pby



无法直接获得积分的解析表达式 ,

因此用的样本

布为建议分布。假

MH 方法来抽取i

时刻的







:1 1,,

bn m 。基本思想是，选取一

个适当的建议分布，通常选正态分

态为



1 t

b  ，从建议分布



1 2

ibb



 中随机抽取一个潜在的转移值*

b ，从

均匀分布



0, 1 U机抽取一个数如果

，则接受 *

b 作为链在下一时刻的状态

值，即



bb 否则，设

 

1 tt



 ，其中





,, Ni

设抽样链处于第1t状

中随



1 *

ubb









 ；

u ，





1 * ,

,min ,

pb y





 

1, ,

log ,

i t

bii

b pby

Epby

 































并选取 2



的值，以使得整个迭代过程从建议分布中产

生的潜在转移值的接受率位于区间 [0.25,0.34](Gelman,

。

et al. [13] )

Monte Carlo EM 加速算法

1) 选取初值





，令 m= 0 。

抽样

2) 用MH算法生成N个随机



1 ,,

bb ，然

后用似 (Monte Carlo)，记

为

这 N 个样本近计算条件期望







m Y



。

3) 将







ˆ ,



极大化，解出





，使得

 







ˆ ˆ

,max,

mm m

QYQY



 



 。

4) 求解



，使得下式达到最小



1 ln

N k

k bb





 

 。

5) N-R 步：令

 



1 1

mmmm

 

 



 

(5) ，直到存在某个M使得

6) 重复步骤(2)~

 

1 MM













时停止迭代，并取 ˆ



 ，其

中



为指定的充分小的正数。

[ 注 ]：如果迭代收敛，则ˆ



就是



的极大似然

定理 1 ：设



估

计。











 

2 , m



充分靠近 ˆ



，







ˆ 0













ˆ m





。如果正定，且其Hesse矩

阵满对一切 m，当

所得序列

足 Lipschitz 条件。则n充分大时，









收敛到最优解 ˆ



，

收敛速度

对模型讨论个体

在局部影响分析中，讨论了组内加权

扰动、组间加权扰动和随机效应方差的扰动。

并且序列具有二阶

。

证明：见文献 [11]的定理2.1。

3. 影响分析

由于本文是基于纵向数据，因此

删除和组删除。

3.1. 个体删除模型的影响分析

设







c cij

PL y



是删除模型中第 i组第 j个观测值

然函数，相应的 Q

函数为

后所得到的完全数据的惩罚对数似













ˆˆ

ijc ij

QEPLYY









，且假定 ˆ



和



ˆ ij



分别是







和







达到最大值时候

的取值。如果 ˆ



和



ˆ ij



相差很大，则认为第 i组第j

为强实际计算中，如果对每一个



个观测值影响点。



半变系数再生散度混合效应模型的影响分析

都要进行迭代计计算，因此根据 Zhu,

et al. [12] ，用



ˆ ij

算，则量非常大



的一步近似



ˆ ij



来减少计算量：











ˆˆ ˆˆˆˆ

ij ij



 







其中



(1)







ˆˆ ˆ

ij ij



 





 



 ，

 

ˆ ˆˆ



 





 。

又由于



ˆˆ





用 Q

不能定量地表达影响，仿照 Cook

距离，函数构造广义 Cook 距离；

的小











ˆˆ ˆˆˆˆ

ij ij ij

CD Q

 

 



根据 (1)式，可得到一步近似公式：





 









1 ˆˆˆˆˆˆ

ij ij ij

CD QQQ



 







。

在纵向数据模型中，同一组中的观测值通常有

据对于模型的影

组删除模型，同样可以推导出一步近似公式：



3.2. 组删除模型的影响分析

相

同的协变量，因此有必要研究一组数

响。对于











ˆ ˆˆˆ ˆˆ



 





 (2)

其中



















1 ˆ

ˆ ˆˆ

n ico

cij

EPLY Y











 











根据 (2) ，可得到广义Cook距离的一步近似公式：





 







1 ˆ ˆˆˆ

( T

i ii

CD QQQ



 







。

设是定义在开区间

3.3. 局部影响分析



1 ,, T

 

 域 q

R 上

的向量。令



似然函数

为扰动模型的完 Y



。假定存在 0

全惩罚对数



使得







ooo

PLYPL

 







和







PL Y







cccc

PL Y



对所有的



都成立。设 ˆ



和







分别使得 Q 函数





ˆˆ

cco

QEPLYY

 





和





ˆˆ

QEPLY

















 达到最大

值。以上的条件期望均是对条件分布



pYY



求积分。

根据 Zhu, et al.[12]，构造模型的 Q 函数距离：













ˆˆ

LQQ







 ，

其二阶近似为：







1 ˆ

IIT T

LhQ



 h





 

 。

下文，将分别研究三种加权扰动情况：组内加权

动，随机效应方差的扰动。

3.3.1. 组内加权扰动

数据

扰动，组间加权扰

在不考虑数据结构的情况下，在所有观测数据中

找强影响点或异常点，比较常用的方法是给每个

加权。令





1 21

,,, ,, T

nmn

 

 为扰动向量，

当 0



1 n



时，模型为无扰动模型，其中 1

n 为所有元素

为 1的n



维向量，则组内加权扰动的似然函数可表

示为：



 

og;;

m i

cijijij ij

PLYa ydyu

 



bbnK

































 





 

通过对上式求导我们有：



 



cccc

ij ij ij

ij ij

PL YPL Y

PL Y



 

 



 

 



  



PL Y









 



其中：







24 2

, 1

;

, 11

2 ;

, 0

cc T

ij ij

ij ijij

ij ij

PL Y vd

PL Y Bw xd

PL Y d ay

PL Y



 



 





 







 









 

 

 





由此可得到







11 ,, mnm

 

  。

3.3.2. 组间加权扰动

在组间数据中找强影响数据组或异常数据组，比

是给每个数据组加权。令较常用的方法

半变系数再生散度混合效应模型的影响分析



1 ,, T

 

 扰动向量，当为无扰

动模



0 1,1,,1 T



 

型，则组间加权扰动的似然函数可表示为：



 

og ;

11 1

m i

i iijij

m TT

a yu















;

22 2

cc j

PLYd y

bbnK

 







 







 

通过对上式求导我们有：





cccc

i ii

cccc

PL YPL Y

PLY PLY



 

 

 

 

 



  









 



其中：



2 1

24 2

, 1

;

, 11

2 ;

, 0

n i

cc T

ijij

T j

n iT

ij ijij

T j

n iij

i ij

PL Y vd

Y Bw xd

PL Y d ay

PL Y



 









 









 











 





 



 







由此可得到。

3.3.3. 随机效应方差的扰动

在模型中，随机效应是从

2 PL



1 ,, m

 

  



0, N





差阵中随机效

1,, m则组

中随机抽

一步研究应的扰动影

间加权扰动

的似

取的，为了进协方

响，假设



Var bi







然函数可表示为：



 

,l ;

m i

cij ij

m T

PLd yu

 















og;

22 2

c ij

ii i

Ya y



 























通过对上式求导我们有：



cccc

i ii

cccc

PL YPL Y

PLY PLY



 

 

 



 



  



















 

, 0

, 1

PL Y

PLY bb













 



 



 



 



 



 



 





由此可得到





1 ,, m

 

  。

4. 实例分析

结合药物血浆渗透数据 (Davidian and Gilti-

nan [14] )。来说明本文给出的统计诊断方法的可操作性

和有效性。

人自愿者，在 8内通过 11 次静脉

药物，测得每位病人血液中药物浓度

。本节用半变系数再生散度混合效应模型拟合该

数据

对 6个病小时

注射相同剂量的

数据

。假设 ij i

yb 服从单参数I型极值分布，则ij i

yb 的

概率密度函数为

其中：











exp exp b yy





ijiijij ijij

其中：





ij ijijiji

xxb

 



 ，则



;

ij ij



满足单位

偏差度函数及再生散度定义的条件。

实际计算中，对于变系数部分，选取固定的节

点且阶数为 3 阶，并通

5 个

过 GCV的方法得到光滑参数

的估计 4

1.34210 n





 ，然后用Monte Carlo EM加速

算法得到诊断统计量的值。下面列出广义 Cook 距离

响的图形：

博成，林金官

第 3 组数据的影响最大，因为

它包 3 号点，同时第1组数据的影

响也

义

Coo 一致，从另一个角度证明了广义

Coo

和局部影

1) 个体删除模型

从图 1中可以发现第 1号和第 23 号点的影响比

较大，其中第 23 号点是数据中数值最大的点。且以上

结果与韦，解锋昌 [15] 的结果一致。

2) 组删除模型

从图 2 中可以发现

含了强影响点第 2

较大，因为包含了强影响点第1号点。

3) 局部影响分析

从图 3、图 4和图 5 可以发现，局部影响分析和广

k 距离的结果基本

k 距离的有效性。

半变系数再生散度混合效应模型的影响分析

Figure 1. Generalized Cook distance of individual delete model

图 1.个体删除模型的广义Cook距离

Figure 2. Generalized Cook distance of group delete model

图 2.组删除模型的广义Cook距离

Figure 4. Between group disturbance

图 4.组间扰动

Figure 5. Rando disturbance

图 5.随机效应方差的扰动

5. 致谢

感谢审稿的宝贵意见以及编辑的帮助。

参考文献 (References)

[1] R. D. Cook, S. Weisberg. Residual and influence in regression.

New York: Chapman and Hall, 1982.

[2] G. Seber, C. J. Wild. Nonlinear Regression. New York: Wiley,

1989.

[3] 姜荣,邵明江,钱伟民.半参数非线性模型中的t-型估计和影

响分析 [J]. 华东师范大学学报 ( 自然科学版 ), 2011, 3: 1-11.

[4] R. J. Beckman, C. J. Nachtsheim and R. D. Cook. Diagnostics

for mixed-model analysis of variance. Technometrics, 1987, 29:

413-

[5] 张浩析[J].应

用数学学报 , 2

] B. Jorgensen. The theory of dispersion models. London: Chap-

m effect variance

(4)

426.

, 朱仲义. 半参数广义线性混合效应模型的影响分

Figure 3. Within group disturbance

图 3.组内扰动007, 30(4): 773-756.

半变系数再生散度混合效应模型的影响分析

d Hall, 1997

, 韦博成 .非线性再生散度模型[M].北京:科学

社 , 2007.

[8] tive mixed

models by using smoothing splines. Journal of the Royal S

tical Society, Series B, 1999, 61(2): 381-400.

[9] C. E. Mood algorithm for general-

nal of the American Statistical

Association,170.

single-index models. Journal of the American Statistical

n, 2002, 97(460): 1042-1054. man an

[7] 唐年胜出版

X. H. Lin, D. W. Zh ang. Infere nce in generalized additatis-[13]A. Gelman, G. O. Roberts and W. R. Gilks. Efficient metropolis

jumping rules. Bayesian statistics 5, Oxford: Oxford University

Press, 1995.

cCulloch. Maximum likelih

ized linear mixed models. Jour

1997, 92(437): 162- [14]

[10] 罗季. Monte Carlo EM加速算法[J].应用概率统计, 2008, 24(3):

311-318.

[11] Y. Yu, D. Ruppert. Penalized spline estimation for partially

linear

Associatio

[12] H. T. Zhu, S. Y. Lee, B. C. Wei and J. Zhu. Case-deletion meas-

ures for models with incomplete data. Biometrika, 2001, 88(3):

727-737.

M. Davison, D. M. Giltinan . Nonlinea r models for repeated meas-

urement data. London: Chapman and Hall, 1995.

[15] 韦博成,林金官,解锋昌.统计诊断[M].北京:高等教育出版

社 , 2009.