A Semiparametric Marginalized Model for Longitudinal Data with Informative Dropout

Abstract

We propose a marginalized joint-modeling approach for marginal inference on the association between longitudinal responses and covariates when longitudinal measurements are subject to informative dropouts. The proposed model is motivated by the idea of linking longitudinal responses and dropout times by latent variables while focusing on marginal inferences. We develop a simple inference procedure based on a series of estimating equations, and the resulting estimators are consistent and asymptotically normal with a sandwich-type covariance matrix ready to be estimated by the usual plug-in rule. The performance of our approach is evaluated through simulations and illustrated with a renal disease data application.

1. Introduction

Longitudinal studies often encounter data attrition because subjects drop out before the designated study end. Both statistical analysis and practical interpretation of longitudinal data can be complicated by dropouts. For example, in the Modification of Diet in Renal Disease (MDRD) study [1, 2], one main interest was to investigate the efficacy of interventions of blood pressure control and diet modification on patients with impaired renal functions. The primary outcome was glomerular filtration rate (GFR), which measured filtering capacity of kidneys, and was repeatedly measured over the study period. However, some patients could leave the study prematurely for kidney transplant or dialysis, which precluded further GFR measurements. This resulted in a dropout mechanism that could relate to patients’ kidney function and correlate with their GFR values. Other patients were followed to the end of the study or dropped out due to independent reasons. Thus, statistical analysis of longitudinal GFR needs to take into consideration the presence of mixed types of informative and independent dropouts.

Many statistical models and inference approaches have been proposed to accommodate the nonignorable missingness into modeling longitudinal data (see reviews [3–8]). According to the target of inference and the interpretation of model parameters, existing methods can be classified into three categories: subject-specific inference, event-conditioning inference, and marginal inference. First, a widely used modeling strategy for longitudinal data with informative dropouts is to specify their joint distribution via shared or correlated latent variables. Under such model assumptions, the longitudinal parameters have a conditional, subject-specific interpretation (e.g., [9–11]). But the interpretation of longitudinal parameters usually changes with the number and characteristics of latent variables assumed, for example, a single random intercept versus a random intercept plus a random slope.

Second, event-conditioning approaches have also been widely used when the target of inference is within subgroups of patients with particular dropout patterns or when the dropout can potentially change the material characteristic of the longitudinal process (e.g., death). The inference is usually conducted conditioning on the dropout pattern or on the occurrence of the dropout event. Thus, model parameters have an event-conditioning subpopulation-averaged interpretation, for example, pattern-mixture models for the group expectation of each dropout pattern [3, 12]; treatment effects among survivors [13]; gender and age effects in mortal cohort [14]. Because the interpretation of such models is made by conditioning on a future event, event-conditioning approaches may be natural in a retrospective setting but may not be directly useful for the evaluation of treatment efficacy prospectively.

Lastly, when the research objective is to study covariate effects at population level in a dropout-free situation, marginal models address this concern directly. When data are without missing or missing completely at random (using Rubin's definition on missingness [15]), the estimation of model parameters can be carried out by the generalized estimating equation (GEE) approach assuming a “working” correlation matrix [16]. When dropouts are missing at random, the inverse probability-weighted GEE methods are commonly used [17, 18]. In the presence of informative dropouts, the class of selection models that were originally proposed to adjust selection bias in econometrics [19] have been widely used for the marginal analysis of longitudinal data [20–22]. Recently, the marginalized transition model [23] and marginalized pattern-mixture model [24] were proposed for binary longitudinal data with finite nonignorable nonresponse patterns. These marginalized approaches provide a powerful tool for studying the marginal association between longitudinal outcomes and covariates while incorporating nonignorable nonresponses.

In this paper, we shall adopt the idea of shared latent variables to account for the dependence between longitudinal responses and informative dropouts while focusing on marginal inference for the longitudinal responses. Here dropouts can occur on a continuous time scale. We develop an effective estimation procedure built on a series of asymptotically unbiased estimating equations with light computational burden. The resulting estimators for longitudinal parameters are shown to be consistent and asymptotically normal, with a sandwich-type variance-covariance matrix that can be estimated by the usual plug-in rule.

The remainder of the paper is organized as follows. In Section 2, we introduce the notation and the proposed semiparametric marginalized model. In Section 3, a simple estimating equation-based procedure is first proposed for the situation with pure informative dropouts and is extended to a more general situation where there is a mixture of random dropouts and informative dropouts. Asymptotic properties of resulting estimators are also studied. Simulation studies and an application to a renal disease data set are given in Section 4. Some remarks are discussed in Section 5. All the technical details are provided in the appendix.

2. Notation and Model Specification

2.1. Data Notation

Consider that a longitudinal study follows 𝑛 subjects over time period [0,𝜏]. For the ith subject 𝑖=1,…,𝑛, the complete-data consist of {𝑌_𝑖𝑗,𝑋_𝑖𝑗,𝑡_𝑖𝑗,𝑗=1,…,𝑛_𝑖}, where 𝑌_𝑖𝑗 is the value of response at the jth observation time 𝑡_𝑖𝑗 and 𝑋_𝑖𝑗 is a 𝑝×1 vector of covariates associated with response 𝑌_𝑖𝑗. Note that 𝑋_𝑖𝑗 includes baseline covariates that are separately denoted by 𝑍_𝑖 and potential time-dependent covariates. Let 𝑇_𝑖 denote the informative dropout time and 𝐶_𝑖 denote the random censoring time that is independent of (𝑌_𝑖𝑗,𝑇_𝑖) given the covariates. In practice, we observe (𝑇^∗𝑖,𝛿_𝑖), where 𝑇^∗𝑖=min(𝑇_𝑖,𝐶_𝑖) and 𝛿_𝑖=𝐼(𝑇_𝑖≤𝐶_𝑖) taking the value of 1 if the informative dropout time is observed and 0 otherwise. Throughout the paper, let 𝐼(⋅) denote the indicator function. Due to the dropout, longitudinal responses and covariates can only be observed at 𝑡_𝑖𝑗≤𝑇^∗𝑖. Hence, the observed data are {𝑌_𝑖𝑗,𝑋_𝑖𝑗,𝑡_𝑖𝑗,𝑇^∗𝑖,𝛿_𝑖,𝑖=1,…,𝑛,𝑗=1,…,𝑚_𝑖}, where 𝑚_𝑖=∑^𝑛𝑖_𝑗=1𝐼(𝑡_𝑖𝑗≤𝑇^∗𝑖).

2.2. Semiparametric Marginalized Latent Variable Model

We first introduce the composition of our proposed model and then discuss the model motivation and interpretation. The first component is a marginal generalized linear model for longitudinal responses 𝑌_𝑖𝑗’s:𝑔E𝑌_𝑖𝑗∣𝑋_𝑖𝑗≜𝑔𝜇_𝑖𝑗=𝛽^𝑀𝑋_𝑖𝑗,(2.1) where 𝑔(⋅) is a known link function and 𝜇_𝑖𝑗 denotes the marginal expectation. The second component is a linear transformation model for the informative dropout time 𝑇_𝑖:𝐻𝑇_𝑖=−𝜃^𝑍_𝑖+𝜂_𝑖,(2.2) where 𝐻(⋅) is an unspecified monotone transformation function and 𝜂_𝑖 is assumed to follow a known continuous distribution 𝐹(⋅) that is independent of 𝑍_𝑖. The last component is a conditional mean model characterizing the dependence between longitudinal responses and informative dropouts:𝑔E𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂_𝑖≜𝑔𝜇_𝑖𝑗𝜂_𝑖=Δ_𝑖𝑗+𝛼^𝑏_𝑖𝑗𝜂_𝑖,(2.3) where the latent random effects 𝐛_𝑖(𝜂_𝑖)={𝑏^𝑖1(𝜂_𝑖),…,𝑏^𝑖𝑛_𝑖(𝜂_𝑖)} are investigator-specified functions of 𝜂_𝑖 and covariates, and Δ_𝑖𝑗 is an implicit parameter whose value is determined by the integral equation matching the conditional mean model (2.3) with the corresponding marginal model (2.1), that is,𝑔⁻¹𝛽^𝑀𝑋_𝑖𝑗=E𝑌_𝑖𝑗∣𝑋_𝑖𝑗=EE𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂_𝑖=𝑔⁻¹Δ_𝑖𝑗+𝛼^𝑏_𝑖𝑗(𝜂)𝑑𝐹(𝜂).(2.4)

The marginal mean model (2.1) directly specifies the marginal relationship between the responses and covariates, and 𝛽_𝑀 is the 𝑝×1 marginal regression parameters of main interest. Next, the semiparametric linear transformation model (2.2) is chosen to provide a flexible survival model for the informative dropout time while it still can be easily incorporated into model (2.3) for the dependence of the longitudinal responses and informative dropouts. Model (2.2) includes the proportional hazards model [25], the proportional odds model [26], and the Box-Cox transformation model as special cases and has been studied intensively in survival analysis literature [27–29]. In addition, as we present in Section 3, the explicit assumption on the error distribution in (2.2) can facilitate the “marginalization” procedure for parameter estimation.

The conditional mean model (2.3) is motivated by the construction of the marginalized random-effects model [30, 31]. As a motivating example, we consider a continuous Gaussian process following a simple random-effects model, 𝑌_𝑖𝑗=𝛽^𝑋_𝑖𝑗+𝑏_0𝑖+𝑏_1𝑖𝑡_𝑖𝑗+𝜀_𝑖𝑗, where (𝑏_𝑖0,𝑏_𝑖1) are the random intercept and slope, and error terms 𝜀_𝑖𝑗,𝑗=1,…,𝑛_𝑖 are assumed to follow 𝑁(𝟎,Σ_𝑖) but independent of (𝑏_𝑖0,𝑏_𝑖1) or 𝜂_𝑖. Note that 𝜀_𝑖𝑗's can still exhibit temporal dependence in addition to what has been accounted by the random effects, that is, Σ_𝑖 with AR(1) covariance structure. Furthermore, as in joint modeling approaches via latent variables, the joint distribution of (𝑏_0𝑖,𝑏_1𝑖,𝜂_𝑖) is assumed to be 𝑁𝟎,^{Σ𝑏𝐶,𝐶}_′,1. It is easy to see that the conditional mean has the expression as model (2.3),E𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂_𝑖=Δ_𝑖𝑗+𝐶₁𝜂_𝑖+𝐶₂𝜂_𝑖𝑡_𝑖𝑗.(2.5) We use model (2.3) primarily as a parsimonious model for the dependence structure between the longitudinal responses and informative dropout times. However, note that although model (2.3) takes a similar form as the marginalized random-effects model, it does not intend to fully specify the joint distribution of the repeated measurements since model (2.3) only specifies the conditional mean function and there is no conditional independence assumed.

Note that Δ_𝑖𝑗 is the solution of (2.4), and thus its value implicitly depends on 𝛽_𝑀, 𝛼, the formulation of 𝑏_𝑖𝑗(𝜂_𝑖) and the distribution of 𝜂_𝑖. The specification of 𝑏_𝑖𝑗(𝜂_𝑖) reflects investigator's assumptions on the dependence structure among the longitudinal responses and their association with the dropout times. It is well known that the dependence assumptions between longitudinal measurements and informative dropouts are usually unverifiable from the observed data, but it intrinsically affects the inference about 𝛽_𝑀. Thus, a sensitivity analysis under various assumptions is always warranted. It is clear that the sensitivity analysis can be easily conducted within the framework of model (2.3). For example, the sensitivity analysis may start with a large model in the specification for 𝛼^𝑏_𝑖𝑗(𝜂_𝑖), for example, 𝛼₁𝜂_𝑖+𝛼₂𝜂_𝑖𝑡_𝑖𝑗+𝛼^3𝜂_𝑖𝑋_𝑖𝑗, and then examine the statistical significance of estimates for 𝛼’s to further simplify the model. As shown in the next subsection, complex structure can be imposed on 𝑏_𝑖𝑗(𝜂_𝑖) without introducing much extra computation. Lastly, we note that the marginal interpretation of longitudinal parameters in model (2.1) is invariant under different specifications of the conditional mean model (2.3).

3. Estimation and Asymptotic Properties

3.1. Conditional Generalized Estimating Equation

First, assume that 𝜂_𝑖 is known. We construct a “conditional” generalized estimating equation for ℬ=(𝛽^𝑀,𝛼^)^. More specifically, the estimating function 𝑈(ℬ) is specified as^𝑛_𝑖=1𝐀^𝑖⋅𝜂_𝑖𝐖_𝑖𝐘_𝑖⋅−𝝁_𝑖⋅𝜂_𝑖=^𝑛_𝑖=1𝐀^𝑖⋅𝜂_𝑖𝐖_𝑖𝐘_𝑖⋅−𝑔⁻¹𝚫_𝑖⋅+𝐛_𝑖⋅𝜂_𝑖𝛼,(3.1) where 𝐘_𝑖⋅ denotes the 𝑚_𝑖×1 vector of observed responses of subject 𝑖; Δ_𝑖⋅=(Δ_𝑖1,…,Δ_𝑖𝑚𝑖)^; 𝐛_𝑖⋅(𝜂_𝑖)={𝑏_𝑖1(𝜂_𝑖),…,𝑏_𝑖𝑚𝑖(𝜂_𝑖)}^; 𝝁_𝑖⋅(𝜂_𝑖)={𝜇_𝑖1(𝜂_𝑖),…,𝜇_𝑖𝑚𝑖(𝜂_𝑖)}^; 𝐀_𝑖⋅(𝜂_𝑖)=[𝜕𝝁_𝑖⋅(𝜂_𝑖)/𝜕𝛽^𝑀,𝜕𝝁_𝑖⋅(𝜂_𝑖)/𝜕𝛼^]; 𝐖_𝑖 is a 𝑚_𝑖×𝑚_𝑖 weight matrix.

It is easy to see that 𝑈(ℬ) has mean zero at the true parameter values ℬ₀ under model (2.3). Note that the vector of marginal parameters 𝛽_𝑀 is implicitly present in 𝑈(ℬ) with Δ_𝑖⋅ through the constrain equation (2.4). Thus, the Jacobian matrix 𝐀_𝑖⋅(𝜂_𝑖) needs to be derived using both the constrain (2.4) and models (2.1) and (2.3), which is different from the ordinary GEE. More specifically, entries of the Jacobian matrix 𝐀_𝑖⋅(𝜂_𝑖) are given by 𝜕𝜇_𝑖𝑗𝜂_𝑖𝜕𝛽_𝑀𝜕𝜇_𝑖𝑗=𝑋_𝑖𝑗Υ_𝑖𝑗EΥ_𝑖𝑗𝜂_𝑖⁻¹Υ_𝑖𝑗𝜂_𝑖,𝜂_𝑖𝜕𝛼=𝑏^𝑖𝑗𝜂_𝑖−E𝑏_𝑖𝑗𝜂_𝑖Υ_𝑖𝑗𝜂_𝑖EΥ_𝑖𝑗𝜂_𝑖⁻¹Υ^𝑖𝑗𝜂_𝑖,(3.2) where Υ_𝑖𝑗=1/𝑔̇{𝑔⁻¹(𝛽^𝑀𝑋_𝑖𝑗)}, Υ_𝑖𝑗(𝜂_𝑖)=1/𝑔̇[𝑔⁻¹{Δ_𝑖𝑗+𝛼′𝑏_𝑖𝑗(𝜂_𝑖)}], and we use 𝑎(̇𝑥) to denote the derivative of a function 𝑎(𝑥) throughout this paper. In particular, we have Υ_𝑖𝑗=𝜇_𝑖𝑗(1−𝜇_𝑖𝑗) and Υ_𝑖𝑗(𝜂_𝑖)=𝜇_𝑖𝑗(𝜂_𝑖){1−𝜇_𝑖𝑗(𝜂_𝑖)} under the logit-link function for binary data; Υ_𝑖𝑗=𝜇_𝑖𝑗 and Υ_𝑖𝑗(𝜂_𝑖)=𝜇_𝑖𝑗(𝜂_𝑖) under the log-link function for count data. Thus, under these canonical link functions, Υ_𝑖𝑗=Var(𝑌_𝑖𝑗∣𝑋_𝑖𝑗) and Υ_𝑖𝑗(𝜂_𝑖)=Var(𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂_𝑖) are the marginal variance and conditional variance of the responses, respectively. In addition, these formulations also facilitate our selection of the weight matrix 𝐖_𝑖. For example, for binary longitudinal data with logit-link function, we can choose a weight matrix as 𝐖^𝑖−1=diag{Var(𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂_𝑖),𝑗=1,…,𝑚_𝑖}.

It is clear that the implementation of the estimating function (3.1) requires the knowledge of 𝜂_𝑖, which is an unknown quantity and has to be estimated first. The estimation of the semiparametric linear transformation model (2.2) has been studied by many authors [27–29]. In particular, Chen et al. [27] proposed a class of martingale-based estimating equations, ^𝑛_𝑖=1^∞0𝑍_𝑖𝑑𝑁_𝑖(𝑡)−𝐼𝑇^∗𝑖≥𝑡𝑑Λ𝜃^𝑍_𝑖+𝐻(𝑡)=0,_𝑛_𝑖=1𝑑𝑁_𝑖(𝑡)−𝐼𝑇^∗𝑖≥𝑡𝑑Λ𝜃^𝑍_𝑖+𝐻(𝑡)=0,∀𝑡≥0,(3.3) where 𝑁_𝑖(𝑡)=𝐼(𝑇^∗𝑖≤𝑡,𝛿_𝑖=1). Then an iterative algorithm can be carried out to solve 𝜃 and 𝐻 simultaneously. We estimate (𝜃,𝐻) using the approach of Chen et al. [27] and shall denote the estimates as Θ≜(𝜃̂,𝐻).

3.2. Estimation Procedure for Pure Informative Dropouts

We first consider the situation of pure informative dropouts, that is, 𝛿_𝑖≡1. Define 𝜂̂_𝑖=𝐻(𝑇_𝑖)+𝜃̂^𝑍_𝑖 and replace 𝜂_𝑖’s in (3.1) with their estimated counterparts 𝜂̂_𝑖’s. Denote the resulting estimating function by 𝑈(ℬ;Θ) and define the estimator of ℬ as the solution to 𝑈(ℬ;Θ)=0. The estimation of ℬ entails an iteration between solving nonlinear equations for Δ_𝑖𝑗 and updating a Newton-Ralphson equation for ℬ. More specifically, given the current estimated value of ℬ^(𝐽) at the Jth step, we first estimate Δ^{(𝐽)𝑖𝑗} from𝑔⁻¹𝛽^{(𝐽)′𝑀}𝑋_𝑖𝑗=𝑔⁻¹Δ^𝑖𝑗+𝛼^(𝐽)𝑏_𝑖𝑗(𝜂)𝑑𝐹(𝜂),(3.4) and then update the parameters ℬ by ℬ^(𝐽+1)=ℬ^(𝐽)+^𝑛_𝑖=1𝐀^{𝑖(𝐽)}𝜂̂_𝑖𝐖^𝑖(𝐽)𝐀^𝑖(𝐽)𝜂̂_𝑖⁻¹^𝑛_𝑖=1𝐀^{𝑖(𝐽)}𝜂̂_𝑖𝐖^𝑖(𝐽)𝐘^𝑖⋅−𝝁^{(𝐽)𝑖⋅}𝜂̂_𝑖,(3.5) where 𝐀^𝑖(𝐽)(𝜂̂_𝑖), 𝐖^𝑖(𝐽), and 𝝁^{(𝐽)𝑖⋅}(𝜂̂_𝑖) are evaluated at the current parameter values ℬ^(𝐽) and Δ^{(𝐽)𝑖𝑗}. The algorithm is iterated until it converges. Because 𝜂_𝑖 is assumed to follow an explicit parametric distribution 𝐹, it greatly simplifies the marginalization procedure (3.4). We propose to use the Gaussian-quadrature approach [32] to numerically evaluate (3.4) and 𝐀^𝑖(𝐽)(𝜂̂_𝑖). Since the integrand of (3.4) is monotonic in Δ^{(𝐽)𝑖𝑗} and so is the whole integral, it is easy to calculate a large number of Δ^{(𝐽)𝑖𝑗},𝑖=1,…,𝑛,𝑗=1,…,𝑚_𝑖, in all iterative steps. Moreover, the numerical integration is only upon the one-dimensional space of 𝜂_𝑖 and requires light computation even with complex structure assumed on 𝑏_𝑖𝑗(𝜂_𝑖). The proposed iterative algorithm has been implemented using “R” codes, which are available from the authors upon request.

3.3. Estimation Procedure for Mixed Types of Dropouts

We generalize the proposed estimation function (3.1) to accommodate the situation where there are mixed informative dropouts and random censoring. More specifically, the modified estimating equation is given by𝑈^∗(ℬ;Θ)=^𝑛_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖^∗𝑖𝐘_𝑖⋅−𝝁^∗𝑖⋅𝜂^∗𝑖,𝛿_𝑖=0,(3.6) where 𝜂^∗𝑖=𝐻(𝑇^∗𝑖)+𝜃^𝑍_𝑖; the jth component of 𝝁^∗𝑖⋅(𝜂^∗𝑖,𝛿_𝑖) is 𝜇^∗𝑖𝑗𝜂^∗𝑖,𝛿_𝑖=𝛿_𝑖𝑔⁻¹Δ_𝑖𝑗+𝛼^𝑏_𝑖𝑗𝜂^∗𝑖+1−𝛿_𝑖E𝑔⁻¹Δ_𝑖𝑗+𝛼^𝑏_𝑖𝑗(𝜂)∣𝜂≥𝜂^∗𝑖,(3.7) and the Jacobian matrix 𝐀^∗𝑖⋅(𝜂^∗𝑖,𝛿_𝑖)=[𝜕𝝁^∗𝑖⋅(𝜂^∗𝑖,𝛿_𝑖)/𝜕𝛽^𝑀,𝜕𝝁^∗𝑖⋅(𝜂^∗𝑖,𝛿_𝑖)/𝜕𝛼^]. When 𝛿_𝑖=1, the ith component of 𝑈^∗(ℬ;Θ) is the same as the one in (3.1). For 𝛿_𝑖=0, the entries of 𝐀^∗𝑖⋅(𝜂^∗𝑖,0) are given by 𝜕𝜇^∗𝑖𝑗𝜂^∗𝑖,0𝜕𝛽_𝑀𝜕𝜇_∗𝑖𝑗=𝑋_𝑖𝑗Υ_𝑖𝑗EΥ_𝑖𝑗𝜂_𝑖⁻¹EΥ_𝑖𝑗(𝜂)∣𝜂≥𝜂^∗𝑖,𝜂^∗𝑖,0𝜕𝛼=E𝑏^𝑖𝑗(𝜂)Υ_𝑖𝑗(𝜂)∣𝜂≥𝜂^∗𝑖−E𝑏_𝑖𝑗𝜂_𝑖Υ_𝑖𝑗𝜂_𝑖EΥ_𝑖𝑗𝜂_𝑖⁻¹EΥ_𝑖𝑗(𝜂)∣𝜂≥𝜂^∗𝑖.(3.8) In addition, the entries of the weight matrix can be changed to Var{𝑌_𝑖𝑗∣𝑋_𝑖𝑗,𝜂^∗𝑖,𝛿_𝑖} accordingly. Conditional expectations of various functions given 𝜂≥𝜂^∗𝑖 are computed using the Gaussian-quadrature method. Let 𝜂̂^∗𝑖=𝐻(𝑇^∗𝑖)+𝜃̂^𝑍_𝑖 and replace 𝜂^∗𝑖's in (3.6) with their estimated counterparts 𝜂̂^∗𝑖’s. Denote the resulting estimating function by 𝑈^∗(ℬ;Θ). Then the estimator ℬ of ℬ can be obtained from the equation 𝑈^∗(ℬ;Θ)=0 using the same iterative algorithm described in the previous subsection.

3.4. Asymptotic Properties of ℬ

In this subsection, we establish the asymptotic properties of ℬ. Towards this end, we need the following assumptions. (C1)The covariates 𝑋_𝑖𝑗’s are bounded with probability 1. (C2)The true parameter values ℬ₀ and 𝜃₀ belong to the interior of a known compact set, and the true transformation function 𝐻₀ has a continuous and positive derivative. (C3)Let Λ(⋅) denote the cumulative hazard function of 𝜂_𝑖. Define 𝜆(𝑡)=Λ̇(𝑡) and 𝜓(𝑡)=𝜆̇(𝑡)/𝜆(𝑡). Then 𝜆(⋅) is positive, 𝜓(⋅) is continuous, and lim_{𝑡→−∞}𝜆(𝑡)=0=lim_{𝑡→−∞}𝜓(𝑡). (C4)𝜏 is finite and satisfies 𝑃(𝑇>𝜏)>0 and 𝑃(𝐶=𝜏)>0. (C5)The matrix Ω≡E{𝐀^∗1⋅(𝜂^∗1,𝛿₁)𝐖₁𝐀^∗1⋅(𝜂^∗1,𝛿₁)} is positive finite, and the number of repeated measurements 𝑚_𝑖≪𝑁.

The regularity conditions (C1)–(C4) are also used by Chen et al. [27] to derive the consistency and asymptotically normality of the estimators Θ. Condition (C5) is needed to establish the consistency and asymptotic normality of ℬ, which is given in the following theorem.

Theorem 3.1. Under conditions (C1)–(C5), with probability 1,  |ℬ−ℬ₀|→0. In addition, one has, as 𝑛→∞, √𝑛ℬ−ℬ₀⟶^𝐷𝑁0,Ω⁻¹𝑉Ω⁻¹.(3.9)

The definition of 𝑉 and a sketch of the proof for Theorem 3.1 are given in the appendix. The asymptotic variance-covariance matrix can be consistently estimated by its empirical counterpart Ω⁻¹𝑉Ω⁻¹, which can be easily obtained using the usual plug-in rule.

4. Numerical Studies

4.1. Simulations

We conducted a series of simulation studies to evaluate the finite-sample performance of our proposed approach. Consider a binary longitudinal process with the marginal probability of success as 𝑔⁻¹(𝛽₀+𝛽₁𝑡+𝛽₂𝑍), where 𝑔 was the logit-link function; observations occurred at 𝐭_𝑖={𝑡_𝑖𝑗=𝑗,𝑗=0,1,…,5}; 𝑍 was generated from a Bernoulli distribution with the success probability of 0.5, and (𝛽₀,𝛽₁,𝛽₂)=(−1.5,0.3,1). The informative dropout time 𝑇_𝑖 was generated from a linear transformation model 𝐻(𝑇_𝑖)=−𝜃𝑍_𝑖+𝜂_𝑖, where 𝜃=−0.5. We considered three distributions for 𝜂_𝑖: the standard normal distribution (N.), the extreme value distribution (E.), and the logistic distribution (L.), corresponding to the normal transformation model, the Cox proportional hazard model, and the proportional odds model, respectively, for the informative dropout time. We then generated the binary response 𝑌_𝑖𝑗 independently from a Bernoulli distribution with the success probability of 𝑔⁻¹{Δ_𝑖𝑗+𝛼𝑏_𝑖𝑗(𝜂_𝑖)}, where Δ_𝑖𝑗 was calculated to match the marginal mean value as in (2.4) and 𝛼 indicated the level of dependence. We considered several combinations to specify the dependence between longitudinal outcomes and informative dropout times. More specifically, when 𝛼=0, there was no informative dropouts; when 𝛼=0.5 (or 0.25) and 𝑏_𝑖𝑗(𝜂_𝑖)=𝜂_𝑖, the dependence existed and was linear in the latent variable 𝜂_𝑖; when 𝛼=0.25 (or 0.5) and 𝑏_𝑖𝑗(𝜂_𝑖)=𝜂_𝑖𝑡_𝑖𝑗, the dependence was present through an interaction between the latent variable and the observation time.

For each scenario, we considered samples of size 100 and 200 and conducted 500 runs of simulations. The Gaussian-quadrature approximation was calculated using 50 grid points. We first considered the situation of pure informative dropouts and generated the dropout time 𝑇_𝑖 from the transformation model with 𝐻₀(𝑡)=2{arctan(𝑡)+𝜋/2}. Under the assumptions of 𝜂_𝑖 following the normal, the extreme value, and the logistic distributions, the average numbers of repeated measurements were 3.91, 3.37, and 3.94, respectively. The estimation results on 𝛽₁ and 𝛼 are summarized in Table 1. The proposed estimators are unbiased under all simulated scenarios, and the Wald-type 95% confidence intervals all have reasonable empirical coverage probabilities. The performance of the proposed method is consistent with different distributional assumptions of 𝜂_𝑖 and different specifications of the dependence structure, and the results improve as the sample size increases.

Table 1 

Simulation results for pure informative dropouts.

Next, we consider the situation where there are mixed informative dropouts and random censoring. For simplicity, let 𝐶_𝑖 be an administrative censoring at the end of the study, that is, 𝜏=6. The informative dropout time 𝑇_𝑖 was generated from the transformation model with 𝐻(𝑡)=log(𝑡)−1 and 𝜂_𝑖 followed the standard normal distribution. This yielded the proportion of informative dropouts of 69.6% and the average number of repeated measurements about 4. Other settings were kept the same as in the previous simulations. The simulation results are presented in Table 2. Again, the proposed approach gives unbiased parameter estimates and reasonable coverage probabilities under all the scenarios. For comparison, we also implemented the ordinary GEE method [16]. When the informative dropout is absent, that is, 𝛼=0, the GEE method yields consistent parameter estimates of 𝛽₁ as expected. But when there is informative dropout (𝛼≠0), the performance of the GEE method deteriorates quickly as the magnitude of the dependence between the longitudinal data and informative dropout increases.

Table 2 

Simulation results for mixed types of dropouts.

Last, we conducted sensitivity analysis for the proposed approach and our simulations consisted of two parts. First, as discussed in Section 2, to better characterize the dependence structure between longitudinal responses and informative dropouts, we would suggest to start with a large model in the specification for 𝛼^𝑏_𝑖𝑗(𝜂_𝑖) and then examine the statistical significance of estimates for 𝛼’s to further simplify the model. We simulated data from a simple model with either 𝛼^𝑏_𝑖𝑗(𝜂_𝑖)=𝛼₁𝜂_𝑖 or 𝛼^𝑏_𝑖𝑗(𝜂_𝑖)=𝛼₂𝜂_𝑖𝑡_𝑖𝑗, and then applied the proposed approach by assuming a bigger model (2.3) as Δ_𝑖𝑗+𝛼₁𝜂_𝑖+𝛼₂𝜂_𝑖𝑡_𝑖𝑗. The simulation results are summarized in the top panel in Table 3. The proposed method can reasonably well estimate all the parameters, and in particular, could correctly indicate the unnecessary zero term. Second, we simulated data from Δ_𝑖𝑗+𝛼₁𝜂_𝑖+𝛼₂𝜂_𝑖𝑡_𝑖𝑗 but fitted misspecified models that omitted some terms. The results are summarized in the lower panel of Table 3. It is evident that misspecified models lead to biased estimates for both longitudinal regression coefficients and dependence parameters.

Table 3 

Sensitivity analysis for misspecified models under mixed types of dropouts.

4.2. Application to Renal Disease Data from MDRD Study

Here we considered a subgroup of 129 patients with low-protein diet in MDRD study B, among whom, 62 patients were randomized to the group of normal-blood-pressure control and 67 patients were randomized to the group of low-blood-pressure control. Besides the randomized intervention, other covariates included time in study (time), baseline disease progression status (Prog), baseline blood pressure (Bp), and log-transformed baseline urine protein level (log.Pt). There were 52 (40.3%) patients left the study prematurely for kidney transplant or dialysis and were treated as informative dropouts.

We applied the proposed approach to estimate the marginal effects of covariates on GFR values. To account for the possible informative dropouts, we assumed that the dependence term 𝛼^𝑏_𝑖𝑗(𝜂_𝑖) had a form of 𝛼₁𝜂_𝑖+𝛼₂𝜂_𝑖𝑡_𝑖𝑗, analogous to the joint modeling approach with latent random intercept and random slope used in Schluchter et al. [33]. We considered the situations of 𝜂_𝑖 following the standard normal, the extreme value, or the standard logistic distributions. Because the outcome was a continuous variable, we used the identity link function.

Our results are presented in Table 4 and compared with the results from the ordinary GEE [16] with an independent working correlation matrix. More specifically, the slope estimates from the proposed approach indicate a much faster decreasing rate of GFR (e.g., time Est = −0.27, SE = 0.03, under the normality assumption for 𝜂_𝑖) than the result from the ordinary GEE method (time Est = −0.14, SE = 0.03). A possible explanation is that those patients remaining under observation usually have better kidney functions and thus higher GFR values. The ordinary GEE approach that treats the observed patients as random representatives of the population tends to underestimate the degressive trend of GFR.

Table 4 

Estimates of regression coefficients for the MDRD study.

The estimates for the intervention on blood pressure control show positive effect of the low-blood-pressure control on the longitudinal GFR development. Although the results are not statistically significant, the estimates from the proposed method (e.g., Intervention Est = 0.82, SE = 1.07, under the normality assumption for 𝜂_𝑖) are about twice large of the values from the ordinary GEE method (Intervention Est = 0.35, SE = 0.90). Moreover, for the proposed approach, the results under different distributional assumptions for 𝜂_𝑖 are quite similar. The estimates of the dependency parameters (𝛼₁,𝛼₂) are positive and statistically significant. This indicates that higher GFR values are positively associated with longer dropout times in the study. In addition, our proposed approach shows that the baseline urine protein level is significantly associated with the longitudinal GFR development, but the ordinary GEE method does not show such significance. The results obtained using our proposed method are also consistent with those reported in Schluchter et al. [33].

5. Discussion

In this paper, we propose a semiparametric marginalized model for marginal inference of the relationship between longitudinal responses and covariates in the presence of informative dropouts. The regression parameters represent the covariate effects on the population level. The proposed estimators are expected to be insensitive to misspecification of the latent variable distribution [31], which is desirable pertaining to the sensitivity analysis on unverifiable assumptions for the informative dropouts. In practice, the choice between marginal models and other types of joint modeling approaches should be determined by study objective.

To estimate the regression parameters in the proposed marginalized model, we proposed a class of simple conditional generalized estimating equations and demonstrated its computational convenience. In general, a likelihood-based approach can be used to achieve more efficient inference and is also of great interest. For example, a marginalized random effects model [30, 31] can be used for the longitudinal process and a frailty model [34] can be used for the dropout time. Furthermore, latent variables (𝑏_𝑖𝑗,𝜂_𝑖) can be modeled by a copula distribution or non-Gaussian distributions [35]. The likelihood-based methods enjoy the high efficiency and facilitate the implementation of classical model selection procedures, such as AIC or BIC; however, intensive computations are often involved when the dimension of random effects is high.

Appendix

A. Proof of Theorem 3.1

Under the conditions (C1)–(C4), Chen et al. [27] established the consistency of 𝜃̂ and the uniform consistency of a transformed 𝐻(⋅), that is, sup_{𝑡∈[0,𝜏]}|exp{𝐻(𝑡)}−exp{𝐻₀(𝑡)}|=𝑜_𝑝(1). We first derive the consistency of the proposed estimator ℬ, which is the solution of the equation 𝑈^∗(ℬ;Θ)=0. Note that −(1/𝑛)(𝜕𝑈^∗(ℬ;Θ)/𝜕ℬ)=(1/𝑛)∑_𝑛𝐀^{∗𝑖⋅}(𝜂̂^∗𝑖,𝛿_𝑖)𝐖_𝑖𝐀^∗𝑖⋅(𝜂̂^∗𝑖,𝛿_𝑖) is a positive definite matrix, and by the law of large numbers and the consistency of 𝜃̂ and 𝐻, it converges uniformly to a deterministic positive definite matrix Ω(ℬ) over a compact set of ℬ. In addition, we have that (1/𝑛)𝑈^∗(ℬ;Θ) converges uniformly to a deterministic function 𝑢^∗(ℬ;Θ₀) satisfying 𝑢^∗(ℬ₀;Θ₀)=0 and −(𝜕𝑢^∗(ℬ;Θ₀)/𝜕ℬ)∣_ℬ=ℬ0=Ω(ℬ₀)=Ω. Thus, the estimating equation 𝑈^∗(ℬ;Θ)=0 exists a unique solution ℬ. Since ℬ₀ is the unique solution of 𝑢^∗(ℬ₀;Θ₀)=0, the consistency of ℬ easily follows.

To prove the asymptotic normality, by the Taylor expansion, we have√𝑛ℬ−ℬ₀=⎧⎪⎨⎪⎩−1𝑛𝜕𝑈^∗ℬ^∗;Θ𝜕ℬ⎫⎪⎬⎪⎭⁻¹1√𝑛𝑈^∗ℬ⁰;Θ,(A.1) where ℬ^∗ lies between ℬ and ℬ₀. Since ℬ is consistent and −(1/𝑛)(𝜕𝑈^∗(ℬ;Θ)/𝜕ℬ) converges uniformly to Ω(ℬ), we have that −(1/𝑛)(𝜕𝑈^∗(ℬ^∗;Θ)/𝜕ℬ) converges to Ω. Furthermore, the Taylor expansion of 𝑈^∗(ℬ₀;Θ) around Θ₀ gives1√𝑛𝑈^∗ℬ⁰;Θ=1√𝑛^𝑛_𝑖=1−𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝐘_𝑖⋅−𝝁^∗𝑖⋅𝜂^∗𝑖,𝛿_𝑖_𝑛_𝑖=1=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝝁^∗𝑖⋅𝜂̂^∗𝑖,𝛿_𝑖−𝝁^∗𝑖⋅𝜂^∗𝑖,𝛿_𝑖√𝑛𝑈^∗ℬ₀;Θ₀−^𝑛_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝜕𝝁^∗𝑖⋅(𝜂,⋅)×𝜃̂−𝜃₀𝜕𝜂||||_{(𝜂∗𝑖,𝛿𝑖)}^𝑍_𝑖+𝑂_𝑖+𝐻𝑇^∗𝑖;𝜃₀−𝐻₀𝑇^∗𝑖+𝑜_𝑝(1),(A.2) where 𝑂_𝑖=𝜕𝐻(𝑇^∗𝑖;𝜃)/𝜕𝜃∣_𝜃=𝜃0 and its asymptotic representation can be found in the appendix of Chen et al. [27]. The asymptotic representations of 𝜃̂ and 𝐻 have also been derived by Chen et al. [27],√𝑛𝜃̂−𝜃=Σ^∗−11√𝑛^𝑛_𝑖=1^𝜏0Λ𝑍_𝑖−𝜇_𝑏(𝑡)𝑑𝑀_𝑖(𝑡)+𝑜_𝑝(1),_∗𝐻𝑡;𝜃⁰−Λ^∗𝐻₀(𝑡)=𝑛^𝑛−1_𝑖=1^𝑡0𝜆^∗𝐻₀(𝑠)𝐵₂(𝑠)𝑑𝑀^𝑖(𝑠)+𝑜_𝑝𝑛^−1/2,(A.3) where 𝑁_𝑖(𝑡)=𝛿_𝑖𝐼(𝑇^∗𝑖≤𝑡) and 𝜉_𝑖(𝑡)=𝐼(𝑇^∗𝑖≥𝑡) are the counting and at-risk processes, respectively, 𝑀_𝑖(𝑡)=𝑁_𝑖(𝑡)−∫^𝑡0𝜉_𝑖(𝑠)𝑑Λ{𝜃^0𝑍_𝑖+𝐻₀(𝑠)} is the mean-zero martingale process, and the definitions of Σ_∗ and the functions Λ^∗(⋅), 𝜆^∗(⋅), 𝐵₂(⋅), and 𝜇_𝑏(⋅) are given in the appendix of Chen et al. [27].

Plugging these terms back to the expansion of (1/√𝑛)𝑈^∗(ℬ₀;Θ) specified in (A.2), some rearrangement yields that it is equal to 1√𝑛^𝑛_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝐘_𝑖⋅−𝝁^∗𝑖⋅𝜂^∗𝑖,𝛿_𝑖−^𝜏0𝑄₁𝑍_𝑖−𝜇_𝑏(𝑡)+𝑄₂(𝑡)𝑑𝑀_𝑖(𝑡)+𝑜_𝑝(1),(A.4) where 𝑄₁=lim_𝑛→∞𝑛^𝑛−1_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝜕𝝁^∗𝑖⋅(𝜂,⋅)𝜕𝜂||||_{(𝜂∗𝑖,𝛿𝑖)}𝑍^𝑖𝑄₂+𝑂_𝑖^Σ^∗−1,(𝑡)=lim_𝑛→∞𝑛^𝑛−1_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝜕𝝁^∗𝑖⋅(𝜂,⋅)𝜕𝜂||||_{(𝜂∗𝑖,𝛿𝑖)}𝐵𝑡,𝑇^∗𝑖𝜉_𝑖(𝑡)𝐵₂(𝑡).(A.5) The definition of 𝐵(𝑡,𝑠) can be found in Chen et al. [27].

Hence, (1/√𝑛)𝑈^∗(ℬ₀;Θ) has been written as a standardized summation of independent terms with mean zero. By the central limit theorem, it is asymptotically equivalent to a multivariate Gaussian variable with zero mean and covariance matrix 𝑉, which is the limit of 1𝑛^𝑛_𝑖=1𝐀^{∗𝑖⋅}𝜂^∗𝑖,𝛿_𝑖𝐖_𝑖𝐘_𝑖⋅−𝝁^∗𝑖⋅𝜂^∗𝑖,𝛿_𝑖−^𝜏0𝑄₁𝑍_𝑖−𝜇_𝑏(𝑡)+𝑄₂(𝑡)𝑑𝑀_𝑖(𝑡)^⊗2.(A.6) From (A.1), it is easy to see that the estimator ℬ is asymptotically normal with mean zero and the variance-covariance matrix Ω⁻¹𝑉Ω⁻¹, which can be consistently estimated by its empirical counterpart  Ω⁻¹𝑉Ω⁻¹ using the usual plug-in rule.

Acknowledgments

The authors would like to thank the editor and the referees for their instructive comments. This research was supported partially by NIH grants RO1 CA-140632 and RO3 CA-153083.