Local sensitivity to non-ignorability in joint models

Abstract

This article deals with the analysis of sensitivity to the non-ignorability of the dropout process in joint models (JMs). We investigate the behaviour of the maximum likelihood estimates for the longitudinal process in a neighbourhood of ignorability through the Index of Local Sensitivity to Non Ignorability (ISNI). Some concerns may arise since the ISNI is an absolute measure of change in parameter estimates induced by departures from the MAR assumption; for this reason, we introduce a relative index based on the ratio between the ISNI and a measure of its variability under the MAR assumption, highlighting potential interpretation and drawbacks of this approach. The local sensitivity of the JM and the performance of the relative index are discussed in a simulation study, by varying the number of repeated measurements per individual and the random effect covariance structure. The approach is also applied to a benchmark dataset on Primary Biliary Cirrhosis (PBC).

Keywords

Joint Models dropout ignorability sensitivity survival analysis

1 Introduction

Longitudinal studies often record two types of outcomes: repeated measurements for a response of interest and realizations of a time-to-event process describing individual participation to the study. An example can be derived from AIDS studies, where interest lies in the longitudinal evolution of markers such as the CD4 cells count, as well as in the time-to-AIDS or time-to-death. The longitudinal response may be directly related to the event process; after an event has occurred, longitudinal measurements may no longer be collected or considered relevant, therefore inducing dropout.

In many cases, the dropout process may depend on the unobserved values of the longitudinal response, corresponding to a not random missingness (Little and Rubin, 2002). Literature has so far concentrated on three modelling frameworks for the analysis of a longitudinal outcome subject to a (nonrandom) dropout process: Selection Models (SeMs), see e.g. (Diggle et al. 1994), Pattern Mixture Models (PMMs), (Little, 1993), and Shared Parameter Models (SPMs), (Follman and Wu, 1995). Even though these models can be quite flexible in the specification of the dropout mechanism, not too much trust must be placed on inferences due to potential sensitivity of parameter estimates to (unverifiable) modelling assumptions, see Little (1995), Rizopoulos et al. (2008a), Rizopoulos et al. (2008b) among others. Therefore, we need to investigate sensitivity of model parameter estimates to assumptions on the dropout mechanism; the majority of the proposals have focused on Selection and Pattern Mixture Models, see e.g. Troxel et al. (2004), Ma et al. (2005), Molenberghs (2005). Not too much have been done for the class of Shared Parameter Models, with the only exception of the recent work by Creemers et al. (2010).

In the present article, we consider a specific class of SPMs, referred to as joint models (JM), see Wulfsohn and Tsiatis (1997), where the longitudinal and the dropout processes may share a set of time-dependent covariates, as well as the corresponding fixed and random effects. In this framework, the time that a subject spends in the study is described by a proportional hazard model (Cox, 1972). Since the longitudinal evolution of the response of interest is measured with error, we assume that the risk of dropout is influenced by its (error-free) expected value, as well as by some baseline, time-invariant, covariates. As a sensitivity measure, we extend the Index of Local Sensitivity (ISNI) proposed by Troxel et al. (2004) and Ma et al. (2005), to joint models. The goal is to measure how much the maximum likelihood estimates are influenced by assumptions about the dropout mechanism. The main advantage of this approach is that it is based on quantities that can be calculated without currently fitting the full joint model.

By adapting the standard ISNI formulation to JMs, some concerns arise since ISNI is an absolute measure of change in parameter estimates induced by departures from the MAR assumption, see Troxel et al. (2004). We introduce a relative index based on the ratio between the ISNI and its standard error under MAR, and highlight potential interpretation and drawbacks of this formulation. To analyse the index behaviour in empirical applications and study the sensitivity of JM to non-ignorability, we review a benchmark dataset on primary biliary cirrhosis (PBC), see Murtaugh et al. (2002), where 312 patients have been considered to test for treatment effect on survival after adjusting for the longitudinal bilirubin levels. We also discuss the performance of absolute and relative indices through a simulation study where the number of subjects, the association structure between the longitudinal and the survival processes and the random effect covariance structure vary. This would help us analyze the proposed index behaviour when the model is correctly specified and when it is partially/globally misspecified.

The rest of the article is organized as follows: Section 2 describes the JM framework in details. Section 3 gives an insight on sensitivity analysis and introduces the ISNI formulation in JMs, while Section 4 entails the definition of the relative index. The sensitivity of JM and the empirical behaviour of the indices are discussed in Sections 5 for the PBC dataset and in Section 6 through a simulation study. Finally, Section 7 gives some concluding remarks and outlines potential developments.

2 The joint model

Let Y_i(t_ij) represent a longitudinal continuous response recorded at occasion j = 1, . . ., r_i for the ith subject (i = 1, . . ., n) and let y_i(t_ij) be the corresponding observed outcome. Throughout the article, we will assume that Y_i(t_ij) follows a linear mixed model of the form:

Y_{i} (t_{i j}) = β^{⊤} X_{i} (t_{i j}) + b_{i}^{⊤} Z_{i} (t_{i j}) + ϵ_{i} (t_{i j}),

(2.1)

where $ϵ_{i} (t_{i j}) ~ N (0, σ^{2})$ is the measurement error, $b_{i} ~ MVN (0, D)$ is a set of random coefficients, while X_i(t_ij) and Z_i(t_ij) are the rows of design matrices corresponding to fixed and random effects for subject i at time t_ij, respectively. We consider the case when X_i(t_ij) contains time-invariant covariates, the response measurement occasion time t_ij, j = 1, . . ., r_i, as well as the interaction of t_ijand some other covariates, while Z_i(t_ij) is a subset of X_i(t_ij), along the lines of Wulfsohn and Tsiatis (1997). The random terms b_i and e_i are assumed to be independent, since the first is time invariant and is shared by all subject i responses, while the latter represents unstructured, time varying, random variation from the ‘true’ individual signal. Often, longitudinal outcomes are not observed for the whole study period since individuals may leave the study before its end, potentially as a response to a secondary event. In this case, the occurrence of the event at a given time may induce dropout in the longitudinal outcome since no longitudinal measurements are collected afterwards. If the dropout is non random, i.e., if the probability of dropout, conditional on the observed individual characteristics, depends, either directly or indirectly, on the unobserved longitudinal responses, standard estimation procedures may lead to inconsistent parameter estimates. Formally, let $T_{i} = \min (T_{i}^{*}, C_{i})$ denote the observed failure time for the ith individual, taken as the minimum between the true event time $T_{i}^{*}$ and the censoring time C_i, and δ_i be the event indicator defined as $I (T_{i}^{*} \leq C_{i})$ . The longitudinal response y_i(t_ij) is therefore observed for $t_{i j} < T_{i}$ and missing for $t_{i j} \geq T_{i} .$ The observed failure time T_i may represent either a true event time or a censoring time. A clear distinction should be made between dropouts due to censoring events, where longitudinal responses are no longer available (e.g., as for death events with $T_{i}^{*}$ > C_i and δ_i = 0), and dropouts with secondary events, where longitudinal responses could have been registered should the subject have been participating to the study ( $T_{i}^{*} \leq C_{i}$ and δ_i = 1). We assume that the longitudinal process is associated with $T_{i}^{*},$ i.e., with the true event time, but is independent of C_i; in this sense, the censoring is assumed to be random. As far as the longitudinal process is concerned, t_ij represents the (potentially continuous) time scale for each observation, $t_{i j} \in (0, T_{i})$ . However, in the analysed examples, the longitudinal outcome is repeatedly observed in predefined occasions (usually the scheduled visits j = 1, . . ., r_i), while the dropout process includes the measurement of the actual event (or censoring) time.

Our aim is at investigating the effect of potentially informative dropouts on parameter estimates for the longitudinal model. JMs represent an appealing solution to account for the association between the longitudinal and the survival processes. The literature on JMs for longitudinal data dates back at least to Faucett and Thomas (1996) and Wulfsohn and Tsiatis (1997); since then several authors, see among others Rizopoulos et al. (2009) and Rizopoulos and Ghosh (2011), have proposed extensions and different estimation methods. A nice overview of existing methods is given by Tsiatis and Davidian (2004). In this framework, the repeated measurements and the time-to-event are assumed to share a set of time-varying covariates as well as the corresponding fixed and (subject-specific) random coefficients, the latter inducing dependence in the longitudinal profile as well as between the two processes. More specifically, the risk of experiencing the event at time T_i, described through a proportional hazard model, Cox, 1972, depends on the expected (error-free) value of the longitudinal outcome at the same time, treated as a time-dependent covariate with regression parameter α.

The joint model, in the formulation of Wulfsohn and Tsiatis (1997), was originally proposed to account for the effect of time dependent covariates measured with error on the survival process. With this purpose, the authors postulated a linear mixed model for the time dependent covariate and introduced the error-free part (i.e., the expected value of the longitudinal covariates) in the hazard function for the event of interest. Over the years, this model has became popular also in the context of longitudinal data with dropout, see for instance Tsiatis and Davidian (2004) and Rizopoulos (2010), to mention just a few, as a generalization of shared parameter models, Follmann and Wu (1995).

The model can be expressed by a set of two equations, one for the longitudinal and the other one for the survival part:

(\begin{array}{l} Y_{i} (t_{i j}) = m_{i} (t_{i j}) + ε_{i} (t_{i j}), j = 1, . . ., r_{i}, i = 1, . . ., n \\ h (t_{i j} | M_{i} (t_{i j}), W_{i}) = h_{0} (t_{i j}) \exp {γ^{⊤} W_{i} + α m_{i} (t_{i j})}, \end{array}

(2.2)

where $m_{i} (t_{i j}) = β^{⊤} X_{i} (t_{i j}) + b_{i}^{⊤} Z_{i} (t_{i j})$ is the expected value of the longitudinal process at the j-th occasion for the ith individual, $M_{i} (t_{i j}) = {m_{i} (u), 0 \leq u \leq t_{i j}}$ represents the history of m_i(.) up to time t_ij, W_i is a row vector of baseline covariates and h₀(t_ij) denotes the baseline risk function. As it can be easily noticed, model (2.2) can be seen as a time-dependent survival model, with time varying covariate Y_i(t_ij) introduced without the measurement error e_i (t_ij). Thus, we have two potential sources of non-ignorability: when α ࣔ 0 the two processes are not independent, since they share the same set of random coefficients. Furthermore, fixed effects β appear in both submodels and (fixed) parameter distinctiveness does not hold. The non-ignorability parameter α can be interpreted as a regression parameter in a proportional hazard model and quantifies the effect of the expected value of the longitudinal process on the hazard of the event as far as other baseline covariates are kept fixed. A significant α estimate denotes a dependence between the expected longitudinal response and the dropout process and may suggest a not random missing data mechanism.

Model (2.2) allows the dropout process to depend on covariates (the same or different than those affecting the longitudinal process) in the definition of the relative risk model, by including the covariate vector W_i. Moreover, not only individual-specific heterogeneity, but the overall longitudinal expected value of the response at a fixed time point influences the hazard of dropout at the same time. A similar procedure is adopted in selection models, where the dropout probability at a fixed time point is assumed to be influenced by the current (but unobserved) response at the same time; the JM approach, on the other hand, may be more robust, since the response in the survival model is error-free.

A relevant issue concerns the meaning of a (near) null estimate for α. This leads to a missing completely at random (MCAR) model, (Little and Rubin, 2002), that is, to the assumption that the dropout mechanism, conditional on available covariates, does not depend on the longitudinal response, either observed or missing. In fact when α = 0 model parameters in the two submodels are distinct and the joint probability for the dropout and the longitudinal process can be factorized as follows:

p (T_{i}, δ_{i}, Y_{i} (t_{i j})) = p (T_{i}, δ_{i}) p (Y_{i} (t_{i j})) .

The same factorization holds for the log-likelihood function with respect to the fixed effects in the longitudinal process. As we adopt a maximum likelihood approach to parameter estimation, we know that parameters estimates derived from maximizing the likelihood of the longitudinal process, that is p(Y_i(t_ij)), yield maximum likelihood estimates that are valid both under MCAR and under MAR assumptions, i.e., under the hypothesis that the dropout mechanism depends on the observed responses only. No setup may describe, within this modelling formulation, a proper MAR mechanism. However, what we are actually comparing is the results of the mixed model fitted alone to the data to the results of the joint model. In that respect and because the linear mixed model alone produces valid results under both MCAR and MAR, it makes sense to talk about MAR in our settings. When α ! 0, on the other hand, the dropout process depends on both the observed and the unobserved response through m_i (T_i) and parameter distinctiveness does not hold. In these cases, we are in presence of non-ignorabile dropout process.

The basic assumption of JMs is that of conditional independence, expressed by the following equations:

p (T_{i}, δ_{i}, y_{i} | b_{i}; θ) = p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ)

p (y_{i} | b_{i}; θ) = \prod_{j} p {y_{i} (t_{i j}) | b_{i}; θ} .

That is, conditionally on the random effects, the repeated measurements for the i-th individual are independent, and the same is true for the longitudinal and the survival processes. The random coefficients allow for within-individual dependence in the longitudinal process and for dependence between the longitudinal and the survival processes.

2.1 Score Vector and Hessian Matrix

In this section, we present the general formulation of the log-likelihood and corresponding score vector and Hessian matrix in the case of JMs, needed for the sensitivity analysis described in Section 3.

By indicating with θ the vector of all model parameters, the log-likelihood function for a JM can be written as follows:

ℓ (θ) = ℓ (θ | T_{i}, δ_{i}, y_{i}) = \sum_{i} \log \int_{b_{i}} p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ) d b_{i},

(2.3)

where $p (y_{i} | b_{i}; θ) = \prod_{j = 1}^{r_{i}} p (y_{i} (t_{i j}) | b_{i}; θ),$ and

p (T_{i}, δ_{i} | b_{i}; θ) = p {(T_{i} | b_{i}; θ)}^{δ_{i}} S {(T_{i} | M_{i} (T_{i}), W_{i}; θ)}^{1 - δ_{i}}

= h {(T_{i} | M_{i} (T_{i}), W_{i}; θ)}^{δ_{i}} S (T_{i} | M_{i} (T_{i}), W_{i}; θ) .

If a joint model parameterization is adopted, the survival function at time T_i is given by

S (T_{i} | M_{i} (T_{i}), W_{i}; θ) = \exp {- \int_{0}^{T_{i}} h_{0} (s) \exp {γ^{⊤} W_{i} + α m_{i} (s)} d s} .

(2.4)

The integral in (2.4) is solved using a Gauss-Kronrod numerical rule, see Kronrod (1964).

The score function, see e.g., Rizopoulos et al. (2009), takes the form:

s (θ) = \frac{\partial ℓ (θ)}{\partial θ} = \sum_{i} s_{i} (θ) = \sum_{i} \frac{\partial}{\partial θ^{⊤}} \log \int p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ) d b_{i}

= \sum_{i} \frac{1}{p (T_{i}, δ_{i}, y_{i}; θ)} \frac{\partial}{\partial θ^{⊤}} \int p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ) d b_{i}

= \sum_{i} \frac{1}{p (T_{i}, δ_{i}, y_{i}; θ)} \int \frac{\partial}{\partial θ^{⊤}} p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ) d b_{i}

= \sum_{i} \int [\frac{\partial}{\partial θ^{⊤}} \log {p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ)}]

\times \frac{p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i}; θ)}{p (T_{i}, δ_{i}, y_{i}; θ)} d b_{i}

= \sum_{i} \int q (θ, b_{i}) p (b_{i} | T_{i}, δ_{i}, y_{i}; θ) d b_{i},

(2.5)

where

q (θ, b_{i}) = \frac{\partial}{\partial θ^{⊤}} \log p (T_{i}, δ_{i}, y_{i}, b_{i}; θ)

= \frac{\partial}{\partial θ^{⊤}} {\log p (T_{i}, δ_{i} | b_{i}; θ + \log p (y_{i} | b_{i}; θ) + \log p (b_{i}; θ)} .

For an index pair (u, v), indicating the (u, v)-th component of the parameter vector θ, and the i-th individual, the Hessian matrix takes the following

H_{i} (θ_{u}, θ_{v}^{⊤}) = \frac{\partial s_{i} (θ_{u})}{\partial θ_{v}^{T}} = \frac{\partial}{\partial θ_{v}} \int q (θ_{u}, b_{i}) p (b_{i} | T_{i}, δ_{i}, y_{i}; θ_{u}) d b_{i}

= \int \frac{\partial q (θ_{u}, b_{i})}{\partial θ_{v}^{⊤}} p (b_{i} | T_{i}, δ_{i}, y_{i}; θ) d b_{i} + I_{i},

(2.6)

where

I_{i} = \int q (θ_{u}, b_{i}) {\frac{\partial \log p (b_{i} | T_{i}, δ_{i}, y_{i}; θ_{u})}{\partial θ_{v}^{⊤}}}^{⊤} p (b_{i} | T_{i}, δ_{i}, y_{i}; θ_{u}) d b_{i}

\begin{array}{l} = \int q (θ_{u}; b_{i}) {\frac{\partial \log p (T_{i}, δ_{i} | b_{i}; θ_{u}) + \log p (y_{i} | b_{i}; θ_{u}) + \log p (b_{i}; θ_{u})}{\partial θ_{v}^{⊤}} \\ - {\frac{\partial \log p (T_{i}, δ_{i}, y_{i}; θ_{u})}{\partial θ_{v}^{⊤}}}}^{⊤} p (b_{i} | T_{i}, δ_{i}, y_{i}; θ_{u}) d b_{i} \end{array}

= \int q (θ_{u}, b_{i}) {q (θ_{v}, b_{i}) - s_{i} (θ_{v})}^{⊤} p (b_{i} | T_{i}, δ_{i}, y_{i}; θ_{u}) d b_{i} .

From a computational point of view, to calculate the score vector and the Hessian matrix, we apply Gauss–Hermite quadrature since the integrals with respect to b_i in Equations (2.5)–(2.6) do not have a closed form.

3 Sensitivity analysis

In the previous section, we have claimed that a joint model corresponds to a MNAR mechanism. To illustrate this, let us suppose that T_i denote the event (e.g., dropout) time for the i-th individual and let us indicate by $y_{i}^{o}$ and $y_{i}^{m}$ the set of observed (before T_i) and missing (at and after T_i) longitudinal measurements for the i-th subject. The missing data mechanism, i.e., the conditional distribution of the dropout process given the complete longitudinal data vector $(y_{i}^{o}, y_{i}^{m}),$ is given by:

p (T_{i} | y_{i}^{o}, y_{i}^{m}) = \frac{\int p (T_{i} | b_{i}) p (y_{i}^{o}, y_{i}^{m} | b_{i}) p (b_{i}) d b_{i}}{\int p (y_{i}^{o}, y_{i}^{m} | b_{i}) p (b_{i}) d b_{i}}

= \int p (T_{i} | b_{i}) p (b_{i} | y_{i}^{o}, y_{i}^{m}) d b_{i},

(3.1)

which depends on $y_{i}^{m}$ through the posterior distribution of the random effects. As noted before, statistical literature has so far focused on the development of sensitivity analysis tools for Selection and Pattern Mixture Models, see e.g., Troxel et al. (2004), Molenberghs (2005) and Ma et al. (2005), while there seems to be a lack of research in the context of joint models. For general shared parameter models, Creemers et al. (2010) performed an analysis around the missing at random model by considering a sensitivity parameter that is forced to assume different (but fixed) values around zero. The aim of the present article is, rather, to investigate local sensitivity of JMs by assessing the behaviour of parameter estimates in the longitudinal submodel in a neighbourhood of the MAR hypothesis (α = 0). As stressed before, in a joint model the longitudinal and the dropout processes are linked by, at least, a common random effect structure. For this reason, the majority of the literature regarding sensitivity analysis in this context has focused on the choice of the random effect distribution. The general conclusion that has emerged is that shared parameter models are rather robust to misspecification of this distribution, in particular as the number of repeated measurements increases, see for instance Song et al. (2002), Rizopoulos et al. (2008a) and Rizopoulos et al. (2008b). While the random effect distribution represents an interesting issue, we will mainly focus on the influence of the adopted modelling structure on inferences about the fixed effects in the longitudinal submodel.

3.1 ISNI in JMs

In this section, we review the general formulation of ISNI, proposed by Troxel et al. (2004) and Ma et al. (2005), and further developed by Xie (2008) for longitudinal non-Gaussian responses and give a specific formulation in the case of JMs. A more detailed derivation can be found in Appendix A. This tool aims at investigating the local sensitivity of parameter estimates to departures from the MAR assumption, i.e., how much maximum likelihood parameter estimates for the longitudinal model are influenced by the hypothesis about ignorability of the dropout mechanism. Sensitivity could be also analyzed for the fixed effects in the time-to-event model; however, the focus here is on modelling the longitudinal response accounting for potential informative missing data. For this reason, the longitudinal model plays a primary role and sensitivity will be discussed for fixed effects in that model only. Ma et al. (2005) give an extended overview of the ISNI introduced by Troxel et al. (2004), and carry through this analysis for longitudinal data. In particular, they consider a selection model where the probability of dropout at a given time point is a function of the response measured at the same time and at the previous one. Our case is different and considers a special case of selection models, where the dropout and the longitudinal processes do not share the current response, but only the expected, error-free, part. This particular model has been mentioned in literature as ‘joint model’, see for instance Rizopoulos (2010). In many practicals cases, for instance when the actual response at the dropout time is not available, or it is measured with error, this modelling approach is particularly appealing. Moreover, as it has been pointed out for instance by Rizopoulos and Ghosh (2011), this approach is convenient from a survival analysis standpoint, when one would account for the presence of time dependent covariates measured with error; from a longitudinal standpoint, when the interest is in obtaining unbiased longitudinal parameter estimates; finally, when the focus is on both the processes. Ma et al. (2005) convincingly illustrated the nice features of ISNI for selection models. Motivated by this, the aim of the present article is to investigate how the ISNI can be adapted to work with shared parameter model and in particular with the joint model, which is one of the three main frameworks for dealing with dropout.

Let $\hat{θ} (α) = (\hat{β} (α), \hat{σ} (α), \hat{D} (α), \hat{γ} (α))$ denote the vector of ML estimates for the longitudinal and the survival process parameters in the JM, corresponding to a given value of α. On the other hand, let ${\hat{θ}}_{0} = ({\hat{β}}_{0}, {\hat{σ}}_{0}, {\hat{D}}_{0}, {\hat{γ}}_{0})$ be the corresponding ML estimates under the MAR model, where α = 0. The ISNI measures the rate of change of $\hat{θ} (α)$ from ${\hat{θ}}_{0},$ for a unit displacement of α from 0; it is based on the derivative of $\hat{θ} (α)$ with respect to α, evaluated at $θ = {\hat{θ}}_{0}$ and α = 0. To derive the index formulation, the likelihood function is expanded around $({\hat{θ}}_{0}, α = 0)$ as follows:

L (θ, α) ≅ L (θ_{0}, α = 0) + [{(θ - {\hat{θ}}_{0})}^{⊤}, α] s (θ, α) |_{θ = θ_{0}, α = 0}

+ \frac{1}{2} [{(θ - {\hat{θ}}_{0})}^{⊤}, α] \frac{\partial s (θ, α)}{\partial θ^{⊤}} |_{θ = {\hat{θ}}_{0}, α = 0} {[{(θ - {\hat{θ}}_{0})}^{⊤}, α]}^{⊤},

(3.2)

where $s (θ, α) | θ = \hat{θ} 0, α = 0$ and $\frac{\partial s (θ, α)}{\partial θ^{⊤}} |_{θ = {\hat{θ}}_{0}, α = 0}$ represent the score vector and the Hessian matrix for the parameter vector evaluated at $θ = {\hat{θ}}_{0}$ and a = 0, i.e., at the ML estimates for the MAR model.

For the sake of simplicity, we will indicate these quantities as $s (θ) |_{α = 0}$ and $\frac{\partial s (θ)}{\partial θ^{⊤}} |_{α = 0} .$ In JMs, the score vector and the Hessian matrix are described in Equations (2.5) and (2.6). For the Hessian matrix, in particular, the (u, v)-th component is given by

H (θ_{u}, θ_{v}^{⊤}) |_{α = 0} = \frac{\partial s (θ_{u})}{\partial θ_{v}^{⊤}} |_{α = 0} .

(3.3)

The expansion in (3.2) allows us to write the ML estimate of θ as a function of α:

\hat{θ} (α) ≅ {\hat{θ}}_{0} + \frac{\partial \hat{θ} (α)}{\partial α} |_{α = 0} α .

(3.4)

The index of sensitivity to nonignorability is defined as the derivative of $\hat{θ} (α)$ with respect to α, evaluated at α = 0, see Troxel et al. (2004) and Ma et al. (2005).

Since our interest lies in the study of the sensitivity for the longitudinal parameter vector β, we report the corresponding specific index as follows:

ISNI = \frac{\partial \hat{β} (α)}{\partial α} |_{α = 0} = - {(H (β β^{⊤}) |_{α = 0})}^{- 1} H (β α) |_{α = 0} .

(3.5)

The main advantage of this approach is that it actually does not require to fit the MNAR model, but only to compute the corresponding Hessian matrix at α = 0. Further, it may be noted that (3.4) and (3.5) yield to the following incremental limit:

ISNI = \lim_{α \to 0} \frac{1}{α} (\hat{β} (α) - {\hat{β}}_{0}),

(3.6)

which represents the ratio of the distance between the MNAR and the MAR estimates to the value of the nonignorability parameter as α → 0.

The approximation in (3.2) is a second order Taylor’s expansion of the likelihood similar to the one adopted by Siannis et al. (2005) to study the dependence between survival and the censoring mechanisms. In that case, the approximation is used to build confidence intervals for the parameter of interest and holds only for small values of the dependence parameter. In the present context, Equation (3.2) is used to derive the index of sensitivity to non-ignorability, which is the derivative of the maximum likelihood estimator for a → 0. The interest is on changes of the parameter estimate when we move away from a MAR (or MCAR) assumption, i.e., when α changes from 0. For this reason, the assumption of a small dependence parameter is not necessary. The ISNI definition is based on the Hessian matrix for the JM as α → 0, and it does not actually depend on the value of α, but only on the available data and the sample size. The estimates for α obtained by analyzed data give an idea of the strength of potential dependence between the response and the survival time distribution. For example, if the ISNI suggests that a large displacement in the longitudinal model parameter estimates may occur when α moves away from zero and the α parameter is large, this may imply that the adopted modelling structure could produce parameter estimates that are far from those obtained under MAR assumption and that care is needed to interpret JM results. While an estimate α ≠ 0 is not enough to state dependence, since the estimates depends on the adopted modelling formulation, a large value for the ISNI points out that model parameter estimates may be heavily influenced (in a sensitivity perspective) by the adopted assumptions about the missing data mechanism.

This interpretation is confirmed by Ma et al. (2005), where the ISNI is seen as a measure of the extent to which small departures from ignorability affect the MLE parameters. Along the article, we report the α estimates as an additional tool to quantify the impact that the assumptions on the dropout process suggested in the JM formulation have on results for the analyzed data.

4 A relative index formulation

As it can be noticed, Equation (3.5) represents an absolute measure of parameter change and, therefore, cannot be straight to interpret. To allow an easier interpretation, Troxel et al., (2004) consider the ratio of the ISNI to the standard error of the corresponding MAR estimate rescaled for the variability of the response(v_y), stating sensitivity when this ratio is greater than 1. This is referred to as ‘sensitivity transformation’ (1/c, where $c = s e (\hat{θ} 0) σ_{y} /ISNI$ ) and it represents a scale-free measure to avoid dependence of the ISNI on the scale of Y. The standard error of the MAR estimate is a sensible solution in selection and pattern mixture models where the parameter meaning does not change when we move from the MAR hypothesis. However, in JMs, parameter non-distinctiveness is a major issue and should we use the standard error of the MAR estimate, it would not be considered at all. In this case, the ISNI may vary also for sampling data variability and model approximation, which could be present also when α = 0. Therefore, we need to assess whether the observed ISNI value is indicating a substantial departure of the MNAR estimates from the MAR ones, or rather it is simply due to the sampling variability in the index or to numerical/model approximations. Therefore, the relative formulation suggested by Troxel et al. (2004), which has been proven to work well in SeMs and PMMs, should be extended to deal with the new modelling framework.

For this purpose, we propose to provide a variability estimate based on parametric bootstrap. Let us consider a dataset with longitudinal outcome $y_{i} ~ N (m_{i}, σ),$ m_i = {m_i(t_ij), J = 1, . . ., r_i}, and survival time T_i following a PH model with regression parameter γ and Weibull baseline hazard h₀(T_i); on this data, we consider a MAR approach and fit two separate longitudinal and survival models (i.e., under model (2.2) with α = 0). We draw β samples from the estimated longitudinal and survival distributions with parameters fixed at the ML estimates $\hat{m} i, \hat{σ}, \hat{γ}$ and $\hat{ξ}$ . For each simulated sample, we compute the ISNI and estimate the corresponding standard error as:

{\overset{\land}{se}}_{B} (ISNI) = \sqrt{\frac{1}{B - 1} \sum_{q = 1}^{B} {(I \overset{\land}{SN} I_{q} - \bar{ISNI})}^{2}},

(4.1)

where $I \overset{\land}{SN} I_{q}$ is the estimated ISNI in the q-th sample and $\bar{ISNI}$ is the corresponding sample mean. Therefore, we propose to use the following relative definition

i s n i = \frac{ISNI}{{\overset{\land}{se}}_{B} (ISNI)},

(4.2)

stating sensitivity when | isni |≥1.

The assessment of the sampling variability when the MAR assumption is true, could be of interest for different reasons: first, to study the precision of the ISNI estimate; second, it could lead to a valid alternative for a relative index of sensitivity. This should help us separate the ‘geometric’ feature of the ISNI (representing potential sensitivity of model parameter) from numerical changes due to the association structure as described by the JM.

As it can be easily observed by looking at Equation (2.2), the joint model structure is based on the assumption that some effects are shared by the longitudinal and the survival processes; therefore, we could observe sensitivity to nonignorable missingness even for an approximately null value for the ‘true’ α and for a modest amount of dropouts. Adopting the JM formulation, the β are computed through the joint likelihood, i.e., the score function for β depends, by definition, on the longitudinal and the survival information; the corresponding meaning obviously changes when compared to the estimate obtained from the longitudinal model under the MAR hypothesis. Therefore, we could not be able to distinguish between the sensitivity due to the dependence structure from the one due to effective changes in the meaning of the parameter vector when α moves away from zero and to numerical approximation. Even if expression (3.6) does not depend on the value of α, since the ISNI is defined as the incremental limit of the ML estimate, one would expect that, when the true model is MAR (or MCAR), i.e., when α = 0, parameter estimates are still sensitive to non-ignorability. In fact, some variability of the index could be also present when α is close to zero, for example due to parameter non-distinctness or numerical/model approximations. Therefore, we need to assess whether the observed ISNI value could be due to a strong departure of the MNAR estimates from the MAR ones, or rather to the sampling variability in the index.

Previous considerations suggest that an estimate of the ISNI variability calculated under the MAR assumption could be an appealing measure of the index variability, and the relative formulation $^B(ISNI)$ may account for issues due to parameter non-separability and model specification.

5 Primary Biliary Cirrhosis data

This example comes from the primary biliary cirrhosis (PBC) data collected by the Mayo Clinic from 1974 to 1984, see Murtaugh et al. (2002). PBC is a fatal liver desease characterized by inflammatory destruction of the small bile ducts within the liver, which may lead to cirrhosis of the liver. In this study, 312 patients were considered; 158 were randomly assigned to receive D-penicillamine and 154 placebo. By the end of the study, 140 patients (45%) died, 143 (46%) were alive and 9% were transplanted. We are interested in testing for treatment effect on survival after adjusting for the longitudinal bilirubin levels. The Kaplan-Meier estimate for the time to death is depicted in Figure 1.

By looking at the K-M estimate, the treatment does not seem to significantly affect the survival time; this may be due to the fact that the treatment effect changes with time. Furthermore, patients did not return to the study centers at prespecified time points to provide serum bilirubin measurements and thus we observe great variability between their visiting patterns. In particular, patients have on average 6.23 visits (standard deviation 3.77 visits), resulting in a total of 1945 observations.

Figure 1

Source: Authors’ own.

Taking advantage of the randomization set-up of the study, we fit a joint model where we only correct for treatment. Given that the distribution of the observed bilirubin serum shows substantial skeweness, we consider its natural logarithm. Here, we model the longitudinal dependence of the the bilirubin serum and the survival of enrolled patients by employing the following model for the bivariate (longitudinal and survival) process:

{\begin{cases} \log (y_{i} (t_{i j})) = (β_{0} + b_{i 0}) + (β_{1} + b_{i 1}) x_{i 1} + β_{2} x_{i 2} + β_{3} [x_{i 1} \times x_{i 2}] \\ h (T_{i}) = h_{0} (T_{i}) \exp {γ_{0} + γ_{1 i} w_{i 1} + α m_{i} (T_{i})} \end{cases},

(5.1)

where x_i₁ represents the visit time, x_i₂ the treatment, and x_i₁ × x_i₂ the interaction between treatment and time, while w_i₁ is the baseline treatment. It can be noticed that the parameter vector could be interpreted in a causal framework; in fact, the study is randomized and potential confounding factors are controlled by individual-specific, time-constant, random effects in b_i, shared by both the longitudinal and the survival processes. The aim is at exploring sensitivity of parameter estimates in the longitudinal model to the assumption about non-ignorability of the survival process. The results of sensitivity analysis are shown in Table 1, where β₀ and β represent the parameter estimates for the MAR model and the JM, respectively.

While a quite significant change in parameter estimates is experienced by the time and the treatment marginal effects and the application of the JM suggests a non-ignorable dropout mechanism ( $\hat{α} = 1.212$ with $p < 0.0001$ ), the ISNI and the $ISNI/ σ_{y} se (\hat{β} 0)$ values do highlight a negligible sensitivity for small departure from the MAR assumption. The proposed relative formulation, instead, highlights some, potentially interesting, changes in the intercept and time parameter estimates and a slight change in the treatment effect.

It is interesting to note that with an estimate $\hat{α} = 1.212$ , the JM analysis suggests a non-ignorable dropout mechanism, but this does not mean that the α estimate and the ISNI are related. In fact, a high value of the ISNI highlights that the model formulation is sensitive to even small changes in the α parameter, i.e., for small departures from MAR.

As it is shown in Table 2, it is evident that bootstrap estimate of the ISNI standard error is not of the same order of the rescaled standard error of the MAR estimate. Therefore, this variability measure could account for the ‘geometric’ feature of the ISNI (representing potential sensitivity of model parameter estimates due to the shape of the likelihood surface in a neighbourhood of α = 0) and numerical changes due to sampling variability in the analyzed data and to approximation used in the JM estimation algorithm.

Table 1

Absolute and relative ISNI for the PBC data set. $\hat{α}$ = 1.212 with p < 0.0001

	${\hat{β}}_{0}$	$\hat{β}$	ISNI	$ISNI/ σ_{y} se ({\hat{β}}_{0})$	$^$
Intercept	0.563	0.662	0.391	0.031	1.050
Time	0.180	0.192	–0.045	–0.001	–1.011
Treatment	–0.133	–0.021	–0.334	–0.038	–0.467
Treat*Time	–0.004	0.005	–0.025	–0.001	–0.006

Source: Authors’ own.

Table 2

Standard errors of the MAR parameter estimates (rescaled for the response standard deviation σ_y) and the ISNI for the PBC data set. $\hat{α}$ = 1.212 with p < 0.0001

	$se ({\hat{β}}_{0})$	$^(ISNI)$
Intercept	0.086	0.372
Time	0.018	0.044
Treatment	0.119	0.715
Treat*Time	0.025	4.490

Source: Authors’ own.

6 Simulation study

To investigate the empirical behaviour of the ISNI and the proposed relative formulation in joint models, we performed the following simulation study.

6.1 Study design

We simulate a longitudinal response y_i(t_ij) from a Gaussian distribution with mean $m_{i} (t_{i j}) = (β_{0} + b_{i 0}) + (β_{1} + b_{i 1}) x_{1 i} + β_{2} x_{2 i} + β_{3} [x_{1 i} \times x_{2 i}]$ according to model (5.1) adopted for the PBC dataset. For the dropout process, we consider the hazard function

h_{i} (T_{i}) = h_{0} (T_{i}) \exp {γ w_{i} + α m_{i} (T_{i})},

(6.1)

In this study, x_1i is a sequence from 0 to K (maximum observation time) and represents the longitudinal measurement occasions, x_2i is a Bernoulli variable with parameter 0.5, representing a treatment variable and [x_1i × x_2i] is the interaction between time and treatment. On the other hand, w_i = x_2i. Individual censoring times have been randomly drawn from an exponential distribution with mean chosen to result in about 50% censoring.

We consider N = 1000 samples of size n = 312 observed for K = 5 or K = 15 occasions. We simulate different scenarios by varying the longitudinal parameter values, the α values and the random effect covariance matrix. For the latter, we consider the following matrices:

D_{1} = (\begin{array}{l} 1.066 0.073 \\ 0.073 0.029 \end{array}),

and

D_{2} = (\begin{array}{l} 0.5 0.01 \\ 0.01 0.5 \end{array}) .

While D₁ is the estimated random effect covariance matrix for the PBC dataset, D₂ describes a homoscedastic random effect covariance with a positive association between b_i₁ and b_i₂. The response variance and the survival covariate effect w_i, as well as the baseline hazard h₀(T_i), are fixed to the estimates for the PBC dataset, i.e., σ² = 0.158 and γ = 0.07.

On each simulated dataset, we compute the ISNI in (2.3), the sensitivity transformation $1 / c = ISNI/ σ_{y} se (β_{0})$ proposed by Troxel et al. (2004) and the proposed relative isni in (4.2) for the fixed effects in the longitudinal process.

Further, we perform a sensitivity analysis to model misspecification. We consider three potential cases of misspecification and fix α = 0:

The observed data are simulated according to the joint model

{\begin{array}{l} y_{i} (t_{i j}) = (β_{0} + b_{i 0}) + (β_{1} + b_{i 1}) x_{1 i} + β_{2} x_{2 i} + β_{3} [x_{1 i} \times x_{2 i}] \\ h_{i} (T_{i}) = h_{0} (T_{i}) \exp {γ w_{i} + α y_{i} (t_{i j}) \end{array}

and the absolute and relative ISNIs are computed for parameter estimates for model (6.1);

The longitudinal response is drawn from a skew-normal distribution, Azzalini (1985), with shape parameter equal to 5, while the ISNIs are computed for the JM with Gaussian response;

The survival time is assumed to follow an exponential distribution with rate 0.5, while the sensitivity indexes are computed for the JM with survival time following a proportional hazard model.

The aim is to assess robustness of the ISNI to model specification; this is quite a major point, since the ISNI is known to measure departures from the MAR estimates only when the fitted model is the true one.

The simulation results are organized as follows. In Section 6.2 we explore the sensitivity under the JM parametrization. Section 6.3 gives some interpretational remarks, while Section 6.4 describes the ISNI behaviour when the model is misspecified.

6.2 Local sensitivity for JMs

Tables 3–5 show the median of the ISNI, the Troxel sensitivity transformation 1/c and the proposed relative $i s n i = ISNI/ \overset{\land}{se} (I \overset{\land}{SN} I)$ for different datasets. It can be noticed that both the ISNI and the isni suggest some sensitivity of the JM as α departs from zero. On the other hand, this behaviour is less evident by looking at the sensitivity transformation 1/c.

Table 3
Simulation results. Median over 1000 datasets of ISNI, Troxel sensitivity transformation 1/c and the proposed relative $i s n i = ISNI/ \overset{\land}{se} (l \overset{\land}{SN} I)$ for different values of non-ignorability parameter α, random effects covariance matrices D₁ and D₂, number of repeated measurements K and true longitudinal parameter β_True = (–0.7, 0.2, –0.1, –0.05).

D₁
D₂

β_True ${\hat{β}}_{Mar}$ ISNI 1/c isni ${\hat{β}}_{Mar}$ ISNI 1/c isni

K = 5

α = 0 Int –0.70 –0.70 –0.14 –0.03 –0.39 –0.66 –0.05 –0.01 –0.43

Time 0.20 0.20 0.08 0.00 0.24 0.21 0.04 0.00 0.39

Drug –0.10 –0.10 0.00 0.00 0.01 –0.10 –0.00 0.00 –0.00

TimeDrug –0.05 –0.00 –0.00 0.00 –0.01 –0.00 –0.00 0.00 –0.01

K = 15

α = 0 Int –0.70 –0.71 0.30 0.06 0.34 –0.69 –0.23 –0.04 –0.00

Time 0.20 0.20 0.68 0.03 0.83 0.18 0.12 0.02 0.00

Drug –0.10 –0.09 –0.02 –0.01 –0.02 –0.10 0.01 0.00 0.00

TimeDrug –0.05 –0.00 0.03 0.00 0.02 0.01 –0.00 –0.00 0.00

K = 5

α = 1.2 Int –0.70 –0.67 –0.05 –0.01 –0.02 –0.74 1.73 0.29 0.90

Time 0.20 0.52 0.83 0.11 0.50 0.00 3.01 0.42 1.54

Drug –0.10 –0.08 –0.00 0.00 –0.00 –0.12 –0.10 –0.02 –0.06

TimeDrug –0.05 –0.00 –0.01 –0.00 –0.01 –0.01 0.11 0.02 0.04

K = 15

α = 1.2 Int –0.70 –0.71 –0.24 –0.05 –0.84 –0.68 –0.06 –0.01 –0.10

Time 0.20 0.20 0.10 0.00 0.35 0.20 0.02 0.00 0.03

Drug –0.10 –0.10 0.01 0.00 0.03 –0.08 0.01 0.00 0.01

TimeDrug –0.05 –0.00 –0.00 0.00 –0.01 –0.00 –0.00 0.00 –0.00

K = 5

α = –1.2 Int –0.70 –0.66 –0.05 –0.01 –0.08 –0.73 0.60 0.10 0.76

Time 0.20 0.52 0.32 0.04 0.57 0.00 1.26 0.17 1.67

Drug –0.10 –0.12 –0.02 –0.00 –0.03 –0.10 –0.01 –0.00 –0.03

TimeDrug –0.05 –0.02 0.00 0.00 0.00 –0.01 0.02 0.00 0.03

K = 15

α = –1.2 Int –0.70 –0.70 –0.19 –0.04 –0.25 –0.70 –0.30 –0.04 –0.20

Time 0.20 0.20 0.73 0.03 1.07 0.19 0.13 0.02 0.08

Drug –0.10 –0.11 –0.01 –0.00 –0.00 –0.11 0.01 0.00 0.01

TimeDrug –0.05 –0.00 0.02 0.00 0.02 –0.02 –0.00 –0.00 –0.00

		D₁	D₂
α = 0	Int	–0.70	–0.70	–0.14	–0.03	–0.39	–0.66	–0.05	–0.01	–0.43
	Time	0.20	0.20	0.08	0.00	0.24	0.21	0.04	0.00	0.39
	Drug	–0.10	–0.10	0.00	0.00	0.01	–0.10	–0.00	0.00	–0.00
	Time*Drug	–0.05	–0.00	–0.00	0.00	–0.01	–0.00	–0.00	0.00	–0.01
K = 15
α = 0	Int	–0.70	–0.71	0.30	0.06	0.34	–0.69	–0.23	–0.04	–0.00
	Time	0.20	0.20	0.68	0.03	0.83	0.18	0.12	0.02	0.00
	Drug	–0.10	–0.09	–0.02	–0.01	–0.02	–0.10	0.01	0.00	0.00
	Time*Drug	–0.05	–0.00	0.03	0.00	0.02	0.01	–0.00	–0.00	0.00
	K = 5
α = 1.2	Int	–0.70	–0.67	–0.05	–0.01	–0.02	–0.74	1.73	0.29	0.90
	Time	0.20	0.52	0.83	0.11	0.50	0.00	3.01	0.42	1.54
	Drug	–0.10	–0.08	–0.00	0.00	–0.00	–0.12	–0.10	–0.02	–0.06
	Time*Drug	–0.05	–0.00	–0.01	–0.00	–0.01	–0.01	0.11	0.02	0.04
	K = 15
α = 1.2	Int	–0.70	–0.71	–0.24	–0.05	–0.84	–0.68	–0.06	–0.01	–0.10
	Time	0.20	0.20	0.10	0.00	0.35	0.20	0.02	0.00	0.03
	Drug	–0.10	–0.10	0.01	0.00	0.03	–0.08	0.01	0.00	0.01
	Time*Drug	–0.05	–0.00	–0.00	0.00	–0.01	–0.00	–0.00	0.00	–0.00
	K = 5
α = –1.2	Int	–0.70	–0.66	–0.05	–0.01	–0.08	–0.73	0.60	0.10	0.76
	Time	0.20	0.52	0.32	0.04	0.57	0.00	1.26	0.17	1.67
	Drug	–0.10	–0.12	–0.02	–0.00	–0.03	–0.10	–0.01	–0.00	–0.03
	Time*Drug	–0.05	–0.02	0.00	0.00	0.00	–0.01	0.02	0.00	0.03
	K = 15
α = –1.2	Int	–0.70	–0.70	–0.19	–0.04	–0.25	–0.70	–0.30	–0.04	–0.20
	Time	0.20	0.20	0.73	0.03	1.07	0.19	0.13	0.02	0.08
	Drug	–0.10	–0.11	–0.01	–0.00	–0.00	–0.11	0.01	0.00	0.01
	Time*Drug	–0.05	–0.00	0.02	0.00	0.02	–0.02	–0.00	–0.00	–0.00

Source: Authors’ own.

Table 4

Simulation results. Median over 1000 datasets of ISNI, Troxel sensitivity transformation 1/c and the proposed relative $i s n i = ISNI/ \overset{\land}{se} (l \overset{\land}{SN} I)$ for different values of non-ignorability parameter α, random effects covariance matrices D₁ and D₂, number of repeated measurements K and true longitudinal parameter β_True = (0.7, –0.2, –0.1, –0.05).

		D₁					D₂
		β_True	${\hat{β}}_{Mar}$	ISNI	1/c	isni	${\hat{β}}_{Mar}$	ISNI	1/c	isni
K = 5
α = 0	Int	0.700	0.716	–0.131	–0.029	–1.788	0.733	–0.024	–0.005	–0.026
	Time	–0.200	–0.197	0.074	0.003	1.018	–0.200	0.224	0.010	0.530
	Drug	–0.100	–0.114	0.004	0.001	0.032	–0.115	–0.003	–0.001	–0.001
	Time*Drug	–0.050	–0.002	–0.001	0.000	–0.007	–0.002	0.018	0.001	0.026
K = 15
α = 0	Int	0.700	0.689	–0.038	–0.008	–0.024	0.699	–0.241	–0.040	–0.100
	Time	–0.200	–0.199	0.278	0.013	0.554	–0.192	0.129	0.021	0.054
	Drug	–0.100	–0.101	0.006	0.002	0.014	–0.102	0.005	0.001	0.002
	Time*Drug	–0.050	0.001	–0.013	–0.001	–0.022	–0.017	–0.003	–0.001	–0.001
K = 5
α = 1.2	Int	0.700	0.723	0.139	0.025	0.114	0.666	0.432	0.074	0.307
	Time	–0.200	0.087	1.807	0.237	1.068	–0.433	1.657	0.221	1.347
	Drug	–0.100	–0.105	–0.017	–0.005	–0.013	–0.077	–0.024	–0.006	–0.027
	Time*Drug	–0.050	–0.016	0.095	0.017	0.062	–0.011	–0.079	–0.014	–0.054
K = 15
α = 1.2	Int	0.700	0.717	–0.199	–0.043	–0.495	0.720	–0.147	–0.033	–0.274
	Time	–0.200	–0.198	0.080	0.003	0.200	–0.201	0.177	0.007	0.316
	Drug	–0.100	–0.094	0.001	0.001	0.001	–0.098	–0.022	–0.007	–0.041
	Time*Drug	–0.050	–0.007	–0.002	0.000	–0.002	–0.002	0.020	0.001	0.028
K = 5
α = –1.2	Int	0.700	0.739	0.064	0.011	0.194	0.678	0.030	0.005	0.032
	Time	–0.200	0.071	0.811	0.109	1.527	–0.431	0.673	0.087	1.152
	Drug	–0.100	–0.104	–0.023	–0.006	–0.037	–0.090	0.007	0.002	0.014
	Time*Drug	–0.050	0.009	0.062	0.011	0.108	–0.003	–0.055	–0.010	–0.081
K = 15
α = –1.2	Int	0.700	0.684	–0.208	–0.045	–0.673	0.694	–0.386	–0.059	–0.007
	Time	–0.200	–0.198	0.166	0.007	0.628	–0.211	0.156	0.023	0.003
	Drug	–0.100	–0.066	0.004	0.001	0.035	–0.099	0.033	0.007	0.001
	Time*Drug	–0.050	0.000	–0.005	0.000	–0.019	–0.001	–0.011	–0.002	0.000

Source: Authors’ own.

Table 5

		D₁					D₂
		β_True	${\hat{β}}_{Mar}$	ISNI	1/c	isni	${\hat{β}}_{Mar}$	ISNI	1/c	isni
K = 5
α = 0	Int	0.700	–0.704	–0.143	–0.029	–0.393	–0.663	–0.046	–0.009	–0.428
	Time	0.200	0.201	0.083	0.004	0.240	0.211	0.036	0.002	0.391
	Drug	0.100	–0.098	0.004	0.001	0.005	–0.104	–0.001	0.000	–0.004
	Time*Drug	–0.050	–0.004	–0.004	0.000	–0.012	–0.003	–0.001	0.000	–0.007
K = 15
α = 0	Int	0.700	–0.713	0.295	0.061	0.342	–0.689	–0.229	–0.037	–0.001
	Time	0.200	0.204	0.679	0.032	0.831	0.185	0.120	0.018	0.001
	Drug	0.100	–0.090	–0.022	–0.007	–0.024	–0.097	0.009	0.002	0.000
	Time*Drug	–0.050	–0.004	0.027	0.002	0.017	0.011	–0.004	–0.001	0.000
K = 5
α = 1.2	Int	0.700	–0.666	–0.048	–0.009	–0.019	–0.736	1.726	0.295	0.901
	Time	0.200	0.524	0.827	0.109	0.504	0.003	3.015	0.419	1.536
	Drug	0.100	–0.083	–0.001	0.000	–0.003	–0.116	–0.100	–0.025	–0.059
	Time*Drug	–0.050	–0.001	–0.011	–0.002	–0.005	–0.013	0.110	0.021	0.043
K = 15
α = 1.2	Int	0.700	–0.712	–0.243	–0.053	–0.839	–0.681	–0.059	–0.013	–0.099
	Time	0.200	0.199	0.100	0.004	0.353	0.204	0.018	0.001	0.032
	Drug	0.100	–0.096	0.005	0.002	0.029	–0.076	0.005	0.001	0.009
	Time*Drug	–0.050	–0.004	–0.003	0.000	–0.008	–0.001	–0.002	0.000	–0.004
K = 5
α = –1.2	Int	0.700	–0.660	–0.051	–0.009	–0.080	–0.728	0.602	0.097	0.764
	Time	0.200	0.519	0.323	0.041	0.569	0.001	1.258	0.170	1.667
	Drug	0.100	–0.116	–0.017	–0.004	–0.025	–0.103	–0.014	–0.003	–0.028
	Time*Drug	–0.050	–0.019	0.001	0.000	0.001	–0.009	0.017	0.003	0.029
K = 15
α = –1.2	Int	0.700	–0.702	–0.194	–0.042	–0.251	–0.697	–0.300	–0.043	–0.202
	Time	0.200	0.203	0.727	0.033	1.066	0.193	0.126	0.018	0.082
	Drug	0.100	–0.110	–0.006	–0.002	–0.005	–0.105	0.011	0.003	0.009
	Time*Drug	–0.050	–0.003	0.017	0.001	0.019	–0.016	–0.003	–0.001	–0.002

Source: Authors’ own.

Tables 3–5 highlight sensitivity for the intercept and the time effects. The observed sensitivity of the drug effect is negligible with ISNI values approaching to zero. This may be due to the fact that the drug effect is time-invariant. Moreover, it is worth noticing that the sensitivity of JMs decreases as the number of repeated measurements increases and that the random effect covariance matrix does not play a substantial role in parameter sensitivity. These empirical evidences could suggest that, if the JM is the ‘true’ model and as the number of repeated measures increases, the longitudinal fixed effect estimates are not substantially influenced by the presence of MNAR mechanisms.

6.3 Interpretational remarks

As pointed out in Section 6.2, the ‘sensitivity transformation’ of Troxel et al. (2004) suggests a low sensitivity even when the non-ignorability parameter is far from zero. As it can be evinced by Equation (3.6), the ISNI variability takes into account both the variability of the MAR parameter estimates and the variability due to the dependence structure. Rather, 1/c considers the MAR parameter variation and may underestimate the effective sensitivity. The absolute ISNI, on the other hand, suggests sensitivity when α ≠ 0, but it might be misleading since it is an absolute measure and the corresponding values cannot be directly interpreted. In the simulation study, the proposed relative index (isni) leads to a clearer interpretation, since it provides a direct comparison of the ISNI and a measure of its variability in a neighbourhood of the MAR hypothesis. Thus, we compare the ISNI with its natural variation, calculated under α = 0; sensitivity is therefore evident for a strong curvature of the log-likelihood function or for a strong departure from MAR. Rather than a theoretical measure, isni is an empirical relative index, which is strongly linked to the analyzed data.

6.4 Local sensitivity under misspecification

The ISNI is known to depend on the parametric structure of the model. In this Section we show the sensitivity behaviour of JMs when the model is misspecified.

Table 6 suggests that the joint model structure and the longitudinal distribution do not seem to influence the JM sensitivity to non-ignorability. The survival distribution, on the other hand, may influence the indexes value. Therefore, it is suggested to check the assumptions on the survival model before handling a sensitivity analysis with the ISNI.

Table 6
Simulation study for model misspecification. Median over 1000 datasets of ISNI, Troxel sensitivity transformation 1/c and the proposed relative $i s n i = ISNI/ \overset{\land}{se} (l \overset{\land}{SN} I)$ for source of misspecification relative to model structure (a), the longitudinal distribution (b) and the survival distribution (c), α = 0, β_True = (0.7, 0.2, 0.1, –0.05) and random effects covariance matrix D₁.

β_True ISNI 1/c isni

Miss. (a) Int –0.663 –0.046 –0.009 –0.428

Time 0.211 0.036 0.002 0.391

Drug –0.104 –0.001 0.000 –0.004

TimeDrug
–0.003
–0.001
0.000
–0.007

Miss. (b) Int –0.713 0.295 0.061 0.342

Time 0.204 0.679 0.032 0.831

Drug –0.090 –0.022 –0.007 –0.024

TimeDrug
–0.004
0.027
0.002
0.017

Miss. (c) Int –0.704 4.936 1.106 5.948

Time 0.195 2.389 0.345 2.853

Drug –0.098 0.236 0.070 0.245

Time*Drug –0.007 0.170 0.033 0.232

		β_True	ISNI	1/c	isni
Miss. (a)	Int	–0.663	–0.046	–0.009	–0.428
	Time	0.211	0.036	0.002	0.391
	Drug	–0.104	–0.001	0.000	–0.004
	Time*Drug	–0.003	–0.001	0.000	–0.007
Miss. (b)	Int	–0.713	0.295	0.061	0.342
	Time	0.204	0.679	0.032	0.831
	Drug	–0.090	–0.022	–0.007	–0.024
	Time*Drug	–0.004	0.027	0.002	0.017
Miss. (c)	Int	–0.704	4.936	1.106	5.948
	Time	0.195	2.389	0.345	2.853
	Drug	–0.098	0.236	0.070	0.245
	Time*Drug	–0.007	0.170	0.033	0.232

Source: Authors’ own.

7 Discussion

In this article, we have extended the notion of ISNI, see Troxel et al. (2004) and Ma et al. (2005), to the class of joint models, Wulfsohn and Tsiatis (1997). We have focused on changes registered in fixed parameters for the longitudinal process due to potential misspecification in the missing data mechanism and in the dependence structure.

We have proposed a relative index of local sensitivity, given by the ratio of the ISNI values to the corresponding standard error estimated via parametric bootstrap. In this way, we have provided an estimate of the sampling variability of the index when the MAR hypothesis is true, which could help define a further relative measure of parameter sensitivity to departures from the MAR assumption. This approach is based on the empirical evidence that, also when the MAR hypothesis is true, the ISNI may take non null values; therefore, this absolute measure of displacement could be compared to a measure of its sampling variability, which is, in the context of joint models, likely linked to parameter non-separability. The proposed relative ISNI (isni) seems to lead to a clearer interpretation of parameters sensitivity, while the classical ISNI and sensitivity transformation may be misleading, at least in some cases, when a joint model is adopted.

The data application and the simulation study have shown a slight sensitivity of model parameter estimates when the departure from MAR is not too large, and a decreasing effect of the assumptions about ignorability for the missing data mechanism as the length of the individual sequences increases. The random effect covariance structure does not seem to play a substantial role. In addition, the sensitivity with respect to model misspecification seems not to be meaningful, with the only exception of the misspecification for the survival distribution.

In conclusion, as further illustrated in a benchmark data example, the performed sensitivity analysis based on the proposed relative formulation can be helpful to properly account for the peculiarities of joint models.

Footnotes

Calculating the ISNI for the shared parameter models

The score vectors with respect to the longitudinal fixed effects and the association parameter for the joint model assume the form

s (β) = \sum_{i = 1}^{n} δ_{i} α X_{i}^{⊤} (T_{i}) - \int_{b_{i}} \int_{0}^{T_{i}} h_{0} (s) α X_{i}^{⊤} (s) \exp {γ^{⊤} W_{i} + α m_{i} (s)} d s

+ X_{i} {(t_{i j})}^{⊤} σ^{- 2} (y_{i} - m_{i} (t)) p (b_{i} | T_{i}, δ_{i}, y_{i}; θ) d b_{i},

and

s (α) = \sum_{i = 1}^{n} δ_{i} m_{i} (T_{i}) - \int_{b_{i}} \int_{0}^{T_{i}} h_{0} (s) m_{i} (s) \exp {γ^{⊤} W_{i} + α m_{i} (s)} d s

p (b_{i} | T_{i}, δ_{i}, y_{i}; θ) d b_{i} .

The posterior distribution of the random effects is given by

p (b_{i} | T_{i}, δ_{i}, y_{i}; θ) = \frac{p (T_{i}, δ_{i} | b_{i}; θ) p (y_{i} | b_{i}; θ) p (b_{i})}{p (T_{i}, δ_{i}, y_{i})} .

We derive the ISNI for the longitudinal fixed covariate effects as follows:

ISNI = - {(σ^{- 2} X_{i} {(t_{i j})}^{⊤} X_{i} (t_{i j}) - \sum_{i = 1}^{n} \int_{0}^{T_{i}} h_{0} (s) α^{2} X_{i} {(s)}^{⊤} X_{i} (s) μ_{1, b_{i}} (s) d s + I_{1})}^{- 1}

\times \sum_{i = 1}^{n} δ_{i} X_{i} (t_{i j}) - \int_{0}^{T_{i}} [h_{0} (s) X_{i} (s) μ_{1, b_{i}} (s)] [1 + α m_{i} (s)] d s + I_{2},

where $μ_{1, b_{i}} (t_{i}) = E_{b_{i} | T_{i}, δ_{i}, y_{i}} [\exp {γ^{⊤} W_{i} + α m_{i} (T_{i})}],$

i_{1} = \int_{b_{i}} q (β, b_{i}) {q (β, b_{i}) - s_{i} (β)}^{⊤} p (b_{i} | t_{i}, δ_{i}, y_{i}),

and

q (β, b_{i}) = α δ_{i} X_{i} (t_{i j}) - \int_{0}^{T_{i}} h_{0} (s) α X_{i} (s) \exp {γ^{⊤} W_{i} + α m_{i} (s)} d s

+ X_{i} {(t_{i j})}^{⊤} σ^{- 2} (y_{i} - m_{i} (t)),

while $I_{2} = \int_{b_{i}} q (β, b_{i}) {q (α, b_{i}) - s_{i} (α)}^{⊤} p (b_{i} | T_{i}, δ_{i}, y_{i})$

and

q (α, b_{i}) = δ_{i} m_{i} (T_{i}) - \int_{0}^{T_{i}} h_{0} (s) m_{i} (s) \exp {γ^{⊤} W_{i} + α m_{i} (s)} d s .

References

Azzalini

(1985) A class of distributions which includes the normal ones. Scandinavian Journal of Statistics, 12, 171–78.

Cox

(1972) Regression models and life tables (with discussion). Journal of the Royal Statistical Society, Series B, 34, 187–220.

Creemers

Hens

Aerts

Molenberghs

Verbeke

Kenward

(2010) A sensitivity analysis for shared-parameter models for incomplete longitudinal data. Biometrical Journal, 52, 111–25.

Diggle

Heagerty

Liang

Zeger

(1994) Analysis of Longitudinal Data. Oxford: Oxford University Press.

Faucett

Thomas

(1996) Simultaneously modelling censored survival data and repeatedly measured covariates. Statistics in Medicine, 15, 1663–1685.

Follmann

(1995) An approximate generalized linear model with random effects for informative missing data. Biometrics, 51, 151–68.

Kronrod

(1964) Doklady akad (in russian). Nauk SSSR, 154, 283–86.

Little

RJA

(1993) Pattern-mixture models for multivariate incomplete data. Journal of the American Statistical Association, 88, 125–34.

Little

RJA

(1995) Modelling the dropout mechanism in repeated-measures studies. Journal of American Statical Association, 90, 1112–21.

10.

Little

RJA

Rubin

(2002) Statistical analysis with missing data. New York: John Wiley & Sons.

11.

Troxel

Heitjan

(2005) An index of local sensitivity to nonognorabile drop-out in longitudinal modelling. Statistics in Medicine, 24, 2129–50.

12.

Molenberghs

Verbeke

(2005) Models for discrete longitudinal data. New York: Springer.

13.

Murtaugh

Dickinson

Van Dam

Malincho

Grambsch

Langworthy

Gips

(2002) Primary billiary cirrhosis: Prediction of short term survival based on repeated patient visit. Hepatology, 20, 126–34.

14.

Rizopoulos

(2010) JM: An r package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software, 35, 9.

15.

Rizopoulos

Ghosh

(2011) A bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event. Statistics in Medicine, 30, 1366–80.

16.

Rizopoulos

Lesaffre

(2009) Fully exponential laplace approximations for the joint modelling of survival and longitudinal data. Journal of the Royal Statistical Society, Series B, 71, 637–54.

17.

Rizopoulos

Verbeke

Molenberghs

(2008) Shared parameter models under random effects misspecification. Biometrika, 95, 63–74.

18.

Rizopoulos

Vanrenterghem

(2008) A two-part joint model for the analysis of survival and longitudinal binary data with excess zeros. Biometrics, 64, 611–19.

19.

Siannis

(2005) Sensitivity analysis for informative censoring in parametric survival models. Biostatistics, 6, 77–91.

20.

Song

Davidian

Tsiatis

(2002) A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. Biometrics, 58, 742–53.

21.

Troxel

Heitjan

(2004). Shared parameter models for the joint analysis of longitudinal data and event times. Statistica Sinica, 14, 1221–37.

22.

Tsiatis

Davidian

(2004) An overview of joint modeling of longitudinal and time-to-event data. Statistica Sinica, 14, 793–818.

23.

Wulfsohn

Tsiatis

(1997) A joint model for survival and longitudinal data measured with error. Biometrics, 53, 330–39.

24.

Xie

(2008) A local sensitivity analysis approach to longitudinal non-gaussian data with non-ignorable dropout. Statistics in Medicine, 27, 3155–77.