Generalized linear mixed joint model for longitudinal and survival outcomes

Abstract

Longitudinal studies often entail categorical outcomes as primary responses. When dropout occurs, non-ignorability is frequently accounted for through shared parameter models (SPMs). In this context, several extensions from Gaussian to non-Gaussian longitudinal processes have been proposed. In this paper, we formulate an approach for non-Gaussian longitudinal outcomes in the framework of joint models. As an extension of SPMs, based on shared latent effects, we assume that the history of the response up to current time may have an influence on the risk of dropout. This history is represented by the current, expected, value of the response. Since the time a subject spends in the study is continuous, we parametrize the dropout process through a proportional hazard model. The resulting model is referred to as Generalized Linear Mixed Joint Model (GLMJM). To estimate model parameters, we adopt a maximum likelihood approach via the EM algorithm. In this context, the maximization of the observed data log-likelihood requires numerical integration over the random effect posterior distribution, which is usually not straightforward; under the assumption of Gaussian random effects, we compare Gauss-Hermite and Pseudo-Adaptive Gaussian quadrature rules. We investigate in a simulation study the behaviour of parameter estimates in the case of Poisson and Binomial longitudinal responses, and apply the GLMJM to a benchmark dataset.

Appendices

Appendix A: Score vectors and Hessian matrices for the Poisson and the Binomial joint model

As highlighted in Sect. 2, for parameter estimation in the GLMJMs, the EM algorithm requires the computation of the observed data score vector and Hessian matrix for the complete joint model. The score can assume the form:

$$ \mathcal{S}({\boldsymbol{\theta}})= \sum _i \int\omega({\boldsymbol{\theta}},\mathbf{b}_i)p(\mathbf{b}_i \mid T_i, \delta_i,y_i;{\boldsymbol{\theta}})d\mathbf{b}_i, $$

(12)

where ω(θ,b _i )=∂{logp(T _i ,δ _i |b _i ;θ)+logp(y _i |b _i ;θ)+logp(b _i ;θ)}/∂θ ^⊤ denotes the complete data score vector. Therefore, it can be considered as the expected value of the complete data score vector with respect to the posterior distribution of the random effects. In the M-step, at the q-th iteration, equation (12) is solved with respect to θ, with p(b _i |T _i ,δ _i ,y _i ;θ) calculated by fixing θ=θ ^(q−1).

Since this work has focused on the longitudinal parameter estimation under non-ignorability, in this section we report the score vector and the Hessian matrices for the specific longitudinal processes considered (Poisson and Binomial distributions).

For the Poisson longitudinal distribution, the score vector for the fixed effects in the longitudinal model can be written as

while the Hessian matrix assumes the form

For the Binomial case, on the other hand, the score vector and the Hessian matrix are given by, respectively:

and

where

$$\mathbf{A}_i = \frac{ [ 2 \mathbf{x}_i \exp \{{\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \} ] [ 1 - 2 \mathbf{x}_i \exp \{ {\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \} - 1/2 \mathbf{x}_i \exp \{ {\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \}^2 ]}{ [ 1+ \exp \{ {\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \} ]^4}, $$

while

$$\mathbf{B}_i = \frac{\exp \{{\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \}}{ [1+\exp \{ {\boldsymbol{\beta}}^{\mathsf{T}} \mathbf{x}_i +\mathbf{b}_i^{\mathsf{T}}\mathbf{z}_i \} ]^2}. $$

Appendix B: The pseudo-adaptive Gaussian quadrature

As it has been pointed out in Sect. 3, the computation of the observed data log-likelihood requires the use of quadrature rules, such as Gaussian quadrature. In this Appendix, we describe an alternative numerical integration approach based on the Pseudo adaptive Gauss-Hermite rule, proposed by Rizopoulos (2012) for standard JMs.

The standard Gauss-Hermite rule (GH) approximates the integral in (3) by a weighted sum of integrand evaluations at pre-specified abscissas, see e.g. Press et al. (2007). The quality of the approximation improves with the number of quadrature points. However, as the GH requires the evaluation of the integrand over the Cartesian product of the abscissas corresponding to each dimension in the random effect vector, the computational effort increases exponentially with the random effect dimension. Another limitation of the GH approach is that the corresponding approximation depends on the location of the quadrature points with respect to the main mass of the integrand, that may be far from each other. For instance, the random effect distribution may be far from Gaussianity.

To solve these problems, an adaptive quadrature approach may be used, that appropriately centres and scales the integrand at each iteration of the optimization algorithm. However, this kind of approach implies a substantial computational burden due to the location of the mode and the computation of the Hessian matrix of the random effect posterior distribution at each iteration.

To decrease the computational effort, we propose to use a pseudo-adaptive quadrature rule. This rule is based on the Bayesian Central Limit (BCL) theorem, see Cox and Hinkley (1974); the posterior distribution of the random effects can be approximated by a multivariate Gaussian distribution as the number of repeated measurements increases. The logarithm of this distribution can in fact be written as

(13)

It is clear from (13) that as n _i increases, the predominant term becomes the longitudinal density logp(y _i |b _i ;β). According to the BCL and under general regularity conditions, as n _i →∞, we have

$$p( \mathbf{b}_i; T_i, \delta_i, \mathbf{y}_i; {\boldsymbol{\beta}}) \stackrel {p}{\rightarrow}\mathcal{N}\bigl( \tilde{\mathbf{b}}_i, \tilde{\mathbf{H}}_i^{-1}\bigr), $$

where $\tilde{\mathbf{b}}_{i} = \text{arg\,max}_{b}\{\log p(\mathbf{y}_{i} | \mathbf{b}_{i}; {\boldsymbol{\beta}})\}$ and $\tilde{\mathbf{H}}_{i} = -\partial^{2} \* \log p(\mathbf{y}_{i} | \mathbf{b}_{i}; {\boldsymbol{\beta}}) / \partial\mathbf{b} \partial\mathbf{b}^{\mathsf{T}} |_{\mathbf{b} = \tilde{\mathbf{b}}}$. Practically, as n _i increases, it is sufficient to re-center and re-scale the integrand for each subject using only the information coming from the longitudinal outcome model.

Applying this procedure to the joint model, we propose to first fit a generalized linear mixed model, extract the empirical Bayes estimates $\tilde{\mathbf{b}}_{i} = \text{arg\,max}_{b}\{\log p(\mathbf{y}_{i} | \mathbf{b}_{i}; \tilde{\mathbf{\theta}}_{y})\}$ and the corresponding covariance matrix $\tilde{\mathbf{H}}_{i}^{-1} = (\mathbf{B}_{i}\mathbf{B}_{i}^{T})^{-1}$; successively, we apply the following transformation to compute the score vector:

where q is the number of random effects, $\sum_{t_{1}, \ldots, t_{q}}$ is a shorthand for $\sum_{t_{1}=1}^{K}, \ldots, \sum_{t_{q}=1}^{K}$, K denoting the common number of quadrature points, b _t =(b _t1,…,b _tq ) are the abscissas with corresponding weights π _t , $\tilde {r}_{t} = \tilde{\mathbf{b}}_{i} + \sqrt{2} \mathbf{B}_{i}^{-1} \mathbf{b}_{t}$ are the re-scaled abscissas and $\tilde{\mathbf{\theta}}_{i}$ denotes the maximum likelihood estimates for the longitudinal model.

This procedure is very similar to the GH rule, but it is implemented only once, at the start of the EM algorithm, leading to faster estimation.