## Abstract

Random effects models are commonly used to analyze longitudinal categorical data. Marginalized random effects models are a class of models that permit direct estimation of marginal mean parameters and characterize serial correlation for longitudinal categorical data via random effects (Heagerty, 1999). Marginally specified logistic-normal models for longitudinal binary data. *Biometrics***55**, 688–698; Lee and Daniels, 2008. Marginalized models for longitudinal ordinal data with application to quality of life studies. *Statistics in Medicine***27**, 4359–4380). In this paper, we propose a Kronecker product (KP) covariance structure to capture the correlation between processes at a given time *and* the correlation within a process over time (serial correlation) for bivariate longitudinal ordinal data. For the latter, we consider a more general class of models than standard (first-order) autoregressive correlation models, by re-parameterizing the correlation matrix using partial autocorrelations (Daniels and Pourahmadi, 2009). Modeling covariance matrices via partial autocorrelations. *Journal of Multivariate Analysis***100**, 2352–2363). We assess the reasonableness of the KP structure with a score test. A maximum marginal likelihood estimation method is proposed utilizing a quasi-Newton algorithm with quasi-Monte Carlo integration of the random effects. We examine the effects of demographic factors on metabolic syndrome and C-reactive protein using the proposed models.

## Introduction

Data collected over time from the same subjects are called longitudinal data. Within-subject measurements (over time) typically exhibit serial correlation that must be taken into account for proper inference on covariate effects (Fitzmaurice and Laird, 1993).

Marginalized likelihood-based models are a class of models that have been developed for the analysis of longitudinal categorical data (Heagerty, 1999, 2002; Lee and Daniels, 2007, 2008; Schildcrout and Heagerty, 2007; Lee *and others*, 2009, 2011; Lee and Mercante, 2010). In these models, the correlation structure is modeled via random effects or a Markov correlation structure (or both), while the population averaged response is directly modeled as a function of covariates, which induces restrictions on the correlation model. In this paper, we focus on marginalized models with random effects.

Marginalized models using the random effects are called marginalized random effects models (MREMs). In MREMs, serial correlation is often modeled using a first-order autoregressive correlation structure (AR(1)). However, this assumption is restrictive. Here, we consider a more flexible class of models for serial autocorrelation by re-parameterizing the correlation matrix using partial autocorrelations (Daniels and Pourahmadi, 2009), which has not been used to address serial correlation in these marginalized models.

Although the literature on multivariate longitudinal ordinal data are somewhat limited, a number of likelihood-based models have been proposed. Random effects models were proposed for joint modeling of clustered ordinal responses (Gueorguieva, 2001) and longitudinal bivariate ordinal responses (Todem *and others*, 2007). Similar models were proposed for item response theory models for multiple ordinal outcomes using random effects in Liu and Hedeker (2006), and item response models were also extended to accommodate longitudinal ordinal data under a missing at random (MAR) assumption (Liu, 2008). However, the above approaches were developed for generalized linear mixed models in which conditional (or subject-specific) covariate effects are directly modeled and little work has been done on likelihood-based marginal models such as marginalized models, which directly model marginal (or population-averaged) covariate effects. In this paper, we extend MREMs (Heagerty, 1999; Lee and Daniels, 2008) to accommodate multivariate longitudinal ordinal data using random effects to explain the serial correlation within responses over time and the correlation between responses at a given time. In particular, the *multivariate* longitudinal component of the covariance structure is modeled using a Kronecker product structure; to assess the feasibility of this assumption, we propose a score test and examine its type I error and power via simulations. Related models (but with no testing framework and no flexible modeling of the serial correlation matrix) were proposed for multivariate longitudinal binary data in Lee *and others* (2009) and in Ilk and Daniels (2007); the latter uses common random effects and Markov dependence structure.

We use the model proposed here to make inferences on data from the Korean Genomic Epidemiology Study (KoGES). The KoGES started in 2001 with a cohort of middle-aged Korean adults aged 39–69 years. Participants were examined every 2 years for up to 8 years to monitor the development of metabolic syndrome. Metabolic syndrome is a condition associated with the risk for developing type 2 diabetes, coronary heart disease, and other diseases related to plaque buildup in artery walls (e.g. stroke and peripheral vascular disease) (Expert Panel on Detection, Evaluation, and Treatment of High Blood Cholesterol in Adults, 2001). Metabolic syndrome is defined as three or more of the following five disorders: abdominal obesity (waist circumference >90 cm in men or >80 cm in women), high blood pressure (systolic BP levels >130 mmHg or diastolic BP levels >85 mmHg), high impaired fasting glucose (>110 mg/dl), hypertriglyceridemia (>150 mg/dl), and low high-density lipoprotein cholesterol (<40 mg/dl in men or <50 mg/dl in women). The primary outcome measures were the five metabolic syndrome indicators. All of these attributes are associated with increased levels of C-reactive protein (CRP) (Ridker *and others*, 2003). CRP is a protein produced in the liver during episodes of acute inflammation or infection. High levels of CRP serves as a general indicator of acute inflammation. In practice, CRP levels >2.5 mg/l are regarded as a problem. Also, recent research has suggested that chronic low-grade inflammation may contribute to the pathogenesis of metabolic syndrome (Ye *and others*, 2007; Devaraj *and others*, 2009). This motivates our interest in modeling these jointly. In addition, it is of primary interest how demographic factors affect metabolic syndrome and CRP at the population level, which motivates our development of marginalized models in this setting (bivariate, ordinal, longitudinal data).

The paper is organized as follows. In Section 2, we propose an MREM for multivariate longitudinal ordinal data. In Section 3, we conduct a simulation study to examine bias and coverage of estimators of the marginal mean parameters in the setting of ignorable missingness as well as the operating characteristics of the proposed score test. In Section 4, we analyze data from the longitudinal study on metabolic syndrome. Finally, conclusions and extensions are provided in Section 5.

## Model

We consider two longitudinal ordinal responses. The extension to more than two responses is (in principle) straightforward. Denote the vector of longitudinal *K*-category ordinal responses for the *i*th subject by $$\textbf {y}_i=(\textbf {y}_{i1}^{\mathrm {T}},\ldots ,\textbf {y}_{in_i}^{\rm T})^{\rm T}$$ where $$\textbf {y}_{it}^{\mathrm {T}}=(y_{it1},y_{it2})$$ are two ordinal responses at time *t*. We assume that *y*_{itj}, *t*=1,…,*n*_{i}, *j*=1,2, are conditionally independent given **b**_{i⋅j}=(*b*_{i1j},…,*b*_{inij}), where **b**_{i⋅j} is the vector of random effects for the *j*th ordinal response for subject *i*. Also *y*_{it1} and *y*_{it2} are conditionally independent given the vector of all the random effects for subject *i*, **b**_{i}=(*b*_{i11},*b*_{i12},…,*b*_{ini1},*b*_{ini2}), and the response vectors on different subjects are independent. Let **x**_{it} denote the vector of covariates for **y**_{it⋅}. We specify the marginal mean models as

*β*

_{j0k}:

*j*=1,2 are intercepts that satisfy the monotonicity property (strictly increasing in

*k*for a fixed

*j*) and

*β*_{j}:

*j*=1,2 are (

*p*×1)-dimensional unknown marginal mean parameters. Models (2.1) and (2.2) are cumulative logit models for ordinal data (McCullagh, 1980).

For the correlation models, we specify

_{it1}and Δ

_{it2}are intercepts (more in Section 2.2),

**0**is a 2

*n*

_{i}×1 vector of zeros, and

*Ω*_{i}is a 2

*n*

_{i}×2

*n*

_{i}matrix. Although correlations of within-subject measurements in

*Ω*_{i}may not be of primary interest, they must be taken into account to make proper inferences. Models (2.3) and (2.4) are models for association that allow the marginal models in (2.1) and (2.2) to be preserved. In (2.3) and (2.4), we need to ensure the monotonicity of

**Δ**(in

*k*) to guarantee the validity of the models.

For the multivariate marginalized models given by (2.1)–(2.4), if *β*_{j01}<⋯<*β*_{j0K−1}, then Δ_{itj1}<⋯<Δ_{itjK−1} for *j*=1,2.

The proof is similar to that of Theorem 1 in Lee and Daniels (2008).

### Covariance structure for random effects

In multivariate longitudinal data analysis, two types of correlations among responses must be considered: correlations between responses at the same time and correlations within and between responses at different times. To parsimoniously account for this correlation in the covariance matrix, we use a Kronecker product (KP) structure. The covariance matrix (*Ω*_{i}) can be written as the KP of the covariance matrix accounting for the serial correlation within a response over time (**R**_{i}) and covariance between responses at the same time (**Σ**), that is, *Ω*_{i}=**R**_{i}⊗**Σ**, where ⊗ denotes the KP (Galecki, 1994; Chaganty and Naik, 2002). We consider the following forms for **R**_{i} and **Σ**:

**R**

_{i}and

**Σ**are not unique; for

*c*≠0,

*c*

**R**

_{i}⊗(1/

*c*)

**Σ**=

**R**

_{i}⊗

**Σ**. However, this non-identifiability can be fixed by specifying

**R**

_{i}as a correlation matrix (Galecki, 1994) as we have here. The parameters in

*Ω*_{i}provide a measure of random variation of bivariate responses at the same time ($$\sigma _{1}^2$$, $$\sigma _{2}^2$$, and

*σ*

_{12}) and over time (

*σ*

_{12}

*ρ*

_{i,t,t′}, $$\sigma _1^2\rho _{i,t,t'}$$, and $$\sigma _2^2\rho _{i,t,t'}$$).

#### Flexible models for the serial correlation matrix, *R*_{i}

The positive-definiteness of the matrix **R**_{i} is a major obstacle for constructing parsimonious models. Most of the approaches to put structure on a covariance matrix (e.g. using the Modified Choleski Decomposition, Pourahmadi, 1999) do not extend in a convenient manner to a correlation matrix. As a result, a relatively simple structure for **R**_{i} such as AR(1) is often assumed (Heagerty, 1999; Lee and Daniels, 2008). However, it is often too strong an assumption. To allow more flexible structures, we re-parameterize the correlation matrix **R**_{i} in terms of the lag-1 correlations *ρ*_{i,j,j+1} and the partial autocorrelations (Joe, 2006, Daniels and Pourahmadi, 2009), *π*_{i,j,k}=*cor*(*Y*_{ij},*Y*_{ik}|*Y*_{il},*j*<*l*<*k*), the partial correlations between *Y*_{ij} and *Y*_{ik} adjusted for the intervening variables, {*j*+1,…,*k*−1}. These can be organized in a partial autocorrelation matrix *Π*_{i} with elements *π*_{i,jk} in the *jk*th position and 1’s on the main diagonal, where *π*_{i,j,j+1}=*ρ*_{i,j,j+1} for *j*=1,…,*n*_{i}−1. We note that the matrix *Π*_{i} is not required to be positive-definite and its off-diagonal entries are free to independently vary in the interval (−1,1).

The partial correlations are computed using the expression

**R**

_{2}(

*j*,

*k*) is the correlation matrix corresponding to the components (

*j*+1,…,

*j*+

*k*−1).

To model the elements of *Π*_{i}, we transform the off-diagonal elements of the matrix to the entire real line using Fisher’s *Z*-transformation as the link function, which is given by

**is a (**

*γ**q*×1)-dimensional unknown parameter and

**w**

_{i,jk}is a unit-specific covariate vector. The choice of

**w**

_{i,jk}allows the correlation structure to vary by unit-specific covariates and also allows for flexible parsimonious structure for the correlation matrix

**R**

_{i}. For example, we might specify

**w**

_{i,jk}=(

*I*(|

*k*−

*j*|=1),

*I*(|

*k*−

*j*|=1)*age) corresponding to a banded Toeplitz structure of the autocorrelation matrix

*Π*_{i}. The corresponding correlation matrix,

**R**

_{i}has an AR(1) structure with a lag-1 correlation parameter that varies (linearly) with the covariate, age. More general structures are given in Sections 3 and 4. A discussion of priors for

*γ*in Bayesian partial autocorrelation models can be found in Wang and Daniels (2013).

#### Score test for the covariance structure

The KP covariance structure makes an implicit assumption that, for both response variables, the longitudinal correlation structure is the same and that the covariance between both response variables does not depend on time and remains constant for all time points. To check the reasonableness of this assumption jointly with the structure of the serial correlation matrix (**R**_{i}), we propose a score test (Fahrmeir and Tutz, 2000) for the following hypothesis:

*Ω*_{i}is a function of

**and**

*σ***in (2.6) and (2.8), but we suppress this dependence in what follows (in addition to the subscript**

*γ**i*). Let

**be the vector of lower triangular elements of**

*ω***. Suppose that $$(\hat {\boldsymbol {\beta }}(\hat {\boldsymbol {\omega }}_0), {\hat {\boldsymbol {\omega }}}_0)$$ are the maximum likelihood estimates of**

*Ω***and**

*β***under the null hypothesis**

*ω**H*

_{0}. Then

*H*

_{0}, where $$\textbf {A}({\hat {\boldsymbol {\omega }}_0})$$ is the submatrix of the inverse of the observed information matrix for

**=(**

*θ***,**

*β***) corresponding to**

*ω***,**

*ω**d*is the number of constraints imposed by the null hypothesis, and

**R**

_{i}obviously impacts the number of constraints

*d*. For further details and an assessment of the power of each departure from

*H*

_{0}, see Section 3.2.

### Connection between marginal and conditional models

The marginal and conditional probabilities in (2.1)–(2.4) are related as follows:

*j*=1,2, and

*f*(⋅) is a univariate normal distribution with mean 0 and variance Var(

*b*

_{itj}). Therefore, the intercepts Δ

_{itjk}are a function of

*β*_{j}and Var(

*b*

_{itj}). Note that the integral in (2.9) is one-dimensional. Details for computing Δ

_{itjk}are given in supplementary material available at Biostatistics online.

### Maximum likelihood estimation

Let ** θ**=(

**,**

*β***,**

*σ***). The likelihood function is the integral over random effects of a product of multinominals,**

*γ**j*=1,2 and

*ϕ*(⋅) is a multivariate normal density with mean vector 0 and variance–covariance matrix

*Ω*_{i}. The marginalized likelihood in (2.10) is not available in closed form. We use a quasi-Monte Carlo (QMC) method (Niederreiter, 1992; Caflisch, 1998; Judd, 1998) to integrate out the random effects. QMC uses uniformly distributed deterministic sequences (Niederreiter, 1992). The QMC method has been shown to perform better than other numerical integration methods such as Gauss–Hermite and Monte Carlo methods for high-dimensional and highly correlated data in related settings (Gonzalez

*and others*, 2006).

Maximizing the log-likelihood with respect to *θ* yields the likelihood equations

*K*−1)+2

*p*+

*q*+3)-dimensional likelihood equations are given in supplementary material available at Biostatistics online.

The matrix of second derivatives of the observed data log-likelihood has a complex form. To avoid computing them directly, we use a quasi-Newton method to solve the likelihood equations,

**I**

_{e}(

**) is an empirical and consistent estimator of the information matrix at iteration**

*θ**g*,

## Simulations

### Examination of bias and coverage probabilities

We conducted several simulation studies to examine the bias and coverage for estimators of the marginal mean parameters in the setting of mis-specification of the correlation model under ignorable (MAR) missingness.

We simulated bivariate longitudinal ordinal data under the marginalized model proposed in Section 2 with two covariates, time and group (two levels). The marginal cumulative probabilities were specified as

*j*=1,2 where

*β*_{1}=(

*β*

_{101},

*β*

_{102},

*β*

_{103},

*β*

_{11},

*β*

_{12},

*β*

_{13})=(−0.8,0.5,1.6,−0.4,−0.1,0.7),

*β*_{2}=(

*β*

_{201},

*β*

_{202},

*β*

_{203},

*β*

_{21},

*β*

_{22},

*β*

_{23})=(−0.6,0.8,1.5,1.0,−0.2,0.9),

*t*=1,…,6, time

_{it}=(

*t*−1)/10, and group

_{i}equals 0 or 1 with an approximately equal sample size per group. The conditional probabilities in (2.3) and (2.4) were specified with $$\boldsymbol {\sigma }=(\sigma _{1}^2,\sigma _{21},\sigma _{2}^2)=(1.45,1.04,1.01)$$,

**=(0.5,0.1), and $$\textbf {w}_{i,jk}^{\mathrm {T}}=(I\{|j-k|=1\},I\{|j-k|=2\})$$.**

*γ*We generated 200 simulated datasets, each with a sample size of 300 (for details, see supplementary material available at Biostatistics online). We fit four models. Model 1 is the true model; Model 2 is a special case of Model 1 with $$\boldsymbol {\Omega }=\textbf {R}\otimes \hbox {diag}(\sigma _1^2,\sigma _2^2)$$; Model 3 is a special case of Model 1 with **b**_{it}=**b**_{i0}∼*N*(0,**Σ**) (ignoring the temporal correlation); Model 4 is a special case of Model 1 with first-order structure in **R** having $$\textbf {w}_{i,jk}^{\mathrm {T}}=(I\{|j-k|=1\})$$. Note that Models 2–4 all mis-specify the correlation structure in different ways.

For missingness, we specified the following MAR dropout model:

Table 1 presents the point estimates, absolute values of biases, and the coverage probabilities (of 95% confidence intervals) for the marginal mean parameters. We saw considerable bias in the Models 2 and 3, with the percentage relative bias as large as 13.75% (*β*_{11} for Model 2), 30.75% (*β*_{11} for Model 3), 11.40% (*β*_{21} for Model 3), and 18.50% (*β*_{22} for Model 2) for the coefficients of time and group, respectively. In addition, we point out the biases for time (8.00% for *β*_{21}) and the time by group interaction (6.67% for *β*_{23}) parameters in Model 4; this is not surprising since the mis-specification in Model 4 is in the longitudinal correlation matrix *R*. The coverage probabilities for all coefficients in Models 2 and 3 were similar and several were much less than the desired 95%.

Model 1 | Model 2 | Model 3 | Model 4 | ||||||
---|---|---|---|---|---|---|---|---|---|

Para. | Truth | Mean | AB | Mean | AB | Mean | AB | Mean | AB |

β_{101} | −0.8 | −0.80 | 0.00 | −0.83 | 0.03 | −0.80 | 0.00 | −0.80 | 0.00 |

(0.954) | (0.921) | (0.933) | (0.941) | ||||||

β_{102} | 0.5 | 0.50 | 0.00 | 0.48 | 0.02 | 0.50 | 0.00 | 0.50 | 0.00 |

(0.951) | (0.937) | (0.939) | (0.951) | ||||||

β_{103} | 1.6 | 1.60 | 0.00 | 1.60 | 0.00 | 1.59 | 0.01 | 1.60 | 0.00 |

(0.958) | (0.889) | (0.939) | (0.953) | ||||||

β_{11} | −0.4 | −0.40 | 0.00 | −0.35 | 0.05 | −0.52 | 0.12 | −0.41 | 0.01 |

(0.979) | (0.952) | (0.894) | (0.963) | ||||||

β_{12} | −0.1 | −0.108 | 0.008 | −0.080 | 0.020 | −0.090 | 0.010 | −0.091 | 0.009 |

(0.953) | (0.952) | (0.951) | (0.955) | ||||||

β_{13} | 0.7 | 0.74 | 0.04 | 0.73 | 0.03 | 0.78 | 0.08 | 0.75 | 0.05 |

(0.959) | (0.955) | (0.889) | (0.951) | ||||||

β_{201} | −0.6 | −0.59 | 0.01 | −0.61 | 0.01 | −0.60 | 0.00 | −0.62 | 0.02 |

(0.969) | (0.968) | (0.933) | (0.944) | ||||||

β_{202} | 0.8 | 0.80 | 0.00 | 0.80 | 0.00 | 0.80 | 0.00 | 0.80 | 0.00 |

(0.956) | (0.928) | (0.900) | (0.944) | ||||||

β_{203} | 1.5 | 1.50 | 0.00 | 1.48 | 0.02 | 1.50 | 0.00 | 1.50 | 0.00 |

(0.958) | (0.948) | (0.939) | (0.952) | ||||||

β_{21} | 1.0 | 0.99 | 0.01 | 1.03 | 0.03 | 0.89 | 0.11 | 1.08 | 0.08 |

(0.964) | (0.944) | (0.917) | (0.963) | ||||||

β_{22} | −0.2 | −0.22 | 0.02 | −0.16 | 0.04 | −0.21 | 0.01 | −0.18 | 0.02 |

(0.964) | (0.948) | (0.939) | (0.944) | ||||||

β_{23} | 0.9 | 0.93 | 0.03 | 0.84 | 0.06 | 0.98 | 0.08 | 0.86 | 0.04 |

(0.968) | (0.944) | (0.917) | (0.943) |

Model 1 | Model 2 | Model 3 | Model 4 | ||||||
---|---|---|---|---|---|---|---|---|---|

Para. | Truth | Mean | AB | Mean | AB | Mean | AB | Mean | AB |

β_{101} | −0.8 | −0.80 | 0.00 | −0.83 | 0.03 | −0.80 | 0.00 | −0.80 | 0.00 |

(0.954) | (0.921) | (0.933) | (0.941) | ||||||

β_{102} | 0.5 | 0.50 | 0.00 | 0.48 | 0.02 | 0.50 | 0.00 | 0.50 | 0.00 |

(0.951) | (0.937) | (0.939) | (0.951) | ||||||

β_{103} | 1.6 | 1.60 | 0.00 | 1.60 | 0.00 | 1.59 | 0.01 | 1.60 | 0.00 |

(0.958) | (0.889) | (0.939) | (0.953) | ||||||

β_{11} | −0.4 | −0.40 | 0.00 | −0.35 | 0.05 | −0.52 | 0.12 | −0.41 | 0.01 |

(0.979) | (0.952) | (0.894) | (0.963) | ||||||

β_{12} | −0.1 | −0.108 | 0.008 | −0.080 | 0.020 | −0.090 | 0.010 | −0.091 | 0.009 |

(0.953) | (0.952) | (0.951) | (0.955) | ||||||

β_{13} | 0.7 | 0.74 | 0.04 | 0.73 | 0.03 | 0.78 | 0.08 | 0.75 | 0.05 |

(0.959) | (0.955) | (0.889) | (0.951) | ||||||

β_{201} | −0.6 | −0.59 | 0.01 | −0.61 | 0.01 | −0.60 | 0.00 | −0.62 | 0.02 |

(0.969) | (0.968) | (0.933) | (0.944) | ||||||

β_{202} | 0.8 | 0.80 | 0.00 | 0.80 | 0.00 | 0.80 | 0.00 | 0.80 | 0.00 |

(0.956) | (0.928) | (0.900) | (0.944) | ||||||

β_{203} | 1.5 | 1.50 | 0.00 | 1.48 | 0.02 | 1.50 | 0.00 | 1.50 | 0.00 |

(0.958) | (0.948) | (0.939) | (0.952) | ||||||

β_{21} | 1.0 | 0.99 | 0.01 | 1.03 | 0.03 | 0.89 | 0.11 | 1.08 | 0.08 |

(0.964) | (0.944) | (0.917) | (0.963) | ||||||

β_{22} | −0.2 | −0.22 | 0.02 | −0.16 | 0.04 | −0.21 | 0.01 | −0.18 | 0.02 |

(0.964) | (0.948) | (0.939) | (0.944) | ||||||

β_{23} | 0.9 | 0.93 | 0.03 | 0.84 | 0.06 | 0.98 | 0.08 | 0.86 | 0.04 |

(0.968) | (0.944) | (0.917) | (0.943) |

Displayed are the average regression coefficient estimates and the absolute values of biases (AB). The coverage probabilities (CP) of 95% confidence intervals are in the parentheses.

Overall, the simulation results demonstrate the large biases that can occur in the marginal mean parameters when the correlation is mis-specified in the presence of MAR missingness.

### Examination of score test for the covariance structure

We also conducted a simulation study to examine the performance of the score test for the covariance structure of the random effects covariance matrix in Section 2.1.2. The type I error and the power were assessed for several cases.

We first simulated bivariate longitudinal ordinal data using the same marginal and conditional probabilities as in Section 3.1 except for the structure of the partial autocorrelation matrix and time dimension (here, *T*=3). We used first-order structure using ** γ**=(0.5) and

**w**

_{i,jk}=(

*I*{|

*j*−

*k*|=1}). We assessed Type I error rates for sample sizes of 1000 and 2000 based on 500 replicated datasets. The Type I error rates were 0.10 and 0.06, respectively. This indicates that the score test for KP structure may need a fairly large sample size for the asymptotic results to hold. However, the test avoids the need to fit a model with a more complex covariance structure to assess its reasonableness for the data.

We conducted further simulations to assess power. We simulated bivariate longitudinal ordinal data using the same marginal and conditional probabilities as in Section 3.1 and with a KP structure for the random effects covariance matrix, but now with the correlation matrix **R**_{i} having second-order structure using ** γ**=(0.5,0.4) and

**w**

_{i,jk}=(

*I*{|

*j*−

*k*|=1},

*I*{|

*j*−

*k*|=2}). Then we fit the marginalized model with first-order structure for

**R**

_{i}. The powers for sample sizes of 1000 and 2000 based on 200 replicated datasets were high, 0.98 and 0.99, respectively.

We also assessed the power of the score test when the KP structure was violated, but **R**_{i} was correctly specified. In particular, we assumed the following non-KP structure for ** Ω** in (2.5):

*α*

_{1}=0.7,

*α*

_{2}=0.4,

*α*

_{12}=0.2;

*τ*

_{1}=2.0,

*τ*

_{2}=1.0. Then we fit the marginalized model assuming the KP structure (with first-order structure for

**R**

_{i}). The powers for sample sizes of 1000 and 2000 based on 200 replicated datasets were 0.54 and 0.90, respectively.

Based on these limited simulations, it appears that the test may have more power to detect departures of the structure of the correlation matrix within a KP structure than departures from the KP structure itself. In addition, a fairly large sample size is needed to obtain the correct type I error.

## Analysis of metabolic syndrome data

### Description of data

We use the models here to make an inference on the population-averaged effects of demographic factors on metabolic syndrome and the level of CRP as well as examine the association between these two outcomes. We consider only the last three waves of this study as CRP was only collected on a handful of subjects in the first wave.

The two response variables, metabolic syndrome and CRP, were categorized on a three-point ordinal scale (Post *and others*, 2011). The metabolic syndrome variable takes the value 1 when the number of the five disorders present is zero or 1; the variable takes value 2 when two disorders are present; the variable takes value 3 when three or more disorders are present (indicating a diagnosis of metabolic syndrome) (Post *and others*, 2011). For CRP, the three values correspond to a level less than 1, between 1 and 3, and greater than or equal to 3 (Reilly *and others*, 2003). The proportions of metabolic syndrome and CRP over the three visits are summarized in Table 2.

Visit | ||||||
---|---|---|---|---|---|---|

1 | 2 | 3 | ||||

Meta | CRP | Meta | CRP | Meta | CRP | |

1 | 234 (0.124) | 788 (0.673) | 420 (0.250) | 1084 (0.646) | 282 (0.199) | 884 (0.624) |

2 | 1237 (0.656) | 220 (0.188) | 952 (0.567) | 308 (0.183) | 849 (0.600) | 303 (0.214) |

3 | 416 (0.220) | 163 (0.139) | 307 (0.183) | 287 (0.171) | 285 (0.201) | 229 (0.162) |

Total | 1887 | 1171 | 1679 | 1679 | 1416 | 1416 |

Visit | ||||||
---|---|---|---|---|---|---|

1 | 2 | 3 | ||||

Meta | CRP | Meta | CRP | Meta | CRP | |

1 | 234 (0.124) | 788 (0.673) | 420 (0.250) | 1084 (0.646) | 282 (0.199) | 884 (0.624) |

2 | 1237 (0.656) | 220 (0.188) | 952 (0.567) | 308 (0.183) | 849 (0.600) | 303 (0.214) |

3 | 416 (0.220) | 163 (0.139) | 307 (0.183) | 287 (0.171) | 285 (0.201) | 229 (0.162) |

Total | 1887 | 1171 | 1679 | 1679 | 1416 | 1416 |

Proportions are given in parentheses.

The demographic factors of interest included sex (1=male; 0=female), age, alcohol consumption (Drink1=1 if drinking in the past, 0 otherwise; Drink2=1 if drinking currently, 0 otherwise), and smoking status (Smoke1=1 if smoking in the past, 0 otherwise; Smoke2=1 if smoking currently, 0 otherwise).

There was missingness in the data. The total number of observations of CRP was smaller at the first visit than at the second and third visits. Although there were some subjects with one of the paired responses missing, we used all data available. In the analysis, we assume the missingness is ignorable.

### Computations

The quasi-Newton algorithm is not trivial computationally due to the need to obtain estimates and derivatives of Δ_{it} using Gauss–Hermite quadrature for all subjects and at all times within each quasi-Newton step. Each quasi-Newton step (in which all the Δ_{it} need to be computed) took about 1 min on a Pentium with a 1.6 GHz processor for the OMREM with a quasi-MC sample size of 10 000 and a 40-point Gauss–Hermite quadrature. Note that increasing the quasi-MC sample size to 20 000 resulted in an insignificant difference in the maximized log-likelihood values. By using good initial values obtained by fitting an independent cumulative logit model in standard software minimized the number of iterations until convergence. For example, in our analysis below, we obtained convergence in 40 iterations of the quasi-Newton step (with a convergence criterion of 10^{−5}).

### Models fit

We fit seven models under an assumption of ignorable missingness. Model 0 was two independent cumulative logit models (assuming all the responses were independent). The other models were based on the specification of *z*(*π*) in (2.8) using *w*_{i,jk} and Σ_{2} (Table 3). Model 3 has a random effects covariance matrix with first-order structure. Models 4, 5, and 6 allow the entire random effects covariance matrix with first-order structure to differ by sex, drinking type, and smoking type, respectively. Models 1 and 2 are special cases of Model 3 ignoring the correlation between two responses and the temporal correlation, respectively. We did not consider partial autocorrelation models with non-zero higher lag partial autocorrelations due to convergence problems with only three (time) waves.

Model | $$z(\pi)=\mathbf {w}_{i,jk}^{\mathrm {T}}\boldsymbol {\gamma }$$ in R_{i} | Σ |
---|---|---|

1 | γ_{1}I{|j−k|=1} | $$\hbox {diag}(\sigma _1^2,\sigma _2^2)$$ |

2 | γ=0 | Unstructured |

3 | γ_{1}I{|j−k|=1} | Unstructured |

4 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Sex | Unstructured |

5 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Drink1+γ_{3}I{|j−k|=1} Drink2 | Unstructured |

6 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Smoke1+γ_{3}I{|j−k|=1} Smoke2 | Unstructured |

Model | $$z(\pi)=\mathbf {w}_{i,jk}^{\mathrm {T}}\boldsymbol {\gamma }$$ in R_{i} | Σ |
---|---|---|

1 | γ_{1}I{|j−k|=1} | $$\hbox {diag}(\sigma _1^2,\sigma _2^2)$$ |

2 | γ=0 | Unstructured |

3 | γ_{1}I{|j−k|=1} | Unstructured |

4 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Sex | Unstructured |

5 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Drink1+γ_{3}I{|j−k|=1} Drink2 | Unstructured |

6 | γ_{1}I{|j−k|=1}+γ_{2}I{|j−k|=1} Smoke1+γ_{3}I{|j−k|=1} Smoke2 | Unstructured |

Table 4 gives maximum likelihood estimates, maximized log-likelihood values, and Akaike Information Criterion (AICs) for all seven models. Comparison of deviances for Models 0 versus 1, 1 versus 3, 2 versus 3, 3 versus 4, 3 versus 5, and 3 versus 6 which were nested, yielded △*D*_{01}=1447.61, *p*-value = 0.00 on 3 degrees of freedom, △*D*_{13}=111.274, *p*-value = 0.00 on 1 degrees of freedom, △*D*_{23}=7.808, *p*-value = 0.005 on 1 degrees of freedom, △*D*_{34}=7.100, *p*-value = 0.008 on 1 degrees of freedom, △*D*_{35}=7.122, *p*-value = 0.028 on 2 degrees of freedom, and △*D*_{36}=6.694, *p*-value = 0.035 on 2 degrees of freedom. These comparisons indicated that the Models 4–6 fit better than simpler models (Models 0–3). Comparison of AIC suggested Model 4 fit slightly better than Models 5 and 6. Thus, we focus on Model 4 for further analysis.

Model 0 | Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | Model 6 | |
---|---|---|---|---|---|---|---|

Response=Metabolic syndrome | |||||||

Int1 | 0.306 (0.201) | −0.285 (0.276) | −0.267 (0.267) | −0.263 (0.267) | −0.259 (0.267) | −0.260 (0.267) | −0.261 (0.267) |

Int2 | 3.214* (0.208) | 2.625* (0.284) | 2.642* (0.275) | 2.646* (0.275) | 2.651* (0.275) | 2.650* (0.275) | 2.650* (0.275) |

Sex | 0.459* (0.092) | 0.468* (0.115) | 0.468* (0.113) | 0.468* (0.113) | 0.468* (0.113) | 0.469* (0.113) | 0.469* (0.113) |

Age | −0.037* (0.004) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) |

Drink1 | 0.017 (0.133) | 0.072 (0.124) | 0.073 (0.122) | 0.074 (0.122) | 0.075 (0.122) | 0.071 (0.122) | 0.076 (0.122) |

Drink2 | 0.160* (0.067) | 0.127 (0.075) | 0.122 (0.074) | 0.122 (0.074) | 0.123 (0.074) | 0.121 (0.074) | 0.122 (0.074) |

Smoke1 | −0.243* (0.106) | −0.277* (0.117) | −0.282* (0.116) | −0.282* (0.116) | −0.282* (0.116) | −0.282* (0.116) | −0.283* (0.116) |

Smoke2 | −0.628 (0.107) | −0.562* (0.124) | −0.545* (0.123) | −0.546* (0.123) | −0.546* (0.123) | −0.546* (0.123) | −0.547* (0.123) |

Response=CRP | |||||||

Int1 | 2.351* (0.229) | 2.357* (0.286) | 2.147* (0.275) | 2.150* (0.275) | 2.153* (0.275) | 2.152* (0.275) | 2.152* (0.275) |

Int2 | 3.433* (0.233) | 3.435* (0.289) | 3.227* (0.279) | 3.231* (0.279) | 3.234* (0.279) | 3.233* (0.279) | 3.233* (0.279) |

Sex | −0.325* (0.100) | −0.337* (0.120) | −0.336* (0.117) | −0.337* (0.117) | −0.337* (0.117) | −0.337* (0.117) | −0.337* (0.117) |

Age | −0.032 (0.004) | −0.032* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028 (0.005) |

Drink1 | 0.289 (0.153) | 0.336* (0.149) | 0.350* (0.146) | 0.350* (0.146) | 0.351* (0.146) | 0.349* (0.147) | 0.351* (0.146) |

Drink2 | 0.121 (0.075) | 0.109* (0.085) | 0.110* (0.083) | 0.111* (0.083) | 0.111* (0.083) | 0.110* (0.083) | 0.111* (0.083) |

Smoke1 | 0.213 (0.116) | 0.186 (0.130) | 0.183 (0.127) | 0.184 (0.128) | 0.184 (0.128) | 0.184 (0.128) | 0.183 (0.128) |

Smoke2 | −0.292 (0.117) | −0.237 (0.132) | −0.223 (0.131) | −0.223 (0.131) | −0.223 (0.131) | −0.223 (0.131) | −0.224* (0.131) |

Σ^{1/2} | |||||||

s_{11} | 2.410* (0.092) | 2.359* (0.088) | 2.382* (0.090) | 2.379* (0.090) | 2.389* (0.091) | 2.386* (0.090) | |

s_{21} | 0.737* (0.070) | 0.740* (0.071) | 0.740* (0.071) | 0.741* (0.071) | 0.741* (0.071) | ||

s_{22} | 1.832* (0.082) | 1.634* (0.083) | 1.645* (0.086) | 1.643* (0.086) | 1.648* (0.086) | 1.647* (0.086) | |

R_{i} (γ) | |||||||

Int | 2.597* (0.132) | 2.966* (0.133) | 2.804* (0.186) | 2.625* (0.199) | 2.736* (0.180) | ||

Sex | 0.535* (0.256) | ||||||

Drink1 | 1.980* (0.456) | ||||||

Drink2 | 0.356 (0.269) | ||||||

Smoke1 | 1.116* (0.285) | ||||||

Smoke2 | 0.025 (0.329) | ||||||

Loglik. | −8349.404 | −7625.599 | −7573.866 | −7569.962 | −7566.412 | −7566.401 | −7566.615 |

AIC | 16762.807 | 15289.198 | 15185.732 | 15179.924 | 15174.822 | 15176.802 | 15177.230 |

Model 0 | Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | Model 6 | |
---|---|---|---|---|---|---|---|

Response=Metabolic syndrome | |||||||

Int1 | 0.306 (0.201) | −0.285 (0.276) | −0.267 (0.267) | −0.263 (0.267) | −0.259 (0.267) | −0.260 (0.267) | −0.261 (0.267) |

Int2 | 3.214* (0.208) | 2.625* (0.284) | 2.642* (0.275) | 2.646* (0.275) | 2.651* (0.275) | 2.650* (0.275) | 2.650* (0.275) |

Sex | 0.459* (0.092) | 0.468* (0.115) | 0.468* (0.113) | 0.468* (0.113) | 0.468* (0.113) | 0.469* (0.113) | 0.469* (0.113) |

Age | −0.037* (0.004) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) | −0.026* (0.005) |

Drink1 | 0.017 (0.133) | 0.072 (0.124) | 0.073 (0.122) | 0.074 (0.122) | 0.075 (0.122) | 0.071 (0.122) | 0.076 (0.122) |

Drink2 | 0.160* (0.067) | 0.127 (0.075) | 0.122 (0.074) | 0.122 (0.074) | 0.123 (0.074) | 0.121 (0.074) | 0.122 (0.074) |

Smoke1 | −0.243* (0.106) | −0.277* (0.117) | −0.282* (0.116) | −0.282* (0.116) | −0.282* (0.116) | −0.282* (0.116) | −0.283* (0.116) |

Smoke2 | −0.628 (0.107) | −0.562* (0.124) | −0.545* (0.123) | −0.546* (0.123) | −0.546* (0.123) | −0.546* (0.123) | −0.547* (0.123) |

Response=CRP | |||||||

Int1 | 2.351* (0.229) | 2.357* (0.286) | 2.147* (0.275) | 2.150* (0.275) | 2.153* (0.275) | 2.152* (0.275) | 2.152* (0.275) |

Int2 | 3.433* (0.233) | 3.435* (0.289) | 3.227* (0.279) | 3.231* (0.279) | 3.234* (0.279) | 3.233* (0.279) | 3.233* (0.279) |

Sex | −0.325* (0.100) | −0.337* (0.120) | −0.336* (0.117) | −0.337* (0.117) | −0.337* (0.117) | −0.337* (0.117) | −0.337* (0.117) |

Age | −0.032 (0.004) | −0.032* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028* (0.005) | −0.028 (0.005) |

Drink1 | 0.289 (0.153) | 0.336* (0.149) | 0.350* (0.146) | 0.350* (0.146) | 0.351* (0.146) | 0.349* (0.147) | 0.351* (0.146) |

Drink2 | 0.121 (0.075) | 0.109* (0.085) | 0.110* (0.083) | 0.111* (0.083) | 0.111* (0.083) | 0.110* (0.083) | 0.111* (0.083) |

Smoke1 | 0.213 (0.116) | 0.186 (0.130) | 0.183 (0.127) | 0.184 (0.128) | 0.184 (0.128) | 0.184 (0.128) | 0.183 (0.128) |

Smoke2 | −0.292 (0.117) | −0.237 (0.132) | −0.223 (0.131) | −0.223 (0.131) | −0.223 (0.131) | −0.223 (0.131) | −0.224* (0.131) |

Σ^{1/2} | |||||||

s_{11} | 2.410* (0.092) | 2.359* (0.088) | 2.382* (0.090) | 2.379* (0.090) | 2.389* (0.091) | 2.386* (0.090) | |

s_{21} | 0.737* (0.070) | 0.740* (0.071) | 0.740* (0.071) | 0.741* (0.071) | 0.741* (0.071) | ||

s_{22} | 1.832* (0.082) | 1.634* (0.083) | 1.645* (0.086) | 1.643* (0.086) | 1.648* (0.086) | 1.647* (0.086) | |

R_{i} (γ) | |||||||

Int | 2.597* (0.132) | 2.966* (0.133) | 2.804* (0.186) | 2.625* (0.199) | 2.736* (0.180) | ||

Sex | 0.535* (0.256) | ||||||

Drink1 | 1.980* (0.456) | ||||||

Drink2 | 0.356 (0.269) | ||||||

Smoke1 | 1.116* (0.285) | ||||||

Smoke2 | 0.025 (0.329) | ||||||

Loglik. | −8349.404 | −7625.599 | −7573.866 | −7569.962 | −7566.412 | −7566.401 | −7566.615 |

AIC | 16762.807 | 15289.198 | 15185.732 | 15179.924 | 15174.822 | 15176.802 | 15177.230 |

Standard errors are given in the parentheses. *s*_{11}, *s*_{21}, and *s*_{22} are elements of Σ^{1/2}.

*Significance with 95% confidence level.

The top part of Table 4 gives estimates of coefficients and associated standard errors for metabolic syndrome,

*β*

_{101},

*β*

_{102})=(−0.259,2.651). The effects of gender (SEX: $$\hat \beta =0.468$$, SE = 0.113), age (AGE: $$\hat \beta =-0.026$$, SE = 0.005), and past or current smoking (SMOKE1: $$\hat \beta =-0.282$$, SE = 0.116; SMOKE2: $$\hat \beta =-0.546$$, SE = 0.123) were significant. This indicates that the log odds of the sum of five metabolic syndrome indicators below any given level (instead of above it) was higher for males than for females, decreased as the participant’s age increases, and was smaller for smokers than for non-smokers. Also, there was a significant difference between the two smoking types’ marginal cumulative probabilities of metabolic syndrome ($$\hat \beta _{{\mathrm {smoke}}1}-\hat \beta _{{\rm smoke}2}=0.264$$, SE = 0.094).

The next part of Table 4 displays the estimates of coefficients and standard errors for CRP,

*β*

_{201},

*β*

_{202})=(2.153,3.234). The regression coefficients for participant’s sex ($$\hat \beta =-0.337$$, SE = 0.117), age ($$\hat \beta =-0.028$$, SE = 0.005), and past ($$\hat \beta =0.351$$, SE = 0.146) or current drinking ($$\hat \beta =0.111$$, SE = 0.083) were significant. This indicates that the log odds of CRP below any given level (instead of above it) was lower for males than for females, decreased as participant’s age increases, and was higher for drinkers than for nondrinkers. There was not a significant difference between smokers and nonsmokers. However, the difference of the regression coefficients for the two smoking types (past versus present) was significant ($$\hat \beta _{{\mathrm {smoke}}1}-\hat \beta _{{\rm smoke}2}=0.407$$, SE = 0.118).

We jointly tested the KP and serial correlation matrix structure using the score test proposed in Section 2.1.2. The score test statistic was 23.5 on 16 degrees of freedom. The corresponding *p*-value was 0.10. Therefore, we conclude that there was not strong evidence against the KP structure with a first-order serial correlation matrix that depends on sex.

The estimates of the serial correlation parameters *γ*_{1} and *γ*_{2} were significant ($$\hat \gamma _1=2.804$$, SE = 0.186; $$\hat \gamma _2=0.535$$, SE = 0.256) and implied lag-1 correlations of 0.993 and 0.997 for females and males, respectively. The estimated correlation of *b*_{it1} and *b*_{it2} was moderate (0.410) and statistically significant (covariance = 1.770, SE = 0.217, *p*-value <0.0001). This indicates that metabolic syndrome and CRP were positively correlated (as expected *a priori*).

### Comparisons between models

There were numerous differences between the simplest models (Models 0 and 1) and the other models, including the smaller estimated coefficients for SMOKE2 for metabolic syndrome in Models 0 and 1 versus the other models, and the smaller estimated slopes for Age in CRP in Models 0 and 1 versus the other models. In addition, there were large differences in the intercepts in Models 0 and 1 versus the other models. The standard errors of the estimated marginal mean parameters in Model 1 were larger than those in Models 2–6, which indicates the increased efficiency of Models 2–6 versus Model 1. The smaller standard errors for Model 0 were from the incorrect assumption of complete independence which (typically) results in too small standard errors and considerable bias.

## Conclusion

We have developed marginalized models for multivariate longitudinal ordinal outcomes. These directly model marginal probabilities as a function of covariates while accounting for the longitudinal correlation of responses and the correlation of responses at the same time via random effects. These two types of correlation are taken into account using a covariance matrix with a KP structure and serial correlation captured by flexible models using partial autocorrelations; the overall covariance structure can be checked using a score test.

Simulations emphasized the importance of modeling the correlation correctly in the presence of ignorable missingness and demonstrate the operating characteristics of the score test.

For the KoGES study, we found that metabolic syndrome is more prevalent in females, while higher levels of CRP are more common in males. Smoking in the past and the present increases the chance of metabolic syndrome, but does not correspond to higher levels of CRP.

The proposed models in this paper could be extended to a more general random effect covariance matrix, including allowing the correlation in **Σ** to be a function of covariates. We also plan on doing further exploration of the score test asymptotics and power. Finally, we are also working on further improving the computation time and developing an R package for these models.

## Supplementary material

Supplementary material is available online at http://biostatistics.oxfordjournals.org.

## Funding

This project was supported by NIH grant CA85295 and Basic Science Research Program through the National Research Foundation of Korea (KRF) funded by the Ministry of Education, Science and Technology (2012-004002).

## Acknowledgments

We would like to thank Dr Jungbok Lee of Korea University for providing the data and the latter for their help in data collection and clarifying some issues with data. *Conflict of Interest*: None declared.