 Statistical and Application 统计学与应用, 2012, 1, 37-43 http://dx.doi.org/10.12677/sa.2012.12008 Published Online December 2012 (http://www.hanspub.org/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. 19th, 2012; revised: Nov. 26th, 2012; accepted: Dec. 13th, 2012 Abstract: This paper proposes several case-deletion as well as local 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 model as missing data and estimate unknown parameters by acceleration of Monte Carlo EM algorithm. On the basis of the Q-function 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; Acceleration 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]研究了非线性再生散度随机效应 Copyright © 2012 Hanspub 37  半变系数再生散度混合效应模型的影响分析 模型，讨论了该模型的几何结构、渐近性质和统计诊 断等问题。 变系数模型在生物医学、公共卫生、经济、农业、 制造业、道路安全等众多领域的数据分析中有广泛的 应用。本文研究的半变系数再生散度混合效应模型是 变系数模型的推广，对于此类模型的参数和非参数的 估计关键在于条件期望的计算，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次的观察值 关于 随机效应 的条件密度为： ij y i b 22 2 1 ,,,; exp; 2 ~0, ij iijijij yb ij i i TTT ijijijijij i pyb aydyu bN uxwv zb 其中 ,,,,, 1,,;1,, ijij ijijiji x wvyzimjn 表示 n个 独立的观察数据点， ij u为未知函数， 是位置参数， 2 是散度参数， 为已知函数， 为已知单 位偏差度函数， ;a ;d 为联系函数。根据唐年胜和韦博 成[7]不妨设假设 ij 为典则联系，即 ij 。 uu 2.2. 非参数函数的 P-样条估计 对于未知单变量函数 ，本文采用 P-样条 估 计。根据 Yu, et al.[11]，假设： 01 1 Kl l iililri r wwww r 为K个样条节点， 且为整数。Yu, et al. 1l [11]详细研究了节点的选择方法，对于光滑函数(取 2l 并固定选取 5~10 个节点)。通常情况下，取预测 等分位点为节点。如果函数有不连续点，则在 其附近要有一个节点；如果函数有限多极大值和极小 值，则需要取 10个以上的节点。设样条系数为 01 ,,,T lK ，样条基为： 变量的 其中 1 K rr 1 ,, ,T iiiK w ww 1, , , ii Bw wl ll 则函数 的样条 T ii wBw 。函数为 阵的形式， 将上述向量结合写成矩 12 ,,, T m xx x, 12 ,,, T m www w , X Vv , 12 ,,, T TT T m Yyyy 12 ,,,T m v v, 12 ag, ,,dim Z zzz，其中 12 ,,, ii i in i y y 12 Bw Bw ,i z类似的定 ,,,T m Bw ,,, T iii xwv 和 yy Bw 义。 T wBw ，则模型可写成如下矩阵形式： 22 1 ; exp;ay dyu 2 ,,,2 ~0, yb TT pyb bN uBXV Zb 由Yu, et al.[11]，上述模型的惩罚对数似然函数为： 22 1 ,, ,,, ,, T PL YL YnK 2 0 oo oo 2 12 2 11 1 ,, , log, ,, 1 exp d 2 oo n mi ij i yb ij i ij T iii LY pyb bbb 其中： 表示观测到的数据集， o Y0 是光滑参数，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 算法对模型的参数和非参数进行估计。具体的做法是 Copyright © 2012 Hanspub 38  半变系数再生散度混合效应模型的影响分析 将随机效应 i b看作缺失数据 m Y，并用 , com YYY表示完全数据，则完全数据的惩罚对数似 然函数为： 2 , , cc n Y 2 2 11 1 1 , 1 log ;; 2 11 1 ln 22 2 mi ij ij ij ij mTT ii i PL aydy u bbnK 利用 EM 算法求解。标准的 EM 算法包含 E步和 M 步，给定初值 0 ； 1 , step:Max cc o mm EQ MQ 其中 ，m表示EM 算法中迭 的 次数。 step:mm EPLYY 2 ,,, T TT T E步是对分布 代 , ii pby 求条件期望，而 M步 是求解 m 使得 m达到最大。在相对较弱的 条件下，通过反复的迭 Q 代， m 能收敛到 的极大似 然估计 。根据以上的 算法，可以得到每次迭代 计算 EM 的极大似然估计方程为 ： ,0 m o Y 条件分布 cc m PL Y QE , ii pby 无法直接获得积分的解析表达式, 因此用 的样本 布为建议分布。假 MH 方法来抽取 i b 时刻的 :1 1,, n i bn m。基本思想是，选取一 个适当的建议分布，通常选正态分 态为 1t i b，从建议分布 12 , t ibb Nb 中随机抽取一个潜在的转移值 * i b，从 均匀分布 0, 1U机抽取一个数 如果 ，则接受 * i b作为链在下一时刻的状态 值，即 t ii bb否则，设 1tt ii bb ，其中 ,,Ni 设抽样链处于第 1t状 中随 1* , t ii ubb * ； u， * 1*, ,min , i i pb y b 1 2 1, , log , i t it ii bii T ii bpb y Epby bb 并选取 2 b 的值，以使得整个迭代过程从建议分布中产 生的潜在转移值的接受率位于区间[0.25,0.34](Gelman, 。 et al.[13]) Monte Carlo EM 加速算法 1) 选取初值 0 ，令 m = 0。 抽样 2) 用MH算法生成 N个随机 1,, N bb，然 后用似(Monte Carlo)，记 为 这N个样本近 计算条件期望 , mY 。 3) 将 ˆ Q ˆ, m QY 极大化，解出 m E M ，使得 ˆˆ ,max, mm m EM QYQY 。 4) 求解 ，使得下式达到最小 1 1 1ln 2 Nk ii kbb N 。 5) N-R步：令 11 mmmm QQ (5)，直到存在某个 M使得 . 6) 重复步骤(2)~ 1MM 时停止迭代，并取ˆ M ，其 中 为指定的充分小的正数。 [注]：如果迭代收敛，则 ˆ 就是 的极大似然 定理 1：设 估 计。 m QC 2,m 充分靠近 ˆ ， ˆ0 m Q ˆm Q 。如果 正定，且其 Hesse 矩 阵满对一切 m，当 所得序列 足Lipschitz 条件。则 n充分大时， m 收敛到最优解ˆ ， 收敛速度 对模型讨论个体 在局部影响分析中，讨论了组内加权 扰动、组间加权扰动和随机效应方差的扰动。 并且序列具有二阶 。 证明：见文献[11]的定理2.1。 3. 影响分析 由于本文是基于纵向数据，因此 删除和组删除。 3.1. 个体删除模型的影响分析 设 ccij PL y 是删除模型中第 i组第 j个观测值 然函数，相应的 Q 函数为 后所得到的完全数据的惩罚对数似 ˆˆ , co ijc ij QEPLYY ，且假定 ˆ 和 ˆij 分别是 ˆ Q 和 ˆ ij Q 达到最大值时候 的取值。如果ˆ 和 ˆij 相差很大，则认为第i组第 j 为强 实际计算中，如果对每一个 ij 个观测值 影响点。 Copyright © 2012 Hanspub 39  半变系数再生散度混合效应模型的影响分析 都要进行迭代计计算，因此根据 Zhu, et al.[12]，用 ˆij 算，则 量非常大 的一步近似 1 ˆij 来减少计算量： 1 1 ˆˆ ˆˆˆˆ ij ij QQ 其中 (1) ˆ ˆˆ ˆ ij ij QQ ， 2 ˆ ˆˆˆ T QQ 。 又由于 ˆˆ ij 用Q 不能定量地表达影响，仿照 Cook 距离，函数构造广义 Cook 距离； 的小 ˆˆ ˆˆˆˆ T ij ij ij CD Q 根据(1)式，可得到一步近似公式： 1 1ˆˆ ˆˆˆˆ T ij ij ij CD QQQ 。 在纵向数据模型中，同一组中的观测值通常有 据对于模型的影 组删除模型，同样可以推导出一步近似公式： 3.2. 组删除模型的影响分析 相 同的协变量，因此有必要研究一组数 响。对于 1 1 ˆˆˆˆ ˆˆ ii QQ (2) 其中 ˆ 1ˆ ˆˆˆ ˆ , ii nico cij j QQ EPLY Y 根据(2)，可得到广义 Cook 距离的一步近似公式： 1 1ˆˆˆˆ (T iii CD QQQ 。 设是定义在开区间 3.3. 局部影响分析 1,,T q 域q R上 的向量。令 , c PL 似然函数 c 为扰动模型的完Y 。假定存在 0 全惩罚对数 使得 0 ooo PLYPL , o Y 和 PL Y 0 , cccc PL Y 对所有的 都成立。设 ˆ 和 ˆ 分别使得 Q函数 ˆˆ , cco QEPLYY 和 ˆˆ o Y ,, , QEPLY cc 达到最大 值。以上的条件期望均是对条件分布 ˆ , mo pYY 求积分。 根据 Zhu, et al.[12]，构造模型的 Q函数距离： ˆˆ 2 Q LQQ ， 其二阶近似为： 1ˆ IIT T Q LhQ h 。 下文，将分别研究三种加权扰动情况：组内加权 动，随机效应方差的扰动。 3.3.1. 组内加权扰动 数据 扰动，组间加权扰 在不考虑数据结构的情况下，在所有观测数据中 找强影响点或异常点，比较常用的方法是给每个 加权。令 1 21 1 ,,, ,,T nmn m 为扰动向量， 11 当0 1n 时，模型为无扰动模型，其中1 1 n为所有元素 为1的n 维向量，则组内加权扰动的似然函数可表 示为： 2 2 11 og;; 2 11 ln 22 mi cijijij ij ij m PLYa ydyu 1 1 1 ,l 1 2 n c TT ii i bbnK 通过对上式求导我们有： 22 2 2 ,, ,, , , cccc TT ij ij ij T cc T ij ij PL YPL Y PL Y 2 , cc PL Y 其中： 2 2 2 2 2 2 24 2 2 2 ,1 2 ,1 2 ; ,11 2; ,0 cc T ij ij T ij T cc ij ijij T ij ij cc ij ij ij cc T ij PL Yvd PL YBw xd ay PL Yday PL Y 由此可得到 11 ,,mnm 。 3.3.2. 组间加权扰动 在组间数据中找强影响数据组或异常数据组，比 是给每个数据组加权。令 较常用的方法 Copyright © 2012 Hanspub 40  半变系数再生散度混合效应模型的影响分析 1,, T m 扰动向量，当 为无扰 动模 01,1,,1 T 型，则组间加权扰动的似然函数可表示为： 2 2 11 1 1 og ; 11 1 n mi i iijij j mTT a yu 1 ,l ; 2 ln 22 2 cc j i ii i PLYd y bbnK 通过对上式求导我们有： 22 22 2 ,, ,, ,, , cccc TT iii T cccc T ii PL YPL Y PLY PLY 其中： 2 21 21 2 2 24 2 2 1 2 ,1 2 ,1 2 ; ,11 2; ,0 ni cc T ijij Tj i niT cc ij ijij Tj i niij cc ij j iij cc T i PL Yvd YBw xd ay PL Yday PL Y 由此可得到 。 3.3.3. 随机效应方差的扰动 在模型中，随机效应 是从 2PL 1,,m i b 0,N 差阵中随机效 1,,m则组 中随机抽 一步研究 应的扰动影 间加权扰动 的似 取的，为了进 协方 响，假设 , ii Var bi 然函数可表示为： 2 2 1 ,l ; 11 ln n mi cij ij mT PLd yu bb 11 1 1 og; 2 1 22 2 c ij ij T ii i i Ya y nK 通过对上式求导我们有： 22 22 2 ,, ,, ,, , cccc TT iii T cccc T ii PL YPL Y PLY PLY 2 2 2 2 2 11 ,0 ,0 ,0 ,1 2 cc T i cc T i cc i cc T ii TT i PL Y PL Y PL Y PLY bb 由此可得到 1,,m 。 4. 实例分析 结合药物血浆渗透数据(Davidian and Gilti- nan[14])。来说明本文给出的统计诊断方法的可操作性 和有效性。 人自愿者，在8内通过 11 次静脉 药物，测得每位病人血液中药物浓度 。本节用半变系数再生散度混合效应模型拟合该 数据 对6个病 小时 注射相同剂量的 数据 。假设 ij i yb服从单参数 I型极值分布，则ij i yb的 概率密度函数为 其中： exp expb yy ijiijij ijij 其中： py ij ijijiji x xxb ，则 ; ij ij dy 满足单位 偏差度函数及再生散度定义的条件。 实际计算中，对于变系数部分，选取固定的节 点且阶数为3阶，并通 5个 过GCV 的方法得到光滑参数 的估计 4 1.34210n ，然后用 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距离的有效性。 Copyright © 2012 Hanspub 41  半变系数再生散度混合效应模型的影响分析 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. [6 Copyright © 2012 Hanspub 42  半变系数再生散度混合效应模型的影响分析 Copyright © 2012 Hanspub 43 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. Inference 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. |