总结
我们在Tweedie指数分散模型分布的基础上,引入了一类新的广义线性混合模型,适用于广泛的离散、连续和混合数据。使用随机效应的最佳线性无偏预测,我们获得了Godambe意义下回归参数的最佳估计函数,从而为整个类提供了有效的公共拟合算法。虽然允许全参数推断,但我们的主要结果仅取决于未观察到的随机效应的一阶矩和二阶矩假设。此外,我们还获得了回归参数和离散参数的一致估计。我们通过分析癫痫数据和蛋糕烘焙数据来说明该方法。通过仿真和渐近证明,这表明了该方法在分析聚集非正态数据时的有用性。
1.简介
已经为聚类数据的广义线性混合模型开发了各种方法(例如,参见Breslow和Clayton(1993)还有Lee和Nelder(1996))但定义可能性的高维积分的难解性阻碍了计算效率和统计一致性方法的发展。密集型计算方法,如贝叶斯方法(Zeger和Karim,1991)蒙特卡洛电磁法(McCulloch,1997)已提出;但这种方法对于具有复杂随机效应结构的海量数据集通常是不切实际的(Skrondal和Rabe-Hesketh,2004).
被广泛采用的更简单的方法涉及Henderson线性混合模型方程的推广(1975),通过区分回归参数和随机效应的“联合”可能性获得(Sholl,1991;Breslow和Clayton,1993; 沃尔芬格和奥康奈尔,1993; 麦吉尔克里斯,1994; Lee和Nelder,1996). 由此产生的随机效应预测因子,尽管在分析上是可处理的,但在线性混合模型中是最佳线性无偏预测因子(BLUP)的不完美类似物,应该更恰当地称为条件模式预测因子,因为它们与给定数据的随机效应的条件分布模式相对应。
对于非正态分布,这些条件模式预测通常既不是线性的,也不是无偏的,而是基于给定数据的随机效应条件分布的某种正态近似(Draper,1996)通常将随机效应分布限制为正态或共轭分布(Breslow和Clayton,1993; Lee和Nelder,1996). 条件模式方法处理大量随机效应的能力也受到操作大型矩阵的必要性的限制。
更重要的是,条件模式预测通常会导致有偏估计方程和不一致的参数估计,尤其是在簇数增加的情况下(Jiang,1998; 哈等人。,2001). 为了获得一致的参数估计,Jiang(1998)介绍了计量经济学文献中用于处理广义线性混合模型的模拟矩方法,但该方法既没有计算效率,也没有统计效率。
本文基于Tweedie指数分散模型分布的柔性类(Jörgensen,1987,1997). 我们的模型显示了方差分量分解,它巧妙地反映了聚类数据的层次结构。我们得到了Godambe意义下回归参数的最优估计方程(1976)基于给定数据的随机效应BLUP。我们还证明了聚类数随样本大小增加的情况下的渐近性质。我们的方法仅依赖随机效应的第一和第二矩。与条件模式方法相比,该方法计算简单,效率更高,不需要求大矩阵的逆。
在简要讨论了条件期望之后与第节中的条件模式2,第节介绍了Tweedie随机效应模型三.在节中4和5我们分别讨论了随机效应预测和模型估计。章节中给出了示例和模拟6和7.
2.随机效应预测
随机效应的预测在随机效应模型中起着重要作用(Robinson,1991). 为了强调我们的方法与现有方法之间的差异,我们现在将条件期望和BLUP与条件模式预测进行比较。
我们首先引入部分观测联合对数似然的概念。让Y(Y)是回应和U型是未观察到的随机效应。我们假设Y(Y)鉴于U型以及U型是参数化的。部分观察到的“联合”原木定义为
哪里我(β,α;Y(Y)|U型)是的对数条件密度Y(Y)鉴于U型带回归参数β和色散参数α、和我(γ;U型)是的对数密度U型带参数γ.
对于正规理论线性混合模型,对于固定(β,γ),解决方案方程的导数我(β,α,γ;年,u个)/∂u个=0符合条件期望和BLUPU型鉴于Y(Y)(亨德森,1975). 预测器因此,在广义线性混合模型的背景下,作为对U型鉴于Y(Y)在EM算法的E步骤中(Dempster等人。,1977). 然而,实际上是条件模式预测器U型鉴于Y(Y),自
其中第二项不依赖于U型然而,条件模式不一定是非对称条件分布的条件平均值的良好近似;评估其偏离条件平均值也不容易。
3.Tweedie混合模型
在这一节中,我们介绍了一类基于Tweedie指数色散分布的广义线性模型,该模型具有两层嵌套随机效应,1947,1984). 具体来说,随机变量Y(Y)据说遵循Tweedie指数色散分布Tw对(μ,τ2)如果它的密度是这样的
其中的显式表达式c(c)对(年;τ2)由约根森(1987,1997)但在我们接下来的章节中,这些都是无关紧要的。这个Tweedie家族也被称为power–variance家族,因为E类(Y(Y))=μ和var(Y(Y))=τ2μ对.
Tweedie家族实际上涵盖了广泛的知名发行版。使用通常的符号,我们有Tw1(μ,1)=泊松(μ),吐温2(μ,τ2)=伽马(μ,τ2)和Tw三(μ,τ2)=逆-高斯(μ,τ2). 正态和某些复合泊松分布对应于Tweedie分布对=0和1<对分别小于2。有关Tweedie家族的更多详细信息,请参见Jörgensen(1987,1997).
3.1. 模型
在本节中,我们考虑三级层次模型,其中每个模型由以下部分组成米独立簇索引我。在每个集群内我,有J型我由索引的相关子簇(我,j个). 此外,在每个子集群中(我,j个),有n个ij公司相关观察结果。
让Y(Y)ijk公司代表k个子星团中的第个观测(我,j个)具有协变量x个ijk公司。我们假设存在特定于簇和子簇的随机效应U型我和U型ij公司对于群集我和子集群(我,j个)分别是。让Y(Y)=(Y(Y)111,…,Y(Y)11n个11,…,Y(Y)百万焦米1,…,Y(Y)百万焦米n个百万焦米)T型表示观测矢量。让U型=(U型*,U型**)表示随机效应的向量,其中U型*和U型**代表(U型1,…,U型米)T型和(U型11,…,U型百万焦米)T型分别是。我们通过以下三个假设来描述我们的模型。
假设1。集群随机效应U型1,…,U型米是正的、独立的且与
1
假设2。给定集群级随机效应U型*=u个*=(u个1,…,u个米),亚簇级随机效应U型11,…,U型百万焦米为正且条件独立,其中U型ij公司,给定U型*=u个*,取决于u个我仅,具有
2
参数ν2可能取决于我和j个解释边际异质性;参见第节5.2.
假设3。给定随机效应U型=u个,的组件Y(Y)是条件独立的Y(Y)ijk公司,给定U型=u个,取决于u个ij公司只有,这是
三
哪里.
我们对随机效应的假设1和2具有某些独立性,并且仅基于前两个矩。这是可取的,因为潜在的随机机制通常并不完全已知。然而,我们的假设1和2确实涵盖了广泛的参数随机效应,包括众所周知的伽马分布、逆高斯分布和对数正态分布(Ma,1999); 这种参数解释通常有助于模型检查和模型解释。根据Tweedie族,可以指定一大类随机效应分布,如下所示:U型我~Tw第页(1,σ2)和Tweedie分布Tw第页,吐温q个和Tw对分别称为簇、子簇和观测水平分布。为了避免非正面随机效应,第页2和q个需要2个。观测水平分布的条件平均值为E类(Y(Y)ijk公司|U型)=μijk公司U型ij公司; 因此,在log-link函数下,回归参数和随机效应都是线性的。此外,公约E类(U型我)=E类(U型ij公司)=1由这个乘法模型隐含。因此,我们E类(Y(Y)ijk公司)=E类{E类(Y(Y)ijk公司|U型)}=μijk公司(U型ij公司)=μijk公司.Lee和Nelder(1996)指出,这最后一个等式避免了区分人口平均数和特定主题模型的需要(Liang和Zeger,1986).
3.2. 协方差结构和模型分解
使用Kronecker表示法δ(秒,我)为1,如果秒=我否则,协方差结构可以表示如下。协方差结构反映了集群内的内部依赖性:
4
由方程给出的观测值的协方差(4)显示方差分量结构。除伽马混合模型外,响应的相关矩阵通常取决于回归参数。请注意,我们的伽马混合模型意味着与标准固定效应情况一样的恒定变异系数。
为了将我们的模型假设与其他假设进行比较,我们给出了以下内容乘法分解对于我们的模型:
哪里V(V)ij公司=U型ij公司/U型我在文献中,随机效应U型我和V(V)ij公司通常被认为是独立的。也可以很容易地验证cov(U型我,V(V)ij公司)在我们的Tweedie混合模型中=0(Ma,1999). 这种乘法分解将我们的模型与莫顿考虑的乘法随机效应模型联系起来(1987)还有费斯和哈里斯(1991). 除了乘法分解,该模型的响应还包括以下内容加性分解具有三个不相关但相关的组件:
5
右侧的三个分量分别对应于假设3、2和1中给出的(条件)分布的残差。通过下一节描述的协方差结构,可以很容易地证明这三个残差是不相关的。具有单一随机效应的Tweedie混合模型是通过设置ν2=0和J型我=1代表所有我.
4.随机效应的最佳线性无偏预测
让协方差矩阵的逆Y(Y)用var表示−1(Y(Y)). 与Ma一样等人。(2003),我们通过以下正统BLUP预测随机效应U型鉴于Y(Y):
6
这是线性无偏预测值U型鉴于Y(Y)使随机效应之间的均方距离最小化U型以及它们在线性函数类中的预测器Y(Y)(哈维,1981; 约根森等人。,1996b) ●●●●。
应用导出的协方差矩阵逆的显式表达式附录A到等式(6),集群随机效应预测值可以表示为
哪里
用于固定(我,j个). 类似地,子簇随机效应预测值可以表示为
这些预测因素和当回答是肯定的时候,它们显然是非否定的。
随机效应分量之间均方距离的显式表达式U型其预测因子如下:
7
和
8
与Ma一样等人。(2003)在“小离散渐近性”(Jörgensen,1987)大样本渐近性;见马(1999).
5.模型估算
5.1. 参数估计
在已知离散参数的情况下,首先考虑回归参数的估计。区分Tweedie混合模型的部分观测“联合”对数似然数据和随机效应β生成部分观察到的“关节”得分函数。用BLUP代替随机效应,我们得到回归参数的无偏估计函数β:
9
这相当于使用BLUP来近似条件期望U型鉴于Y(Y)在EM算法的E步骤中。
与Ma一样等人。(2003),此的全局矩阵表达式ψ(β)灵敏度矩阵之间的简单关系S公司(β)=E类{∂ψ(β)/∂βT型},可变性矩阵V(V)(β)=E类{ψ(β)ψT型(β)}和Godambe信息矩阵J型(β)=S公司(β)V(V)(β)−1 S公司(β)T型可以获得。由于接头密度(Y(Y),U型)在回归参数中形成指数族β,马的证明等人。(2003)适用于ψ(β)在方程式中(9)只要两者的得分函数都是线性的U型和Y(Y)然而,这种线性是清楚的,因为这个分数函数实际上是方程的右侧(9)具有被替换为U型ij公司因此,我们得出以下结果。
定理1。对于Tweedie混合模型,预测得分函数ψ(β),灵敏度矩阵S公司(β)和可变性矩阵V(V)(β)可以表示为:
通常,第三个语句对其他估算函数无效。Tweedie混合模型的这种关系类似于Fisher信息矩阵的性质。我们的定理1也适用于Jörgensen使用的估计方程等人。(1996a、,b条,c(c))在增长曲线模型和状态空间模型的背景下,因为我们的理由成立,而不诉诸于潜在过程之间的依赖性假设。因此,使用该定理可以大大简化其哥达姆信息矩阵的计算。
由于定理1适用于每个簇,因此方程中的估计函数(9)可以表示为
这个估计函数可以很容易地被证明是最优的,因为它达到了估计量的最小渐近协方差β在所有线性估计函数的类中Y(Y); 参见Crowder(1986,1987)和Ma(1999). 当Y(Y)只能指定到协方差,这种形式的估计方程似乎只在前两个矩的假设下采用了模型的最大效率。因此,文献中通过不同的方法,如拟似然法(McCullagh,1983),广义估计方程(Liang和Zeger,1986)、边缘准似似群(Breslow和Clayton,1993)和多级建模(Goldstein,1995); 然而,与这些边际建模方法不同,我们的方法明确预测了随机效应。此外,我们的模型具有明确的方差分量结构,其中相关矩阵取决于回归和离散参数。
我们现在假设色散参数未知。受广义线性模型的启发,我们采用了以下方法调整后的皮尔逊估值器(约根森等人。,1996b),即通过其偏差校正调整的Pearson估计器。
基于Tweedie混合模型的矩结构,我们可以构造色散参数的无偏估计σ2,ν2和ρ2如下(Ma,1999):
10
11
和
哪里c(c)我和c(c)ij公司在方程式中定义(7)和(8)分别是。和剩余最大似然估计量一样,我们也可以考虑小样本自由度校正。例如,直觉上,我们可以替换米通过米减去估计的回归参数数量。模拟研究(Ma,1999; 哈和李,2005)已经证明了这一点特别的小样本校正可能会带来更好的性能。
5.2. 渐近性质和边际异质性
通过出租ξ表示(σ2,ν2,ρ2)T型和更换在方程式中(10), (11)和(12)由(σ2,ν2,ρ2)T型很明显,这三个方程可以重写为无偏估计方程的总和:与一起我们得到了一组无偏估计方程。标准渐近理论在估计方程中的应用(ψ(β,ξ)T型,φ(β,ξ)T型)T型=0,我们得到根的序列,,与一致(β\升降箱3.6pt(磅)\脚本样式{{\rm{T}}},ξ\升降箱3.6pt(磅)\脚本样式{{\rm{T}})T型渐近正态和渐近平均(βT型,ξT型)T型和其哥达姆信息矩阵的逆给出的渐近方差(βT型,ξT型)T型因此仍然是β,而现在必须在哥达姆信息矩阵的基础上为(βT型,ξT型)T型,作为米→∞.
要获得的渐近方差β在这个扩展的参数空间下,我们使用了Jörgensen和Knudsen的结果(2004)关于具有扰动参数不敏感性质的估计方程。他们已经证明了β在此扩展参数下,空间再次由J型−1(β)就好像色散参数与我们的估计方程一样ψ(β,ξ)是ξ不敏感,即。E类{∂ψ(β,ξ)/∂ξ}=0然而,注意数量Y(Y)−E类(Y(Y))不涉及色散参数ξ,这个ξ-我们的Tweedie混合模型(Ma,1999). 没有ξ-不敏感,渐近方差对于离散参数未知的情况,通常涉及回归参数和离散参数的变异性和灵敏度矩阵。如果没有这个ξ-由于不敏感,这种幼稚的插件估计器的使用通常需要更加谨慎,因为它忽略了因估计离散参数的需要而产生的额外可变性。
在到目前为止讨论的模型中,假设分散参数在集群中保持不变,即边际方差仅以系统方式随边际均值变化。然而,在生物学研究中,由于遗传、环境或其他未知因素,集群间的边际异质性往往超过平均变化所能解释的范围。为了解释响应的边际分布中的这种不规则异质性,我们可以允许特定于簇的分散参数和不在方程中跨簇求平均(11)和(12). Thej个-比分散参数也可以同样适用于平衡设计J型1=…=J型米=J型例如纵向数据。估计值对于大量集群是一致的,而对于不是。
5.3. 计算程序和残差分析
我们采取了以下措施牛顿记分算法由约根森介绍等人。(1996a)求解估算方程ψ(β)=0。因此,生成的算法给出了以下更新值β:
Godambe信息矩阵J型(β)=−S公司(β)这里所起的作用类似于Fisher评分算法中的Fisher信息矩阵。
自S公司(β)=−V(V)(β),通过插入var的表达式,可以得到灵敏度矩阵的精确表达式−1(Y(Y)我)英寸附录A转化为定理1中可变性矩阵的表达式:
其中第一项与假设独立观测的模型的Fisher信息一致,第二项和第三项对应于由聚类和亚聚类随机效应引起的变化。正如马所言等人。(2003)回归参数估计的渐近方差β如果两个模型的估计值相同,则基于Tweedie混合模型的结果大于基于标准Tweedi模型的结果,且没有随机影响。
回归参数估计的初始值被视为回归参数估计值,这些参数估计值是从假设独立响应的标准广义线性模型技术中获得的。我们对和作为集群内响应的平均值我除以所有响应的平均值和子集群内响应的平均数(我,j个)分别除以所有响应的平均值。初始弥散参数估计值由Pearson估计值通过方程计算得出(10)–(12),但省略了偏差校正项。然后,该算法在通过Newton计分算法更新回归参数估计值、通过BLUP更新随机效应预测值和通过调整的Pearson估计值更新离散参数估计值之间进行迭代。
响应和随机效应的残差分析是我们方法的一个重要组成部分。将响应的边际残差分解为方程中的三个不相关残差(5)对应于模型中给出的三级分布假设。因此,我们可以通过估计残差来检查1、2和3级的分布假设和分别是。所有这三个估计残差均为零均值,可以标准化为单位方差。残差分析的基本思想是使用标准化残差图,其方法与广义线性模型中的方法大致相同。
6.数据分析
通过对癫痫数据(Thall和Vail,1990)和蛋糕烘焙数据(Cochran和Cox,1957)分别在本节中列出。
6.1. 计数数据
在Ma等人。(2003),我们对一个特殊模式的泊松混合模型的方法被应用于一个大规模数据集的生存分析,其中大约5000万个观测值被分为两个级别。为了便于将BLUP方法与其他方法进行比较,我们考虑了Thall和Vail提供的众所周知的癫痫发作计数数据(1990). 在这项临床试验中,59名癫痫患者被随机分为新的抗癫痫药物progabet(Trt=1)或安慰剂(Trt=0)。试验开始时可用的基线数据包括随机分组前8周内癫痫发作次数和患者年龄。多变量反应变量包括四次就诊前两周内的癫痫发作计数。初步分析表明,在第四次访问期间,计数要低得多,因此Breslow和Clayton(1993)引入了线性趋势协变量访视(编码为(−0.3,−0.1,0.1,0.3))。
我们对这些数据的重新分析基于泊松混合模型。我们取基线发作计数(base)的四分之一的对数、年龄(age)的对数、Trt、Visit和相互作用项base。Trt作为协变量。由于每次访问都没有重复,因此省略了第三个索引。我们表示年ij公司患者的癫痫发作次数我访问中j个; 因此随机效应u个我和u个ij公司分别处于患者和就诊级别。模型在假设1-3下规定如下:
哪里j个=1,2,3,4和我=1,…,59. 受试者内部的聚集效应和观察水平的过度分散表现为嵌套随机效应u个我和u个ij公司分别是。
我们考虑四个具有不同色散参数结构的模型其中模型(ν2=0),型号(ν2),型号()和型号()假设分别,即模型(ν2=0)是仅具有簇级随机效应的泊松混合模型。型号(ν2),型号()和型号()是具有两级随机效应的模型,对应于不同的子簇级随机效应结构。我们拟合了这四个模型,并通过残差分析进行了模型检验。模型的正常绘图(ν2),型号()和型号()如图所示1(a) –1(c) ,图1(d) –1(f)和图1(g) –1(i) 分别是。三行中的图对应于标准化的聚类、子聚类和观测水平残差。明显偏离常态;然而,对于随机效应广义线性模型,残差的正态性似乎从未在理论上得到证明,尽管它仍然被用于进行常规正态图,作为文献中模型检查的一部分(Lee和Nelder,1996). 然而,适度偏离正常状态可能是
![癫痫数据的1级、2级和3级残差的正态图:(a)模型(ν2),1级;(b) 模型(ν2),2级;(c) 模型(ν2),3级;(d) 模型(vj2),1级;(e) 型号(vj2),级别2;(f) 模型(vj2),3级;(g) 模型(vi2),1级;(h) 模型(vi2),2级;(i) 模型(vi2),3级](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/jrsssb/69/4/10.1111_j.1467-9868.2007.00603.x/1/m_rssb_69_4_625_f1.gif?Expires=1722480092&Signature=YNlO3PaA42YqAnjOLk0R79IRLjrf9pFIuRIjY2I50ccy6alzekzWPhaIgNWEkDxJP4w3ofMUonm~VCBJbl4BCRtKr4b7Qo24zrbDNyTauUVppNXUugDr-QzyJR2In0xmYlmw3TURtRigzN-2GtzonhRDX2eMoJVZ9ybfLvaSCUYtbmEbGhKGdjC1h4pMQ4l9pXWnASyHuAr~bsNggIOOuWOV5fHPCSCjmn9lr6FoI37dTjt4bbe8I4EUkigk7lAuHuaeHJeIoqVsRAEV0uSlWwc6nWkyC9Yh46dn0ecFSGHOZWebH8HSMSFaHMJ7MTjZ4D9sOYt6s09D6j49Z0WXOg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
图1
癫痫数据的1、2和3级残差的正常曲线图:(a)模型(ν2),1级;(b) 型号(ν2),2级;(c) 型号(ν2),3级;(d) 型号(),1级;(e) 型号(),2级;(f) 型号(),级别3;(g) 型号(),1级;(h) 型号(),2级;(i) 型号(),级别3
预计为一小组离散数据。事实上,适度偏离常态可能并不意味着偏离模型假设。我们将在第节中进一步讨论这个问题7.
回归和离散参数估计值如表所示1.型号(ν2=0),型号(ν2)和型号()以及它们的回归参数估计值和相应的标准误差与以前的研究(Thall和Vail,1990; 布雷斯洛和克莱顿,1993; Lee和Nelder,1996). 处理和基线之间的相互作用保留在这三个模型中,因为近似对-此交互作用的值略大于0.05。作为Thall和Vail(1990)这意味着治疗组的预测平均癫痫发作率高于或低于安慰剂组,相应地,基线计数是否超过临界阈值(减去Trt回归参数与交互作用项回归系数的比值)。因此,对于高癫痫发作率的患者,这种新药程序显然是禁忌的。应谨慎对待这一建议的禁忌症,因为它是基于受试者之间分散参数保持不变的模型,即边际方差仅以系统的方式随边际平均值变化。
参数. | 以下模型的结果:. |
---|
. | GLM公司†. | 模型(ν2=0). | 模型(ν2). | 模型. | 模型. |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | −2.80 | 0.41 | −1.35 | 1.22 | −1.35 | 1.22 | −1.37 | 1.16 | −1.11 | 0.72 |
底座 | 0.95 | 0.04 | 0.88 | 0.14 | 0.88 | 0.14 | 0.87 | 0.13 | 0.94 | 0.07 |
Trt公司 | −1.34 | 0.16 | −0.89 | 0.41 | −0.89 | 0.42 | −0.91 | 0.40 | −0.49 | 0.28 |
基础。Trt公司 | 0.56 | 0.06 | 0.34 | 0.21 | 0.34 | 0.21 | 0.35 | 0.20 | 0.03 | 0.13 |
年龄 | 0.90 | 0.12 | 0.51 | 0.36 | 0.51 | 0.36 | 0.52 | 0.34 | 0.36 | 0.21 |
访问 | −0.08 | 0.10 | −0.22 | 0.10 | −0.22 | 0.22 | −0.28 | 0.23 | −0.35 | 0.18 |
σ2 | | | 0.28 | | 0.19 | | 0.16 | | 0.007 | |
| | | | | 0.36 | | 0.33 | | 0.03 | |
| | | | | | | 0.32 | | 0.08 | |
| | | | | | | 0.67 | | 0.54 | |
| | | | | | | 0.25 | | 12.9 | |
参数. | 以下模型的结果:. |
---|
. | GLM公司†. | 模型(ν2=0). | 模型(ν2). | 模型. | 模型. |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | −2.80 | 0.41 | −1.35 | 1.22 | −1.35 | 1.22 | −1.37 | 1.16 | −1.11 | 0.72 |
底座 | 0.95 | 0.04 | 0.88 | 0.14 | 0.88 | 0.14 | 0.87 | 0.13 | 0.94 | 0.07 |
晶体管 | −1.34 | 0.16 | −0.89 | 0.41 | −0.89 | 0.42 | −0.91 | 0.40 | −0.49 | 0.28 |
基础。晶体管 | 0.56 | 0.06 | 0.34 | 0.21 | 0.34 | 0.21 | 0.35 | 0.20 | 0.03 | 0.13 |
年龄 | 0.90 | 0.12 | 0.51 | 0.36 | 0.51 | 0.36 | 0.52 | 0.34 | 0.36 | 0.21 |
访问 | −0.08 | 0.10 | −0.22 | 0.10 | −0.22 | 0.22 | −0.28 | 0.23 | −0.35 | 0.18 |
σ2 | | | 0.28 | | 0.19 | | 0.16 | | 0.007 | |
| | | | | 0.36 | | 0.33 | | 0.03 | |
| | | | | | | 0.32 | | 0.08 | |
| | | | | | | 0.67 | | 0.54 | |
| | | | | | | 0.25 | | 12.9 | |
参数. | 以下模型的结果:. |
---|
. | GLM公司†. | 模型(ν2=0). | 模型(ν2). | 模型. | 模型. |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | −2.80 | 0.41 | −1.35 | 1.22 | −1.35 | 1.22 | −1.37 | 1.16 | −1.11 | 0.72 |
底座 | 0.95 | 0.04 | 0.88 | 0.14 | 0.88 | 0.14 | 0.87 | 0.13 | 0.94 | 0.07 |
Trt公司 | −1.34 | 0.16 | −0.89 | 0.41 | −0.89 | 0.42 | −0.91 | 0.40 | −0.49 | 0.28 |
基础。Trt公司 | 0.56 | 0.06 | 0.34 | 0.21 | 0.34 | 0.21 | 0.35 | 0.20 | 0.03 | 0.13 |
年龄 | 0.90 | 0.12 | 0.51 | 0.36 | 0.51 | 0.36 | 0.52 | 0.34 | 0.36 | 0.21 |
访问 | −0.08 | 0.10 | −0.22 | 0.10 | −0.22 | 0.22 | −0.28 | 0.23 | −0.35 | 0.18 |
σ2 | | | 0.28 | | 0.19 | | 0.16 | | 0.007 | |
| | | | | 0.36 | | 0.33 | | 0.03 | |
| | | | | | | 0.32 | | 0.08 | |
| | | | | | | 0.67 | | 0.54 | |
| | | | | | | 0.25 | | 12.9 | |
参数. | 以下模型的结果:. |
---|
. | GLM公司†. | 模型(ν2=0). | 模型(ν2). | 模型. | 模型. |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | −2.80 | 0.41 | −1.35 | 1.22 | −1.35 | 1.22 | −1.37 | 1.16 | −1.11 | 0.72 |
底座 | 0.95 | 0.04 | 0.88 | 0.14 | 0.88 | 0.14 | 0.87 | 0.13 | 0.94 | 0.07 |
Trt公司 | −1.34 | 0.16 | −0.89 | 0.41 | −0.89 | 0.42 | −0.91 | 0.40 | −0.49 | 0.28 |
基础。Trt公司 | 0.56 | 0.06 | 0.34 | 0.21 | 0.34 | 0.21 | 0.35 | 0.20 | 0.03 | 0.13 |
年龄 | 0.90 | 0.12 | 0.51 | 0.36 | 0.51 | 0.36 | 0.52 | 0.34 | 0.36 | 0.21 |
访问 | −0.08 | 0.10 | −0.22 | 0.10 | −0.22 | 0.22 | −0.28 | 0.23 | −0.35 | 0.18 |
σ2 | | | 0.28 | | 0.19 | | 0.16 | | 0.007 | |
| | | | | 0.36 | | 0.33 | | 0.03 | |
| | | | | | | 0.32 | | 0.08 | |
| | | | | | | 0.67 | | 0.54 | |
| | | | | | | 0.25 | | 12.9 | |
对数据的初步分析显示,患者之间存在大量不规则的边缘异质性。考虑到模型中患者的受试者特定分散参数()治疗与基线之间的相互作用变得明显不显著,治疗效果仍具有统计学意义。因此,在适当考虑了患者之间的边缘异质性后,新药普罗格贝对高癫痫发作率患者的明显禁忌症消失了。这个模型似乎有效地捕捉到了这种不规则的异质性。最小值、第一四分位数、第三四分位数和最大值表中显示了1标签为对于j个=1,2,3,4. 估计s范围为0.03~12.9,分别对应于213例和225例患者。型号()在检测异常值方面似乎也非常敏感。的散点图图中的s。2确定135、227、112、207和225名这些患者符合Thall和Vail报告的异常值(1990)布雷斯洛和克莱顿(1993)他们通过一些绘图来识别这些异常值。
6.2. 连续数据
科克伦和考克斯(1957)介绍了一项巧克力蛋糕烘焙实验的数据。测试了三个配方,每个配方重复15次,总共得到45批蛋糕混合物。每批蛋糕分为六块,在六个不同的温度下烘焙。在烘烤之后,测量断裂角作为实验的响应。费斯和哈里斯(1991)发现基于方差分析的残差图揭示了批次间的异质性。这也可以从批次响应的方框图中看出(Ma,1999).
科克伦和考克斯(1957)利用方差分析对蛋糕烘焙数据进行分析。方差分析模型的残差法向图显示曲率(Ma,1999). 费斯和哈里斯(1991)通过使用乘法模型重新分析蛋糕烘焙数据,他们发现变异系数恒定的假设得到了有力支持。因此,我们根据伽马混合模型重新分析蛋糕烘焙数据(ν2),其中随机效应的两个级别分别为批处理和观察级别。温度和配方被视为因素。我们根据Wald测试筛选这些因素。我们发现,配方效应以及温度和配方之间的任何相互作用似乎可以忽略不计,但温度效应具有高度统计意义。事实上,根据温度绘制的断裂角方框图显示,随着温度的升高,断裂角呈增加趋势(Ma,1999).
科克伦和考克斯(1957)还有费斯和哈里斯(1991)得出了类似的结论,但他们发现配方效果在0.05的统计显著性水平上显著。簇、子簇和观测水平残差的正态图除了簇水平残差(Ma,1999). 为了进行系统比较,我们在表中显示了方差、伽马和伽马混合模型分析的估计值和相应的标准误差,无交互作用2方差分析模型的参数估计值与标准广义线性模型和伽马混合模型的估计值有很大不同;然而,对于这种平衡设计,最后两个模型的回归参数再次几乎相同。在考虑了随机影响后,配方因素的统计显著性降低,而温度因素的统计明显性提高。
参数. | 以下模型的结果:. |
---|
. | 方差分析. | 广义线性. | 模型(ν2). |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | 32.122 | 0.474 | 3.466 | 0.015 | 3.466 | 0.030 |
R1级 | −0.739 | 0.581 | −0.023 | 0.018 | −0.023 | 0.037 |
R2级 | −0.261 | 0.335 | −0.008 | 0.010 | −0.008 | 0.021 |
T1类 | 0.989 | 0.821 | 0.034 | 0.026 | 0.034 | 0.014 |
T2段 | 0.819 | 0.474 | 0.028 | 0.015 | 0.028 | 0.008 |
T3航站楼 | 0.598 | 0.335 | 0.020 | 0.010 | 0.020 | 0.006 |
T4类 | 1.092 | 0.260 | 0.033 | 0.008 | 0.033 | 0.005 |
T5类 | 0.647 | 0.212 | 0.020 | 0.007 | 0.020 | 0.004 |
σ2 | | | | | 0.038 | |
ν2 | | | | | 0.019 | |
ρ2 | 0.030 | | 0.059 | | 0.00025 | |
参数. | 以下模型的结果:. |
---|
. | 方差分析. | 广义线性. | 模型(ν2). |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | 32.122 | 0.474 | 3.466 | 0.015 | 3.466 | 0.030 |
R1级 | −0.739 | 0.581 | −0.023 | 0.018 | −0.023 | 0.037 |
R2级 | −0.261 | 0.335 | −0.008 | 0.010 | −0.008 | 0.021 |
T1类 | 0.989 | 0.821 | 0.034 | 0.026 | 0.034 | 0.014 |
T2段 | 0.819 | 0.474 | 0.028 | 0.015 | 0.028 | 0.008 |
T3航站楼 | 0.598 | 0.335 | 0.020 | 0.010 | 0.020 | 0.006 |
T4类 | 1.092 | 0.260 | 0.033 | 0.008 | 0.033 | 0.005 |
T5类 | 0.647 | 0.212 | 0.020 | 0.007 | 0.020 | 0.004 |
σ2 | | | | | 0.038 | |
ν2 | | | | | 0.019 | |
ρ2 | 0.030 | | 0.059 | | 0.00025 | |
参数. | 以下模型的结果:. |
---|
. | 方差分析. | 广义线性. | 模型(ν2). |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | 32.122 | 0.474 | 3.466 | 0.015 | 3.466 | 0.030 |
R1级 | −0.739 | 0.581 | −0.023 | 0.018 | −0.023 | 0.037 |
R2级 | −0.261 | 0.335 | −0.008 | 0.010 | −0.008 | 0.021 |
T1类 | 0.989 | 0.821 | 0.034 | 0.026 | 0.034 | 0.014 |
T2段 | 0.819 | 0.474 | 0.028 | 0.015 | 0.028 | 0.008 |
T3航站楼 | 0.598 | 0.335 | 0.020 | 0.010 | 0.020 | 0.006 |
T4类 | 1.092 | 0.260 | 0.033 | 0.008 | 0.033 | 0.005 |
T5类 | 0.647 | 0.212 | 0.020 | 0.007 | 0.020 | 0.004 |
σ2 | | | | | 0.038 | |
ν2 | | | | | 0.019 | |
ρ2 | 0.030 | | 0.059 | | 0.00025 | |
参数. | 以下模型的结果:. |
---|
. | 方差分析. | 广义线性. | 模型(ν2). |
---|
. | 估算. | 标准误差. | 估算. | 标准误差. | 估算. | 标准误差. |
---|
常量 | 32.122 | 0.474 | 3.466 | 0.015 | 3.466 | 0.030 |
R1级 | −0.739 | 0.581 | −0.023 | 0.018 | −0.023 | 0.037 |
R2级 | −0.261 | 0.335 | −0.008 | 0.010 | −0.008 | 0.021 |
T1类 | 0.989 | 0.821 | 0.034 | 0.026 | 0.034 | 0.014 |
T2段 | 0.819 | 0.474 | 0.028 | 0.015 | 0.028 | 0.008 |
T3航站楼 | 0.598 | 0.335 | 0.020 | 0.010 | 0.020 | 0.006 |
T4类 | 1.092 | 0.260 | 0.033 | 0.008 | 0.033 | 0.005 |
T5类 | 0.647 | 0.212 | 0.020 | 0.007 | 0.020 | 0.004 |
σ2 | | | | | 0.038 | |
ν2 | | | | | 0.019 | |
ρ2 | 0.030 | | 0.059 | | 0.00025 | |
7.仿真
在本节中,我们在仿真研究的基础上评估BLUP方法的性能。我们将重点关注泊松混合模型,因为这些模型在我们方法中的计数、二项式和生存数据分析中发挥着重要作用(Ma等人。, 2003). 为了研究BLUP方法在实际条件下的性能,我们通过使用泊松模型和在假设1-3下生成的嵌套伽马随机效应,通过模拟将癫痫数据集“复制”200次。我们取协变量Constant、Base、Trt、Age、Visit和Base。癫痫数据的Trt作为模拟的协变量。表中列出了真实模型参数三作为“真值”β0,β1,…,β5,σ2和ν2.
. | . | 以下参数和真值的结果:. | . |
---|
. | . | . | β0,−1.3. | β1, 0.88. | β2,−0.88. | β三, 0.50. | β4,−0.23. | β5, 0.34. | σ2, 0.24. | ν2, 0.44. |
---|
偏差† | GLMм | | −0.08 | −0.01 | 0.06 | 0.02 | 0.01 | −0.03 | | |
| BLUP公司 | | −0.10 | −0.01 | −0.02 | 0.02 | 0 | 0.01 | 0.01 | −0.02 |
标准误差 | GLM公司 | 模拟§ | 1.66 | 0.18 | 0.55 | 0.48 | 0.27 | 0.28 | | |
| | 估计§§ | 0.41 | 0.04 | 0.16 | 0.12 | 0.10 | 0.06 | | |
| BLUP公司 | 模拟 | 1.45 | 0.14 | 0.40 | 0.42 | 0.21 | 0.21 | 0.07 | 0.07 |
| | 估计 | 1.20 | 0.14 | 0.40 | 0.35 | 0.22 | 0.21 | | |
| | | 以下参数和随机效应的结果: |
| | | β0 | β1 | β2 | β三 | β4 | β5 | u个1 | u个11 |
覆盖范围计数* | GLM公司 | 模拟 | 190 | 190 | 190 | 191 | 191 | 188 | | |
| | 估计 | 75 | 71 | 81 | 75 | 96 | 61 | | |
| BLUP公司 | 模拟 | 190 | 192 | 192 | 191 | 190 | 192 | 183 | 199 |
| | 估计 | 175 | 188 | 191 | 176 | 194 | 188 | 169 | 189 |
. | . | 以下参数和真值的结果:. | . |
---|
. | . | . | β0,−1.3. | β1, 0.88. | β2,−0.88. | β三, 0.50. | β4,−0.23. | β5, 0.34. | σ2, 0.24. | ν2, 0.44. |
---|
偏差† | GLMм | | −0.08 | −0.01 | 0.06 | 0.02 | 0.01 | −0.03 | | |
| BLUP公司 | | −0.10 | −0.01 | −0.02 | 0.02 | 0 | 0.01 | 0.01 | −0.02 |
标准误差 | GLM公司 | 模拟§ | 1.66 | 0.18 | 0.55 | 0.48 | 0.27 | 0.28 | | |
| | 估计§§ | 0.41 | 0.04 | 0.16 | 0.12 | 0.10 | 0.06 | | |
| BLUP公司 | 模拟 | 1.45 | 0.14 | 0.40 | 0.42 | 0.21 | 0.21 | 0.07 | 0.07 |
| | 估计 | 1.20 | 0.14 | 0.40 | 0.35 | 0.22 | 0.21 | | |
| | | 以下参数和随机效应的结果: |
| | | β0 | β1 | β2 | β三 | β4 | β5 | u个1 | u个11 |
覆盖范围计数* | GLM公司 | 模拟 | 190 | 190 | 190 | 191 | 191 | 188 | | |
| | 估计 | 75 | 71 | 81 | 75 | 96 | 61 | | |
| BLUP公司 | 模拟 | 190 | 192 | 192 | 191 | 190 | 192 | 183 | 199 |
| | 估计 | 175 | 188 | 191 | 176 | 194 | 188 | 169 | 189 |
. | . | 以下参数和真值的结果:. | . |
---|
. | . | . | β0,−1.3. | β1, 0.88. | β2,−0.88. | β三, 0.50. | β4,−0.23. | β5, 0.34. | σ2, 0.24. | ν2, 0.44. |
---|
偏差† | GLMм | | −0.08 | −0.01 | 0.06 | 0.02 | 0.01 | −0.03 | | |
| BLUP公司 | | −0.10 | −0.01 | −0.02 | 0.02 | 0 | 0.01 | 0.01 | −0.02 |
标准误差 | GLM公司 | 模拟§ | 1.66 | 0.18 | 0.55 | 0.48 | 0.27 | 0.28 | | |
| | 估计§§ | 0.41 | 0.04 | 0.16 | 0.12 | 0.10 | 0.06 | | |
| BLUP公司 | 模拟 | 1.45 | 0.14 | 0.40 | 0.42 | 0.21 | 0.21 | 0.07 | 0.07 |
| | 估计 | 1.20 | 0.14 | 0.40 | 0.35 | 0.22 | 0.21 | | |
| | | 以下参数和随机效应的结果: |
| | | β0 | β1 | β2 | β三 | β4 | β5 | u个1 | u个11 |
覆盖范围计数* | GLM公司 | 模拟 | 190 | 190 | 190 | 191 | 191 | 188 | | |
| | 估计 | 75 | 71 | 81 | 75 | 96 | 61 | | |
| BLUP公司 | 模拟 | 190 | 192 | 192 | 191 | 190 | 192 | 183 | 199 |
| | 估计 | 175 | 188 | 191 | 176 | 194 | 188 | 169 | 189 |
. | . | 以下参数和真值的结果:. | . |
---|
. | . | . | β0,−1.3. | β1, 0.88. | β2,−0.88. | β三, 0.50. | β4,−0.23. | β5, 0.34. | σ2, 0.24. | ν2, 0.44. |
---|
偏差† | GLMм | | −0.08 | −0.01 | 0.06 | 0.02 | 0.01 | −0.03 | | |
| BLUP公司 | | −0.10 | −0.01 | −0.02 | 0.02 | 0 | 0.01 | 0.01 | −0.02 |
标准误差 | GLM公司 | 模拟§ | 1.66 | 0.18 | 0.55 | 0.48 | 0.27 | 0.28 | | |
| | 估计§§ | 0.41 | 0.04 | 0.16 | 0.12 | 0.10 | 0.06 | | |
| 虚张声势 | 模拟 | 1.45 | 0.14 | 0.40 | 0.42 | 0.21 | 0.21 | 0.07 | 0.07 |
| | 估计 | 1.20 | 0.14 | 0.40 | 0.35 | 0.22 | 0.21 | | |
| | | 以下参数和随机效应的结果: |
| | | β0 | β1 | β2 | β三 | β4 | β5 | u个1 | u个11 |
覆盖范围计数* | GLM公司 | 模拟 | 190 | 190 | 190 | 191 | 191 | 188 | | |
| | 估计 | 75 | 71 | 81 | 75 | 96 | 61 | | |
| BLUP公司 | 模拟 | 190 | 192 | 192 | 191 | 190 | 192 | 183 | 199 |
| | 估计 | 175 | 188 | 191 | 176 | 194 | 188 | 169 | 189 |
表中显示了我们模拟的汇总统计数据三对于广义线性模型和BLUP方法,回归参数的偏差相对较小。然而,色散参数σ2和ν2分别被稍微高估和低估。这与Breslow和Clayton报告的模拟研究不同(1993)色散参数总是被低估。
200多个模拟估计值的样本标准误差和200个估计标准误差的平均值分别称为模拟标准误差和估计标准误差(Lin和Breslow,1996)。显然,广义线性模型标准误差估计值具有严重的负偏差。相反,传统的BLUP方法略微低估了年龄回归参数的标准误差,而截距的标准误差更大。由于偏差的平方可以忽略不计,所以均方误差几乎与模拟标准误差的平方相同;因此,表中省略了均方误差三.
在模拟和估计的标准误差的基础上,我们构建了参数的模拟和估计95%置信区间,以及随机效应的95%预测区间(u个1和u个11)假设估计量和预测量的正态性。我们在表中列出三真实值被95%置信区间或预测区间覆盖的次数计数。正如预期的那样,基于广义线性模型估计的95%置信区间表现非常差。相反,对于Base、Trt、Visit和Base,基于正统BLUP方法的回归参数估计95%置信区间的覆盖率约为95%。Trt,但截距和年龄仅为87.5%左右。95%预测区间的覆盖率低至84.5%u个1高达94.5%u个11.
为了检查参数估计值的正态性,我们对模拟中的每个回归和离散参数估计值进行了正态绘图;见马(1999). 即使对于色散参数,偏离正态性的情况似乎也不太严重。相反,即使是模拟数据(Ma,1999). 显然,在评估模型假设时,检查观测数据和模拟数据之间各种曲线图的相似性比检查正态性更有效。
8.讨论
跟随罗宾逊(1991),我们已经表明,BLUP对于预测非正态数据设置中的嵌套随机效应也是一件“好事”,前提是联合得分函数在这两个方面都是线性的Y(Y)和U型,Tweedie混合型号也是如此。该方法仅依赖随机效应的前两个矩,因此它不仅计算效率高,而且对随机效应分布的错误指定具有鲁棒性(Ha和Lee,2005).
Tweedie混合模型同时考虑了聚类数据的分布形状和内部相关性。响应协方差矩阵结构的方差分量分解不仅很容易解释,而且使模型拟合比条件模式方法简单得多。BLUP方法的渐近合理性不依赖于随机效应的正态性。相反,条件模式方法的渐近证明在很大程度上依赖于随机效应的近似正态性或“正确尺度上随机效应的正确变换正态性”(Draper,1996).
我们的BLUP方法得出的回归参数估计方程与通过其他边际建模方法得出的“人口平均”平均推断的回归参数估算方程相同;然而,BLUP方法明确地纳入了随机效应,并通过使用相同的模型合并了“特定主题”和“人口平均”推论(Lee和Nelder,1996). 我们的gamma混合模型为广义估计方程提供了一个参数示例,其中相关结构不依赖于回归参数,而我们模型的边际协方差结构通常与其他方法考虑的结构不同。
在开发随机效应建模方法时,模型检查通常被忽略(参见Zeger和Karim等(1991)布雷斯洛和克莱顿(1993)和Booth等人。(2003))或通过残差的正态性检查(Lee和Nelder,1996); 然而,在文献中还没有发现残差正态性的正当性。从我们的模拟研究中可以看出,适度偏离常态并不一定表示偏离模型假设。
Tweedie广义线性混合模型的正统BLUP方法产生了一种灵活高效的拟合算法,能够处理大量数据集,如Ma在生存分析中的应用所示等人。(2003). 在这里,我们通过分析聚集计数和连续数据来说明该方法。此外,它还可以应用于正概率分量为零的正连续数据,例如汽车保险数据(Jörgensen和Souza,1994; 斯迈思和约根森,2002)或客户支出数据。该方法对二项式混合数据(Ma,1999)以及空间和时空随机效应广义线性模型的适应性将在其他地方讨论。
致谢
这项工作得到了加拿大自然科学与工程研究委员会的资助。作者非常感谢John Petkau博士和Charles McCulloch博士对我们的工作提出的建设性意见,这些意见大大改进了这些结果的表述。作者还感谢联合编辑和裁判的宝贵意见。
工具书类
附录A:协方差矩阵的逆
我们推导了响应协方差矩阵的逆,因为该逆在BLUP方法推导的所有阶段都起着关键作用。这里提出的反演公式普遍适用,而不仅仅适用于BLUP方法。
由于该协方差矩阵相对于簇是块对角的,因此可以反转这些块。让响应向量对应于我第个簇表示为哪里Y(Y)ij公司=(Y(Y)ij公司1,…,Y(Y)ijn公司ij公司)T型响应协方差矩阵的显式表达式可以通过注意该协方差矩阵模式化为
哪里A类ij公司=ρ2诊断{E类(Y(Y)ij公司)对}+ν2(Y(Y)ij公司) (Y(Y)ij公司)T型首先,var的逆(Y(Y)我)可以用通过反转var(Y(Y)我)通过以下著名的谢尔曼-莫里森-伍德伯里公式
哪里一和b条是向量和A类是可逆的。将此公式再次应用于A类ij公司,我们有
哪里
现在可以递归地获得这个逆的显式表达式。
©2007皇家统计学会