Two-part models with stochastic processes for modelling longitudinal semicontinuous data: Computationally efficient inference and modelling the overall marginal mean

Abstract

Several researchers have described two-part models with patient-specific stochastic processes for analysing longitudinal semicontinuous data. In theory, such models can offer greater flexibility than the standard two-part model with patient-specific random effects. However, in practice, the high dimensional integrations involved in the marginal likelihood (i.e. integrated over the stochastic processes) significantly complicates model fitting. Thus, non-standard computationally intensive procedures based on simulating the marginal likelihood have so far only been proposed. In this paper, we describe an efficient method of implementation by demonstrating how the high dimensional integrations involved in the marginal likelihood can be computed efficiently. Specifically, by using a property of the multivariate normal distribution and the standard marginal cumulative distribution function identity, we transform the marginal likelihood so that the high dimensional integrations are contained in the cumulative distribution function of a multivariate normal distribution, which can then be efficiently evaluated. Hence, maximum likelihood estimation can be used to obtain parameter estimates and asymptotic standard errors (from the observed information matrix) of model parameters. We describe our proposed efficient implementation procedure for the standard two-part model parameterisation and when it is of interest to directly model the overall marginal mean. The methodology is applied on a psoriatic arthritis data set concerning functional disability.

1 Introduction

Semicontinuous data arise when the outcome is a mixture of true zeros and continuously distributed positive values. ¹ Some examples in the literature have included average daily alcohol consumption, ¹ hospital lengths of stay ² and medical expenditures.^3,4 In these situations, and more generally, it is natural to view the outcome as a result of two processes, the first determines if the outcome is zero, and if not the second determines the positive value. Two-part models are therefore convenient for the analysis of semicontinuous data and have been used extensively. Recently, Smith et al.^3,4 considered the interesting notion of reparameterising the mean of the positive values in terms of the overall mean, which is arguably a more justified target of inference (see Tom et al. ⁵ and the references therein). We also consider this notion with respect to the overall marginal mean in our framework.

Two-part marginal models and two-part mixed models have both been proposed for the analysis of longitudinal semicontinuous data. The first is motivated by obtaining population-based inference and have been constructed using generalized estimating equations. ⁶ The second is more convenient when patient-specific inference is of interest and are constructed by incorporating correlated patient-specific random effects in both parts of the model. ⁷ This paper focuses on the two-part mixed modelling approach, although considerations are provided on how population-based inference can be obtained.

In some situations, correlated patient-specific random effects models will not provide an adequate fit to the data. This may especially be the case when the lengths of follow-up are relatively long. Here, it may be less plausible to assume that patients can only have consistently high or low outcomes throughout their entire follow-up. In terms of the correlation structure, it may not be reasonable to assume constant correlation between outcomes from the same patient regardless of their gap times (which is induced by patient-specific random effects). Flexible two-part models that allow for random changes in the trajectory through serially correlated stochastic processes may then be more plausible and these have been proposed in the literature. Albert and Shen ⁸ and Ghosh and Albert ⁹ proposed two-part mixed models that consisted of correlated Gaussian processes and random walks (in addition to correlated patient-specific random effects), respectively, in both parts of the model. Albert and Shen ⁸ demonstrated, through their application and a simulation study, that overall conditional means may suffer from bias if serial correlation (which is not captured by patient-specific random effects) is present but ignored. It is also worth noting, both models incorporating stochastic processes provided considerable improvements of fit to their data.

A main drawback of fitting models with stochastic processes is the computationally intensive nature of the model fitting procedure. The primary difficulty results from the following feature: if a patient has m_i observations, then a model consisting of correlated stochastic processes in each part of the model will require $2m_i$ integrations to evaluate the marginal likelihood contribution from that patient (assuming, as is usual, the stochastic processes are realised at the observation times). For manageable values of m_i, Albert and Shen ⁸ and Ghosh and Albert ⁹ have developed methods based on a Monte Carlo Expectation Maximization algorithm and Markov chain Monte Carlo, respectively, to evaluate the marginal likelihood. Both of these procedures can be computationally intensive, with the former also requiring standard errors of parameter estimates to be computed by bootstrap. The primary aim of this paper is to demonstrate, using a property of the multivariate normal distribution and the standard marginal cumulative distribution function identity, how a marginal likelihood can be obtained in terms of the cumulative distribution function of a multivariate normal distribution. Implicitly, because it is possible to efficiently evaluate the cumulative distribution function of a multivariate normal distribution, maximum likelihood estimation can be used to obtain parameter estimates and (asymptotic) standard errors (from the observed information matrix) of model parameters.

The rest of this paper is organised as follows. In Section 2, the motivating application concerning functional disability in psoriatic arthritis is introduced. Section 3 describes the flexible two-part modelling framework of Albert and Shen ⁸ and Ghosh and Albert ⁹ (including additional comments regarding implementation). Section 4 proposes an efficient maximum likelihood estimation procedure for the models in Section 3. Section 5 applies the methodology in Section 4 to the data described in Section 2. While retaining the flexibility of using stochastic processes models and the practicality of the proposed efficient implementation procedure, Section 6 extends the modelling framework of Section 3 to allow for the direct modelling of the overall marginal mean. Finally, concluding remarks are made in Section 7.

2 Functional disability in psoriatic arthritis

Psoriatic arthritis (PsA) is an inflammatory arthritis associated with the skin condition psoriasis. Because of both skin and joint involvement of the disease, PsA can result in patients having severe physical functional disability. The dominant measure of functional disability in PsA, as well as in many other disease areas, ¹⁰ is the self-reported health assessment questionnaire (HAQ). This produces an essentially continuous measure^11–14 between zero, representing no disability, and three, representing severe disability.

The HAQ scores of 698 patients observed longitudinally at the University of Toronto PsA clinic were considered for this analysis. Figure 1 shows the frequencies of HAQ scores from these patients. From Figure 1, it is evident that a large proportion of zeros exist in this data set (1526/4811 = 0.32). The clumping at zero, together with the continuous distributed outcomes for the non-zero values, suggests that the HAQ score can be viewed as a semicontinuous outcome. Su et al.^12,13 considered two-part models with patient-specific random effects for analysing an earlier version of this PsA data set. In this paper, we relax the assumption of constant patient-specific random effects to patient-specific stochastic processes and consider the extent to which they improve understanding of the disability process. This includes making easy interpretable inference on the overall marginal mean HAQ scores, a concept that has not been considered before with stochastic processes models (see Section 6 for more details). On average, patients had 6.89 clinic visits (ranging from 2 to 20) with mean inter-visit and follow-up times of 1 year and 5 months (standard deviation (SD) of 1 year and 1 month) and 8 years and 3 months (SD of 5 years and 10 months), respectively.

Figure 1.

Frequencies of HAQ scores in our data.

3 Model

Let Y_ij ( $i=1,\dots ,N$ ) denote the semicontinuous response from patient i at time t_ij ( $j=1,\dots ,m_i$ ), where t_ij represents the time of the jth observation from patient i. Because of true zeros, it is natural to decompose the response into

$$U_{ij}=\{\begin{array}{cc}1& :Y_{ij}>0\\ 0& :Y_{ij}=0\end{array}$$

and $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ , where $g(·)$ is a monotonic function such that $g(0)=0$ and $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ is positive and approximately Gaussian with constant variance ${\sigma }^2$ . For convenience, the model for U_ij is referred to as the binary component, while the model for $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ is referred to as the continuous component.

We now describe the flexible modelling framework. Let $X_{\mathrm{ij}}$ and $Z_{\mathrm{ij}}$ be column vectors of covariates that influence the probability of Y_ij > 0 and the mean of $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ , respectively. Then conditional on correlated patient-specific random effects $(B_i^r=b_i^r,C_i^r=c_i^r)$ and correlated stochastic processes $\{B_i^s(t_{\mathrm{ij}})=b_{\mathrm{ij}}^s,C_i^s(t_{\mathrm{ij}})=c_{\mathrm{ij}}^s\}$ , where the random effects are assumed independent of the stochastic processes, we model U_ij as Bernoulli with response probability

$$\mathbb{P}(U_{\mathrm{ij}}=1|b_{\mathrm{ij}}^s,b_i^r)=\Phi (X_{\mathrm{ij}}^{\top }\beta +b_{\mathrm{ij}}^s+b_i^r)$$

(1)

where $\Phi (·)$ is the cumulative distribution function of a standard Gaussian distribution (i.e. probit model), and $[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0;c_{\mathrm{ij}}^s,c_i^r]$ as Gaussian with mean $Z_{\mathrm{ij}}^{\top }\gamma +c_{\mathrm{ij}}^s+c_i^r$ and constant variance ${\sigma }^2$ (i.e. linear mixed effect model on $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ ). Here, $\beta$ and $\gamma$ are column vectors of regression coefficients. The patient-specific random effects $(B_i^r,C_i^r)$ allow patients to have a consistently high or low probability of having disability and a consistently high or low mean for the non-zero HAQ scores across time. While the patient-specific stochastic processes $\{B_i^s(t_{i1}),\dots ,B_i^s(t_{{\mathrm{im}}_i}),C_i^s(t_{i1}),\dots ,C_i^s(t_{{\mathrm{im}}_i})\}$ can capture serial correlation and non-predictable changes in unobserved heterogeneity. ⁹

We assume $\{B_i^r,C_i^r\}$ follows a bivariate normal distribution with mean vector zero and

$$\text{Var}(B_i^r)={\sigma }_b^2,\text{Var}(C_i^r)={\sigma }_c^2,\text{Cov}(B_i^r,C_i^r)={\sigma }_b{\sigma }_c\rho$$

(2)

where ${\sigma }_b^2$ and ${\sigma }_c^2$ are variance parameters and ρ is the correlation between $B_i^r$ and $C_i^r$ . Furthermore, we consider two classes of stochastic processes for $\{B_i^s(t),C_i^s(t)\}$ that are subsequently described. For convenience, define $B_i(t)=B_i^r+B_i^s(t)$ and $C_i(t)=C_i^r+C_i^s(t)$ , i.e. the patient-specific random effects $B_i^r$ and $C_i^r$ are absorbed into the stochastic processes $B_i(t)$ and $C_i(t)$ , respectively, and let the covariance matrix of $\{B_i(t_{i1}),\dots ,B_i(t_{{\mathrm{im}}_i}),C_i(t_{i1}),\dots ,C_i(t_{{\mathrm{im}}_i})\}$ be

$${\Sigma }_i=\left(\begin{array}{cc}{\Sigma }_{ib}& {\Sigma }_{ibc}\\ {\Sigma }_{ibc}& {\Sigma }_{ic}\end{array}\right).$$

3.1 Correlated Gaussian processes

The first and most general model that we consider is defined when $\{B_i^s(t),C_i^s(t)\}$ are correlated stationary Gaussian processes. That is the model proposed by Albert and Shen ⁸

$$\begin{array}{l}\text{Cov}\{B_i^s(t_{ij}),B_i^s(t_{ik})\}={\sigma }_{gb}^2{\rho }_{gb}^{|t_{ij}-t_{ik}|}, \text{Cov}\{C_i^s(t_{ij}),C_i^s(t_{ik})\}={\sigma }_{gc}^2{\rho }_{gc}^{|t_{ij}-t_{ik}|}\\ \text{Cov}\{B_i^s(t_{ij}),C_i^s(t_{ik})\}={\sigma }_{gb}{\sigma }_{gc}{\rho }_g{\rho }_{gbc}^{|t_{ij}-t_{ik}|}\end{array}$$

(3)

where ${\sigma }_{\mathrm{gb}}^2$ and ${\sigma }_{\mathrm{gc}}^2$ are variance parameters, ρ_g is the correlation between the Gaussian processes at the same time point, and ρ_gb, ρ_gc, ρ_gbc are the degradation parameters governing the serial correlation within and between processes, respectively. Following Albert and Shen, ⁸ the processes $B_i(t)$ and $C_i(t)$ are taken to be exchangeable Ornstein-Uhlenbeck (EOU) processes, and the model containing these processes is called the general model, i.e. (1–3). Some special cases of the general model are

• Shared EOU process model when $c_{\mathrm{ij}}=\theta b_{\mathrm{ij}}$ ,

• Correlated OU processes model when ${\sigma }_b^2\equiv {\sigma }_c^2\equiv \rho \equiv 0$ ,

• Shared OU process model when $c_{\mathrm{ij}}^s=\theta b_{\mathrm{ij}}^s$ and ${\sigma }_b^2\equiv {\sigma }_c^2\equiv \rho \equiv 0$ ,

• Correlated random effects model when ${\sigma }_{\mathrm{gb}}^2\equiv {\sigma }_{\mathrm{gc}}^2\equiv {\rho }_g\equiv {\rho }_{\mathrm{gb}}\equiv {\rho }_{\mathrm{gc}}\equiv {\rho }_{\mathrm{gbc}}\equiv 0$ ,

• Shared random effect model when $c_i^r=\theta b_i^r$ and ${\sigma }_{\mathrm{gb}}^2\equiv {\sigma }_{\mathrm{gc}}^2\equiv {\rho }_g\equiv {\rho }_{\mathrm{gb}}\equiv {\rho }_{\mathrm{gc}}\equiv {\rho }_{\mathrm{gbc}}\equiv 0$ ,

where θ is a parameter to be estimated.

3.1.1 Remarks on ρ_gb, ρ_gc and ρ_gbc

Although the general model is very flexible, it will not always be mathematically valid. Let the covariance matrices ${\Sigma }_{\mathrm{ib}}^s,{\Sigma }_{\mathrm{ic}}^s$ and ${\Sigma }_{\mathrm{ibc}}^s$ have (j, k)th entry ${\sigma }_{\mathrm{gb}}^2{\rho }_{\mathrm{gb}}^{|t_{\mathrm{ij}}-t_{\mathrm{ik}}|},{\sigma }_{\mathrm{gc}}^2{\rho }_{\mathrm{gc}}^{|t_{\mathrm{ij}}-t_{\mathrm{ik}}|}$ and ${\sigma }_{\mathrm{gb}}{\sigma }_{\mathrm{gc}}{\rho }_g{\rho }_{\mathrm{gbc}}^{|t_{\mathrm{ij}}-t_{\mathrm{ik}}|}$ respectively, i.e. described by equation (3). If ρ_gb, ρ_gc and ρ_gbc are unconstrained (as specified by Albert and Shen ⁸ ), the matrix ${\Sigma }_i^s$ where

$${\Sigma }_i^s=\left(\begin{array}{cc}{\Sigma }_{ib}^s& {\Sigma }_{ibc}^s\\ {\Sigma }_{ibc}^s& {\Sigma }_{ic}^s\end{array}\right)$$

will not in general be a valid covariance matrix since ${\Sigma }_i^s$ , although symmetric, is not constrained to be positive semi-definite and therefore $B_i^s(t)$ and $C_i^s(t)$ will not necessarily form a jointly Gaussian process. The primary difficulty results when ρ_g (the correlation between $B_i^s(t)$ and $C_i^s(t)$ at each time t) is close to one because the processes $B_i^s(t)$ and $C_i^s(t)$ are similar and therefore it will not be plausible for them to degrade at vastly different rates (i.e. for ρ_gb, ρ_gc and ρ_gbc to be vastly different). A reasonable approximation in this situation would be to constrain the degradation and cross degradation parameters to be same, specifically ${\rho }_{\mathrm{gb}}\equiv {\rho }_{\mathrm{gc}}\equiv {\rho }_{\mathrm{gbc}}\equiv {\rho }_{g1}$ . This constraint would then enforce ${\Sigma }_i^s$ to be a valid covariance matrix since the Schur component ${\Sigma }_{\mathrm{ib}}^s-{\Sigma }_{\mathrm{ibc}}^s({\Sigma }_{\mathrm{ic}}^s){-}^1{\Sigma }_{\mathrm{ibc}}^s\equiv {\sigma }_{\mathrm{gb}}^2(1-{\rho }_g^2){\Sigma }_i^s({\rho }_{g1})$ , where ${\Sigma }_i^s({\rho }_{g1})$ has (j, k)th entry ${\rho }_{g1}^{|t_{\mathrm{ij}}-t_{\mathrm{ik}}|}$ , is constrained to be positive semi-definite. The resulting correlation structure would then be

$$\begin{array}{l}\text{Cor}\{B_i^s(t_{ij}),B_i^s(t_{ik})\}=\text{Cor}\{C_i^s(t_{ij}),C_i^s(t_{ik})\}={\rho }_{g1}^{|t_{ij}-t_{ik}|}\\ \text{Cor}\{B_i^s(t_{ij}),C_i^s(t_{ik})\}={\rho }_g{\rho }_{g1}^{|t_{ij}-t_{ik}|}\end{array}$$

(4)

In the motivating application, ρ_g was estimated close to one. Slight deviations from the correlation structure described by equation (4) (for example ${\rho }_{\mathrm{gb}}\equiv {\rho }_{\mathrm{gc}}\equiv {\rho }_{g1}$ and ${\rho }_{\mathrm{gbc}}={\rho }_{g1}\delta$ where $\delta \in (0,1)$ ) resulted in non-positive semi-definite matrices for various ${\Sigma }_i$ , and therefore the model fitting procedure was problematic. Note that a further simplification would be to constrain ${\rho }_g=1$ (in addition to ${\rho }_{g1}={\rho }_{\mathrm{gb}}$ ), this would result in the shared EOU process model. If, however, ρ_g takes a smaller value, and therefore the two Gaussian processes are less correlated, it would then be more plausible for the Gaussian processes to degrade at different rates. Hence, having unconstrained degradations parameters will likely be less problematic.

For completeness, note that

$$\begin{array}{l}{\hspace{0.17em}\left\{{\Sigma }_{ib}\right\}}_{jk}={\sigma }_b^2+{\sigma }_{gb}^2{\rho }_{gb}^{|t_{ij}-t_{ik}|}\\ {\hspace{0.17em}\left\{{\Sigma }_{ic}\right\}}_{jk}={\sigma }_c^2+{\sigma }_{gc}^2{\rho }_{gc}^{|t_{ij}-t_{ik}|}\\ {\left\{{\Sigma }_{ibc}\right\}}_{jk}={\sigma }_b{\sigma }_c\rho +{\sigma }_{gb}{\sigma }_{gc}{\rho }_g{\rho }_{gbc}^{|t_{ij}-t_{ik}|}\end{array}$$

(5)

3.2 Correlated random walks

The second model structure that we consider is defined when $\{B_i^s(t),C_i^s(t)\}$ are correlated continuous-time random walks. That is the model proposed by Ghosh and Albert. ⁹ Specifically, define sequentially $\{B_i(t_{\mathrm{ij}+1}),C_i(t_{\mathrm{ij}+1})\}|\{B_i(t_{\mathrm{ij}})=b_{\mathrm{ij}},C_i(t_{\mathrm{ij}})=c_{\mathrm{ij}}\}$ to be bivariate normal with mean $(b_{\mathrm{ij}},c_{\mathrm{ij}})$ and covariance matrix

$$\left(\begin{array}{cc}{\sigma }_{wb}^2(t_{ij+1}-t_{ij})& {\sigma }_{wb}{\sigma }_{wc}{\rho }_w(t_{ij+1}-t_{ij})\\ {\sigma }_{wb}{\sigma }_{wc}{\rho }_w(t_{ij+1}-t_{ij})& {\sigma }_{wc}^2(t_{ij+1}-t_{ij})\end{array}\right)$$

In addition $(b_{i1},c_{i1})=(b_i^r,c_i^r)$ are initiated at realisations of the patient-specific random effects. Here, ${\sigma }_{\mathrm{wb}}^2,{\sigma }_{\mathrm{wc}}^2$ and ρ_w are variance and correlation parameters that quantify serial correlation (both within and across processes). This model will be denoted by a correlated random walks (CRW) model and it contains as special cases.

Shared random walk model when $c_{\mathrm{ij}}=\theta b_{\mathrm{ij}}$ ,

Correlated random effects model when ${\sigma }_{\mathrm{wb}}^2\equiv {\sigma }_{\mathrm{wc}}^2\equiv {\rho }_w\equiv 0$ ,

Shared random effect model when $c_i^r=\theta b_i^r$ and ${\sigma }_{\mathrm{wb}}^2\equiv {\sigma }_{\mathrm{wc}}^2\equiv {\rho }_w\equiv 0$ .

Although the CRW model is less flexible than the general model, it has the advantage, from its sequential construction, of always being well defined even when the parameters are unconstrained (apart from the usual constraint that correlation parameters have modulus less than or equal to unity). Moreover

$$\begin{array}{l}{\hspace{0.17em}\left\{{\Sigma }_{ib}\right\}}_{jk}={\sigma }_b^2+{\sigma }_{wb}^2\{\text{min}(t_{ij},t_{ik})-t_{i1}\}\\ {\hspace{0.17em}\left\{{\Sigma }_{ic}\right\}}_{jk}={\sigma }_c^2+{\sigma }_{wc}^2\{\text{min}(t_{ij},t_{ik})-t_{i1}\}\\ {\left\{{\Sigma }_{ibc}\right\}}_{jk}={\sigma }_b{\sigma }_c\rho +{\sigma }_{wb}{\sigma }_{wc}{\rho }_w\{\text{min}(t_{ij},t_{ik})-t_{i1}\}\end{array}$$

(6)

4 Efficient maximum likelihood estimation procedure for stochastic processes models

This section describes our efficient maximum likelihood estimation procedure for the flexible models described in Section 3. Firstly, in Section 4.1, we describe a generic likelihood function for all of the described models. The multivariate normal identity that can be used to evaluate certain multi-dimensional integrals in terms of a multivariate normal cumulative distribution function is introduced in Section 4.2. Finally, in Section 4.3, we outline how to apply the multivariate normal identity in Section 4.2 to the generic likelihood function in Section 4.1, thus culminating in a computationally efficient likelihood. For completeness, we also provide computational simplifications for correlated stochastic processes models in the appendix.

4.1 Likelihoods

For ease of exposition, we describe the likelihood contribution from patient i. The likelihood can then be obtained by taking the product of all likelihood contributions from each patient. Firstly, we consider models that contain two (correlated) stochastic processes. For these models, the likelihood contribution from patient i is

$$\begin{array}{l}{L}_i(\Theta )={\displaystyle {\int }_{{\text{b}}_{\text{i}}}{\displaystyle {\int }_{{\text{c}}_{\text{i}}}\left[{\displaystyle \prod _{j=1}^{m_i}\Phi }{({\text{X}}_{ij}^{\top }\text{β}+b_{ij})}^{u_{ij}}{\{1-\Phi ({\text{X}}_{ij}^{\top }\text{β}+b_{ij})\}}^{1-u_{ij}}\right]}}\\ \times \left[{\displaystyle \prod _{j=1}^{m_i}{\left[\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left\{-\frac{{(g(y_{ij})-{\text{Z}}_{ij}^{\top }\text{γ}-c_{ij})}^2}{2{\sigma }^2}\right\}\right]}^{u_{ij}}}\right]{\varphi }^{(2m_i)}({\text{b}}_i,{\text{c}}_i;0,{\Sigma }_i)d{\text{b}}_{i}d{\text{c}}_i\end{array}$$

(7)

where $\Theta$ is a vector comprising all of the unknown parameters, $b_i=(b_{i1},\dots ,b_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top }$ , $c_i=(c_{i1},\dots ,c_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top },{\phi }^{(m)}(.;\mu ,\Sigma )$ is an m dimensional multivariate normal density with mean vector $\mu$ and covariance matrix $\Sigma$ , and ${\Sigma }_i$ is defined by either equation (5) or equation (6). Similarly, for models containing a single stochastic process (i.e. shared process models), the likelihood contribution from patient i is

$$\begin{array}{l}{L}_i(\Theta )={\displaystyle {\int }_{{\text{b}}_{\text{i}}}\left[{\displaystyle \prod _{j=1}^{m_i}\Phi }{({\text{X}}_{ij}^{\top }\text{β}+b_{ij})}^{u_{ij}}{\{1-\Phi ({\text{X}}_{ij}^{\top }\text{β}+b_{ij})\}}^{1-u_{ij}}\right]}\\ \times \left[{\displaystyle \prod _{j=1}^{m_i}{\left[\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left\{-\frac{{(g(y_{ij})-{\text{Z}}_{ij}^{\top }\text{γ}-\theta b_{ij})}^2}{2{\sigma }^2}\right\}\right]}^{u_{ij}}}\right]{\varphi }^{(m_i)}({\text{b}}_i;0,{\Sigma }_{ib})d{\text{b}}_i\end{array}$$

(8)

where ${\Sigma }_{\mathrm{ib}}$ can again be obtained from equation (5) or equation (6). We now define our generic likelihood contribution from patient i which encompasses all of the described models. Throughout we apply the following notation: $0$ and $1$ are $m_i\times 1$ vectors with all entries being zero and one respectively, $\text{diag}(v)$ is a matrix with diagonal elements v and zero otherwise, and $I_d$ is a d × d identity matrix. We also follow the convention that binary operations with a scalar and vector or matrix argument and unary operations with a vector argument are performed element-wise. In matrix form, we have

$$\begin{array}{l}{L}_i(\Theta )={\displaystyle {\int }_{{\text{l}}_i}{\Phi }^{(m_i)}}\left({\text{A}}_{i1}{\text{μ}}_{ib}+{\text{A}}_{i2}{\text{l}}_i;0,{\text{I}}_{m_i}\right){\left(\frac{1}{\sqrt{2\pi {\sigma }^2}}\right)}^{{\displaystyle \sum _{j=1}^{m_i}{u}_{ij}}}\\ \times \exp \left\{-\frac{{(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic}-{\text{A}}_{i4}{\text{l}}_i)}^{\top }(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic}-{\text{A}}_{i4}{\text{l}}_i)}{2{\sigma }^2}\right\}{\varphi }^{\{\text{dim}({\text{l}}_i)\}}({\text{l}}_i;0,{\Sigma }_{il})d{\text{l}}_i\end{array}$$

(9)

where $y_i=(y_{i1},\dots ,y_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top },{\mu }_{\mathrm{ib}}=X_i\beta ,{\mu }_{\mathrm{ic}}=Z_i\gamma ,X_i=(X_{i1},\dots ,X_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top },Z_i=(Z_{i1},\dots ,Z_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top }$ , and $A_{i1}=\text{diag}(2u_i-1),A_{i3}=\text{diag}(u_i)$ are $m_i\times m_i$ matrices with $u_i=(u_{i1},\dots ,u_{{\mathrm{im}}_i}){\hspace{0.17em}}^{\top }$ . Here ${\Phi }^{(d)}(.;0,\Sigma )$ represents the distribution function of ${\phi }^{(d)}(.;0,\Sigma )$ and $l_i$ is a (to be specified) column vector of random effects. Note that equation (9) has resulted from repeated application of the identity $1-\Phi (x)=\Phi (-x)$ .

The likelihood contribution from patient i, $L_i(\Theta )$ , is then obtained by specifying the vector of random effects $l_i$ and its covariance matrix ${\Sigma }_{\mathrm{il}}$ together with the $m_i\times \text{dim}(l_i)$ matrices $A_{i2}$ and $A_{i4}$ which describe how the random effects act on the binary and continuous components of the model. For (7), $l_i=(b_i,c_i),A_{i2}=(\text{diag}(2u_i-1),\text{diag}(0))$ and $A_{i4}=(\text{diag}(0),\text{diag}(u_i))$ . While for equation (8), $l_i=b_i,A_{i2}=\text{diag}(2u_i-1)$ and $A_{i4}=\text{diag}(\theta u_i)$ . Similarly, for the correlated random effects model, $l_i=(b_i,c_i),A_{i2}=(2u_i-1,0)$ and $A_{i4}=(0,u_i)$ , and for the shared random effect model, l_i = b_i, $A_{i2}=2u_i-1$ and $A_{i4}=\theta u_i$ .

4.2 Multivariate normal identity

In order to evaluate the likelihood described by equation (9), we derive a multivariate normal identity that makes use of a property of the multivariate normal distribution and the standard marginal cumulative distribution function identity. Firstly, suppose that $\omega =({\omega }_1,{\omega }_2){\hspace{0.17em}}^{\top }$ follows a multivariate normal distribution with mean vector $(0,\eta ){\hspace{0.17em}}^{\top }$ where ${\omega }_1$ and $0$ are $k_1\times 1$ vectors and ${\omega }_2$ and $\eta$ are $k_2\times 1$ vectors, respectively. Furthermore, suppose that the covariance matrix of $\omega$ is the $(k_1+k_2)\times (k_1+k_2)$ matrix $\Sigma$ where the first k₁ rows of $\Sigma$ is the $k_1\times (k_1+k_2)$ matrix $({\Sigma }_{22},{\Sigma }_{12}^{\top })$ and the remaining k₂ rows of $\Sigma$ is the $k_2\times (k_1+k_2)$ matrix $({\Sigma }_{12},{\Sigma }_{11})$ respectively. It is a well-known result that ${\phi }^{(k_1+k_2)}(\omega ;(0,\eta ){\hspace{0.17em}}^{\top },\Sigma )={\phi }^{(k_1)}({\omega }_1;{\Sigma }_{12}^{\top }{\Sigma }_{11}^{-1}({\omega }_2-\eta ),{\Sigma }_{22}-{\Sigma }_{12}^{\top }{\Sigma }_{11}^{-1}{\Sigma }_{12}){\phi }^{(k_2)}({\omega }_2;\eta ,{\Sigma }_{11})$ where the right-hand side is the product of the conditional density of ${\omega }_1|{\omega }_2$ and the marginal density of ${\omega }_2$ . By applying the standard marginal cumulative distribution function identity $F_{{\omega }_1}({\omega }_1)={\int }_{{\omega }_2}{F}_{{\omega }_1|{\omega }_2}({\omega }_1|{\omega }_2)f_{{\omega }_2}({\omega }_2)\text{d}{\omega }_2$ where the integrand is based on the right-hand side of the above result, we obtain the multivariate normal identity:

$${\Phi }^{(k_1)}({\omega }_1;0,{\Sigma }_{22})={\int }_{{\omega }_2}{\Phi }^{(k_1)}({\omega }_1-{\Sigma }_{12}^{\top }{\Sigma }_{11}^{-1}({\omega }_2-\eta );0,{\Sigma }_{22}-{\Sigma }_{12}^{\top }{\Sigma }_{11}^{-1}{\Sigma }_{12}){\phi }^{(k_2)}({\omega }_2;\eta ,{\Sigma }_{11})\text{d}{\omega }_2$$

(10)

by noting that the marginal distribution of ${\omega }_1$ is multivariate normal with mean vector $0$ and covariance matrix ${\Sigma }_{22}$ .

Returning to the application, the general idea is to rearrange equation (9) to take the form of the right-hand side of equation (10), and then to use equation (10) to compute the integrations over the random effects in terms of an m_i dimensional normal cumulative distribution function. Because there exists efficient implementations of the multivariate normal cumulative distribution function, this approach will allow for the efficient computation of the generic likelihood. We note that Barrett et al. ¹⁵ used equation (10) to obtain computationally efficient likelihoods of flexible models that jointly consider longitudinal and time to event outcomes. Equation (10) also arises frequently in results concerning the multivariate skew normal distribution.^16–19

4.3 Re-expressing the likelihoods

This section demonstrates how equation (9) (the likelihood contribution from patient i) can be re-expressed. We firstly consider the integrand terms resulting from the continuous component and random effects. That is

$$(\frac{1}{\sqrt{2\pi {\sigma }^2}}{)}^{\sum _{j=1}^{m_i}{u}_{\mathrm{ij}}}\exp \{-\frac{(g(y_i)-A_{i3}{\mu }_{\mathrm{ic}}-A_{i4}{l}_i){\hspace{0.17em}}^{\top }(g(y_i)-A_{i3}{\mu }_{\mathrm{ic}}-A_{i4}{l}_i)}{2{\sigma }^2}\}{\phi }^{\{\text{dim}(l_i)\}}(l_i;0,{\Sigma }_{\mathrm{il}})$$

(11)

By completing the square in $l_i$ (see the appendix for more details), equation (11) can be rearranged as

$$L_{i1}{\phi }^{\{\text{dim}(l_i)\}}(l_i;h_i,H_i^{-1})$$

(12)

where

$$\begin{array}{l}{\text{H}}_i={\text{A}}_{i4}^{\top }{\text{A}}_{i4}/{\sigma }^2+{({\Sigma }_{il})}^{-1}\\ {\text{h}}_i={\text{H}}_i^{-1}{\text{A}}_{i4}^{\top }(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic})/{\sigma }^2\end{array}$$

(13)

and

$$L_{i1}=(\frac{1}{\sqrt{2\pi {\sigma }^2}}{)}^{\sum _{j=1}^{m_i}{u}_{\mathrm{ij}}}\frac{1}{|{\Sigma }_{\mathrm{il}}{H}_i{|}^{1/2}}\exp \{-\frac{(g(y_i)-A_{i3}{\mu }_{\mathrm{ic}}){\hspace{0.17em}}^{\top }(g(y_i)-A_{i3}{\mu }_{\mathrm{ic}})}{2{\sigma }^2}+\frac{h_i^{\top }{H}_i{h}_i}{2}\}$$

(14)

is independent of $l_i$ . Substituting equation (12) into equation (9), we consider the integral (ignoring $L_{i1}$ )

$${\int }_{l_i}{\Phi }^{(m_i)}(A_{i1}{\mu }_{\mathrm{ib}}+A_{i2}{l}_i;0,I_{m_i}){\phi }^{\{\text{dim}(l_i)\}}(l_i;h_i,H_i^{-1})\text{d}{l}_i.$$

(15)

We can re-express the argument and covariance matrix of the multivariate normal distribution function in equation (15) as

$$\begin{array}{l}{\text{A}}_{i1}{\text{μ}}_{ib}+{\text{A}}_{i2}{\text{l}}_i={\text{A}}_{i1}{\text{μ}}_{ib}+{\text{A}}_{i2}{\text{h}}_i-(-{\text{H}}_i^{-1}{\text{A}}_{i2}^{\top }{)}^{\top }{\text{H}}_i({\text{l}}_i-{\text{h}}_i)\\ {\text{I}}_{m_i}={\text{I}}_{m_i}+{\text{A}}_{i2}{\text{H}}_i^{-1}{\text{A}}_{i2}^{\top }-(-{\text{H}}_i^{-1}{\text{A}}_{i2}^{\top }{)}^{\top }{\text{H}}_i(-{\text{H}}_i^{-1}{\text{A}}_{i2}^{\top })\end{array}$$

(16)

Therefore equation (15), after applying the multivariate normal identity (described equation (10)), is equivalent to

$${\Phi }^{(m_i)}(A_{i1}{\mu }_{\mathrm{ib}}+A_{i2}{h}_i;0,I_{m_i}+A_{i2}{H}_i^{-1}{A}_{i2}^{\top })$$

(17)

Based on the above expressions, the likelihood contribution from patient i can now be re-expressed as

$$\begin{array}{l}{L}_i(\Theta )={\Phi }^{(m_i)}\left({\text{A}}_{i1}{\text{μ}}_{ib}+{\text{A}}_{i2}{\text{h}}_i;0,{\text{I}}_{m_i}+{\text{A}}_{i2}{\text{H}}_i^{-1}{\text{A}}_{i2}^{\top }\right){\left(\frac{1}{\sqrt{2\pi {\sigma }^2}}\right)}^{{\displaystyle \sum _{j=1}^{m_i}{u}_{ij}}}\\ \times \frac{1}{|{\Sigma }_{il}{\text{H}}_i{|}^{1/2}}\exp \left\{-\frac{{(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic})}^{\top }(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic})}{2{\sigma }^2}+\frac{{\text{h}}_i^{\top }{\text{H}}_i{\text{h}}_i}{2}\right\}\end{array}$$

(18)

where

$$\begin{array}{l}{\text{A}}_{i1}=\text{diag}(2{\text{u}}_i-1)\\ {\text{A}}_{i3}=\text{diag}({\text{u}}_i)\\ {\text{H}}_i={\text{A}}_{i4}^{\top }{\text{A}}_{i4}/{\sigma }^2+{({\Sigma }_{il})}^{-1}\\ {\text{h}}_i={\text{H}}_i^{-1}{\text{A}}_{i4}^{\top }(g({\text{y}}_i)-{\text{A}}_{i3}{\text{μ}}_{ic})/{\sigma }^2\end{array}$$

(19)

and ${\Sigma }_{\mathrm{il}}={\Sigma }_i$ or ${\Sigma }_{\mathrm{ib}}$ with $A_{i2}$ and $A_{i4}$ defined by the specified model.

From equations (18) and (19), it is now evident that evaluating the integrations involved in $L_i(\Theta )$ reduces to computing the cumulative distribution function of a multivariate normal distribution. This can be performed efficiently, for example by using the R ²⁰ package mnormt. ²¹ The model fitting procedure is then completed by maximizing the log-likelihood, for example by using the BFGS ²² optimization technique, to obtain parameter estimates and asymptotic standard errors (from the observed Fisher information matrix) of model parameters.

5 Application: Patient-specific inference

Using the estimation procedure described in Section 4, we demonstrate how patient-specific inference on the probability of being disabled and the transformed mean HAQ score conditional on disability can be obtained. Specifically, how a unit change in covariate values impacts these quantities for any specific patient. We consider the covariate effects of the number of clinically damaged joints (time-dependent), the number of actively inflamed joints (time-dependent), sex (coded as 1 for males and 0 for females), arthritis duration in years (time-dependent), and age at onset of arthritis in years (standardise). Following Su et al.,^12,13 no transformation was applied to the non-zero HAQ scores, i.e. g(y) = y.

Initially, models with two stochastic processes were fitted to the HAQ data. This resulted in large estimated correlation parameters between the random effects (i.e. $\rho \approx 1$ ) and stochastic processes for both the correlated Gaussian processes and random walks cases (i.e. ρ_g and ${\rho }_w\approx 1$ ). These results therefore suggested a single stochastic process would be sufficient for describing the data. The shared EOU model was then fitted. However, the analysis provided evidence for model over-parameterisation as ${\hat{\sigma }}^2$ appeared to converge at virtually zero and a positive-definite observed Fisher information matrix could not be attained (even when a considerably smaller tolerance level than the default was specified for the computation of multivariate normal probabilities). We therefore considered the shared random walk and OU process models, and for comparative purposes, the shared random effect model. The models containing stochastic processes were fitted using the likelihood described by equations (18) and (19), while the shared random effect model was fitted using numerical integration (since only a single integration per patient is required). The same parameter estimates for the shared random effect model were obtained when equations (18) and (19) were used in the model fitting procedure.

Table 1 presents the results of the fitted models. Across the models, the covariate effects on the mean conditional on disability are seen to be relatively similar as the confidence intervals generally overlap. In addition, the models are in agreement with regard to the association of each covariate apart from arthritis duration. Arthritis duration is statistically significant in the shared random effect model but is not statistically significant in the models that incorporate stochastic processes. It is interesting to note that there are strong agreements regarding the covariate effect of the number of active joints (similar parameter estimates across models and relatively narrow confidence intervals). The models indicate an additional actively inflamed joint will increase the mean HAQ score conditional on disability by approximately 0.21 for any specific patient. For the binary component, the covariate effects are again seen to be relatively similar due to the overlapping confidence intervals. Their interpretation through the direction of association and statistical significance are also consistent across models. The covariate effects from the shared random effect model does, however, consistently demonstrate attenuation to the null when compared to the other models with stochastic processes.

Table 1.

Table displaying patient-specific effects and corresponding 95% Wald intervals on the probability of being disabled and the mean HAQ score conditional on disability.

	Shared random walk	Shared OU process	Shared random effect
Binary component
Damaged joints	0.031 (0.013, 0.049)	0.04 (0.021, 0.059)	0.012 (0.002, 0.022)
Active joints	0.16 (0.13, 0.18)	0.17 (0.15, 0.2)	0.15 (0.13, 0.16)
Sex	−1.72 (−1.9, −1.54)	−2.17 (−2.92, −1.42)	−1.34 (−1.65, −1.02)
Arthritis duration	0.042 (0.025, 0.058)	0.044 (0.022, 0.065)	0.034 (0.023, 0.044)
Age at arthritis onset a	0.56 (0.45, 0.66)	0.68 (0.48, 0.88)	0.45 (0.3, 0.6)
Intercept	1.77 (1.62, 1.92)	1.97 (1.07, 2.87)	1.11 (0.81, 1.41)
Continuous component
Damaged joints	0.0065 (0.0032, 0.0097)	0.0078 (0.0046, 0.011)	0.0033 (0.00086, 0.0058)
Active joints	0.021 (0.019, 0.023)	0.021 (0.019, 0.023)	0.02 (0.018, 0.022)
Sex	−0.29 (−0.34, −0.24)	−0.35 (−0.46, −0.23)	−0.29 (−0.37, −0.21)
Arthritis duration	0.0025 (−0.0008, 0.0058)	0.0035 (−0.001, 0.0081)	0.0067 (0.0041, 0.0093)
Age at arthritis onset a	0.076 (0.048, 0.1)	0.093 (0.046, 0.14)	0.086 (0.048, 0.12)
Intercept	0.62 (0.57, 0.66)	0.59 (0.44, 0.74)	0.63 (0.56, 0.7)
θ	0.2 (0.19, 0.22)	0.19 (0.17, 0.21)	0.27 (0.24, 0.3)
${\sigma }^2$	0.074 (0.069, 0.079)	0.066 (0.06, 0.072)	0.12 (0.11, 0.12)
${\sigma }_b^2$	6.3 (5.64, 7.04)		3.29 (2.64, 4.1)
${\sigma }_{\mathrm{gb}}^2$		10.11 (8.03, 12.74)
${\sigma }_{\mathrm{wb}}^2$	0.58 (0.52, 0.65)
ρgb		0.95 (0.94, 0.96)
Log-likelihood	−3279.08	−3282.11	−3500.48

a

Denotes the standardised version of the covariate.

A generalized likelihood ratio test of ${\sigma }_{\mathrm{wb}}^2=0$ and ${\rho }_{\mathrm{gb}}=1$ produced p values of < 0.001 therefore suggesting preference towards the shared random walk and OU process models respectively when compared to the shared random effect model. Since the shared random walk and OU process models contain the same number of parameters, information criteria, such as AIC, would indicate (weakly) that the shared random walk model is preferable. It is also worth noting that the heterogeneity parameter in the binary component (i.e. ${\sigma }_b^2$ or ${\sigma }_{\mathrm{gb}}^2$ ) is significantly lower in the shared random effect model. For this model, this parameter governs both the heterogeneity and correlation due to repeated measurements and therefore in light of greater unaccounted heterogeneity (compared to the models with stochastic processes), less correlation is expected. ²³ In the continuous component, where ${\sigma }^2$ also accounts for heterogeneity, a smaller difference between the heterogeneity parameters (i.e. ${\theta }^2{\sigma }_b^2$ or ${\theta }^2{\sigma }_{\mathrm{gb}}^2$ ) is seen; in the order of the models displayed in the table (from right to left), the heterogeneity parameters are 0.24, 0.36 and 0.25, respectively.

6 Modelling the overall marginal mean

In many cases, it is of interest to obtain population-based inference in addition/as opposed to patient-specific inference. For example, for strategic public health policy purposes, it would be more clinically meaningful to obtain covariate effects on quantities of interest after averaging over all patients. Currently, the proposed models are parametrised to allow easily interpretable patient-specific covariate effects, those with $B_i(t_{\mathrm{ij}})=b_{\mathrm{ij}}$ and $C_i(t_{\mathrm{ij}})=c_{\mathrm{ij}}$ , to act on the patient-specific mean of the transformed positive values (i.e. $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0,C_i(t_{\mathrm{ij}})=c_{\mathrm{ij}}]$ ) and the patient-specific probability of a having a positive value (i.e. $\mathbb{P}(U_{\mathrm{ij}}=1|B_i(t_{\mathrm{ij}})=b_{\mathrm{ij}})$ ). However, under this parametrisation, it no longer becomes straightforward to obtain easily interpretable population-level covariate effects on the marginal mean of the transformed positive values (the mean of the transformed positive values after averaging over all b_ij and c_ij, i.e. $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0]$ ) since it is a highly non-linear function of the linear predictors in the binary and continuous components. ⁵ Thus, the effect of a single covariate is generally interpreted by fixing other covariates at certain values. ⁸ This problem remains even when population-level covariate effects on the overall marginal mean of the transformed values (i.e. $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})]$ ) are of primary interest, which has strongly been argued as an important target of inference; ²⁴ it is estimated using data from the same patients over time (unlike $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0]$ ) and it is a measure of the undecomposed outcome. We reiterate that in considering the overall marginal mean of the transformed values as a target of inference, we assume that the monotonic transformation function is such that $g(0)=0$ and $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ is positive and approximately Gaussian with constant variance ${\sigma }^2$ .

In order to obtain population-based inference on the overall marginal mean of the transformed values, Smith et al. ⁴ proposed the following model parameterisation

$$\begin{array}{l}\mathbb{P}(U_{ij}=1|B_i^r=b_i^r)=g_1({\text{X}}_{ij}^{\top }\text{β}+b_i^r)\\ \mathbb{E}[g(Y_{ij})|C_i^r=c_i^r]=g_2({\text{Z}}_{ij}^{\top }\text{α}+c_i^r)\end{array}$$

(20)

where $g_1(·)$ and $g_2(·)$ are monotonic link functions and $B_i^r,C_i^r$ are, as before, zero mean bivariate normal patient-specific random effects. Recall that transformation and link functions differ in that transformation functions are applied prior to modelling. In their specific context, Smith et al. ⁴ considered the identity transformation for $g(·)$ but allowed the positive values of Y_ij to follow a log-skew-normal distribution. Under this parametrisation, for a suitably chosen link function such as $g_2^{-1}$ being the identity or log link, it is implicit that easily interpretable covariate effects on the overall marginal mean of $g(Y_{\mathrm{ij}}),\alpha$ , can now be obtained. Smith et al. ⁴ implemented this model by using a Bayesian estimation approach with

$$\mathbb{\text{E}}[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0,B_i^r=b_i^r,C_i^r=c_i^r]=\frac{g_2(Z_{\mathrm{ij}}^{\top }\alpha +c_i^r)}{g_1(X_{\mathrm{ij}}^{\top }\beta +b_i^r)}$$

specified in the likelihoods defined by equation (7) or equation (8). Note that $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0,B_i^r=b_i^r,C_i^r=c_i^r]$ is no longer parametrised to be equivalent to a monotonic function of a linear predictor, as was specified before. While this approach for modelling the overall marginal mean is intuitive, it is clear that the multivariate normal identity in Section 4.2 can no longer be used to compute the integrations over the multi-dimensional random effects in the marginal likelihood. Thus, as mentioned in the introduction, implementation of such models can be computationally challenging, especially for our situation where it would be of interest to consider b_ij and c_ij (i.e. realisations of stochastic processes) instead of $b_i^r$ and $c_i^r$ (i.e. realisations of patient-specific random effects) in equation (20).

We now propose another method which would allow easily interpretable covariate effects to act on the overall marginal mean of $g(Y_{\mathrm{ij}})$ . In contrast, this method facilitates the inclusion of stochastic processes because it retains the proposed efficient implementation procedure described in Section 4. To the best of our knowledge, there are no other methods in the literature that facilitates the practical implementation of stochastic processes models for directly modelling the overall marginal mean.

We first begin by computing the overall marginal mean of $g(Y_{\mathrm{ij}})$ when

$$\begin{array}{l}\mathbb{P}(U_{ij}=1|B_i(t_{ij})=b_{ij})=g_1({\text{X}}_{ij}^{\top }\text{β}+b_{ij})\\ \mathbb{E}[g(Y_{ij})|Y_{ij}>0,C_i(t_{ij})=c_{ij}]={\Delta }_{ij}+c_{ij}\end{array}$$

(21)

where Δ _ij is a function of covariates (at the jth visit from patient i) and regression coefficients only. That is c_ij is now assumed to act additively on the mean of the transformed positive values and not on the overall mean of the transformed values as is the case in equation (20). For models with two processes, the overall marginal mean of $g(Y_{\mathrm{ij}})$ is defined by

$$\mathbb{\text{E}}[g(Y_{\mathrm{ij}})]={\int }_b{\int }_c\Phi (X_{\mathrm{ij}}^{\top }\beta +b_{\mathrm{ij}})({\Delta }_{\mathrm{ij}}+c_{\mathrm{ij}}){\phi }^{(2)}(b_{\mathrm{ij}},c_{\mathrm{ij}};0,{\Sigma }_{\mathrm{ij}})\text{d}{b}_{\mathrm{ij}}\text{d}{c}_{\mathrm{ij}}$$

where

$${\Sigma }_{ij}=\left(\begin{array}{cc}{\sigma }_{bij}^2& {\sigma }_{bij}{\sigma }_{cij}{\rho }_{bcij}\\ {\sigma }_{bij}{\sigma }_{cij}{\rho }_{bcij}& {\sigma }_{cij}^2\end{array}\right)$$

and ${\sigma }_{\mathrm{bij}}^2,{\sigma }_{\mathrm{cij}}^2$ and ρ_bcij are the variances and correlation of $B_i(t_{\mathrm{ij}})$ and $C_i(t_{\mathrm{ij}})$ , respectively. Similarly, the overall marginal mean for a shared process model is given by

$$\mathbb{\text{E}}[g(Y_{\mathrm{ij}})]={\int }_b\Phi (X_{\mathrm{ij}}^{\top }\beta +b_{\mathrm{ij}})({\Delta }_{\mathrm{ij}}+\theta b_{\mathrm{ij}}){\phi }^{(1)}(b_{\mathrm{ij}};0,{\sigma }_{\mathrm{bij}}^2)\text{d}{b}_{\mathrm{ij}}.$$

Conveniently, these integrals can be computed analytically and this results in

$$\Phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}){\Delta }_{\mathrm{ij}}+\frac{{\sigma }_{\mathrm{bij}}{\sigma }_{\mathrm{cij}}{\rho }_{\mathrm{bcij}}}{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}\phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})$$

and

$$\Phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}){\Delta }_{\mathrm{ij}}+\frac{\theta {\sigma }_{\mathrm{bij}}^2}{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}\phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})$$

respectively. The derivation of the first overall marginal mean of $g(Y_{\mathrm{ij}})$ (resulting from models with two processes) can be found in the supplementary material of Tom et al., ⁵ and the second overall marginal mean of $g(Y_{\mathrm{ij}})$ is derived in the appendix. If we specify $\mathbb{\text{E}}[g(Y_{\mathrm{ij}})]=Z_{\mathrm{ij}}^{\top }\alpha$ , we can then reparametrise

$${\Delta }_{\mathrm{ij}}=[Z_{\mathrm{ij}}^{\top }\alpha -\frac{{\sigma }_{\mathrm{bij}}{\sigma }_{\mathrm{cij}}{\rho }_{\mathrm{bcij}}}{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}\phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})]/\Phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})$$

(22)

and

$${\Delta }_{\mathrm{ij}}=[Z_{\mathrm{ij}}^{\top }\alpha -\frac{\theta {\sigma }_{\mathrm{bij}}^2}{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}}\phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})]/\Phi (\frac{X_{\mathrm{ij}}^{\top }\beta }{\sqrt{1+{\sigma }_{\mathrm{bij}}^2}})$$

in the respective models. Thus, as in equation (20), $\alpha$ offers easily interpretable covariate effects of $Z_{\mathrm{ij}}$ on the overall marginal mean of $g(Y_{\mathrm{ij}})$ (by definition). In particular, a unit change in components of $Z_{\mathrm{ij}}$ will increase the overall marginal mean of $g(Y_{\mathrm{ij}})$ by the respective components in $\alpha$ . However, from equations (21) and (22), it is also evident that replacing $Z_{\mathrm{ij}}^{\top }\gamma$ with Δ _ij in Section 4 will still allow the proposed efficient estimation procedure to be applied. It is also possible to reparametrise patient-specific covariate effects $\beta$ in the binary component in terms of population-level covariate effects $\xi$ , specifically $\beta =\xi \sqrt{1+{\sigma }_{\mathrm{bij}}^2}$ , since it can be shown that $\mathbb{P}(U_{\mathrm{ij}}=1)=\Phi (X_{\mathrm{ij}}^{\top }\beta /\sqrt{1+{\sigma }_{\mathrm{bij}}^2})$ . This relationship is easily proved. In the motivating application, this reparametrisation led to a numerically unstable optimization routine, therefore $(\beta ,\alpha )$ was estimated with $\hat{\xi }$ obtained as $\hat{\beta }/\sqrt{1+{\hat{\sigma }}_{\mathrm{bij}}^2}$ and standard errors were calculated using the delta method.

6.1 Population-based inference

Using the parameterisations described in the previous subsection, we demonstrate how population-based inference on the probability of being disabled and the overall marginal mean HAQ score can be obtained. Specifically, on averaging across patients, how a unit change in covariate values impacts these quantities. For illustrative purposes, the same covariates as those considered in the patient-specific case are considered. Note that for generalized linear models, conditional and marginal covariate effects will generally differ unless certain random effects distributions and link functions are chosen. ²⁵

As mentioned, marginal covariate effects on the probability of being disabled, $\xi$ , were obtained from $\beta /\sqrt{1+{\sigma }_{\mathrm{bij}}^2}$ with $\beta$ and ${\sigma }_{\mathrm{bij}}^2$ (the variance of $B_i(t_{\mathrm{ij}})$ ) estimated using the model fitting procedure. The shared OU process and random effect models, where ${\sigma }_{\mathrm{bij}}^2$ does not depend on j (or i) for these models, were considered. For models with random walks, ${\sigma }_{\mathrm{bij}}^2$ varies with j and therefore $\xi$ will have a time-dependent interpretation. For simplicity, these models are not considered. The shared random effect model was fitted using the parameterisation described by equation (20), with $g_1(y)=\Phi (y)$ and $g_2(y)=y$ , and using the parameterisation described by equations (21) and (22), thus the same link functions are used and the inferences (at the population-level) from these models are comparable. These models will be denoted as shared random effect model-overall and -conditional respectively. Note that unlike at the population-level, the patient-specific assumptions from the shared random effect model-overall and -conditional are vastly different. The shared random effect-overall model assumes that the overall patient-specific mean, i.e. $\mathbb{\text{E}}[Y_{\mathrm{ij}}|B_i^r=b_i^r,C_i^r=c_i^r]$ , has a linear form, namely $Z_{\mathrm{ij}}^{\top }\alpha +c_i^r$ . While the shared random effect-conditional model assumes that this quantity takes a particular non-linear form, namely $({\Delta }_{\mathrm{ij}}+c_i^r)\Phi (X_{\mathrm{ij}}^{\top }\beta +b_i^r)$ . As before, the shared random effect models (both -overall and -conditional) were fitted using numerical integration and maximum likelihood estimation under the assumption that the positive values follow a normal distribution with constant variance.

Table 2 presents the results. Population-level covariate effects on the overall marginal mean are seen to be relatively similar across models due to the considerable overlap in confidence intervals. All three models are in strong agreement regarding the population-level covariate effect of the number of active joints. That is, on average, patients with an additional actively inflamed joint have an overall mean HAQ score increased by approximately 0.02. In contrast to the patient-specific case, the population-level covariate effects on the probability of being disabled are now more consistent across models. A generalized likelihood ratio test of ${\rho }_{\mathrm{gb}}=1$ produced a p value of < 0.001 and therefore the shared OU process model is to be preferred over the shared random effect-conditional model. Log-likelihood values also indicate slight preference to the shared random effect-conditional model (−3507.63) over the shared random effect-overall model (−3582.78).

Table 2.

Table displaying population-level effects and corresponding 95% Wald intervals on the probability of being disabled and the overall marginal mean HAQ score.

	Shared OU process	Shared random effect-conditional	Shared random effect-overall
Binary component
Damaged joints	0.012 (0.006, 0.018)	0.0057 (0.0006, 0.011)	0.0066 (0.002, 0.011)
Active joints	0.05 (0.043, 0.058)	0.063 (0.054, 0.071)	0.051 (0.044, 0.058)
Sex	−0.65 (−0.88, −0.43)	−0.64 (−0.8, −0.48)	−0.61 (−0.76, −0.46)
Arthritis duration	0.013 (0.0057, 0.02)	0.014 (0.0095, 0.02)	0.012 (0.0078, 0.017)
Age at arthritis onset a	0.19 (0.14, 0.24)	0.2 (0.13, 0.27)	0.19 (0.12, 0.26)
Intercept	0.59 (0.35, 0.83)	0.56 (0.42, 0.7)	0.68 (0.55, 0.81)
Overall marginal mean
Damaged Joints	0.0064 (0.0036, 0.0094)	0.0029 (0.00054, 0.0053)	0.004 (0.0016, 0.0065)
Active joints	0.02 (0.018, 0.022)	0.021 (0.019, 0.023)	0.022 (0.02, 0.024)
Sex	−0.3 (−0.4, −0.19)	−0.28 (−0.36, −0.21)	−0.31 (−0.39, −0.23)
Arthritis duration	0.0035 (0.00032, 0.0068)	0.006 (0.0038, 0.0082)	0.0054 (0.0031, 0.0077)
Age at arthritis onset a	0.073 (0.049, 0.098)	0.08 (0.048, 0.11)	0.09 (0.052, 0.13)
Intercept	0.62 (0.5, 0.73)	0.61 (0.54, 0.68)	0.59 (0.52, 0.66)
θ	0.18 (0.16, 0.2)	0.27 (0.24, 0.29)	0.22 (0.2, 0.24)
${\sigma }^2$	0.064 (0.057, 0.071)	0.12 (0.11, 0.12)	0.12 (0.11, 0.13)
${\sigma }_b^2$		3.63 (2.92, 4.5)	5.19 (4.26, 6.34)
${\sigma }_{\mathrm{gb}}^2$	11.12 (8.71, 14.19)
ρgb	0.95 (0.94, 0.96)
Log−likelihood	−3277.2	−3507.63	−3582.78

a

Denotes the standardised version of the covariate.

7 Discussion

This paper reconsiders the flexible two-part models of Albert and Shen ⁸ and Ghosh and Albert ⁹ and proposes an efficient method of implementation. Specifically, the problem of integrating over high dimensional random effects is replaced by evaluating the cumulative distribution function of a multivariate normal distribution. This leads to efficient algorithms being employed and results in only an optimization procedure being required for model fitting. Furthermore, while retaining the flexibility of including stochastic processes and the practicality of an efficient model fitting procedure, this paper also provides model parameterisations which allow easily interpretable covariate effects to act on the overall marginal mean. The proposed methodology was applied to a psoriatic arthritis data set with extensive follow-up information.

Through their application and a simulation study, Albert and Shen ⁸ demonstrated that overall conditional means (conditional on realisations of stochastic processes) may suffer from bias if serial correlation is present but a shared random effect model is used instead. Furthermore, as the shared random effect model becomes more misspecified (ρ_gb decreases from one), the degree of bias increases. However, under the same set-up, overall marginal means were less susceptible to bias. In the motivating application, the estimated degradation parameters from the shared OU process models were ${\hat{\rho }}_{\mathrm{gb}}=0.95$ in both applications (Sections 5 and 6.1). The reasonably high estimated correlation may therefore explain why the shared random effect model $({\rho }_{\mathrm{gb}}=1)$ was a reasonable approximation in terms of estimating regression coefficients, although it was substantially the worst fitting model.

Preliminary analyses suggested shared process models were reasonable for our data since $\hat{\rho },{\hat{\rho }}_w$ and ${\hat{\rho }}_g\approx 1$ when the described bivariate processes models were fitted. Although this may not be surprising as both parts of the model are describing the same response process, it is worth noting that the estimated correlation parameter (between processes) can in principle take a value between $(-1,1)$ as evidenced in other works.^7,9 Our preliminary analyses also demonstrated the need for careful evaluation of models fitted as problems with over fitting may arise. This was evident when the estimated random variation parameter was estimated to be virtually zero (i.e. ${\hat{\sigma }}^2\approx 0$ ) and the observed Fisher information matrix was non positive-definite, even when a considerably smaller tolerance level than the default was specified for computing multivariate normal probabilities.

As mentioned in Section 6, the proposed model parameterisations were motivated by making inference on the overall marginal mean. In this regard, covariate effects (both patient-specific and population-level) on the mean of the positive values and its correlation structure were assumed not of interest. If the mean of the positive values is of primary interest, it would be more sensible to directly use equations (18) and (19), as in Section 5, to obtain patient-specific effects or derive a similar parameterisation, as in Section 6, to obtain population-level effects.

A limitation of the current framework is that it is based on the assumption $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ is approximately Gaussian with constant variance ${\sigma }^2$ . Specifically, in situations where $g(·)$ is required to be complex so that this assumption will at least approximately hold, the resulting inferential targets will no longer be intuitively interpretable owing to the complexity of the transformation function. One approach that may weaken the need to assume normality of $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ , particularly when the outcome exhibits a large amount of right skewness (e.g. medical expenditures), would be to make the alternative assumption $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ follows a log-normal distribution. This may allow less complex and hence more interpretable transformation functions to be applied to the outcome without having to strongly violate the assumption on $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ . Under this alternative assumption, we provide details in the supplementary materials of how easily interpretable inference on the overall marginal mean and on the mean of the positive transformed outcomes can be obtained with computationally efficient likelihoods. Similar techniques to those in the supplementary materials can also be used when the assumption $g(Y_{\mathrm{ij}})|Y_{\mathrm{ij}}>0$ follows a log-skew-normal distribution is of interest. Although this comes at the cost of having an increased number of integrations in the marginal likelihood.

Finally, the model described by equations (18) and (19) with possible simplifications described in Appendix 2 is very general. Although it was derived in the context of longitudinal semicontinuous data, it contains the model described by Barrett et al. ¹⁵ for the longitudinal and survival outcomes setting and implicitly provides a model for clustered cross-sectional semicontinuous data, where the index (i, j) specifies the jth outcome from the ith cluster. The multivariate normal identity described in Section 4.2 can also facilitate the fitting of flexible models describing clustered binary data and continuous bounded outcome data. ¹⁴ However, it should be noted that care is required when specifying an appropriate/suitable correlation structure. Particularly, the covariance matrix must be constrained to be symmetric and positive semi-definite otherwise the model fitting procedure will likely be problematic, as was found here. For these alternative situations, the proposed methodology does nevertheless offer a strong basis, especially with regard to implementation, for the developing of flexible models.

Supplemental Material

Supplemental material for Two-part models with stochastic processes for modelling longitudinal semicontinuous data: Computationally efficient inference and modelling the overall marginal mean

Supplemental Material for Two-part models with stochastic processes for modelling longitudinal semicontinuous data: Computationally efficient inference and modelling the overall marginal mean by Sean Yiu and Brian DM Tom in Statistical Methods in Medical Research