|
Article Excerpt
1. INTRODUCTION
Growth curve modeling (Rao 1958) has become a standard tool for examining temporal patterns of repeated measures. Growth curves use various procedures, including random-effects modeling (Laird and Ware 1982; Longford 1993), hierarchical linear modeling (Bryk and Raudenbush 1992), and latent growth modeling (Muthen 1983, 1984, 1996, 1997; Muthen and Curran 1997). All of these related modeling techniques use normal random effects to represent the departure of an individual's latent growth parameters from the population mean growth parameters, which themselves can be modeled as functions of covariate effects, such as intervention effect and interaction with baseline over time (see Muthen et al. 2002). Increasingly, researchers from various disciplines tend to view repeated-measures data from a sample as a composite of several subsets or classes within the population, each following a different developmental pathway. Examples include studies on prostate cancer (Pearson, Morrell, Landis, Carter, and Brant 1994), bladder control (Croudace, Jarvelin, Wadsworth, and Jones 2003), aggressive behavior (Nagin and Tremblay 2001; Muthen et al. 2002), reading achievement (Muthen, Khoo, Francis, and Boscardin 2003), depression (Dew et al. 1997), alcohol use (Muthen and Shedden 1999; Muthen 2000; Muthen and Muthen 2000; Li, Duncan, and Hops 2001; Duncan, Duncan, Strycker, Okut, and Li 2002), and externalizing and internalizing behaviors (Brendgen, Vitaro, Bukowski, Doyle, and Markiewicz 2001). Class memberships of individuals may not be measured directly, but they can be inferred by their developmental courses as well as by their covariates.
In line with these conceptual models, new statistical methods provide a means of identifying these multiple developmental pathways within a population and examining covariate effects on these different population subsets (Muthen and Shedden 1999; Nagin 1999; Muthen et al. 2002). Such growth mixture modeling consists of three levels. Analogous to random-effects modeling (Laird and Ware 1982), level 1 specifies the individual-level observed data as a sum of class-adjusted fixed effects, random effects, and measurement errors specified at each observation time. Level 2 specifies the distribution of the class-specific random effects and covariate effects on the class-specific mean growth trajectories. Level 3 specifies the covariate effects on class membership using the multinomial logistic regression model (McFadden 1974). Detailed formulation of growth mixture models (GrMMs) is given in Section 2.
Most of the existing literature on validation of such mixture modeling relies solely on the number of mixtures but overlooks other parts of the model assumptions. In contrast, this article proposes a comprehensive diagnostic procedure to examine GrMM misspecification in the number of growth classes, mean growth trajectories, and covariance structures. Toward this end, we begin by surveying the literature pertinent to our work.
Model fit can be assessed in different ways, including overall goodness-of-fit indices (e.g., Pearson-type chi-squared test), model comparison measures [e.g., likelihood ratio statistics, Mallows's [C.sub.p] (1973), Akaike information criterion (AIC) (Akaike 1974), and the Bayes information criterion (BIC) (Schwarz 1978)], and residual diagnostics, the subject of this article. Residuals can be effective for detecting specific forms of model misspecification and directing analysts toward a better-fitting model. Such methods have been well established for linear regression models (Cook and Weisberg 1982; Hoaglin, Mosteller, and Tukey 1983), generalized linear models (McCullagh and Nelder 1989), and random-effects models (Waternaux, Laird, and Ware 1989; Bates and Pinheiro 2000). For example, residual diagnostics have been used in probing influential observations (Cook and Weisberg 1982; Beckman, Nachtsheim, and Cook 1987; Lange and Ryan 1989; Bradlow and Zaslavsky 1997), collinearity or omission of important covariates (Hodges 1998; Cook 1998; Atkinson and Riani 2000), and departure from normality (Dempster and Ryan 1985; Lange and Ryan 1989) or variance assumptions (McCullagh and Nelder 1989, pp. 400-401). Residuals can also be used to detect a single outlier or a cluster of outliers (Waternaux et al. 1989; Weiss and Lazaro 1992; Hall, Zeger, and Bandeen-Roche 1996; Langford and Lewis 1998), and provide guidance in fitting those outlying data. Often outliers form a cluster, which in GrMMs can be represented by a new trajectory class. Consequently, our approach relies less on outlier detection and more on checking the adequacy of the mixture components, a perspective similar to that of Lindsay and Roeder (1992).
A central goal of our research is to apply commonly used diagnostic tools to the GrMM setting. But the standard diagnostic techniques mentioned earlier cannot be immediately applied to GrMMs, because of the challenges arising from the unobserved latent class variables. Toward this end, we adopt the pseudoclass technique (Bandeen-Roche, Miglioretti, Zeger, and Rathouz 1997), and modify it to suit our purposes. For a given fit GrMM, our diagnostic method involves three steps. First, for each individual in the sample, we draw one or more random pseudoclasses according to the posterior distribution of class memberships [see (7) in Sec. 3]. Second, we compute from these pseudoclass adjusted residuals for each class (see Sec. 4). Based on our theorem given in Section 3, the empirical distributions of these pseudoclass adjusted residuals are asymptotically indistinguishable from the distributions of residuals formed with actual growth classes. Thus in the third step we examine these residuals using familiar diagnostic procedures. Two alternative methods that we have chosen not to follow are based on bootstrap sampling (see McLachlan and Peel 2000, sec. 6.6) and diagnostics using the most probable class estimate (Biernacki, Celeux, and Govaert 2000), because each requires more intensive computation to evaluate the fit of residuals.
Although the use of a single pseudoclass draw offers a feasible technique for GrMM validation, there remains room for improvement. To improve the asymptotic efficiency of residual estimates, we consider a set of new diagnostic methods drawn on a multiple pseudoclass technique (see the corollary in Sec. 3). Using a numerical approximation to these asymptotic variances (because of their nonclosed-form expressions), in Appendix C we derive a rule of thumb for choosing the minimum number of pseudoclass draws for achieving a targeted relative efficiency. Because filling in the missing latent class variables with multiple independent sets of pseudoclasses can be considered a primitive multiple imputations of missing data (Rubin 1987), it is natural to compare the residual estimates based on these two approaches. We show in (15) that residual estimates due to multiple pseudoclasses and those due to multiple imputations have similar asymptotic variances, but the multiple pseudoclass approach is computationally simpler than multiple imputations.
Three graphical methods are proposed in Section 5 to examine model adequacy using the distributions of pseudoclass adjusted residuals described in Section 4: A time trend plot of residual means for examining the mean growth structure, one involving the quantile-quantile (Q-Q) norm plot technique for examining the adequacy of the number of mixture components, and another involving the empirical Q-Q plot (Chambers, Cleveland, Kleiner, and Tukey 1983) for examining (co)variance restrictions.
To infer model misspecification, the mean departures of these residuals (from the null) and their associated pointwise 95% confidence intervals are calculated based on the formulas derived in Section 3. To further enhance the utility of the three diagnostic methods proposed, we also provide a useful model selection algorithm to build better fitting GrMMs from an initial GrMM. These methods are demonstrated with simulated data in Section 5 and are then applied to longitudinal data from a randomized preventive trial aiming at reducing aggression in elementary school children (see Sec. 6).
We adopt the following definitions and notational conventions in this article. We use class, growth class, trajectory class, or LTC interchangeably to indicate latent trajectory class. We use iid for independent and identically distributed; pdf for probability density function; pmf for probability mass function; [E.sup.x](*) for expectation of "*" with respect to the pdf of x; l(v) for the length of vector v; MLE for maximum likelihood estimator; EB for empirical Bayes; [0.sub.l] for the l X 1 zero vector; [I.sub.l] for the identity matrix of dimension l X l; [[chi square].sub.l] for a chi-squared variable with degree of freedom l; vec(*,...,*) for a vector that stacks all elements in the parentheses; MN(u, V) for the multivariate normal distribution with mean u and covariance V; [phi](u, V) [equivalent to] {(2[pi])[.sup.l.sub.u]|V|}[.sup.-.5]exp{-[1/2]u'[V.sup.-1]u}, where l(u) is the length of u; [PHI](u) for [[integral].sub.-[infinity].sup.u][phi](u', 1)du'; [L.[right arrow]] for "converges in law to as n [right arrow] [infinity]; [c.[right arrow]] for converges completely as n [right arrow] [infinity]; and [a.s.[right arrow]] for converges almost surely to as n [right arrow] [infinity]."
2. THE GROWTH MIXTURE MODEL
We assume that measurement times are the same for all subjects throughout, and we initially assume no missing data. Consider [Y.sub.i] = ([Y.sub.i1],...,[Y.sub.il(y)])' to be the l(y) repeated measures of some continuous outcome for subject i, for i = 1,...,n. Let [C.sub.i] represent the latent trajectory class for subject i, which takes values of 1,...,K. At the first level of the model, let each [Y.sub.i] given [C.sub.i] = k be modeled as a sum of a fixed component [X.sub.ik.sup.f][[alpha].sub.k] (for the class mean), a random component [X.sub.ik.sup.r][[beta].sub.ik] (for within-class variation), and a random error component [e.sub.ik], that is,
[Y.sub.i]|[.sub.[[C.sub.i] = k]] = [X.sub.ik.sup.f][[alpha].sub.k] + [X.sub.ik.sup.r][[beta].sub.ik] + [e.sub.ik], (1)
where [X.sub.ik.sup.f] and [X.sub.ik.sup.r] are the covariates associated with fixed parameters [[alpha].sub.k] and random effects [[beta].sub.ik], and [e.sub.ik] [iid.~] MN([0.sub.l(y)], [[SIGMA].sub.k]).
The second level involves modeling the distribution of [[beta].sub.ik] and explicit specification of [X.sub.ik.sup.f] and [X.sub.ik.sup.r] in terms of the temporal design structure and individual-level covariates. We assume
[[beta].sub.ik] [iid.~] MN([0.sub.l([b.sub.k])], [[PSI].sub.k]), (2)
and write
[X.sub.ik.sup.f] = [T.sub.k.sup.f][Z.sub.ik.sup.f] and [X.sub.ik.sup.r] = [T.sub.k.sup.r][Z.sub.ik.sup.r] (3)
for i = 1,...,n and k = 1,...,K, where [T.sub.k.sup.f] and [T.sub.k.sup.r] are known time measurements associated with [[alpha].sub.k] and [[beta].sub.ik], and [Z.sub.ik.sup.f] and [Z.sub.ik.sup.r] are the associated individual-level covariate matrices.
At the third level, the latent class probabilities are expressed in terms of individual-level covariates [X.sub.i.sup.c] = vec([X.sub.i1.sup.c],...,[X.sub.iK.sup.c]) and unknown fixed parameters [gamma] = vec([[gamma].sub.1],...,[[gamma].sub.K]) using multinomial logistic regression modeling (McFadden 1974). This yields
[[pi].sub.ik]([gamma]) [equivalent to] Pr([C.sub.i] = k|[X.sub.i.sup.c], [gamma]) = [exp([X.sub.ik.sup.c][[gamma].sub.k])]/[[[summation].sub.k'=1.sup.K]exp([X.sub.ik'.sup.c][[gamma].sub.k'])], (4)
where [[gamma].sub.1] is set to [0.sub.[l.sub.[gamma].sub.1]] for model identifiability. We also assume that [C.sub.1],...,[C.sub.n] are mutually independent given ([X.sub.1.sup.c],...,[X.sub.n.sup.c]).
In situations where some [Y.sub.i]'s contain missing values due to an ignorable missing mechanism (Little and Rubin 1987), we replace equations in (2) and (3) with their appropriate marginal structures for the nonmissing components. To summarize, in general, we write (1)-(3) together as
[Y.sub.ik] = [T.sub.ik.sup.f][Z.sub.ik.sup.f][[alpha].sub.k] + [T.sub.ik.sup.r][Z.sub.ik.sup.r][[beta].sub.ik] + [e.sub.ik], (5)
where [e.sub.ik] [approximately] MN([0.sub.l([y.sub.i])], [[SIGMA].sub.ik]) and [[beta].sub.ik] [approximately] MN([0.sub.l([b.sub.k])], [[PSI].sub.k]).
As an example, we construct a K-class GrMM for testing the intervention effect in a longitudinal randomized intervention trial. Let [x.sub.i] denote the indicator of intervention assignment for individual i, where [x.sub.i] = 1 if individual i is assigned to the active intervention condition and [x.sub.i] = for the control condition. Let 0, 1,...,(l(y) - 1) index the time points at which the outcomes are measured. Assume quadratic growth curves for all trajectory classes, where both intercepts and slopes are random but quadratics are fixed. (This last condition of negligible variability in the departure from linearity across subjects appears to fit a number of examples that we have examined, including tracking the course of aggressive behavior, reading achievement, and glycemic control in patients with type 2 diabetes mellitus.) Then the GrMM specification of [[alpha].sub.k], [[beta].sub.k], [[PSI].sub.k], [[SIGMA].sub.k], [T.sub.ik.sup.f], [T.sub.ik.sup.r], [Z.sub.ik.sup.f], [Z.sub.ik.sup.r], [X.sub.ik.sup.c], and [[gamma].sub.k] as defined in (1)-(4) is
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].
In the foregoing expressions, ([[alpha].sub.k1], [[alpha].sub.k2], [[alpha].sub.k3])' is the vector of mean growth factors (intercept, slope, and quadratic) associated with class k for the control group, and the intervention impact on growth factors and class membership can be assessed by ([[alpha].sub.k4], [[alpha].sub.k5], [[alpha].sub.k6])' and [[gamma].sub.k1] (the log-odds of class k relative to class K). For randomized trials, we would typically assume that intervention and control groups have the same baseline distributions for all classes, and thus the population parameters [[alpha].sub.k4] and [[gamma].sub.k1] should both be for k = 1,...,K.
We fit the foregoing GrMMs by computing MLEs for parameters appearing in the marginal likelihood,
L(v; Y, X) = [n.[product].[i=1]]f([Y.sub.i]|[X.sub.i], v), (6)
using the EM algorithm (Dempster, Laird, and Rubin 1977) written in the Mplus program (Muthen and Muthen 2004), where
f([Y.sub.i]|[X.sub.i], v) = [K.summation over (k=1)]{[[pi].sub.ik]([gamma])[phi]([Y.sub.i] - [X.sub.ik.sup.f][[alpha].sub.k], [[OMEGA].sub.ik])},
with v = vec([v.sub.1],...,[v.sub.K]), [v.sub.k] = vec([[alpha].sub.k], [[gamma].sub.k], [[PSI].sub.k], [[SIGMA].sub.k]), [X.sub.i] = vec([X.sub.i.sup.f], [X.sub.i.sup.r], [X.sub.i.sup.c]), [X.sub.i.sup.f] = vec([X.sub.i1.sup.f],...,[X.sub.iK.sup.f]), [X.sub.i.sup.r] = vec([X.sub.i1.sup.r],...,[X.sub.iK.sup.r]), [X.sub.i.sup.c] = vec([X.sub.i1.sup.c],...,[X.sub.iK.sup.c]), and [[OMEGA].sub.ik] = {[X.sup.r]}'[.sub.ik]{[[PSI].sub.k]}{[X.sub.ik.sup.r]} + [[SIGMA].sub.ik].
3. THE PSEUDOCLASS TECHNIQUE
To perform residual diagnostics for GrMM, we first randomly allocate individuals into pseudoclasses as described in Section 3.1, then examine the pseudoclass adjusted residuals with respect to their null distributions. We identify suitable residuals in Section 4 based on the theorem that we derive later in the article. Diagnostic procedures are then proposed in Section 5.
The pseudoclass technique (Bandeen-Roche et al. 1997), originally introduced in the latent class modeling literature, is used here for estimating each individual's class. Each pseudoclass assignment is a random draw based on the estimated conditional probabilities that the individual comes from each class given the data. That is, if [C*.sub.i] denotes the pseudoclass for individual i, then its pmf under (4), (5), and v being evaluated at its MLE [^.v] is
[[pi]*.sub.ik]([^.v]|[Y.sub.i]) [equivalent to] Pr([C*.sub.i] = k|[Y.sub.i], [X.sub.i], [^.v]) = [[phi]([Y.sub.i] - [X.sub.ik.sup.f][^.[alpha].sub.k], [^.[OMEGA].sub.ik])[[pi].sub.ik]([^.[gamma]])]/[[[summation].sub.k'=1.sup.K][phi]([Y.sub.i] - [X.sub.ik'.sup.f][^.[alpha].sub.k'], [^.[sigma].sub.ik'])[[pi].sub.ik']([^.[gamma]])]. (7)
The asymptotic properties associated with a single set of C* = ([C*.sub.1],...,[C*.sub.n])' drawn from (7) have been examined by Bandeen-Roche et al. (1997) in the latent class model framework; data from each pseudoclass behave approximately like samples drawn from the true classes if the model is correctly specified. As shown in the theorem, their result can be extended to our GrMM diagnostics based on a general family of pseudoclass adjusted residuals.
For the theorem that follows, we assume that [Y.sub.1],...,[Y.sub.n] are independently distributed according to (4) and (5). Let Q denote the parameter space for v, and let [^.v] be the MLE of v that maximizes (6). Denote a set of residuals, evaluated at parameter v, by {[[tau].sub.ik](v): 1 [less than or equal to] i [less than or equal to] n, 1 [less than or equal to] k [less than or equal to] K}, where [[tau].sub.ik](v) [equivalent to] [[tau].sub.k]([Y.sub.i], [X.sub.i]; v) such...
|
