Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2703484

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Markov model with shared random effects
- 3. Likelihood functions
- 4. Simulations
- 5. Application: Nun Study
- 6. Higher order Markov chains
- 7. Conclusion and Future Work
- References

Authors

Related links

Comput Stat Data Anal. Author manuscript; available in PMC 2010 July 1.

Published in final edited form as:

Comput Stat Data Anal. 2009 July 1; 53(9): 3334–3343.

doi: 10.1016/j.csda.2009.02.007PMCID: PMC2703484

NIHMSID: NIHMS97412

Lei Yu, PhD,^{1,}^{*} Suzanne L. Tyas, PhD,^{2} David A. Snowdon, PhD,^{3} and Richard J. Kryscio, PhD^{1,}^{4}

See other articles in PMC that cite the published article.

This paper evaluates the effect of ignoring baseline when modeling transitions from intact cognition to dementia with mild cognitive impairment (MCI) and global impairment (GI) as intervening cognitive states. Transitions among states are modeled by a discrete-time Markov chain having three transient (intact cognition, MCI, and GI) and two competing absorbing states (death and dementia). Transition probabilities depend on two covariates, age and the presence/absence of an apolipoprotein E-ε4 allele, through a multinomial logistic model with shared random effects. Results are illustrated with an application to the Nun Study, a cohort of 678 participants 75+ years of age at baseline and followed longitudinally with up to ten cognitive assessments per nun.

In most longitudinal studies on progression of healthy individuals to chronic diseases, such as cancer, AIDS and dementia, the outcome of interest is a series of correlated binary or polytomous responses where these responses are observed at certain time points, sometimes several years apart. Generalized linear mixed models (GLMM) are suggested to account for the dependency among repeated follow-up waves within the same subjects, where unit-specific effects are realizations of some random effects (Stiratelli et al., 1984; Gibbons and Hedeker, 1994; Crouchley, 1995; Skrondal and Rabe-Hesketh, 2004; Salazar et al.,2007).

Salazar (2004; 2007) introduced a multi-state Markov model for longitudinal data with categorical responses. The model maintains the GLMM structure by accounting for conditional effects of covariates given the values of a single shared random/latent effect. In addition, two particular features are presented in his model. First, the dependency among observations on the same subject is addressed by assuming a first-order Markovian structure, which helps to facilitate the expression of the joint distribution of the response vector. Second, the parameterization of the transition probabilities using multinomial logistic regression provides a closed-form expression in the likelihood construction. The model provides a suitable approach to problems of identifying the risk factors associated with the progression of healthy individuals to a chronic disease with death treated as a competing event.

However, Salazar’s model approximates the joint distribution of the response variable using a conditional distribution given the baseline outcome of the response variable. Such an approach could possibly produce a so-called ‘baseline confounding’ problem (Crouchley and Davies, 1999; Ten Have et al., 2002), which might result in biased or inconsistent estimation. The model application in the Nun Study on progression of dementia indicates that among 678 subjects in the cohort, 77 are demented at baseline. These subjects (more than 10 percent) were removed from the analysis since they would contribute nothing to the likelihood when ignoring baseline. It is interesting to see how the model likelihood as well as maximum likelihood estimates (MLEs) will differ once we incorporate baseline information into the model construction.

This paper will focus on addressing this limitation by accommodating the baseline confounding in the Markov model using shared random effects approaches. In the next two sections, the model likelihood construction is discussed in detail. In sections 4 and 5, simulation studies and an application to the Nun Study data discussed in this context by Tyas et al. (2007) are presented. Comparisons are made between the two models with respect to maximum likelihood estimation. Section 6 discusses how this model structure can be modified to accommodate higher orders of the chain and the possibility of testing these chain orders using conventional likelihood ratio tests.

A generalized linear mixed model (GLMM) for a longitudinal analysis is defined as follows, let *i* denote a particular subject under study and *n _{i}* the number of repeated observations for subject

$$\eta (E({y}_{ik}\mid {\mathit{X}}_{ik},{\mathit{\gamma}}_{i}))={\mathit{X}}_{ik}^{\prime}\mathit{\beta}+{\mathit{W}}_{ik}^{\prime}{\mathit{\gamma}}_{i}$$

where *k* = 1,2,…,*n _{i}*. Here

In contrast with the general GLMM, the Markov model introduced by Salazar (2004; 2007) demonstrates two favourable features in modeling longitudinal categorical responses as a multi-state system where series of categorical outcomes are expressed in terms of states, and the onset and progression of these outcomes as transitions between the states.

First, the model relies on a ‘transitional modeling’ (Agresti, 2002) strategy by introducing a multi-state discrete-time Markov chain, which facilitates the expression of the joint distribution function. The natural development of chronic diseases can often be expressed in terms of distinct health stages and the Markov chain is a simple yet powerful tool in describing the progression of healthy individuals through these stages. Assume the Markov property that the conditional distribution *P*(*y _{ik}*|

$$f({y}_{i1},\dots ,{y}_{i{n}_{i}}\mid {y}_{i0})=f({y}_{i1}\mid {y}_{i0})f({y}_{i2}\mid {y}_{i1})\dots f({y}_{i{n}_{i}}\mid {y}_{i{n}_{i}-1})$$

Here each *y _{ik}*,

Second, by applying multinomial logit parameterization, the model provides a closed form in constructing the model likelihood function. For presentation purposes, we assume a finite stochastic system consisting of five transition states with three transient and two competing absorbing states. This corresponds to the five progression stages in the study of dementia (Tyas et al., 2007) and these are (1) intact cognition, (2) mild cognitive impairment (MCI), (3) global impairment (GI), (4) dementia and (5) death. The one-step transition probability matrix could then be presented as below

$$\left[\begin{array}{ccccc}{P}_{11}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{12}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{13}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{14}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{15}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})\\ {P}_{21}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{22}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{23}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{24}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{25}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})\\ {P}_{31}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{32}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{33}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{34}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})& {P}_{35}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

Since
$\sum _{v=1}^{5}}{P}_{sv}(\mathit{X},\mathit{\gamma})=1$ for each row of *s* = 1,2,…,5, a nominal polytomous logistic model for *P _{sv}* can be constructed as

$$log\left(\frac{{P}_{sv}({\mathit{\theta}}_{sv}\mid \mathit{X},\mathit{\gamma})}{{P}_{s1}({\mathit{\theta}}_{s1}\mid \mathit{X},\mathit{\gamma})}\right)={\alpha}_{v}+{\mathit{X}}^{\prime}{\mathit{\beta}}_{v}+{\mathit{\xi}}_{v}(s)+{\mathit{W}}^{\prime}\mathit{\gamma}$$

where *v* = 2,3,4,5. Let ** Θ** represents the set of all unknown parameters (

Following Salazar (2004; 2007) each transition probability can be postulated in the form of

$${P}_{sv}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})=\{\begin{array}{c}\frac{1}{1+{\displaystyle \sum _{h=2}^{5}}exp({\alpha}_{h}+{\mathit{X}}^{\prime}{\mathit{\beta}}_{h}+{\mathit{\xi}}_{h}(s)+{\mathit{W}}^{\prime}\mathit{\gamma})}\\ \frac{exp({\alpha}_{v}+{\mathit{X}}^{\prime}{\mathit{\beta}}_{v}+{\mathit{\xi}}_{v}(s)+{\mathit{W}}^{\prime}\mathit{\gamma})}{1+{\displaystyle \sum _{h=2}^{5}}exp({\alpha}_{h}+{\mathit{X}}^{\prime}{\mathit{\beta}}_{h}+{\mathit{\xi}}_{h}(s)+{\mathit{W}}^{\prime}\mathit{\gamma})}\end{array}$$

The first equation applies for *v* = 1, and the second for *v* = 2,···5.

The estimates produced in Salazar’s multi-state Markov model are based on a likelihood that conditions on the baseline response. The model further assumes that the distribution of the random effects ** γ** is independent of both baseline outcome and covariates. Such an approach could possibly produce a so-called ‘baseline confounding’ problem (Crouchley and Davies, 1999; Ten Have et al., 2002), which might result in biased or inconsistent estimation. To see this, the complete likelihood function for the model is

$${L}_{c}=\int \prod _{i}\prod _{t}f({y}_{it}\mid {y}_{i0},{\mathit{X}}_{i},\mathit{\gamma})dh(\mathit{\gamma}\mid {y}_{i0},{\mathit{X}}_{i})$$

Here $\prod _{t}}f({y}_{it}\mid {y}_{i0},{\mathit{X}}_{i},\mathit{\gamma})$ refers to the product of individual transition probabilities of $\prod _{k=1}^{{n}_{i}}}{P}_{{y}_{ik-1,}{y}_{ik}}(\mathit{\Theta}\mid {\mathit{X}}_{i},\mathit{\gamma})$. Under Salazar’s assumption, the following likelihood is used

$${L}_{s}=\int \prod _{i}\prod _{t}f({y}_{it}\mid {y}_{i0},{\mathit{X}}_{i},\mathit{\gamma})dh(\mathit{\gamma})$$

Crouchley and Davies (1999) argue that if it is possible for the random effects to be independent from the model covariates, the independence of random effects and baseline outcome is difficult to justify. The latent variables which contribute to the random effects are likely to be at least partially responsible for the observed baseline states. This is especially the case for a cohort with heterogeneous baseline where the assumption about the independence between random effect and baseline outcome can not be taken for granted.

Considerable literature has focused on constructing extended likelihood functions to accommodate missing data that are non-ignorable, such as informative drop-out and death in particular (Rubin, 1976; Ten Have et al., 1998; Pulkstenis et al., 1998; Ten Have et al., 2000, Gao, 2004). Sharing of random effects has been a popular approach in this respect. The method incorporates into the likelihood construct both the follow-up response and drop-out response components by assuming that the two share the same random parameters and are conditionally independent given these random effects (Ten Have et al., 1998). In essence, the model likelihood is built up using a separability approach. Suppose ** γ** is a vector of

$${f}^{\u2033}({\mathit{y}}_{\mathit{i}},{z}_{i})=\int {f}^{\prime}({\mathit{y}}_{\mathit{i}}\mid \mathit{\gamma})g({z}_{i}\mid \mathit{\gamma})h(\mathit{\gamma})d\mathit{\gamma}$$

where *f*′ (** y_{i}**|

Ten Have et al. (2002) take this approach in modeling longitudinal binary functional limitation responses. Their model considers both baseline confounding and informative drop-out, in which case the model likelihood consists of three separate components: one for baseline, one for follow-up outcomes and one for time of drop-out. These three pieces are conditionally independent given random effects and their corresponding predictor variables. The key difference between Ten Have’s model and ours is that, on the one hand, the drop-out is of little concern in our case and by omitting the drop-out component it simplifies the model construction. On the other hand, our case involves multinomial responses and the parameterization using polytomous logit under a discrete-time Markov framework is more complicated than the simple Bernoulli approach for binary outcomes.

The inclusion of baseline outcome variable completes the joint distribution function, and the equation now becomes

$$\begin{array}{l}f({y}_{i0},{y}_{i1},{y}_{i2},\dots ,{y}_{i{n}_{i}}\mid \mathit{\gamma})\\ =f({y}_{i1},{y}_{i2},\dots ,{y}_{i{n}_{i}}\mid {y}_{i0},\mathit{\gamma})f({y}_{i0}\mid \mathit{\gamma})\\ =f({y}_{i1}\mid {y}_{i0},\mathit{\gamma})f({y}_{i2}\mid {y}_{i1},\mathit{\gamma})\dots f({y}_{i{n}_{i}}\mid {y}_{i{n}_{i}-1},\mathit{\gamma})f({y}_{i0}\mid \mathit{\gamma})\end{array}$$

Using the previous example of the five-state transition system, let *π _{j}* =

$${\pi}_{j}({\mathit{\varphi}}_{j}\mid {\mathit{X}}_{B},\mathit{\gamma})=\{\begin{array}{c}\frac{1}{1+{\displaystyle \sum _{h=2}^{4}}exp({\tau}_{h}+{\mathit{X}}_{B}{\mathit{\delta}}_{h}+{W}_{B}^{\prime}\mathit{\gamma})}\\ \frac{exp({\tau}_{j}+{\mathit{X}}_{B}{\mathit{\delta}}_{j}+{W}_{B}^{\prime}\mathit{\gamma})}{1+{\displaystyle \sum _{h=2}^{4}}exp({\tau}_{h}+{\mathit{X}}_{B}{\mathit{\delta}}_{h}+{W}_{B}^{\prime}\mathit{\gamma})}\end{array}$$

Again the first equation applies for *j* = 1, and the second for *j* = 2,···4. Here the vector *ϕ** _{j}* (

$$L(\mathit{\Theta},\mathit{\Phi}\mid \mathit{X})=\int \prod _{i}\prod _{k=1}^{{n}_{i}}{P}_{{y}_{ik-1,}{y}_{ik}}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma}){\pi}_{{y}_{i0}}(\mathit{\Phi}\mid {\mathit{X}}_{B},\mathit{\gamma})h(\mathit{\gamma})d\mathit{\gamma}$$

*y _{i}*

Since the last two rows of the transition probability matrix contribute nothing to the likelihood, the range of *s* can be reduced to include only the transient states. The final likelihood function for under this 5-state system becomes

$$\int \prod _{i}\prod _{k=1}^{{n}_{i}}\prod _{\begin{array}{l}s=1\dots 3,\\ v=1\dots 5\end{array}}{({P}_{sv}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma}))}^{{\delta}_{{y}_{ik-1},s}{\delta}_{{y}_{ik},v}}{\pi}_{{y}_{i0}}(\mathit{\Phi}\mid {\mathit{X}}_{B},\mathit{\gamma})h(\mathit{\gamma})d\mathit{\gamma}$$

As an extension, note here that the vector ** γ** is not necessarily random effects per se. It could also be some reparameterization of the random effects such that the model could allow different variance covariance structures of random effects for

Furthermore, the Cholesky decomposition of a positive definite variance covariance matrix can be used to account for the correlation among random effects. To be more specific, suppose ** Σ** is a positive definite variance covariance matrix for the random effect vector

After the random vectors are integrated out, the maximized likelihood estimates (MLEs) can be calculated to make inferences about the parameters of interest. Except under some special assumptions, for example, log-log link function with log-gamma random effects (Pulkstenis, et al., 1998), these integrals have no analytical solutions. The marginal likelihood needs to be resolved using numerical approximation which can be computationally intensive, especially for models with multiple random effects. Several common techniques used to approximate this type of integrations are Laplace method (Gao 2004; Skrondal and Rabe-Hesketh 2004), Binomial approximation (Ten Have and Kunselman 1998; Ten Have et al., 2000; Ten Have et al., 2002), Numerical integration using Gauss-Hermite quadrature or adaptive quadrature (Hedeker and Gibbons, 1994; Skrondal and Rabe-Hesketh, 2004) and Monte Carlo method of importance sampling (Salazar, 2004) or Gibbs sampling (Zeger and Karim, 1991), etc.

Using simulation studies, comparisons are made with respect to parameter estimation between our extended shared random effects model and the model ignoring the baseline. The simulation is set up to have 500 subjects in each iteration and each subject with up to 10 follow-up waves. Depending on one continuous covariate: age and one binary covariate: the presence/absence of an apolipoprotein E-ε4 allele (*APOE-4*), transition probabilities are estimated by multinomial logistic regression. Considering the models under discussion are complex parametric Markov model involving a large number of parameters, one shared random intercept is considered at this time in order to achieve relative fast likelihood convergence. Three cases are examined: the random intercept following a normal distribution with small variance (*σ* = 1) and normal distributions with comparatively larger variances (*σ* = 2 and *σ* = 3).

To demonstrate the impact of the baseline confounding among cohorts with different baseline outcome structure, three separate simulation studies are implemented. The first simulation assumes a single cohort where all the subjects recruited share the same baseline state of intact cognition, regardless of the covariates of interest (*P*(*y _{i}*

The integral is approximated using the Laplace method. We used the dual quasi-Newton algorithm to optimize the log-likelihood functions, and the method is implemented using SAS^{®} NLMIXED procedure. The NLMIXED procedure provides a variety of optimization method which ranges from (1) second derivative methods like Newton Raphson where both gradients and Hessians need to be computed for the optimization, (2) first-derivative methods such as quasi-Newton where gradients are required in finding the optimum, and (3) The no-derivative method such as Nelder-Mead simplex, which only the function value is used in optimizing the underlying likelihood function. The quasi-Newton algorithm is the default optimization algorithm because “it provides an appropriate balance between the speed and stability required for most nonlinear mixed model applications” (SAS online doc). The asymptotic relative bias of the parameter estimates are presented in Table 1, Table 2 and Table 3.

Asymptotic relative bias of parameter estimates for homogeneous cohort (base state: 1=Intact cognition)

Asymptotic relative bias of parameter estimates for heterogeneous cohort, dependent on covariates of interest (base state: 1=Intact cognition)

Asymptotic relative bias of parameter estimates for heterogeneous cohort, with different number of demented at baseline (base state: 1=Intact cognition)

The results for the cohort with homogeneous baseline show that the simulated estimates are almost identical between two models for the parameters associated with *age* and *APOE-4*. At *σ* = 1 for example, the averaged relative biases for the covariate age are both 3.5% and the biases for *APOE-4* positive are 3.0% and 3.2% respectively. The biases for the prior states show a little more fluctuation, but are still quite close.

In contrast, for the heterogeneous cohort where the baseline states depend on the model covariates, the maximum likelihood estimates produced by the two models are quite different. In the case of random intercept variance being 1, the relative biases for the covariate *age* range from −0.7% to 0.5% under our proposed model, while the model ignoring the baseline gives the relative biases ranging from −8.4% to −4.6%. This result indicates that the parameter estimates associated with *age* are underestimated in the model without baseline structure. This is true across different random intercept variance. A similar conclusion can be made with respect to the covariate *APOE-4*, in which case the averaged relative biases from the two models are 2.3% and − 4.2% respectively at *σ* = 1, 3.4% and −18.4% at *σ* = 2, and 5.9% and −35.8% at *σ* = 3. Hence, ignoring the baseline is likely to create a serious downward bias which is likely to increase with *σ* while accounting for the baseline produces a smaller bias.

There are some variations of bias in the maximum likelihood estimation under both models as number of subjects demented at baseline changes. As shown in Table 3, the averaged relative bias under the model with and without the baseline tend to increase as the percentage of subjects demented at baseline gets larger. However such increase is much more conspicuous among those under the model without the baseline. For example, in the case where the cohort has 14% of subjects demented at baseline, the averaged related bias associated with *APOE-4* is −5.1% under the model without baseline versus 0.9% under the model with baseline, while as the percentage of baseline demented subjects increases to 34%, the bias changes to −14.1% versus −4.3% between the models.

The Nun Study data, a longitudinal study of aging and Alzheimer’s disease funded by the National Institute on Aging will be used to illustrate our proposed model. The dataset consists of a cohort of 678 members of the school sisters of Notre Dame religious congregation (Snowdon et al., 1997). Each participant agrees to allow investigators complete access to their convent archives, participate in near-annual assessments of cognitive and physical function and donate their brain at death. 177 participants are excluded from the analysis because of missing covariates or consent withdrawal. One conspicuous feature of the dataset is that over 10 percent (77) of the subjects are diagnosed with dementia at the baseline visit. Instead of removing those subjects from data analysis, the proposed multi-state Markov model helps to accommodate this baseline information into the likelihood by assuming shared random effects.

The first 10 waves of exam results since 1991 are used in this analysis. The transitions are summarized in Table 4. The covariates of interest are: (1) age in years centered at the median of 88 years, and (2) presence of apolipoprotein E-ε4 allele, a well-known risk factor for dementia.

The presence of a shared random effect assures the vector of serial observations on a given subject is correlated. When this is restricted to the transitional likelihood (the likelihood without a baseline), only the second through the last observation in this vector are dependent since the transitional likelihood is conditioned on the first observation. On the other hand the model with the baseline correlates all observations in the vector. The simulation studies in Table 2 show that excluding the baseline likelihood from this shared random effect produces estimates of the beta coefficients that are negatively biased and such bias increases as the variance of the shared random effect increases. Table 5 shows how the corresponding parameter estimates for the transition probabilities are consistently underestimated using real data from the Nun Study. In this application we compared three models: the naïve model where we ignore the random effects, an initial random intercept model without the baseline likelihood, and the shared random intercept model with the baseline likelihood.

Parameter estimates for transitions probabilities in the Nun Study data (base state: 1=Intact Cognition)

The results from Table 5 indicate that estimates under models with and without the baseline likelihood component are different. For example, consider the fixed effect of apolipoprotein E-ε4 allele (APOE4=1). It has been well documented that the presence of APOE4 increases the chance of cognitive impairment. Under the naïve model, the MLEs for APOE4=1 are (0.548, 0.863, 1.199, 0.651), which means that keeping other covariates constant, the odds ratio of having APOE4 present for transitions from intact cognition to MCI is 1.73, intact to global impairment is 2.37, intact to dementia is 3.32 and intact to death is 1.92. In comparison, the odds ratios are 3.14, 4.58, 5.18 and 3.71 under the model without the baseline likelihood component and 5.22, 7.65, 8.58, 6.20 under the model with the baseline component. We can see that although the effect of APOE4 is significant in all three models, the magnitude of odds ratios under the new model is larger.

The model introduced in this paper can be applied to higher order chains. Without loss of generality, in a second order case, the transition probability matrix ** P_{rsv}** has a hierarchical structure

The likelihood for the second order Markov model further breaks down the joint distribution into three likelihood components with shared random effects. The joint distribution can be factorized as the following

$$\begin{array}{l}f({y}_{i0},{y}_{i1},{y}_{i2},\dots ,{y}_{i{n}_{i}}\mid \mathit{\gamma})\\ =\underset{2\text{nd}-\text{order}-\text{Follow}-\text{up}}{\underbrace{f({y}_{i2}\mid {y}_{i1},{y}_{i0},\mathit{\gamma})f({y}_{i3}\mid {y}_{i2},{y}_{i1},\mathit{\gamma})\dots f({y}_{i{n}_{i}}\mid {y}_{i{n}_{i}-1},{y}_{i{n}_{i}-2},\mathit{\gamma})}}\underset{1\text{st}-\text{order}-\text{follow}-\text{up}}{\underbrace{f({y}_{i1}\mid {y}_{i0},\mathit{\gamma})}}\underset{\text{Baseline}}{\underbrace{f({y}_{i0}\mid \mathit{\gamma})}}\end{array}$$

and the likelihood for a particular subject is

$$L(\mathit{\Theta},\mathit{K},\mathit{\Phi}\mid \mathit{X})=\underset{2\text{nd}-\text{order}-\text{follow}-\text{up}}{\underbrace{\int \prod _{k=2}^{{n}_{i}}{P}_{{y}_{ik-2,}{y}_{ik-1,}{y}_{ik}}(\mathit{\Theta}\mid \mathit{X},\mathit{\gamma})}}\underset{1\text{st}-\text{order}-\text{follow}-\text{up}}{\underbrace{{P}_{{y}_{i1},{y}_{i0}}(\mathit{K}\mid \mathit{X},\mathit{\gamma})}}\underset{\text{baseline}}{\underbrace{{\pi}_{{y}_{i0}}(\mathit{\Phi}\mid {\mathit{X}}_{B},\mathit{\gamma})}}\underset{\text{random}\phantom{\rule{0.16667em}{0ex}}\text{effect}}{\underbrace{h(\mathit{\gamma})d\mathit{\gamma}}}$$

The variance covariance structure of the random effect distribution does not have to be the same across the likelihood components. It is also possible for these components to partially share the random effects. Take the random intercept and random slope model for instance: the baseline component shares only the random intercept with the two follow-up components, which share an extra random slope.

In theory this model applies to an arbitrarily higher order Markov chain, while in practice the number of parameters that need to be estimated can add up quickly which, in combination with the numerical integration of random effects, might produce the computational burdens in the likelihood optimization (refer to Table 6 for example).

The advantage of this approach is that; first, it helps to reduce the possible confounding bias that could occur otherwise. This is especially the case when higher order chains are assumed. Referring to the Nun Study data, the flow diagram indicates that there are 77 nuns demented at baseline and after the first follow-up wave, 39 more are demented. In a second order chain scenario, a total of 116 subjects would have to be removed from the analysis, accounting for over 20 percent of the total available data.

Second, the likelihood based approaches facilitate the inferential procedures like common likelihood ratio tests. These can be used to test model fitness, in particular, the hypothesis about the orders of a particular chain. We fit a second and third order Markov models using the Nun Study data. As shown in Table 6, all the fit statistics suggest that the third order model is no better than the second order model, and the likelihood ratio tests as well as the Akaike’s information criterion (AIC) indicate that the second order model might be better than the first order model, while the Bayesian information criterion (BIC) supports the first order model. This mixed result points out that the approximation of the joint distribution *f*(*y _{i}*

Moreover, by modeling the transition among states with a Markov structure, the one-step transition probability matrices constructed based on the parameter estimates provide additional information with respect to the mean time to absorption as well as the odds of absorption in competing absorbing states. This is particularly useful in the study of chronic diseases like Alzheimer’s disease where researchers are interested in the probability of disease onset before dying given a set of risk factors such as age, education, and genetic status.

Subjects in the Nun Study do not share the same baseline state. Among 501 subjects used for analysis, 128 of them had intact cognition at baseline, 249 had MCI, 47 had global impairment and 77 were demented. Although a model without a baseline likelihood component could be considered for a cohort with some homogeneous baseline states^{1}, the diversified states at baseline for subjects in the Nun Study make it important to incorporate the baseline outcomes into the likelihood construction. The proposed multi-state Markov model helps to accommodate this baseline information into the likelihood by assuming shared random effects. Since all the risk factors considered in the application to the Nun Study are the most established risk factors, the comparison of the maximum likelihood estimates shows difference only in magnitude rather than significance, it is feasible however that other risk factors with important but weaker association with cognitive status transitions might be missed without accounting for the baseline information. Moreover, from an epidemiological perspective, without including all baseline information, the prevalence (baseline) cases tend to be different from the incidence (follow-up) cases, which is likely to produce selection bias into the analysis.

The analysis of panel data with categorical outcomes is not a straightforward task. One of the strengths of the approach suggested in this manuscript for constructing a likelihood function for such data is that it assumes a Markov model for transitions among states. This makes it easy to incorporate the baseline status of the individual into the likelihood computation provided we introduce a shared random effect to assure the elements in the entire vector of observations on an individual are correlated. The arithmetic is no more complicated than when the baseline is ignored since we continue to rely on standard statistical software (Procedure NLMIXED in SAS) to fit the expanded likelihood to the data as evidenced by the Nun study data. However, there are some limitations to this approach since this software relies on numerical quadrature techniques. One limitation concerns *k*, the number of states in the process, or *r*, the number of risk factors investigated. As either *k* or *r* increase the number of unknown parameters increases making it difficult to achieve convergence of the likelihood to its maximum. One reason for this is that as *k* increases the possibility of encountering sparse cells in the one step transition matrix increases. The one step transition matrix links the covariate to the transition using a polytomous logistic model and convergence problems arise since that model is sensitive to sparse cells. Similarly as *r* increases the one step matrix gets partitioned according to different combinations of the risk factors again promoting the possibility of encountering sparse cells.

The likelihood construction in this model is based on the first-order Markov assumption, namely that the conditional distribution of the current outcome for a particular subject depends on the previous outcomes only through the most recent one. Whether the data maintains this Markov property directly affects the validity of maximum likelihood estimation. The verification of this assumption is non-ignorable. As in the general GLMM model, conditioning on both measured and unobserved latent variables makes the subject-specific coefficient difficult to interpret, especially when the involved covariates do not vary within individuals. According to Heagerty (1999), these coefficients measure the contrast in covariates when the random effects are held equal, but the random effects are not directly observed. The latent variable assumptions determine what values of random effect are equivalent; the magnitude and interpretation of the fixed effects therefore depend entirely on these assumptions. As a result the model tends to produce biases in regression estimates when the distribution of random effects has been misspecified (Litiere et al., 2007). This raises a new set of issues involving methods of model diagnostics under GLMMs, in particular the analysis of random effects misspecifications, which have not yet been thoroughly explored in the literature.

The work of Richard J. Kryscio was partially supported by a grant from the NIA (AG05144) and by a University of Kentucky Research Professorship.

^{1}To show that model 1 might be used for longitudinal data with homogeneous baseline outcomes, in addition to the simulation result as presented in Table 1, we analyzed the Nun Study data by including only subjects with the same baseline state (MCI). Almost identical MLEs are produced under model 1 and model 2.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Agresti A. Categorical data analysis. 2. John Wiley & Sons Inc.; Hoboken: 2002.
- Crouchley R. A random-effects model for ordered categorical data. J Am Stat Assoc. 1995;90:489–498.
- Crouchley R, Davies RB. A comparison of population-average and random effects models for the analysis of longitudinal count data with baseline information. Journal of the Royal Statistical Society, Series A. 1999;162:331–347.
- Gibbons RD, Hedeker D. Application of random-effects probit regression models. Journal of Consulting and Clinical Psychology. 1994;62:285–296. [PubMed]
- Gao S. A shared random effect parameter approach for longitudinal dementia data with non-ignorable missing data. Statistics in Medicine. 2004;23:211–219. [PMC free article] [PubMed]
- Heagerty PK. Marginally specified logistic-normal models for longitudinal binary data. Biometrics. 1999;55:688–698. [PubMed]
- Litiere S, Alonso A, Molenberghs G. Type I and Type II error under random-effects misspecification in generalized linear mixed model. Biometrics. 2007;63:1038–1044. [PubMed]
- Pulkstenis EP, Ten Have TR, Landis JR. Model for the analysis of binary longitudinal data subject to informative drop-out through remedication. J Am Stat Assoc. 1998;93:438–450.
- Rubin DB. Inference and missing data. Biometrika. 1976;63:581–590.
- Salazar JC, Schmitt FA, Yu L, Mendiondo MM, Kryscio RJ. Shared random effects analysis of multi-state Markov models: application to a longitudinal study of transitions to dementia. Statistics in Medicine. 2007;26:568–580. [PubMed]
- Salazar JC. Multi-state Markov models for longitudinal data. PhD dissertation. Department of Statistics, University of Kentucky; 2004.
- SAS Institute Inc. SAS/STAT User’s Guide, Version 8. SAS Institute Inc.; Cary, NC: 1999.
- Skrondal A, Rabe-Hesketh S. Generalized latent variable modeling: multilevel, longitudinal, and structural equation models. Chapman & Hall/CRC; Boca Raton: 2004.
- Snowdon David A, Greiner Lydia H, Mortimer James A, Riley Kathryn P, Greiner Philip A, Markesbery William R. Brain Infarction and the Clinical Expression of Alzheimer Disease: The Nun Study. JAMA. 1997;277:813–817. [PubMed]
- Stiratelli R, Laird NM, Ware JH. Random-effects models for serial observations with binary response. Biometrics. 1984;40:961–971. [PubMed]
- Ten Have TR, Kunselman A, Pulkstenis EP, Landis JR. Mixed effects logistic regression models for longitudinal binary response data with informative drop-out. Biometrics. 1998;54:367–383. [PubMed]
- Ten Have TR, Miller ME, Reboussin BA, James MK. Mixed effects logistic regression models for longitudinal ordinal functional response data with multiple-cause drop-out from the Longitudinal Study of Aging. Biometrics. 2000;56:279–287. [PubMed]
- Ten Have TR, Reboussin BA, Miller ME, Kunselman A. Mixed effects logistic regression models for multiple longitudinal binary functional limitation response with informative drop-out and confounding by baseline outcomes. Biometrics. 2002;58:137–144. [PubMed]
- Tyas SL, Salazar JC, Snowdon DA, Desrosiers MF, Riley KP, Mendiondo MS, Kryscio RJ. Transitions to mild cognitive impairments, dementia, and death: findings from the Nun Study. Am J Epidemiol. 2007;165:1231–1238. [PMC free article] [PubMed]
- Zeger SL, Karim MR. Generalized linear models with random effects: a Gibbs sampling approach. J Am Stat Assoc. 1991;86:79–86.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |