PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Biom J. Author manuscript; available in PMC 2010 August 1.
Published in final edited form as:
PMCID: PMC2875790
NIHMSID: NIHMS188151

On the Impact of Parametric Assumptions and Robust Alternatives for Longitudinal Data Analysis

Summary

Models for longitudinal data are employed in a wide range of behavioral, biomedical, psychosocial, and health-care related research. One popular model for continuous response is the linear mixed-effects model (LMM). Although simulations by recent studies show that LMM provides reliable estimates under departures from the normality assumption for complete data, the invariable occurrence of missing data in practical studies renders such robustness results less useful when applied to real study data. In this paper, we show by simulated studies that in the presence of missing data estimates of the fixed-effect of LMM are biased under departures from normality. We discuss two robust alternatives, the weighted generalized estimating equations (WGEE) and the augmented WGEE (AWGEE), and compare their performances with LMM using real as well as simulated data. Our simulation results show that both WGEE and AWGEE provide valid inference for skewed non-normal data when missing data follows the missing at random (MAR), the most popular missing data mechanism for real study data.

Keywords: Augmented weighted generalized estimating equations, double robust estimate, missing at random, surrogacy assumption, weighted generalized estimating equations

1 Introduction

Models for longitudinal data are employed in a wide range of behavioral, biomedical, psychosocial, and health-care related research studies. One popular modeling paradigm is the latent variable, or random-effect, based approach for addressing correlated responses arising from such study data (Demidenko, 2004; Fitzmaurice et al. 2004; Raudenbush and Bryk, 2002). For example, for continuous response, the linear mixed-effects model (LMM) is one of the most widely used methods for modeling longitudinal data. By using random-effects to account for correlations across multiple within-subject outcomes, the LMM extends the classic multiple linear regression for modeling cross-sectional data to a general longitudinal data setting.

One major drawback of LMM as well as other random-effects based models is their dependence on distribution assumptions for inference. In recent years, many software packages have implemented robust methods to improve the validity of inference in the presence of departures from assumed parametric models. For example, simulation studies of Maas and Hox (2004; 2005) have found that estimates of the fixed-effects, or population parameters, of LMM are quite robust when model errors depart from the assumed normal distribution. However, these simulation studies have all been carried out under complete data. In most real studies, missing data invariably occurs and the robustness of estimates from LMM requires reassessment in the presence of missing data.

In this paper, we show by simulated data that the robustness of fixed-effect estimates of LMM is compromised under non-random missing data when the error term of the model deviates from the normality assumption. To address this key limitation, we investigate the performance of two robust alternatives, the weighted generalized estimating equations (WGEE) and augmented WGEE (AWGEE), for inference for non-normal longitudinal data. In Section 2, we briefly review the parametric LMM and the distribution-free formulation of this model. In Section 3, we discuss inference by the two dueling paradigms under missing data. In Section 4, we use simulated data to compare the performance of the various approaches for fitting non-normal longitudinal data under non-random missing data. In Section 5, we give our concluding remarks.

2 Models for Longitudinal Data

Over the past three decades, studies in biomedical and behavioral sciences have evolved from simple cross-sectional study designs to modern day longitudinal trials. As longitudinal study designs use subjects as their own controls, they provide a unique opportunity to study changes of outcomes of interest over time, causal effects and disease progression, in addition to providing more power for assessing treatment differences. In this section, we briefly review the two most popular approaches for modeling longitudinal data with continuous response.

2.1 Parametric Linear Mixed-effects Model

First, consider a relatively simple longitudinal study design with only two assessments, i.e., the so-called pre-post design. Let n denote the number of subjects and yit some continuous outcome of interest from the ith subject at time t (= 1, 2). We are interested in modeling the mean response E(yit) of yit over time. One popular approach is the linear mixed-effects model (LMM), which for such a pre-post study design, has the simple form:

yit=β0+β1I{t=2}+bi+εit,biN(0,σb2),εitN(0,σ2),1in,1t2,
(1)

where I{·} denotes a set indicator, and N (μ, σ2) denotes a normal with mean μ and variance σ2. In the LMM above, since E (yit) = β0 + β1I{t=2}, βt represents the (population) mean of yit at pre-(t = 1) and post-treatment (t = 2) and is known as the fixed or population effect. The latent bi in (1) accounts for the variation across the individual subjects around the fixed effect or population mean and is known as the random effect.

Let yi = (yi1,yi2). It then follows from (1) that

V=Var(yi)=(σ2+σb2σb2σb2σ2+σb2)=(σ2+σb2)C2(ρ),
(2)

where ρ=σb2/(σb2+σ2) is known as the intraclass correlation (ICC) and Cm (ρ) denotes a m × m compound symmetry correlation matrix with a correlation coefficient ρ. The LMM for the pre-post study design will be used to illustrate the performance of LMM in Section 4.

More generally, the LMM for a longitudinal data with n assessments has the following form:

yit=xitβ+zitbi+εi,biN(0,b),εiN(0,σ2Im),1in,1tm,
(3)

where xit = (1, xi1t …, xipt) (zit = (1, zi1t …, zilt)) denotes a p × 1 (l × 1) vector of covariates, bi = (bi1,…, bil) a l × 1 vector of normal random-effect, εi = (εi1,, εim) the error term for the model, N (μ, Σ) a multivariate normal with mean μ and variance Σ, Im the m × m identity matrix. Under (3), we have:

Vi=Var(yi|xi,zi)=ZibZi+σ2Im.
(4)

By setting xit=(1,I{t=2}) and zit = 1 in (4), we immediately obtain (2).

Maximum likelihood (ML) is the most popular inference procedure for LMM (Demidenko, 2004; Raudenbush and Bryk, 2002). One major drawback of ML estimate is the dependence on the parametric assumptions; if data does not follow the normality assumptions in (1), model estimates may become biased or inconsistent. In recent years, many packages have adopted the sandwich variance estimate to address this issue (Goldstein, 1995; Rasbash et al., 2000; Raudenbush and Bryk, 2002). In this case, these procedures essentially yield variance estimates equivalent to those from the generalized estimating equations (GEE), which we review next.

2.2 Distribution-free Linear Model

Since the seminal work of Liang and Zeger (1986), the generalized estimating equations (GEE) approach has been widely used as an alternative for modeling longitudinal data. By modeling the marginal mean of the response at each assessment time, GEE eliminates both layers of distribution assumptions for the random effect and error term, thereby providing consistent estimates regardless of data distributions and the complexities of structure of correlated responses.

Within our context, consider the following linear model:

E(yit|xit)=μit=xitβ,1in,1tm.
(5)

In (5), only the marginal mean at each time t is specified, which models the fixed-effect of LMM in (1) at time t. Under GEE, inference is based on the following score-like vector equation:

Wn(β)=i=1nWni(β)=i=1nGi(xi)Si(β)=i=1nGi(xi)(yiμi)=0,
(6)

where μi = (μil, …, μim), Gi (xi) is some matrix function of xi=(xi1,,xim). If the model in (5) is correct, then the GEE in (6) is unbiased, i.e., E [Wn (β)] = 0, regardless of the choice of Gi (xi), ensuring the consistency of estimate of β obtained as the solution to (6) (Liang and Zeger, 1986; Diggle et al. 2002). In most applications, Gi (xi) has the form:

Gi(xi)=DiV1(α),Di=βμi,1in,
(7)

where V (α) denotes a working variance matrix parameterized by α.

The phrase working variance is used to emphasize the fact that V (α) is not necessarily the true variance matrix of yi. For example, the simplest choice is V =diagt (Var(yit)), where diagt (αt) denotes a diagonal matrix with αt on the tth diagonal. In this case, the correlated responses yit are treated as if they are independent. In addition, there is no parameter associated with this particular working independence model. Another popular choice is the exchangeable or uniform compound symmetry correlation matrix, V(α)=diagt(Var(yit))Cm(ρ)diagt(Var(yit)), where Cm(ρ) denotes a m × m matrix correlation matrix with a common correlation ρ for any pair of the component responses of yi.

In addition to β, (6) also depends on α, though we have suppressed this dependence to highlight the fact that (6) is the equation for estimating β. Thus, before proceeding with inference about β, α must be estimated (except for the working independence model). Although the consistency of [beta] does not depend on how α is estimated, judicious choices of the type of estimates of α are required to ensure the asymptotic normality of [beta]. In particular, if [alpha] is n-consistent, [beta] is asymptotically normal with its asymptotic variance Σβ given by (e.g. Liang and Zeger, 1986; Kowalski and Tu, Chap. 4, 2007):

β=B1E(GiSiSiGi)B,B=E(βWni(β,α)).

A consistent estimate of Σβ is obtained by substituting consistent estimates in place of the respective quantities above. For example, we can estimate B by B^=i=1nβWni(β^,α^)/n. Since moment estimates are n-consistent, α is readily estimated for the working independence and exchangeable correlation models (e.g. Liang and Zeger, 1986).

3 Inference under Missing Data

In longitudinal studies, missing data are inevitable; subjects may simply quit or they may not show up at follow-up visits. We characterize the impact of missing data on model estimates through assumptions or missing data mechanisms, which allow us to ignore the multitude of reasons for missing data and focus on addressing their impact on estimation of model parameters. The missing completely at random (MCAR) assumption models a class of missing data that does not affect model estimates when completely ignored. For example, in a treatment study, missing data resulting from patient's relocation or scheduling conflict falls into this category. However, MCAR is not a plausible model when missing data are associated with treatment interventions such as patients' deteriorated or improved health conditions due to treatment. By modeling the occurrence of missing data as a function of observed responses prior to the assessment point, the missing at random (MAR) assumption addresses this class of treatment related or response-dependent missing data.

Within the longitudinal study setting in Section 2, define a missing (or rather, observed) data indicator as:

rit={1ifyitis observed0ifyitis missing,ri=(ri1,,rim),1tm,1in.

We assume no missing data at baseline t = 1 such that ri1 = 1 for all 1 ≤ in. Below, we first briefly review inference for the parametric LMM in (3) and then turn our attention to the distribution-free version in (5).

3.1 Parametric Model

Let yi = (yi1, …, yim) and xi=(xi1,,xim). Let yio and yim denote the observed and missing responses, respectively. Under likelihood based parametric inference, we jointly model the response yi and missing data indicator ri.

The joint density function, f (yi, ri | xi), can be factored into the product of marginal and conditional distributions:

f(yi,ri|xi)=f(yi|xi)f(ri|yi,xi).
(8)

Under MAR, the distribution of ri depends only on the observed response, yio, and thus:

f(ri|yi,xi)=f(ri|yio,yim,xi)=f(ri|yio,xi).
(9)

It follows from (8) and (9) that:

f(yio,ri|xi)=f(yio,yim|xi)f(ri|yio,xi)dyim=f(ri|yio,xi,θy|r)f(yio,yim|xi)dyim=f(yio|xi,θy)f(ri|yio,xi,θy|r).
(10)

If θy and θy|r are assumed disjoint, then following (10) the log-likelihood based on the joint observations (yio,ri) is given by:

l(θ)=i=1nlog(f(yio|xi,θy))+i=1nlog(f(ri|yio|xi,θy|r))=l1(θy)+l2(θy|r).

Thus, inference about θy can simply be based on the log-likelihood l1 (θy).

Most packages provide inference about the parameters of interest θy based on maximizing the likelihood function l1 (θy). Under the model assumptions of LMM, estimates are consistent under both MCAR and MAR. When study data fail to follow the parametric assumptions, maximum likelihood estimates are no longer guaranteed to be consistent. We examine bias from such estimates using simulated data in Section 4.

3.2 Distribution-free Model

3.2.1 Weighted Generalized Estimating Equations

In the presence of missing data, we may apply the GEE in (6) to the observed responses, i.e.,

Wn(β)=i=1nWni(β)=i=1nGi(xi)RiSi=i=1nGi(xi)Δi(yiμi)=0,
(11)

where Δi = diagt(rit). However, the vector estimating equation in (11) is generally biased, i.e., E(Wn(β)) ≠ 0, unless missing data follow the MCAR model. To obtain consistent estimates of β under MAR, we must revise the GEE above.

To illustrate the basic idea for modification, consider the relatively simple, pre-post design, with a homogeneous sample. We are interested in estimating the mean response at pre- and post-assessment, μ = E (yi) = (E (yi1), E (yi2)). By selecting the Gi (xi) according to (7), it follows from (11) that

Wn(μ)=R(α)1(100ri2)[i=1n(yi1μ1yi2μ2)]=0.
(12)

Solving the equations above for μ yields:

μ^=(μ^1,μ^2)=(1ni=1nyi1,1n2i=1nri2yi2),n2=i=1nri2.
(13)

If the missingness of yi2 depends on yi1, it is readily checked that E([mu]2) ≠ μ2, implying that [mu]2 is not a consistent estimate. This is also clear on intuitive grounds. For example, if yi1 and yi2 are positively correlated with higher values of yi1 leading to missing yi2, [mu]2 in (13) will be downwardly biased, since it only averages over the observed yi2 corresponding to lower values of yi1. In treatment studies, this type of response-dependent missingness often occurs if a patient feels that his/her health condition has improved (or deteriorated) during study and decides not to undergo any additional treatment.

Under MAR, the missingness of yi2 only depends on yi1, i.e.,

πi2=P(ri2=1|yi1,yi2)=P(ri2=1|yi1).

This probability πi2 selects which yi2's are to be observed based on the values of yi1. Thus, each ith subject observed at t = 2 represents a subgroup of 1/πi2 subjects with the same baseline value yi1, but unobserved at post-treatment because of the selection process defined by πi2. By augmenting each observed response yi2 at t = 2 with the weight function 1/πi2, we can statistically include the missing responses in the estimation of μ2 by using a weighted GEE (WGEE):

Wn(μ)=R(α)1(100ri2πi2)[i=1n(yi1μ1yi2μ2)]=0.
(14)

It is readily checked that E(Wn(μ)) = 0, enabling (14) to yield a consistent estimate of μ2, μ2^=(i=1nri2πi2yi2)/n (e.g. Kowalski and Tu, Chap. 4, 2007). We can also directly verify this:

μ^2=1ni=1nri2πi2yi2pE(ri2πi2yi2)=E[yi21πi2E(ri2|yi1,yi2)]=E[yi21πi2E(ri2|yi1)]=μ2,

where →p denotes convergence in probability.

By comparing (12) with (14), it is seen that the latter differs from the former only in the definition of Δi. By carrying this modification over to a general setting with m assessments, we obtain the WGEE for inference about β for the distribution-free LMM in (5), which is defined by the same vector equation in (11) except for substituting the following modified Δi:

πit=P(rit=1|xi,yi),Δit=ritπit,Δi=diagt(Δit),1tm,1in.
(15)

It is again readily checked that E [Gi (xi) ΔiSi] = 0, ensuring that the WGEE yields consistent estimates of β (e.g. Robins et al. 1995; Kowalski and Tu, Chap. 4, 2007).

To use WGEE, we must know or have estimates of πit. In some cases, subject dropout is created by study design and πit are known. For example, in some multi-stage trials, patients can only enter the next stage of the study if they satisfy certain criteria such as response to treatment at the previous stage. However, as noted earlier, in most studies, missing data patterns are defined by a host of factors not directly related to study design. We discuss estimation of πit for the general setting after introducing another robust approach for the distribution-free LMM.

3.2.2 Augmented Weighted Generalized Estimating Equations

The WGEE discussed in Section 3.2.1 depends on the model for missing data in (15). In most studies, πit are unknown and must be modeled and estimated. If such a model is misspecified, the WGEE estimate may be inconsistent. In applications, reliable models may also exist for directly relating the missing response to the observed ones and other covariates. The augmented WGEE (AWGEE) is developed to take advantage of this additional source of modeling information to ensure valid inference when the model for πit may be incorrect (Robins et al., 1995; Tsiatis, 2006).

To illustrate, consider again the pre-post study design. Suppose that we can predict yi2 directly based on yi1 using a linear regression:

E(yi2|yi1)=η0+η1yi1,1in.
(16)

Then, we can estimate μ2 without using WGEE by μ2=i=1n(η^0+η^1yi1)/n. This new estimate is consistent if (16) is a correct model, since

μ2=η^0+η^11ni=1nyi1pη0+η1E(yi1)=E[E(yi2|yi1)]=E(yi2)=μ2.

By combining both the prediction model in (16) and the WGEE in (14), we obtain an augmented WGEE to estimate μ as follows:

Wn(μ)=i=1nR(α)1(ΔiSiΔicSi),Δic=ΔiI2,
(17)
Si=(yi1μ1,y^i2μ2),y^i2=η0+η1yi1.

It is readily checked that E[Gi(xi)(ΔiSiΔicSi)]=0 if either (15) or (18) or both are correct. Thus, the AWGEE above yields consistent estimates of μ if at least one of these models is correct. Further, when both models are correct, the AWGEE estimate from (17) may also be more efficient than the WGEE estimate (Robins et al, 1995; Tsiatis, 2006).

The above is readily extended to a more general setting where the prediction model in (16) also involves other baseline covariates. Let ui be a set of baseline variables including yi1 and the prediction model be defined by:

E(yi2|ui)=η0+uiη1,1t2,1in.
(18)

The AWGEE is defined by

Wn(β)=i=1nWni(β)=i=1nGi(xi)(ΔiSiΔicSi)=0,
(19)

where y^i2=η0+uiη1. To ensure consistent estimation for the regression model, we assume a surrogacy-type assumption, [yi2 | ui, xi] = [yi2 | ui], where [yi2 | vi] denotes the conditional distribution of yi2 given vi (e.g. Prentice, 1989; Kowalski and Tu, 2002). Of course, the condition holds if ui includes xi. Under the surrogacy condition, E[Gi(xi)(ΔiSiΔicSi)]=0, if either (15) or (18) or both are correct. Thus, the AWGEE in (15) yields consistent estimates of β if at least one of the missing data models is correct.

Although feasible in principle, it is more complex to implement AWGEE for a general longitudinal study with more than two assessments. For example, for m = 3, we need to consider two missing data patterns when predicting missing yi3: one with observed yi1 and yi2, and the other with observed yi1 only. The number of prediction models grows rapidly as the frequency of assessments increases. Further, it is more intricate to specify the prediction models than models for the missing response probabilities πit.

3.2.3 Estimation of Weight Function and Augmented Term

Under MCAR, ri are independent of xi and yi and πit = P [rit = 1] = πt. In this case, πt are readily estimated by the sample moment: π^t=(i=1nrit)/n(1tm). In many studies, however, πit are dependent of either xi or yi or both. It is difficult to model πit as a function of xi and yi without imposing some additional assumptions regarding the relationship between them. As in the literature, we focus on the MAR mechanism.

As noted earlier, missing data in longitudinal trials often occur as the result of subject dropout due to deteriorated/improved health conditions and other related conditions, exhibiting the so-called monotone missing data pattern (MMDP). The structured patterns under MMDP make it possible to model πit in most studies.

Under MMDP, if yit is observed at time t, then all yis at all earlier times s (< t) are also observed. Let

Hit={xis,yis;1st1},xit=(xi1,,xi1(t1)),yit=(yi1,,yi(t1)),2tm.

The subset Hit contains all observed data prior to time t. Under MAR,

πit=P(rit=1|xi,yi)=P(rit=1|Hit)=P(rit=1|xit,yit).
(20)

Thus, under MMDP and MAR, πit are a function of observed data only, making it possible to estimate these selection probabilities.

Let pit = E (rit = 1 | ri(t−1) = 1, Hit) denote the one-step transition probability of the occurrence of missing data. Then, by invoking MMDP, it is readily checked that

πit=s=2tpis,2tm,1in.
(21)

Thus, we can estimate πit by modeling the pit. Since pit is the probability of a binary response, we model that using logistic regression:

logit(pit(αt))=logit(E(rit=1|ri(t1)=1,Hit))=ξtwit,2tm,
(22)

where wit=(1,xit,yit) and wi=(wi2,,wim). For each t, we can estimate ξt by maximum likelihood or GEE conditional on the observed wit at t − 1 (2 ≤ tm).

For AWGEE, we again consider the pre-post study design. In this case, we can readily estimate η in (18) using GEE, where η=(η0,η1). With an estimate [eta w/ hat], we can predict yit By using [eta w/ hat] and [Xi w/ hat] = ξ^2, we can construct the AWGEE in (19).

3.2.4 Inference for WGEE and AWGEE Estimates

For inference based on WGEE, let Δi, (ξ) be modeled as in (21) and (22). If ξt is estimated by GEE or maximum likelihood, ξ^=(ξ^2,,ξ^m) is the solution to the following vector estimating equation:

i=1nqi(ξ)=i=1n(qi2,,qim)=0,qit=ξt{ri(t1)[ritlog(pit)(1rit)log(1pit)]}.
(23)

Now, let

B=E(βWni),C=E(ξWni),F=E(ξqi(ξ)),

where Wni is defined in (11). Then, as shown in Appendix A, under n-consistency of [alpha], the WGEE estimate [beta] is asymptotically normal with the asymptotic variance given by

β=B1(Var(Wni)+Φ)B,Φ=CF1Var(qi)FCE(WniqiFC)[E(WniqiFC)].
(24)

In (24), Φ accounts for the variability of estimated [Xi w/ hat]. We can estimate Σβ by substituting consistent estimates in place of the respective quantities.

For AWGEE inference, we also need to estimate η for the prediction model in (18). Using GEE, the vector estimating equation is given by:

i=1ngi(η)=i=1nη{ri2[yi2(η0+uiη2)]}=0.
(25)

Let φ^=(φ^1,φ^2) with ϕ1 = ξ and ϕ2 = η. Then, by combining (23) and (25), [phi with hat] can be expressed as the solution to the following joint vector estimating equation:

i=1nsi=i=1n(qi(ϕ1),gi(ϕ2))=0,si=(qi(ϕ1),gi(ϕ2)).

The AWGEE estimate [beta] is also asymptotically normal under n-consistency of [alpha], with the asymptotic variance having the same form as in (24) except for substituting si for qi and redefining C=E(φWni) and F=E(φSi(φ)). Again, we can estimate the asymptotic variance by substituting consistent estimates in place of the respective parameters.

4 Application

We illustrate our considerations with both real and simulated data. We first present an application to data from a longitudinal study in depression research and then investigate the performance of the approach with small to moderate sample sizes by simulation. In all the examples, we set the statistical significance at α = 0.05. All analyses are carried out using a code we have developed for implementing the proposed approach using the R software platform (Free Software Foundation, 1999). This code is available from the author upon request.

4.1 Real Study

In a study on geriatric depression and associated medical comorbidities for old primary care patients, 744 subjects were enrolled from private practices and University-affiliated clinics in general internal medicine, geriatrics, and family medicine in Monroe County, New York (Lyness et al., 2007). All patients age 65 years and older who presented for care on selected days and were capable of giving informed consent. Enrolled subjects underwent semi-structured interviews, administered by trained raters in the subjects' homes or in research offices at the UR Medical Center. The raters' subject interviews included assessments of cognition, functional status, and psychopathology, the latter including the Structured Clinical Interview for DSM-IV (SCID) (Spitzer et al. 1986). Interviews and chart reviews were conducted at study intake, and again at one- and two-year follow-up time points.

In geriatric research, overall functional disability is of particular importance, as it reflects both the mental and physical health conditions of the individual. Primary measures of overall functional status include the Instrumental Activities of Daily Living (IADL), Physical Self-Maintenance Scales (PSMS), Global Assessment of Functioning (GAF), and the Karnofsky Performance Status Scale (KPSS) (Lawton MP and Brody, 1969, Karnofsky DA and Burchenal JH, 1949, Ware JE, Jr. and Sherbourne CD). For illustration purposes, we analyzed the change of IADL from baseline to one-year follow-up, as this measure assesses instrumental activities such as shopping or using the telephone and is particularly popular in geriatric research. Further, we only included the baseline value as a predictor when modeling the missingness of this outcome as well as the outcome itself at the follow-up using the respective logistic (20) and linear (18) models.

Of the 744 enrolled, 468 completed the IADL at the one-year follow-up. Shown in Table 1 are the estimates of the intercept and slope from the fitted logistic regression for modeling the missingness and the linear regression for modeling the outcome of IADL as a function of its baseline value at the follow-up. The baseline IADL was significant in both models, indicating that it did predict the occurrence of missing IADL as well as the outcome itself at the follow-up. Note that the negative sign of the estimate of the coefficient for baseline IADL in the logistic model indicates that the subjects with lower baseline IADL were more likely to come for assessment at the one-year follow-up. As lower IADL is associated with poorer functioning status, the observed sample at the follow-up visit seemed to be biased towards those with more severe overall functional disability at baseline.

Table 1
Estimates of parameters of (1) logistic regression for modeling missingness at one-year follow-up, and (2) linear model for predicting IADL at one-year follow-up for the study on geriatric depression and associated medical comorbidities.

We fit the LMM in (1) and the distribution-free alternative in (5) to examine the change of IADL from baseline to the one-year follow-up. Shown in Table 2 are the estimates of the intercept β0 and slope β1 for the respective models under the different inference procedures. As the estimates of β1 were positive across the board, the mean IADL increased at the follow-up visit, indicating better functioning status for the old primary care patients in this observational study.

Table 2
Estimates of parameters of (1) linear mixed-effects model (ML), and (2) distribution-free linear model (GEE, WGEE and AWGEE) for change of IADL from baseline to one-year follow-up for the study on geriatric depression and medical comorbidities.

However, the magnitude of the estimate of β1 did vary substantially — not only between the models, but also across the different procedures within the same distribution-free linear model. The WGEE and AWGEE yielded quite similar estimates, with the latter AWGEE also providing improved efficiency, as indicated by smaller asymptotic standard errors. The GEE performed poorly, with a whopping 50% downward bias, as compared to its counterparts WGEE and AWGEE estimates. For the between-model comparison, the ML estimate of β1 from the fitted LMM also incurred a downward bias, albeit with a much smaller magnitude relative to the GEE estimate. The downward bias in both cases was consistent with the fact that those assessed at the follow-up visit represented a subgroup with more severe overall functional disability at baseline.

An external file that holds a picture, illustration, etc.
Object name is nihms188151u1.jpg

Shown in the Figure is the normal-based Q-Q plot of the conditional residuals obtained from the estimated fixed and random effects for the fitted LMM model (Nobre and Singer, 2007). The plot indicates clearly that the residuals did not follow a normal distribution, which may explain the difference in the estimates of the fixed effects between the parametric (LMM) and distribution-free models (WGEE and AWGEE).

4.2 Simulation Study

Given the discrepant estimates between LMM and WGEE (AWGEE) for the real study data in 4.2, we conducted a simulation study with a pre-post study design to investigate this issue further. We considered two non-normal distributions for the model error term of the LMM: a rescaled central chi-square with one degree of freedom and a uniform between −1 and 1. Since the results are quite similar, we only discuss and report the results from the chi-square-distributed error. To examine the performance of the models under small, moderate and large samples, we performed the simualation study with three sample sizes: n = 50, 100 and 2,000. All simulations were performed with a Monte Carlo sample of 1,000 using the R software (R Development Core Team, 2007).

We considered the pre- and post-treatment design and simulated the outcome according to the LMM in (1) by setting β0 = β1 = 1 and [set membership]it (t = 1, 2) to follow the rescaled chi-square distribution, (χ12+1)σ2/2, where χ12 denotes a central chi-square with one degree of freedom. We varied σb2 and σ2 to control the within-subject correlation ρ=σb2/(σb2+σ2). We assumed no missing data at baseline t = 1 and simulated the missing response at post-treatment t = 2 under MAR according to the following logistic regression:

logit(πi2)=logit(P(ri2=1|yi1))=ξ0+ξ1yi1,1in.
(26)

We set ξ0 = 0.5 and ξ1 = 1.2 to create about 25% missing response yi2 at t = 2.

Under (1), it is readily checked that regardless of the distributions for bi and [set membership]it (see Appendix)

E(yi2|yi1)=η0+η1yi1,η0=β0(1σb2σb2+σ2)+β1,η1=ρ=σb2σb2+σ2,1in.
(27)

The above was used to predict missing yi2 from yi1 for AWGEE inference about β as discussed in Section 3.2.4. To study the effect of wrong weight function on WGEE and the robustness of AWGEE in such a scenario, we also estimated πi2 under an incorrect model by leaving out yi1 in (26).

We considered the null H0 : β0 = β1, i.e., the mean at post-treatment is twice that at pretreatment. We tested H0 using the Wald statistic, Qn2=nβ^K(K^βK)1Kβ^, which has an asymptotic central χ12 distribution. We estimated the type I error rate α based on the distribution of Qn2 from 1, 000 Monte Carlo (MC) replications, α^=(j=11000I{Qnjq0.95})/1000, where Qnj2 denotes the value of Qn2 from the jth MC replication and q0.95 the 95th percentile of χ12. The maximum likelihood (ML) inference about β = (β0,β1) for LMM was obtained from the LME procedure in R, while the WGEE and AWGEE estimates for the distribution-free alternative in (5) were computed based on the asymptotic results in Section 3.2.

Shown in Table 3 are the estimates of β and associated asymptotic standard errors averaged over 1,000 MC replications obtained from ML, GEE, WGEE and AWGEE for the respective models, with a within-subject correlation ρ = 0.1 (or σb2=2/9 and σ2 = 2). The results confirmed that the baseline mean β0 is consistently estimated by all four procedures. For β1, ML yielded consistent estimates for the normal distributed error, but biased estimates under the rescaled χ12 error. As expected, GEE estimates were biased under both types of error distributions. Note that in the rescaled χ12 error case, the standard errors of these estimates did not increase much, making the upwardly biased estimates yield false significant results in practice.

Table 3
Averaged estimates of β over 1,000 Monte Carlo replications along with asymptotic standard errors (s.e.) and type I error rates α for sample size 50, 100, 2000, with about 25% missing data at post-treatment based on ML for linear mixed-effects ...

Under the correct weight function, WGEE performed well. When the incorrect constant weight function was used, WGEE yielded biased estimates, while the AWGEE estimates remained close to the true value of β1 under the correct prediction model. However, AWGEE did not show any significant gain in efficiency; in fact, the estimate of the slope β1 had a larger standard error under AWGEE than under WGEE across small sample sizes n.

To further investigate the relative efficiency between WGEE and AWGEE, we replicated the above analysis by increasing the within-subject correlation. For example, shown in Tables 4 and and55 are the estimates from one such replicated analysis with ρ = 0.3 (or σb2=1 and σ2 = 2 and 0.6 (or σb2=1.5 and σ2 = 1), respectively. It is seen that all conclusions above remain the same, except for the relative efficiency between WGEE and AWGEE under the correct weight function and prediction model. As ρ increased to 0.3, AWGEE have smaller standard errors than its counterpart WGEE when n = 2000. At ρ = 0.6, not only did AWGEE show smaller standard errors than WGEE across all sample sizes, their differences also widened as compared to those with ρ = 0.1.

Table 4
Averaged estimates of β over 1,000 Monte Carlo replications along with asymptotic standard errors (s.e.) and type I error rates α for sample size 50, 100, 2000, with about 25% missing data at post-treatment based on ML for linear mixed-effects ...
Table 5
Averaged estimates of β over 1,000 Monte Carlo replications along with asymptotic standard errors (s.e.) and type I error rates α for sample size 50, 100, 2000, with about 25% missing data at post-treatment based on ML for linear mixed-effects ...

5 Discussion

We investigated the two primary modeling strategies for longitudinal continuous response, the linear mixed-effects model (LMM) and the distribution-free linear model, with respect to their performance under missing data, and illustrated our considerations using real as well as simulated study data. Our results show that LMM and the GEE procedure for the distribution-free alternative generally yield biased estimates under MAR when the normality assumption for LMM is violated. Further, as indicated by the simulation results, the standard errors of these estimates do not increase to reflect model misspecification, making inference prone to misleading findings. Thus, when modeling longitudinal data, it is important to test the MCAR assumption as discussed in Section 3.2 before applying any of the models and inference procedures considered. If the null of MCAR is rejected, WGEE and/or AWGEE should be considered, unless there is strong evidence to support the use of the alternative linear mixed-effects model.

Our simulation results also indicate that the gain in efficiency by AWGEE over WGEE depends on the magnitude of the within-subject correlation ρ. Within the context of the particular simulation model considered, AWGEE is less efficient than WGEE for small sample size under small ρ such as 0.1. But, AWGEE edged out WGEE to be a more efficient procedure as ρ increased to 0.6.

Acknowledgments

This research was supported in part by NIH grants R01-DA012249, R21-AG023956, UL1 RR024160 and R24-MH071604. We thank Ms. Bliss-Clark for her careful proofreading of the manuscript, and an anonymous, an Associate Editor and Editor Prof. Leonhard Held for their constructive comments that led to a substantial improvement in the presentation of this research.

Appendix

Under (1), we have

ρyb=Corr(yi1,bi)=Cov(yi1,bi)Var(yi1)Var(bi)=σb2σb2+σ2.
(28)

Further, the above holds regardless of the distributions for bi and [set membership]it. Now, consider the linear regression

bi=ϕ0+ϕ1yi1+δi,δi(0,σδ2),1in,

Where (0,σδ2) denotes a distribution with mean 0 and variance σδ2. It follows from (28) and the relationship between linear regression coefficients and correlation that

ϕ1=ρybVar(bi)Var(yi1)=σb2σb2+σ2σb2σb2+σ2=σb2σb2+σ2.

Also, since

0=E(bi)=ϕ0+ϕ1E(yi1)=ϕ0+ϕ1β0,

it follows that φ0=β0σb2/(σb2+σ2).

Under (1) and regardless of the distributions for bi and [set membership]it, we have

E(yi2|yi1)=β0+β1+E(bi|yi1)=β0+β1+ϕ0+ϕ1yi1.

By substituting the expressions of [var phi]0 and [var phi]0 into the above and combining the coefficients, we obtain (27).

References

1. Demidenko E. Mixed Models: Theory and Applications. New York: Wiley; 2004.
2. Diggle PJ, Heagerty P, Liang KY, Zeger SL. Analysis of Longitudinal Data. 2nd. Oxford University Press; 2002.
3. Fitzmaurice GM, Laird NM, Ware JH. Applied Longitudinal Analysis. New York: Wiley; 2004.
4. Goldstein H. Multilevel Statistical Models. Edward Arnold; London; Halsted, New York: 1995.
5. Karnofsky DA, Burchenal JH. The clinical evaluation of chemotherapeutic agents in cancer. In: MacLeod CM, editor. Evaluation of Chemotherapeutic Agents. New York: Columbia; 1949.
6. Kowalski J, Tu XM. A GEE approach to modeling longitudinal data with incompatible data formats and measurement error: application to HIV immune markers. Journal of the Royal Statistical Society, Series C. 2002;51:91–114.
7. Kowalski J, Tu XM. Modern Applied U Statistics. New York: Wiley; 2007.
8. Lawton MP, Brody EM, editors. Gerontologist. Vol. 9. 1969. Assessment of older people; self-maintaining and instrumental activities of daily living; pp. 179–186. [PubMed]
9. Lyness JM, Niculescu A, Tu XM, Reynolds CF, III, Caine ED. The relationship of medical comorbidity to depression in older primary care patients. Psychosomatics. 2007;47:435–439. [PubMed]
10. Maas CJM, Hox JJ. The influence of violations of assumptions on multilevel parameter estimates and their standard errors. Computational Statistics & Data Analysis. 2004;46:427–440.
11. Maas CJM, Hox JJ. Sufficient Sample Sizes for Multilevel Modeling. Methodology. 2005;1:86–92.
12. Nobre JS, Singer JM. Residual analysis for linear mixed models. Biometrical Journal. 2007;49:863–875. [PubMed]
13. Prentice RL. Surrogate endpoints in clinical trials: Definition and operational criteria. Statistics in Medicine. 1989;8:431–440. [PubMed]
14. Rasbash J, Browne W, Goldstein H, Yang M, Plewis I, Healy M, Woodhouse G, Draper D, Langford I, Lewis T. A User's Guide to MLwiN. Multilevel Models Project. University of London; London: 2000.
15. Raudenbush SW, Bryk AS. Hierarchical Linear Models. 2nd. Sage; Thousand Oaks, CA: 2002.
16. R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing; Vienna, Austria: 2007. URL http://www.R-project.org.
17. Robins JM, Rotnitzky A, Zhao LP. Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association. 1995;90:106–121.
18. Spitzer RL, Williams JB, Gibbon M. Structured Clinical Interview for DSM-III-R (SCID) New York: Biometrics Research; 1986.
19. Tsiatis AA. Semiparametric Theory and Missing Data. New York: Springer; 2006.
20. Ware JE, Sherbourne CD. The MOS 36-item short-form health survey (SF-36). I. Conceptual framework and item selection. Medical Care. 1992;30:473–483. [PubMed]
21. Zeger SL, Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics. 1986;42:121–130. [PubMed]