An alternative characterization of MAR in shared parameter models for incomplete longitudinal data and its utilization for sensitivity analysis

Abstract

Dropout is a common complication in longitudinal studies, especially since the distinction between missing not at random (MNAR) and missing at random (MAR) dropout is intractable. Consequently, one starts with an analysis that is valid under MAR and then performs a sensitivity analysis by considering MNAR departures from it. To this end, specific classes of joint models, such as pattern-mixture models (PMMs) and selection models (SeMs), have been proposed. On the contrary, shared-parameter models (SPMs) have received less attention, possibly because they do not embody a characterization of MAR. A few approaches to achieve MAR in SPMs exist, but are difficult to implement in existing software. In this article, we focus on SPMs for incomplete longitudinal and time-to-dropout data and propose an alternative characterization of MAR by exploiting the conditional independence assumption, under which outcome and missingness are independent given a set of random effects. By doing so, the censoring distribution can be utilized to cover a wide range of assumptions for the missing data mechanism on the subject-specific level. This approach offers substantial advantages over its counterparts and can be easily implemented in existing software. More specifically, it offers flexibility over the assumption for the missing data generating mechanism that governs dropout by allowing subject-specific perturbations of the censoring distribution, whereas in PMMs and SeMs dropout is considered MNAR strictly.

1 Introduction

Follow-up studies that include human subjects are known to be highly susceptible to dropout. The extent to which dropout will complicate the subsequent analysis depends on how the missingness is associated with the observed and unobserved longitudinal outcomes. According to the taxonomy introduced by (Rubin (1976)) and (Little and Rubin (2002)), this may happen according to three distinct missing data mechanisms: missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR). Under MCAR, the missingness depends on neither observed nor unobserved outcomes. If the missingness is independent of the unobserved outcomes after conditioning on the observed outcomes, then the missing data mechanism is considered to be MAR. Conversely, under MNAR, conditioning on the observed outcomes cannot disentangle the dependence between the missingness process and the unobserved outcomes.

What separates MNAR from its counterparts in applied practice, though, is the concept of ignorability. That is under the (Bayesian) likelihood framework, one does not need to consider a model for the missingness process and hence ‘ignore’ it. Ignorability holds under the assumption that the process that generates the missing data is MAR and that the missingness and measurement processes depend on independent sets of parameters. While the innate convenience of ignorability popularized MAR models, one needs to consider if such an assumption holds carefully. That is especially the case since the distinction between MAR and MNAR is intractable as for every MNAR model, there exists a MAR counterpart with the same fit to the observed data (Molenberghs et al.(2008)Molenberghs, Beunckens, Sotto, and Kenward). Hence, one should start with an analysis that is valid under MAR and then perform a sensitivity analysis considering plausible MNAR departures from it.

To showcase the importance of such sensitivity analysis, consider the illustrative example in Figure 1 which is based on simulated data. The data were simulated assuming a joint model for longitudinal and survival data. To mimic MAR missingness, censoring was based on the value of previous repeated measurements, and therefore the missing longitudinal trajectories (depicted with dotted grey) consist of the subjects that fulfilled the MAR criterion or experienced the event. We then fitted a model to the complete data, whereas for the incomplete data we fitted an MAR and an MNAR model. All three models resulted in different fits to the data, which consequently lead to different subject-specific predictions as shown in the lower panel of Figure 1. Therefore, sensitivity analysis is necessary to be able to investigate the robustness of the model and the plausibility of the assumptions concerning the mechanism that generated the missing data.

To this purpose, three well-established general frameworks exist, distinguished by their approach in factoring the joint distribution of the outcome and missingness processes to conditionally independent components. In a pattern-mixture model (PMM), the joint distribution is expressed as the product of the marginal distribution of the missingness process and the conditional distribution of the outcome given the missingness process. Contrariwise, in a selection model (SeM), the joint distribution is expressed as the product of the marginal distribution of the outcome and the conditional distribution of the missingness process given the outcome. Alternatively, in a shared-parameter model (SPM), the outcome and the missingness processes are assumed to be driven by a set of unobserved variables, the so-called random effects. The use of latent variables gives rise to the conditional independence assumption under which the outcome and the missingness processes are independent, given the set of random effects.

The SeM framework naturally encompasses MAR, while (Molenberghs et al.(1998)Molenberghs, Michiels, Kenward, and Diggle) provided a characterization of MAR under the PMM setting, specifically for the case of longitudinal data with dropout. As a result, sensitivity analysis is well-developed in the SeM and the PMM frameworks (Molenberghs and Vereke (2005); Molenberghs and Kenward (2007)), whereas the SPM has received less attention. More specifically, (Creemers et al.(2011)Creemers, Hens, Aerts, Molenberghs, Verbeke, and Kenward) provided a characterization of MAR for the SPM case by introducing the general SPM (GSPM) and then proposed a sensitivity analysis approach utilizing this framework (Creemers et al.(2010)Creemers, Hens, Aerts, Molenberghs, Verbeke, and Kenward). Later, (Njagi et al.(2014)Njagi, Molenberghs, Kenward, Verbeke, and Rizopoulos) proposed a characterization of MAR for the case of SPMs for longitudinal and time-to-event data under the GSPM.

In this article, we introduce an alternative characterization of MAR in SPMs for longitudinal and time-to-event data by exploiting the conditional independence assumption and utilizing the censoring distribution. By doing so, we cover a wide range of assumptions for the missing data mechanism on the subject-specific level. Previous work on the topic focused on the GSPM, which includes potentially a broad set of random effects. Here, we focus on the standard SPM (Wu and Carroll (1988); Wu and Bailey (1988); Wu and Bailey (1989); Hogan and Laird (1997); Hogan and Laird (1998)), a particular case of the GSPM, that is most commonly encountered in practice. The benefit of doing so is that our approach is intuitively appealing and can be easily implemented in existing software. Furthermore, it allows for subject-specific perturbations of the censoring distribution, concerning the missing data generating mechanism, whereas in PMMs and SeMs dropout is considered MNAR strictly.

The rest of the article is structured as follows. In Section 2, we propose an alternative characterization of MAR for the SPM model. In Section 3, we discuss its estimation under the Bayesian framework. In Section 4, we illustrate with an application how the SPM can be used for sensitivity analysis with existing software. In Section 5, we present a simulation study to assess the behaviour of the proposed sensitivity analysis under different scenarios with respect to the amount of MAR and MNAR missingness. Discussion follows in Section 6.

Figure 1:

Illustrative Example: (Upper panel) observed longitudinal trajectories (black lines), missing longitudinal trajectories (grey), and fitted models using complete data and using incomplete data under MAR and MNAR. (Lower Panel) observed longitudinal measurements of a new subject and subject-specific predictions based on the three different models

2 A characterization of MAR for the shared parameter model

For each subject $i$ out of a sample of size $N$ from the target population, let ${\mathbf{y}}_i={\left(y_{i1},\dots ,y_{{\mathit{in}}_i}\right)}^{\top }$ denote the vector of size $n_i\times 1$ of planned repeated measurements, with $y_{\mathit{ij}}$ being the value of the longitudinal outcome at time point $t_{\mathit{ij}}$ , $j=1,\dots ,n_i$ , out of a set of $\mathcal{J}$ planned measurement occasions. The vector of planned repeated measurements may further be partitioned into two vectors: ${\mathbf{y}}_i=\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\right)$ , of observed and missing repeated measurements, respectively. Next, let $T_i^{*}$ be the true dropout time, $C_i$ the censoring time and $T_i=\mathit{min}\left(T_i^{*},C_i\right)$ the observed dropout or censoring time. Furthermore, let $\mathbf{\theta }$ and $\mathbf{\psi }=\left\{{\mathbf{\psi }}^{T^{*}},{\mathbf{\psi }}^C\right\}$ be the vectors of parameters for the measurement, dropout and censoring processes respectively. Under these settings, the vectors of observed and missing repeated measurements can be expressed as ${\mathbf{y}}_i^o=\left\{y_{\mathit{ij}};t_{\mathit{ij}}<T_i,\forall j\in \mathcal{J}\right\}$ and ${\mathbf{y}}_i^m=\left\{y_{\mathit{ij}};t_{\mathit{ij}}>T_i,\forall j\in \mathcal{J}\right\}$ . In other words, the vector of observed measurements includes all the observations which were recorded before neither dropout nor censoring occurred. On the contrary, the vector of missing measurements includes the observations which were not recorded due to either dropout or censoring. Note that we regard dropout and censoring as different causes of missingness, a distinction that we will exploit to achieve MAR characterization of SPMs.

What differentiates dropout from censoring in such a setting pertains to the causes of leaving the study. As an example, consider the case of a study, with death as the main outcome and where a longitudinal outcome is also planned to be recorded for a time period. The planned longitudinal outcome measurements will be missing for both the subjects that died during the study but also for subjects that dropped out from the study due to any other reason (e.g., they moved to another country). Depending on the information available and the research setting, one may define death as the only cause of dropout and assume all the other reasons are non-informative and therefore record them as censoring. In the case where no information is available for the cause of dropout, then one would consider only those who completed the study as censored.

Under this framework and with covariate vectors suppressed from the notation, the shared parameter model can be specified as follows:

$$p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m,T_i^{*},C_i;\mathbf{\theta },\mathbf{\psi }\right)=\int p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m,T_i^{*},C_i,{\mathbf{b}}_i;\mathbf{\theta },\mathbf{\psi },\mathbf{D}\right)d{\mathbf{b}}_i,$$

(2.1)

where ${\mathbf{b}}_i\sim \mathcal{N}\left(0,\mathbf{D}\right)$ is a vector of shared random effects assumed to induce the dependence between the processes involved in the joint density of the SPM. That is by conditioning on the random effects, the joint density of the SPM can be decomposed to the following product of conditional densities

$$\int p\left(T_i^{*},\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)p\left(C_i\mid {\mathbf{y}}_i^o,{\mathbf{y}}_i^m,{\mathbf{b}}_i;{\mathbf{\psi }}^C\right)p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,$$

(2.2)

which brings to light the core assumption of conditional independence in SPMs, under which the measurement and the dropout processes are independent conditional on the random effects. Note that the conditional independence assumption commonly does not imply independence of the censoring and the measurement processes. More specifically, in the general case, censoring is allowed to depend on the measurement process. Note that the conditional independence assumption may be extended to allow for independence of the censoring and the measurement processes as well. Finally, the assumption of non-informative censoring, which we discuss later, may be used to relax both these assumptions and allow for independence between the censoring and measurement processes.

Model (2.2) is a conventional SPM for longitudinal and time-to-event data. The term conventional here is in reference to the GSPM framework in the sense that a single common underlying random effects structure is used instead of multiple random effects structures. Another important note for model (2.1) is that by definition it allows for dependence between the dropout process $T_i^{*}$ and the censoring process $C_i$ . When there is no information from the data on the joint distribution of $\left(T_i^{*},C_i\right)$ , the assumption of non-informative censoring may be utilized (Tsiatis (1975)) to further simplify the model. Under both the assumptions of non-informative censoring and conditional independence, (2.2) further simplifies to the following:

$$\int p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)p\left(C_i\mid {\mathbf{y}}_i^o;{\mathbf{\psi }}^C\right)p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,$$

(2.3)

which, in other words, means that the censoring process no longer depends on the random effects and the missing observations. It should be noted that the assumption of non-informative censoring, while commonly applicable, may not always be a reasonable choice depending on the setting under study. Nevertheless, it is a necessary condition for achieving an MAR characterization in SPMs without introducing more latent variables as in (Creemers et al.(2011)Creemers, Hens, Aerts, Molenberghs, Verbeke, and Kenward).

Note that by definition, a subject may either be classified as a dropout or as censored and never both. This means that the decomposition of the joint density as described in (2.3) is, on a subject-specific level, further decomposed to the following:

$$\left\{\begin{array}{cc}\int p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,\hfill & i:\mathit{dropout}\hfill \\ \int p\left(C_i\mid {\mathbf{y}}_i^o;{\mathbf{\psi }}^C\right)p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,\hfill & i:\mathit{censored}.\hfill \end{array}\right.$$

(2.4)

The decomposition in (2.4) reveals that SPMs encompass MAR on the subject-specific level. That is for a subject who dropped out from the study, the model is MNAR because both the dropout, $p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)$ , and measurement, $p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)$ , processes depend on the missing observations through the random effects ${\mathbf{b}}_i$ . For a censored subject, though, the model is MAR since the random effects and the missing observations appear only in factor $p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)$ . Moreover, due to ignorability, the joint density reduces to the density of the observed measurement process $\int p\left({\mathbf{y}}_i^o\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i\right)d{\mathbf{b}}_i$ . Ignorability is a key consequence of MAR in SPMs. It entails that there are not any information to be gained from either the dropout or the censoring processes with respect to the measurement generating process, and therefore MAR can be explored within the SPM framework. This is achieved without broadening the definition of the conventional SPM by introducing additional random effects. Furthermore, the fact that MAR in SPMs is achieved on the subject-specific level allows for the potential of a subject- or cause-specific MAR sensitivity analysis. This is in contrast to the SeM and PMM frameworks which are MNAR or MAR for all the individuals under study. Finally, as we are going to illustrate in the upcoming section, this decomposition allows exploiting existing software for its estimation.

3 Estimation of SPM under MAR

Estimation of the SPM under MAR follows from the conventional joint model for longitudinal and time-to-event data, which exploits the decomposition of the full joint-likelihood function to conditional independent components given the random effects. More specifically, the likelihood contribution of the $i\mathrm{th}$ subject is given by the following:

$$\int p\left(T_i^{*}\mid {\mathbf{y}}_i,{\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)p\left(C_i\mid {\mathbf{y}}_i;{\mathbf{\psi }}^C\right)p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,$$

(3.1)

The form of the conditional densities depends then on the model specification of each component. Let the risk of dropout for subject $i$ depend on a function of the true underlying subject-specific trajectory ${\mathbf{\eta }}_i\left(t\right)$ . Then the instantaneous risk of dropout can be expressed as follows:

$$\begin{array}{cc}{h}_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i,{\mathbf{w}}_i\right)\hfill & =\underset{\Delta t\to 0}{\text{lim}}\frac{1}{\Delta t}\mathit{Pr}\left\{t\le T_i^{*}<t+\Delta t\mid T_i^{*}\ge t,{\mathcal{H}}_i\left(t\right),{\mathbf{w}}_i\right\}\hfill \\ \hfill & =h_0\left(t\right)\text{exp}\left[{\mathbf{\gamma }}^{\top }{\mathbf{w}}_i+f\left\{{\eta }_i\left(t\right),{\mathbf{b}}_i,\mathbf{\alpha }\right\}\right],t>0,\hfill \end{array}$$

(3.2)

where ${\mathcal{H}}_i\left(t\right)=\left\{{\eta }_i\left(s\right),0\le s\le t\right\}$ denotes the history of the true underlying longitudinal process up to time $t$ (and therefore depends on the random effects as well), $h_0\left(·\right)$ is the baseline hazard function and ${\mathbf{w}}_i$ is a vector of baseline covariates with corresponding regression coefficients $\mathbf{\gamma }$ . The association between the features of the longitudinal measurement process, $p\left({\mathbf{y}}_i^o,{\mathbf{y}}_i^m\mid {\mathbf{b}}_i;\mathbf{\theta }\right)$ , and the instantaneous risk for dropout, $h_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i,{\mathbf{w}}_i\right)$ , is quantified by parameter vector $\mathbf{\alpha }$ . Function $f\left(·\right)$ may take various forms but for the remainder of the article we will assume that $f\left({\eta }_i\left(t\right),{\mathbf{b}}_i,\alpha \right)=\alpha {\eta }_i\left(t\right)$ , the so-called current value association, which postulates that the risk of dropout at time $t$ depends only on the current value of the true underlying subject-specific trajectory at the same time ${\mathbf{\eta }}_i\left(t\right)$ . Let ${\mathbf{u}}_i$ be a vector of baseline covariates with corresponding regression coefficients ${\mathbf{\gamma }}_C$ and $g\left({\mathbf{y}}_i^o\right)$ a function of the observed measurements. Then, analogously, the instantaneous risk to be censored can be expressed as follows:

$$\begin{array}{cc}{\lambda }_i\left(t\mid {\mathbf{y}}_i^o,{\mathbf{u}}_i\right)\hfill & =\underset{\Delta t\to 0}{\text{lim}}\frac{1}{\Delta t}\mathit{Pr}\left\{t\le C_i<t+\Delta t\mid C_i\ge t,{\mathbf{y}}_i^o,{\mathbf{u}}_i\right\}\hfill \\ \hfill & =h_0\left(t\right)\text{exp}\left[{\mathbf{\gamma }}_C^{\top }{\mathbf{u}}_i+g\left({\mathbf{y}}_i^o\right)\right],t>0,\hfill \end{array}$$

(3.3)

since based on the assumption of non-informative censoring, it cannot depend on the missing observations and the random effects.

The probability of not dropping out from the study can then be obtained by the following:

$$\begin{array}{cc}{S}_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i,{\mathbf{w}}_i\right)\hfill & =\mathit{Pr}\left(T_i^{*}>t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{w}}_i\right)\hfill \\ \hfill & =\text{exp}\left(-{\int }_0^t{h}_0\left(s\right)\text{exp}\left\{{\mathbf{\gamma }}^{\top }{\mathbf{w}}_i+\alpha {\eta }_i\left(s\right)\right\}\mathit{ds}\right)\hfill \end{array}$$

(3.4)

which means that unlike the instantaneous risk of dropout, the probability of remaining in the study depends on the whole true underlying subject-specific trajectory up to time $t$ . Let ${\mathcal{Y}}_i^o\left(t\right)=\left\{{\mathbf{y}}_i^o\left(l\right),0\le l\le t\right\}$ denote the history of the observed repeated measurements.Then, analogously, the probability of not being censored can be expressed as follows:

$$\begin{array}{cc}{V}_i\left(t\mid {\mathcal{Y}}_i^o\left(t\right),{\mathbf{u}}_i\right)\hfill & =\mathit{Pr}\left(C_i>t\mid {\mathcal{Y}}_i^o\left(t\right),{\mathbf{u}}_i\right)\hfill \\ \hfill & =\text{exp}\left(-{\int }_0^t{\lambda }_0\left(s\right)\text{exp}\left\{{\mathbf{\gamma }}_C^{\top }{\mathbf{u}}_i+g\left({\mathbf{y}}_i^o\right)\right\}\mathit{ds}\right)\hfill \end{array}$$

(3.5)

omitting the conditioning on missing values and the random effects due to the non-informative censoring assumption. Finally, to distinguish between censored subjects and dropouts, let ${\delta }_i$ be an indicator variable which takes the value 1 if a subject dropped out from the study, that is, $T_i^{*}<C_i$ and 0 otherwise.

Under these assumptions and suppressing covariates ${\mathbf{w}}_i$ and ${\mathbf{u}}_i$ from the notation, the conditional density of the dropout process can be expressed as follows:

$$p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)=\left\{\begin{array}{cc}{h}_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)S_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right),\hfill & \hfill {\delta }_i=1\\ S_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right),\hfill & \hfill {\delta }_i=0,\end{array}\right.$$

(3.6)

while the conditional density of the censoring process is given by the following:

$$p\left(C_i\mid {\mathbf{y}}_i^o;{\mathbf{\psi }}^C\right)=\left\{\begin{array}{cc}{\lambda }_i\left(t\right)V_i\left(t\right),\hfill & \hfill {\delta }_i=0\\ V_i\left(t\right),\hfill & \hfill {\delta }_i=1,\end{array}\right.$$

(3.7)

Depending on the dropout and censoring status of the $i\mathrm{th}$ individual, the joint-likelihood in (3.1) can be written as follows:

$$\left\{\begin{array}{cc}{V}_i\left(t\right)\int h_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)S_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,\hfill & {\delta }_i=1\hfill \\ V_i\left(t\right){\lambda }_i\left(t\right)\int S_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,\hfill & {\delta }_i=0.\hfill \end{array}\right.$$

(3.8)

Under (3.8), the model is MNAR for both uncensored and censored subjects. The latter is because even for censored subjects, the term $S_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)$ still appears in the likelihood decomposition. However, for the case that a subject is censored, this term equals to 1 reducing (3.8) to the following:

$$V_i\left(t\right){\lambda }_i\left(t\right)\int p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,$$

(3.9)

which due to ignorability further reduces to the following:

$$\int p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)d{\mathbf{b}}_i,$$

(3.10)

an MAR model for the observed longitudinal measurements. The transition from equation (3.9) to (3.10) holds due to the fact that the instantaneous risk to be censored ${\lambda }_i\left(t\right)$ and the probability of not being censored $V_i\left(t\right)$ are ignorable, since they strictly depend on the history of the observed repeated measurements and not the missing ones (see equations (3.3) and (3.5)). Then when a subject is censored, instead of dropout, he/she is considered MAR allowing for straightforward sensitivity analysis since it can be achieved with existing software just by doing a data manipulation step and changing all the dropout indicators to zero.

As far as the true underlying longitudinal trajectory is concerned, let ${\mathbf{y}}_i\left(t\right)$ be the contaminated, with measurement error ${\mathbf{\epsilon }}_i$ , observed values of the true longitudinal trajectory, ${\mathbf{\eta }}_i$ . Then the longitudinal process can be expressed using a mixed-effects model defined as follows:

$$\left\{\begin{array}{c}{\mathbf{y}}_i\left(t\right)={\mathbf{\eta }}_i+{\mathbf{\epsilon }}_i\left(t\right),\hfill \\ {\mathbf{\eta }}_i\left(t\right)=E\left[{\mathbf{y}}_i\left(t\right)\mid {\mathbf{b}}_i\right]={\mathbf{x}}_i^{\top }\left(t\right)\mathbf{\beta }+{\mathbf{z}}_i^{\top }\left(t\right){\mathbf{b}}_i,\hfill \\ {\mathbf{\epsilon }}_i\sim \mathcal{N}\left(0,{\mathbf{\Sigma }}_i\right),\hfill \\ {\mathbf{b}}_i\sim \mathcal{N}\left(0,\mathbf{D}\right),\hfill \\ {\mathbf{b}}_i\perp {\mathbf{\epsilon }}_i,\hfill \end{array}\right.$$

(3.11)

where the true longitudinal trajectory ${\mathbf{\eta }}_i$ is assumed to be a function of some general time-dependent design vectors ${\mathbf{x}}_i^{\top }\left(t\right)$ and ${\mathbf{z}}_i^{\top }\left(t\right)$ which are associated with a set of fixed-effects $\mathbf{\beta }$ and a set of subject-specific random effects ${\mathbf{b}}_i$ , respectively. The errors are typically assumed to be normally distributed with mean zero and a general variance—covariance matrix ${\mathbf{\Sigma }}_i$ which is commonly specified as $\mathbf{I}\times {\sigma }^2$ . The random effects ${\mathbf{b}}_i$ are also assumed to be normally distributed with mean zero and a general variance—covariance matrix $\mathbf{D}$ , and they are also assumed to be independent from the errors.

The conditional density for the longitudinal responses is then given by the following:

$$\begin{array}{cc}p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)={\left(2\pi {\sigma }^2\right)}^{n_i/2}\text{exp}\left\{-{\mathbf{y}}_i-{\mathbf{x}}_i\mathbf{\beta }-{\mathbf{z}}_i{\mathbf{b}}_i\text{2}/2{\sigma }^2\right\}\hfill & \hfill \end{array}$$

with $\Vert x\Vert ={\displaystyle \sum }_i{\left\{x_i^2\right\}}^{1/2}$ denoting the Euclidean vector norm. The density of the random effects is, then, given by the following:

$$p\left({\mathbf{b}}_i;\mathbf{D}\right)={\left(2\pi \right)}^{-q/2}\mathit{det}{\left(\mathbf{D}\right)}^{-1/2}\mathit{exp}\left(-{\mathbf{b}}_i^{\top }{\mathbf{D}}^{-1}{\mathbf{b}}_i/2\right)$$

(3.12)

Under the Bayesian framework, estimation of the parameters of the SPM is achieved using Markov chain Monte Carlo (MCMC) algorithms. Let $\mathbf{\Theta }=\left\{\mathbf{\beta },{\mathbf{\Sigma }}_i,{\mathbf{\psi }}^{T^{*}},{\mathbf{\psi }}^C,\mathbf{D}\right\}$ be a vector of unknown parameters, then the posterior distribution is proportional to the following:

$$p\left(\mathbf{\Theta },\mathbf{b}\right)\prod _{i}p\left({\mathbf{y}}_i\mid {\mathbf{b}}_i;\mathbf{\theta }\right)p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)p\left(C_i\mid {\mathbf{y}}_i^o;{\mathbf{\psi }}^C\right)p\left({\mathbf{b}}_i;\mathbf{D}\right)p\left(\mathbf{\Theta }\right)$$

(3.13)

The computation of the integrals involved in $p\left(T_i^{*}\mid {\mathbf{b}}_i;{\mathbf{\psi }}^{T^{*}}\right)$ and $p\left(C_i\mid {\mathbf{y}}_i^o;{\mathbf{\psi }}^C\right)$ is achieved using standard Gauss–Kronrod and Gauss–Legendre quadrature rules. For the parameters in $\mathbf{\Theta }$ , standard prior distributions are used. More specifically, for the vector of regression coefficients $\mathbf{\beta }$ of the longitudinal submodel and for the vector of regression coefficients of the dropout submodel $\mathbf{\gamma }$ we use independent univariate diffuse normal priors. For the covariance matrix of the random effects an inverse-Wishart prior is used, while for the variance of the error terms ${\sigma }^2$ an inverse-Gamma prior is used. For the baseline hazard, we use a B-spline approximation expressed as follows:

$$h_0\left(t\right)=\text{exp}\left({\gamma }_{h_0,0}+\sum _{q=1}^Q{\gamma }_{h_0,q}{B}_q\left(t,k\right)\right),$$

(3.14)

with $B_q\left(t,k\right)$ denoting the $q\mathrm{th}$ basis function of a B-spline with knots $k_1,\dots ,k_q$ and ${\mathbf{\gamma }}_{h_0}$ the vector of coefficient corresponding to the spline terms. To achieve the B-spline approximation to the baseline hazard under the Bayesian framework, the following prior is specified:

$$p\left({\mathbf{\gamma }}_{h_0}\mid {\tau }_h\right){\tau }_h^{\rho \left(K\right)/2}\text{exp}\left(-\frac{{\tau }_h}{2}{\mathbf{\gamma }}_{h_0}^{\top }\mathbf{K}{\mathbf{\gamma }}_{h_0}\right)$$

(3.15)

In the above specification, ${\tau }_h$ is a smoothing parameter with a $\mathit{Gamma}\left(1,0.005\right)$ hyper-prior, $\mathbf{K}$ is a penalty matrix and $\rho (\mathbf{K})$ is its rank. For more details, the reader is referred to (Rizopoulos (2016)).

4 The HIV CD4 data: A sensitivity analysis

Figure 2:

Observed longitudinal trajectories of $\sqrt{\mathit{CD}4}$ by treatment arm

As an illustrative case study, we will use data from a randomized clinical trial designed to compare the efficacy and safety of Didanosine (ddI) versus Zalcitabine (ddC) in HIV patients (Abrams et al.(1994)Abrams, Goldman, Launer, Korvick, Neaton, Crane, Grodesky, Wakefield, Muth, Kornegay, Cohn, Harris, Luskin-Hawk, Markowitz, Sampson, Thompson, and Deyton; Goldman et al.(1996)Goldman, Carlin, Crane, Launer, Korvick, Deyton, and D.). In total, 467 advanced HIV patients were included in the study and were randomized to ddI (230; 49.2%) or ddC (237; 50.8%). The primary goal of the trial was to compare survival between the two treatment arms, while a secondary goal of the study was to investigate the association between CD4 cell count (a marker for the strength of the immune system) and the risk of death. For the latter goal, measurements of the CD4 cell count were recorded at scheduled visits at baseline, 2, 6, 12 and 18 months after the start of the study. Figure 2 shows the observed trajectories of the square root CD4 cell count. Until the end of the study, 184 (39.4%) died while only 1 405 (60.1%) out of the 2 335 planned visits were recorded, leading to 930 (39.9%) missing observations. In total, only 24(5.1%) patients were present for all five planned visits, 382(81.8%) dropped out from the study due to death or other reasons, and 61(13.0%) missed at least one planned visit without dropping out though from the remainder of the study. Table 1 shows the number of available measurements per dropout pattern per treatment arm, excluding the 61 cases of intermittent missingness.

Table 1:

Number of observed CD4 measurements per dropout pattern and per treatment arm. ‘OXXXX’ denotes a dropout pattern where only the first planned measurement was obtained (denoted with ‘O’) while the rest were not (denoted with ‘X’), whereas ‘OOOOO’ denotes a pattern for a completer (all planned measurements were obtained)

	ddC		ddI
Dropout Pattern	N	%	N	%
OXXXX	29	14.4	32	15.6
OOXXX	35	17.4	37	18.0
OOOXX	41	20.4	47	22.9
OOOOX	85	42.3	76	37.1
OOOOO	11	5.5	13	6.3
Total	201	100	205	100

We observe that dropout rates are similar between the treatment groups, the closer we look at the beginning of the study, but they slightly start to deviate towards the end of the study. Nevertheless, it is not possible to distinguish if the missing data are missing completely at random, MAR or MNAR. This highlights the importance of being able to explore different scenarios with respect to missing data in the SPM framework via sensitivity analysis.

To model the longitudinal measurements of CD4 cell count, we used a linear mixed-effects model with an interaction between time and ddI in the fixed-effects structure. We also used random intercepts and random slopes to allow for subject-specific deviations. The square root of the CD4 cell count was used in order to meet the normality and homoscedasticity assumptions of the model, which is defined as follows:

$$\sqrt{y_{\mathit{ij}}^o}=\left({\beta }_0+b_{i0}\right)+\left({\beta }_1+b_{i1}\right)t_{\mathit{ij}}+{\beta }_2\mathit{ddI}+{\beta }_3\left(t_{\mathit{ij}}\times \mathit{ddI}\right)+{\epsilon }_{\mathit{ij}},$$

(4.1)

where ${\mathbf{b}}_i\sim \mathcal{N}\left(0,\mathbf{D}\right)$ , independent of ${\epsilon }_{\mathit{ij}}\sim \mathcal{N}\left(0,{\sigma }^2\right)$ .

For the time-to-dropout $T_i$ , we assumed a relative risk model of the form:

$$h_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)=h_0\left(t\right)\text{exp}\left[{\gamma }_1\mathit{ddI}+\mathbf{\alpha }{\eta }_i\left(t\right)\right]$$

(4.2)

where the baseline hazard $h_0\left(t\right)$ was approximated using penalized B-splines with three internal knots placed at quantiles of the observed event times (Eilers and Marx (1996)). The number of knots was kept small due to the low number of time points per subject.

These two sub models define the MNAR joint model for the CD4 cell count and the dropout. We then assumed two different scenarios with respect to the cause of dropout. Under Scenario 1, we assume no information concerning the cause of dropout, which means we treat dropout due to death or any other reason the same. Contrariwise in Scenario 2, we consider as dropouts only the subjects who died during the study. Finally, to achieve MAR, we assume that all the subjects are censored instead of dropping out. As shown in Section 3, this reduces the model to the mixed-effects sub model, which is considered MAR. All models were fitted using R version 3.6.1 and package JMbayes (Rizopoulos (2016)). Table 2 shows the parameter estimates after fitting all three different models.

Table 2:

Posterior means (standard errors) for an MAR and MNAR analyses

	Scenario 1	Scenario 2
Effect	Parameter	MNAR	MNAR	MAR
Intercept	${\beta }_0$	2.4671 (0.0640)	2.4680 (0.0629)	2.4423 (0.0647)
Time	${\beta }_1$	$-$ 0.0655 (0.1701)	$-$ 0.0674 (0.1715)	$-$ 0.0399 (0.0050)
ddI	${\beta }_2$	0.1038 (0.0908)	0.1067 (0.0879)	0.1188 (0.0946)
ddI $\times$ Time	${\beta }_3$	0.0219 (0.2332)	0.0290 (0.2376)	0.0089 (0.0070)

While the results of all the models are close, there are quantitative differences which indicate that the findings might not be stable. This is especially the case when comparing the MAR model with both MNAR models. In Figures 2 and 4, the posterior means of the random effects estimated from the MAR model against the respective estimated random effects from each MNAR model are shown. Similarly, in Figure 5, the posterior means of the random effects estimated from the MNAR model under the 1st Scenario are plotted against the respective estimated random effects from the MNAR model under the 2nd Scenario. These plots give insight into the differences between the three models. We see that especially for the slope components, the differences in the random effects’ estimates are more intense between the MAR model and the MNAR models under both scenarios. These differences become weaker when we compare the MNAR model under the 1st Scenario and the MANR model under the 2nd Scenario, which utilizes the information and distinguishes between dropout due to an event and dropout due to any other reason. This highlights the importance of sensitivity analysis in this context since the differences in the random effects estimates would translate to differences in the subject-specific predictions based on each of the models. The results are stable under different MNAR scenarios, but substantial differences are observed when using an MAR model.

Figure 3:

Scenario 1: Posterior means of the random effects under MNAR and MAR

Figure 4:

Scenario 2: Posterior means of the random effects under MNAR and MAR

Figure 5:

Scenario 2: Posterior means of the random effects under MNAR (1st Scenario) and MNAR (2nd Scenario)

5 Simulation study

We conducted a simulation study to evaluate the performance of the following three models (used in the analysis of the HIV CD4 data):

Model 1: all dropout cases considered MNAR (as the MNAR model used under Scenario 1 MNAR in the analysis of the HIV dataset),

Model 2: all dropout cases considered MAR,

Model 3: dropout cases considered MNAR if dropout due to the event or MAR otherwise (as the MNAR model used under Scenario 1 MNAR in the analysis of the HIV dataset),

under different scenarios for the amount of MAR and MNAR dropouts. More specifically, assuming a total dropout rate of 50%, we considered three different scenarios concerning the amount of dropout type: (a) 5% MAR dropouts versus 45% MNAR dropouts, (b) 25% MAR dropouts versus 25% MNAR dropouts and (c) 45% MAR dropouts versus 5% MNAR dropouts. Our goal is to investigate the behaviour between these models under each scenario.

We assumed 600 subjects and then randomly selected follow-up visits, $t_{\mathit{ij}}$ from a uniform distribution between 0 and 10. We set the total number of measurements per subject to 10, but the final number of measurements may vary depending on whether a subject was a dropout or not. The type of dropout was then determined as follows: if the value of the longitudinal trajectory of a subject exceeded a pre-specified value, then the next time point was selected as the candidate MAR dropout time, which was then compared to the candidate MNAR dropout time. We simulated the MNAR dropout time using a joint model. Depending on which of these time points occurred sooner (if occurred), the subject was classified as MAR dropout, MNAR dropout or completer. The threshold value for MAR varied between the three scenarios in order to attain the target percentages of dropout per type.

For the continuous longitudinal outcome, the data were simulated from a linear mixed-effects model as follows:

$$y_{\mathit{ij}}=\left({\beta }_0+b_{i0}\right)+\left({\beta }_1+b_{i1}\right)t_{\mathit{ij}}+{\beta }_2\mathit{Xcov}+{\beta }_3\left(\mathit{Xcov}\times t_{\mathit{ij}}\right)+{\epsilon }_{\mathit{ij}},$$

(5.1)

where $\mathit{Xcov}$ is a categorical variable with two groups, ${\epsilon }_i\left(t\right)\sim \mathcal{N}\left(0,{\sigma }^2{I}_{n_i}\right)$ and $\mathbf{b}\sim \mathcal{N}\left(0,\mathbf{D}\right)$ . The parameter values used are as follows: ${\beta }_0=2.47$ , ${\beta }_1=-0.067$ , ${\beta }_2=0.107$ , ${\beta }_3=0.029$ , $\sigma =1.431$ and for the variance—covariance matrix of the random effects:

$$\mathbf{D}=\left[\begin{array}{cc}\hfill d_{11}=0.765\hfill & \hfill d_{12}=-0.038\hfill \\ \hfill d_{21}=-0.038\hfill & \hfill d_{22}=4.910\hfill \end{array}\right]$$

The values of the parameters were based on the analysis of the HIV CD4 data. For time-to-dropout, we assumed a relative risk model of the form:

$$h_i\left(t\mid {\mathcal{H}}_i\left(t\right),{\mathbf{b}}_i\right)=h_0\left(t\right)\text{exp}\left[{\gamma }_1\mathit{Xcov}+\mathbf{\alpha }{\eta }_i\left(t\right)\right]$$

(5.2)

where the baseline hazard was simulated from a Weibull distribution $h0\left(t\right)=\xi t^{\xi -1}$ with $\xi =3.765$ .

We simulated 500 datasets per scenario, and we fitted each of the three models on each dataset. To assess the behaviour of the three models, we calculated the Bias and Root Mean Squared Error (RMSE) for both the regression coefficients, and the variance parameter estimates. Table 3 summarizes the results for the regression coefficients, while Table 4 summarizes the results for the variance parameters. The results suggest that when the major cause of missingness is MAR, Models 2 and 3 perform similarly and slightly better than Model 1, as far as the regression coefficients are concerned. As the amount of MNAR dropout increases, the three models seem to have similar performance in the estimation of the regression coefficients with only slight differences.

On the other hand, when looking at the Bias and RMSE for the estimates of variance parameters for the random effects, there are more distinct differences between the models. More specifically, the results suggest that Models 1 and 3 perform better than Model 2. Moreover, for the variances of the random effects, under Model 2, the Bias seems to increase as the balance between dropout types moves from MAR dominance to MNAR dominance, whereas the case is the opposite for Model 3. The differences in the variance parameters’ estimates of the random effects between the models also imply that there are considerable differences in the estimated random effects. The difference in the random effects means that the subject-specific predictions derived from these models will differ.

6 Discussion

In this article, we have proposed an alternative characterization of MAR for the conventional SPM and proposed its application as a sensitivity analysis tool towards MNAR deviations from the MAR assumption. In doing so, we did not broaden the definition of the SPM by adding additional random effects structures and hence retaining the computational feasibility of the model.

Furthermore, we argued how the subject-specific nature of the MAR characterization in SPM comes with the advantage of more flexible comparisons regarding the causes of missingness. This is an important advantage over other MNAR frameworks such as the PMM and SeM since it allows distinguishing groups of subjects as MAR or MNAR depending on the information available for the causes of missingness. This was illustrated in the application to the HIV data, where we were able to consider different MNAR scenarios with respect to the information available on the causes of missingness.

Finally, it should be mentioned that this characterization of MAR in SPMs allows for potential extensions to account for multiple causes of missingness by using for example a multi-state model and/or generalized mixed-effects models for categorical longitudinal outcomes. Dropout is present in most longitudinal studies, and our approach showed how one could perform a sensitivity analysis using standard joint modelling software to explore different causes of missingness. By doing so, we do hope that such sensitivity analyses will become standard procedure and will be routinely reported in longitudinal studies with dropout.

Table 3:

Mean, Bias and RMSE for the regression coefficients from the three models under different scenarios for the amount of dropout caused by each missing data generating mechanism

	Scenario 1: 5% MNAR ’ 45% MAR			Scenario 2: 25% MNAR ’ 25% MAR			Scenario 3: 45% MNAR ’ 5% MAR
Parameter	Mean	Bias	RMSE	Mean	Bias	RMSE	Mean	Bias	RMSE
Model 1
${\beta }_0$	2.4532	$-$ 0.0147	0.0753	2.4667	$-$ 0.0012	0.0715	2.4705	0.0025	0.0711
${\beta }_1$	0.0166	0.0840	0.1477	$-$ 0.0088	0.0585	0.1339	$-$ 0.0129	0.0544	0.1322
${\beta }_2$	$-$ 0.0091	$-$ 0.1158	0.1819	$-$ 0.0109	$-$ 0.1176	0.1827	$-$ 0.0114	$-$ 0.1181	0.1821
${\beta }_3$	0.1056	0.0766	0.1045	0.1045	0.0755	0.1002	0.1048	0.0758	0.0993
Model 2
${\beta }_0$	2.4688	0.0009	0.0738	2.4707	0.0028	0.0716	2.4723	0.0043	0.0712
${\beta }_1$	$-$ 0.0136	0.0538	0.1322	$-$ 0.0127	0.0547	0.1322	$-$ 0.0139	0.0535	0.1319
${\beta }_2$	$-$ 0.0106	$-$ 0.1174	0.1808	$-$ 0.0107	$-$ 0.1174	0.1825	$-$ 0.0116	$-$ 0.1183	0.1823
${\beta }_3$	0.1060	0.0770	0.1043	0.1051	0.0761	0.1002	0.1046	0.0756	0.0994
Model 3
${\beta }_0$	2.4680	0.0001	0.0738	2.4687	0.0007	0.0717	2.4694	0.0014	0.0713
${\beta }_1$	$-$ 0.0129	0.0545	0.1325	$-$ 0.0115	0.0558	0.1328	$-$ 0.0121	0.0553	0.1323
${\beta }_2$	$-$ 0.0107	$-$ 0.1174	0.1810	$-$ 0.0105	$-$ 0.1172	0.1827	$-$ 0.0114	$-$ 0.1181	0.1820
${\beta }_3$	0.1060	0.0770	0.1042	0.1051	0.0761	0.1005	0.1047	0.0757	0.0995

Table 4:

Mean, Bias and RMSE for the variance parameters from the three models under different scenarios for the amount of dropout caused by each missing data generating mechanism

	Scenario 1: 5% MNAR ’ 45% MAR			Scenario 2: 25% MNAR ’ 25% MAR			Scenario 3: 45% MNAR ’ 5% MAR
Parameter	Mean	Bias	RMSE	Mean	Bias	RMSE	Mean	Bias	RMSE
Model 1
$d_{11}$	0.7577	$-$ 0.0074	0.0910	0.7397	$-$ 0.0254	0.0848	0.7394	$-$ 0.0257	0.0824
$d_{12}$	$-$ 0.0514	$-$ 0.0134	0.1206	$-$ 0.0419	$-$ 0.0040	0.1132	$-$ 0.0265	0.0115	0.1107
$d_{22}$	4.7665	$-$ 0.1439	0.3248	4.6891	$-$ 0.2214	0.3543	4.6768	$-$ 0.2337	0.3600
$\sigma$	1.4318	0.0003	0.0182	1.4314	0.0000	0.0160	1.4293	$-$ 0.0022	0.0157
Model 2
$d_{11}$	0.1406	$-$ 0.6245	0.6258	0.1381	$-$ 0.6270	0.6279	0.1381	$-$ 0.6270	0.6278
$d_{12}$	$-$ 0.0222	0.0158	0.0742	$-$ 0.0208	0.0171	0.0708	$-$ 0.0202	0.0178	0.0694
$d_{22}$	5.4293	0.5188	0.8668	5.4355	0.5251	0.8401	5.4484	0.5380	0.8477
$\sigma$	1.4296	$-$ 0.0018	0.0180	1.4302	$-$ 0.0013	0.0159	1.4301	$-$ 0.0014	0.0155
Model 3
$d_{11}$	0.7469	$-$ 0.0182	0.0908	0.7427	$-$ 0.0224	0.0841	0.7392	$-$ 0.0259	0.0829
$d_{12}$	$-$ 0.0271	0.0108	0.1184	$-$ 0.0321	0.0059	0.1133	$-$ 0.0343	0.0037	0.1098
$d_{22}$	4.6621	$-$ 0.2484	0.3817	4.6715	$-$ 0.2389	0.3643	4.6798	$-$ 0.2306	0.3576
$\sigma$	1.4299	$-$ 0.0015	0.0181	1.4300	$-$ 0.0014	0.0160	1.4299	$-$ 0.0016	0.0156