联合模型(Joint Model)中带随机效应协变量的二元时依Copula生存数据与广义线性混合效应纵向数据的极大似然推断：基于MCEM(Monte Carlo Expectation-Maximization)-Metropolis-Hastings-ECME算法

《Computational Statistics & Data Analysis》：Generalized linear mixed-effects joint model for longitudinal and bivariate survival data

　　

本文提出一种基于Monte Carlo期望最大化(Monte Carlo Expectation-Maximization, MCEM)与期望条件最大化(Estimation/Conditional Maximization Either, ECME)及Metropolis-Hastings抽样相结合的极大似然推断框架，用于拟合带随机效应协变量的二元时依copula Clayton生存数据与服从指数族分布的广义线性混合效应模型(Generalized Linear Mixed-effects Model, GLMM)纵向数据的联合模型(Joint Model, JM)。研究人员考虑纵向子模型含可交换随机截距和/或随机斜率bi～Nq(0,Σb)，生存子模型含与bi相关联的时间变协变量及基线危险率经离散化处理的非参数估计，两结局间依赖结构由Clayton copula刻画。由于潜随机效应bi的条件后验分布无解析形式，E步采用Metropolis-Hastings算法从f(bi|yi,Oi,ψ(r))抽取样本并经由重要性加权蒙特卡洛积分近似Q(ψ|ψ(r))；M步采用ECME算法分CM步分别更新固定效应α、离散步长参数γ与βk、离散基线危险跳度λkj及随机效应协方差Σb，其中Σb有闭式解而α、?通过Newton-Raphson迭代，基线危险与copula关联参数ρ经轮廓(profile)似然求解。模拟研究表明该估计量具渐近无偏性、较小均方误及接近名义水平的95%置信区间覆盖率，应用于一肿瘤伴二次复发临床数据集展示其实际效用。标准误由B=200次参数自助法(parametric bootstrap)获得。

论文解读：带随机效应协变量的二元时依Copula生存数据与GLMM纵向数据联合模型的MCEM-MH-ECME极大似然推断

一、研究背景与动机

在随访研究中，受试者常同时记录重复测量的纵向指标（如生化标志物）与多个时间至事件结局（如疾病进展与死亡）。传统做法分别拟合纵向混合效应模型与Cox比例风险模型，忽略了纵向过程与失效时间间的潜在相关性及两失效事件间的相依性，导致参数估计有偏且丧失预测效率。已有联合模型(Joint Model, JM)多针对单生存结局及简单随机截距，且较少同时处理二元或多重时依相依生存结局与广义线性混合效应纵向子模型（含泊松或二项响应）。此外当潜随机效应b_i的后验分布无闭式表达时，标准EM算法难以直接实施。为此研究人员开展本研究，构建含Clayton时依copula二元相依生存子模型与广义线性混合效应模型(Generalized Linear Mixed-effects Model, GLMM)纵向子模型的联合似然框架，设计基于Monte Carlo期望最大化(Monte Carlo Expectation-Maximization, MCEM)、Metropolis-Hastings(MH)抽样及期望条件最大化(Estimation/Conditional Maximization Either, ECME)的数值推断流程以实现极大似然估计(Maximum Likelihood Estimation, MLE)，并验证其有限样本性质。该文发表于《Computational Statistics》。

二、主要关键技术方法

研究人员采用如下技术路线：(1)构造完整数据对数似然?_c(ψ|D,b)=A₁(α,?|b)+A₂(ξ|b)+A₃(Σ_b|b)，其中A₁来自GLMM纵向部分含线性预测η_ij=X_ij^?α+Z_ij^?b_i，A₂来自二元时依Clayton copula生存部分含离散基线累积危险Λ_k(t)跳度λ_kj，A₃为随机效应正态先验；(2)E步：因E{b_i|O_i,ψ^(r)}无解析形式，用对称随机游走Metropolis-Hastings从f(b_i|y_i,O_i,ψ^(r))抽取M个样本{b_i^(m|r)}，以重要性加权蒙特卡洛平均近似Q₁、Q₂、Q₃；(3)M步：用ECME分三条件极大化(CM)步——CM-1用Newton-Raphson更新固定效应α与散度参数?；CM-2对离散基线跳度λ_1j、λ_2j'求偏导置零得闭式剖面估计，同步更新Cox型回归系数γ₁₁,γ₁₂,β₁及γ₂₁,γ₂₂,β₂；CM-3给出Σ_b^(r+1)=n^-1Σ_iE(b_ib_i^?|O_i;ψ^(r))，由MC样本近似；copula参数ρ由轮廓似然E{?l_ρ/?ρ|O;ψ^(r)}=0数值求解；(4)初值由lme4与coxph提供，收敛判据为参数相对变化<10^-4；(5)标准误经B=200次参数自助法(bootstrap)计算95% CI；(6)模拟数据按设定真值生成含n=200,300，删失率(Censoring Rate, CR)=0.3,0.4，Clayton copula Kendall's τ对应ρ=0.3,0.5,0.8，纵向数据取Gaussian(?=0.5)或Poisson(a(?)=1)，随机效应维数q=1或2；(7)实证数据源自某肿瘤患者队列含初次进展与二次复发两时依结局及纵向肿瘤标志物。

三、研究结果

3.2 Implementation of E-step（E步实现）

由于潜随机效应b_i的条件分布无闭型且含指数族项，研究人员采用随机游走Metropolis-Hastings算法以N_q(b_i^(m-1|r),σ_b²Ω_b)为提议分布从f(b_i|y_i,O_i,ψ^(r))抽取样本，接受概率取min{1,f(b_i^*|·)/f(b_i^(m-1|r)|·)}，Ω_b为Fisher信息逆阵加Σ_b^-1，调σ_b²使接受率≈0.25。E{Q₁}中需E(θ_i|O_i,ψ^(r))与E{B(θ_i)|O_i,ψ^(r))，E{Q₃}需E(b_ib_i^?|O_i,ψ^(r))，均由MH样本经重要性加权近似。

3.3 Implementation of M-step（M步实现）

M步沿ECME拆分三部分：Q₃关于Σ_b有显式解Σ_b^(r+1)=(1/n)Σ_iE(b_ib_i^?|O_i;ψ^(r))，MC近似如上。Q₁对(α,?)无闭解，用Newton-Raphson迭代α^(r+1)=α^(r)-[?²Q₁/?α?α^?]^-1·(?Q₁/?α)及类似对?更新。Q₂涉及无限维基线Λ₁(·)、Λ₂(·)，将其约束为仅在观测事件时刻有正跳跃的阶梯函数，导出l_T?1与l_T?2表达式后先对λ_1l令?l_T?1/?λ_1l=0求剖面估计再Newton更新(γ₁₁,γ₁₂,β₁)；同理处理第二生存结局。给定上述估计后代入l_ρ解E{?l_ρ/?ρ|O;ψ^(r)}=0得ρ的估计。

4 Theoretical properties of estimators（估计量的理论性质）

在补充材料所列正则条件下（含参数空间紧性、似然可微、潜变量条件密度有界等），证明由此得到的非参数极大似然估计量(Nonparametric Maximum Likelihood Estimator, NPMLE) ψ?_n在参数全空间上一致相合(Consistent)且√n(ψ?_n-ψ₀)依分布收敛于零均值多元正态分布（渐近正态 Asympotically Normal），其方差可用观察信息矩阵或自助法估计。

5 Simulation studies（模拟研究）

设置Gaussian纵向(q=1,2)与Poisson纵向(q=1,2)，Clayton copula ρ对应Kendall's τ=0.3,0.5,0.8，n=200/300，CR=0.3/0.4。500次重复显示：各参数偏倚(Bias)随n增大趋近0，样本标准差(SSD)与平均估计标准误(ESE)吻合良好，95%经验覆盖率(Coverage Probability, CP)接近名义水平0.95；随机效应维数升高与copula相依增强略增SSD但不破坏一致性；误指定copula族（如用Gaussian copula拟合Clayton生成数据）致ρ及相关联生存参数偏倚增大，证实copula选取敏感性。QQ图显示标准化估计量近似N(0,1)。

6 Real-data analysis（实证数据分析）

将该方法应用于某肿瘤患者队列(n=XXX)，含纵向肿瘤标志物(GLMM拟合)、首次进展与二次复发两时依结局(Clayton copula联结)。结果给出固定效应、随机效应协方差、copula关联参数ρ?及其95%CI，展示方法在实际生物医学联合建模中的可行性。

四、讨论与结论

研究人员成功构建了适用于二元时依copula生存数据与广义线性混合效应纵向数据的联合模型极大似然推断框架，解决了潜随机效应后验无闭式时的EM难计算问题。核心贡献为：(1)引入Metropolis-Hastings近似E步条件期望避免高维积分；(2)用ECME分离处理有限与无限维参数使M步可操作，基线危险以离散跳度剖面估计；(3)给出Σ_b显式MC近似更新及bootstrap标准误。模拟验证估计量具一致性、渐近正态性与合理区间覆盖率。该方法可推广至满足文中条件(A6)之任意copula族（不只Clayton）。局限含MH抽样计算耗时及copula误指定影响。结论为：所提MCEM-MH-ECME算法可有效、稳健地拟合带随机效应协变量之二元时依Copula生存—GLMM纵向联合模型，为多重终点与纵向生物标志物的综合统计推断提供可靠工具。