## 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. Biometrics55, 688–698; Lee and Daniels, 2008. Marginalized models for longitudinal ordinal data with application to quality of life studies. Statistics in Medicine27, 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 Analysis100, 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 ith 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 yitj, t=1,…,ni, j=1,2, are conditionally independent given bij=(bi1j,…,binij), where bij is the vector of random effects for the jth ordinal response for subject i. Also yit1 and yit2 are conditionally independent given the vector of all the random effects for subject i, bi=(bi11,bi12,…,bini1,bini2), and the response vectors on different subjects are independent. Let xit denote the vector of covariates for yit. We specify the marginal mean models as

(2.1)
$$\hbox{logit}\{P(Y_{it1}\leq k\mid \textbf{x}_{it})\}=\beta_{10k}+\textbf{x}_{it}^{\mathrm{T}}\boldsymbol{\beta}_1,\label{eq2.1}$$

(2.2)
$$\hbox{logit}\{P(Y_{it2}\leq k\mid \textbf{x}_{it})\}=\beta_{20k}+\textbf{x}_{it}^{\rm T}\boldsymbol{\beta}_2,\label{eq2.2}$$
where β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

(2.3)
$$\hbox{logit}\{P(Y_{it1} \leq k\mid b_{it1},\textbf{x}_{it})\}=\Delta_{it1k}+b_{it1},\label{eq2.3}$$

(2.4)
$$\hbox{logit}\{P(Y_{it2} \leq k\mid b_{it2},\textbf{x}_{it})\}=\Delta_{it2k}+b_{it2},\label{eq2.4}$$

(2.5)
$$\textbf{b}_{i}\sim \hbox{i.i.d. } N(\textbf{0},\boldsymbol{\Omega}_i) \label{eq2.5},$$
where Δit1 and Δit2 are intercepts (more in Section 2.2), 0 is a 2ni×1 vector of zeros, and Ωi is a 2ni×2ni 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.

Theorem 1

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

Proof

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 (Ri) and covariance between responses at the same time (Σ), that is, Ωi=RiΣ, where ⊗ denotes the KP (Galecki, 1994; Chaganty and Naik, 2002). We consider the following forms for Ri and Σ:

(2.6)
$$\label{eq2.6} \textbf{R}_{i}=\left(\begin{matrix} 1 & \rho_{i12} & \cdots & \rho_{i,1,n_i} \\ \rho_{i12} & 1 & \cdots & \rho_{i,2,n_i} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{i,1,n_i} & \rho_{i,2,n_i} & \cdots & 1 \end{matrix}\right)\quad \hbox{and}\quad \boldsymbol{\Sigma}=\left(\begin{matrix} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_2^2 \end{matrix}\right).$$
A lack of identifiability can result in estimating a covariance matrix with a KP structure. The indeterminacy stems from the fact that Ri and Σ are not unique; for c≠0, cRi⊗(1/c)Σ=RiΣ. However, this non-identifiability can be fixed by specifying Ri 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, Ri

The positive-definiteness of the matrix Ri 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 Ri 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 Ri 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(Yij,Yik|Yil,j<l<k), the partial correlations between Yij and Yik 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 jkth position and 1’s on the main diagonal, where πi,j,j+1=ρi,j,j+1 for j=1,…,ni−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

(2.7)
$$\label{eq2.7} \pi_{i,j,j+k}=\frac{\rho_{i,j,j+k}-\textbf{r}_1^{\mathrm{T}}(j,k)\textbf{R}_2(j,k)^{-1}\textbf{r}_3(j,k)} {[1-\textbf{r}_1^{\rm T}(j,k)\textbf{R}_2(j,k)^{-1}\textbf{r}_1(j,k)]^{1/2} [1-\textbf{r}_3^{\rm T}(j,k)\textbf{R}_2(j,k)^{-1}\textbf{r}_3(j,k)]^{1/2}},$$
where $$\textbf {r}_1^{\mathrm { T}}(j,k)=(\rho _{i,j,j+1},\ldots ,\rho _{i,j,j+k-1})$$, $$\textbf {r}_3^{\mathrm {T}}(j,k)=(\rho _{i,j+k,j+1},\ldots ,\rho _{i,j+k,j+k-1})$$, and R2(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

(2.8)
$$\label{eq2.8} z(\pi_{i,jk}) \equiv \frac{1}{2}\log\left(\frac{1+\pi_{i,jk}}{1-\pi_{i,jk}}\right) =\textbf{w}_{i,jk}^{\mathrm{T}}\boldsymbol{\gamma},$$
where γ is a (q×1)-dimensional unknown parameter and wi,jk is a unit-specific covariate vector. The choice of wi,jk allows the correlation structure to vary by unit-specific covariates and also allows for flexible parsimonious structure for the correlation matrix Ri. For example, we might specify wi,jk=(I(|kj|=1),I(|kj|=1)*age) corresponding to a banded Toeplitz structure of the autocorrelation matrix Πi. The corresponding correlation matrix, Ri 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 (Ri), we propose a score test (Fahrmeir and Tutz, 2000) for the following hypothesis:

$H_0:\boldsymbol{\Omega}_i=\textbf{R}_i\otimes\boldsymbol{\Sigma}\quad \hbox{against } H_a:\boldsymbol{\Omega}_i.$
Note that Ω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 H0. Then
$\textbf{U}^{\rm T}(\hat{\boldsymbol{\beta}}(\hat{\boldsymbol{\omega}}_0), {\hat{\boldsymbol{\omega}}_0}) \textbf{A}({\hat{\boldsymbol{\omega}}_0}) \textbf{U}(\hat{\boldsymbol{\beta}}(\hat{\boldsymbol{\omega}}_0), {\hat{\boldsymbol{\omega}}_0}) \sim \chi_d^2$
asymptotically under H0, 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
$\textbf{U}({\boldsymbol{\beta}}({\boldsymbol{\omega}}),{\boldsymbol{\omega}})=\left(\frac{\partial\log L(\boldsymbol{\theta};y)}{\partial\omega_{11}},\ldots,\frac{\partial\log L(\boldsymbol{\theta};y)}{\partial\omega_{2T,2T}}\right)^{\rm T}.$
The model (2.8) for Ri obviously impacts the number of constraints d. For further details and an assessment of the power of each departure from H0, see Section 3.2.

### Connection between marginal and conditional models

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

(2.9)
$$\label{eq2.9} P_{itjk}^M=\int P_{itjk}^c(b_{itj})f(b_{itj})\,{\mathrm{d}} b_{itj},$$
where $$P_{itjk}^M=P(Y_{itj}\leq k \mid \textbf {x}_{it})$$, $$P_{itjk}^c(b_{itj})=P(Y_{itj}\leq k\mid b_{itj},\textbf {x}_{it})$$ for j=1,2, and f(⋅) is a univariate normal distribution with mean 0 and variance Var(bitj). Therefore, the intercepts Δitjk are a function of βj and Var(bitj). 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,

(2.10)
\begin{align} L(\boldsymbol{\theta} ; y)&=\prod_{i=1}^N \int \prod_{t=1}^{n_i}\prod_{k=1}^{K}(P_{it1k}^c(b_{it1})-P_{it1k-1}^c(b_{it1}))^{y_{it1k-1}}\notag\\ &\quad \times \prod_{k=1}^{K}(P_{it2k}^c(b_{it2})-P_{it2k-1}^c(b_{it2}))^{y_{it2k-1}} \phi(\textbf{b}_{i})\,{\mathrm{d}}\textbf{b}_i,\label{eq2.10} \end{align}
where $$P_{itjk}^c(b_{itj})=P(Y_{itj}\leq k\mid b_{itj},\textbf {x}_{it})$$ for 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

$\sum_{i=1}^N\frac{\partial\log L(\boldsymbol{\theta};\textbf{y}_i)}{\partial \boldsymbol{\theta}}=\sum_{i=1}^N L^{-1}(\boldsymbol{\theta};\textbf{y}_i) \int \frac{\partial L(\boldsymbol{\theta},\textbf{b}_i;\textbf{y}_i)}{\partial\boldsymbol{\theta}} \phi(\textbf{b}_i)\,{\rm d} \textbf{b}_i=0,$
where
\begin{align*} L(\boldsymbol{\theta}; \textbf{y}_i)&=\int \prod_{t=1}^{n_i}\prod_{k=1}^{K}(P_{it1k}^c(b_{it1})-P_{it1k-1}^c(b_{it1}))^{y_{it1k-1}}\\ &\quad \times \prod_{k=1}^{K}(P_{it2k}^c(b_{it2})-P_{it2k-1}^c(b_{it2}))^{y_{it2k-1}} \phi(\textbf{b}_{i})\,{\rm d}\textbf{b}_i \end{align*}
and
(2.11)
$$\label{eq2.11} L(\boldsymbol{\theta},\textbf{b}_i;\textbf{y}_i)=\prod_{t=1}^{n_i}\prod_{k=1}^{K} (P_{it1k}^c(b_{it1})-P_{it1k-1}^c(b_{it1}))^{y_{it1k-1}} \prod_{k=1}^{K}(P_{it2k}^c(b_{it2})-P_{it2k-1}^c(b_{it2}))^{y_{it2k-1}}.$$
The (2(K−1)+2p+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,

$\boldsymbol{\theta}^{(g+1)}=\boldsymbol{\theta}^{(g)}+[\textbf{I}_e(\boldsymbol{\theta}^{(g)};\textbf{y})]^{-1}\frac{\partial\log L}{\partial\boldsymbol{\theta}^{(g)}},$
where Ie(θ) is an empirical and consistent estimator of the information matrix at iteration g,
$\textbf{I}_e(\boldsymbol{\theta} ;\textbf{y})=\sum_{i=1}^N \frac{\partial L(\boldsymbol{\theta};\textbf{y}_i)}{\partial\boldsymbol{\theta}}\frac{\partial L(\boldsymbol{\theta};\textbf{y}_i)}{\partial\boldsymbol{\theta}^{\rm T}}.$
Explicit forms for the terms in the Newton–Raphson algorithm and the derivatives $${\partial \log L}/{\partial \boldsymbol {\theta }}$$ can be found in supplementary material available at Biostatistics online.

## 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

$\hbox{logit}\,P(Y_{itj}\leq k\mid x_{it})=\beta_{j0k}+\beta_{j1} \cdot \hbox{time}_{it}+\beta_{j2} \cdot \hbox{group}_{i}+\beta_{j3} \cdot \hbox{time}_{it}\cdot\hbox{group}_{i},$
for 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, timeit=(t−1)/10, and groupi 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 bit=bi0N(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:

$\hbox{logit}\,P(\hbox{dropout}=t\mid \hbox{dropout}\ge t)=\begin{cases} -2.0+0.2Y_{it-11}+0.2Y_{it-12}, & \hbox{group}_i=0, \\ -2.0+0.1Y_{it-11}+0.1Y_{it-12}, & \hbox{group}_i=1. \end{cases}$
Based on this specification, the observed dropout rates were approximately 50 percent.

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%.

Table 1.

Properties of maximum likelihood estimators

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 wi,jk=(I{|jk|=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 Ri having second-order structure using γ=(0.5,0.4) and wi,jk=(I{|jk|=1},I{|jk|=2}). Then we fit the marginalized model with first-order structure for Ri. 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 Ri was correctly specified. In particular, we assumed the following non-KP structure for Ω in (2.5):

(3.1)
$$\label{eq3.1} {\mathrm{cov}}(b_{it1},b_{it'1})=\tau_1^2\alpha_1^{|t-t'|},\quad {\rm cov}(b_{it2},b_{it'2})=\tau_2^2\alpha_2^{|t-t'|},\quad {\rm cov}(b_{it1},b_{it'2})=\tau_1\tau_2\alpha_{12}^{|t-t'+1|},$$
where α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 Ri). 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.

Table 2.

Summary of ordinal responses of metabolic syndrome and CRP from KoGES Data

Visit

Meta CRP Meta CRP Meta CRP
234 (0.124) 788 (0.673) 420 (0.250) 1084 (0.646) 282 (0.199) 884 (0.624)
1237 (0.656) 220 (0.188) 952 (0.567) 308 (0.183) 849 (0.600) 303 (0.214)
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

Meta CRP Meta CRP Meta CRP
234 (0.124) 788 (0.673) 420 (0.250) 1084 (0.646) 282 (0.199) 884 (0.624)
1237 (0.656) 220 (0.188) 952 (0.567) 308 (0.183) 849 (0.600) 303 (0.214)
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 wi,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.

Table 3.

Description of correlation models in Section 4 in terms of z(π) and Σ

Model $$z(\pi)=\mathbf {w}_{i,jk}^{\mathrm {T}}\boldsymbol {\gamma }$$ in Ri Σ
γ1I{|jk|=1} $$\hbox {diag}(\sigma _1^2,\sigma _2^2)$$
γ=0 Unstructured
γ1I{|jk|=1} Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Sex Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Drink1+γ3I{|jk|=1} Drink2 Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Smoke1+γ3I{|jk|=1} Smoke2 Unstructured
Model $$z(\pi)=\mathbf {w}_{i,jk}^{\mathrm {T}}\boldsymbol {\gamma }$$ in Ri Σ
γ1I{|jk|=1} $$\hbox {diag}(\sigma _1^2,\sigma _2^2)$$
γ=0 Unstructured
γ1I{|jk|=1} Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Sex Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Drink1+γ3I{|jk|=1} Drink2 Unstructured
γ1I{|jk|=1}+γ2I{|jk|=1} Smoke1+γ3I{|jk|=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 △D01=1447.61, p-value = 0.00 on 3 degrees of freedom, △D13=111.274, p-value = 0.00 on 1 degrees of freedom, △D23=7.808, p-value = 0.005 on 1 degrees of freedom, △D34=7.100, p-value = 0.008 on 1 degrees of freedom, △D35=7.122, p-value = 0.028 on 2 degrees of freedom, and △D36=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.

Table 4.

Maximum likelihood estimates for MREMs

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
s11  2.410* (0.092) 2.359* (0.088) 2.382* (0.090) 2.379* (0.090) 2.389* (0.091) 2.386* (0.090)
s21   0.737* (0.070) 0.740* (0.071) 0.740* (0.071) 0.741* (0.071) 0.741* (0.071)
s22  1.832* (0.082) 1.634* (0.083) 1.645* (0.086) 1.643* (0.086) 1.648* (0.086) 1.647* (0.086)
Ri (γ
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
s11  2.410* (0.092) 2.359* (0.088) 2.382* (0.090) 2.379* (0.090) 2.389* (0.091) 2.386* (0.090)
s21   0.737* (0.070) 0.740* (0.071) 0.740* (0.071) 0.741* (0.071) 0.741* (0.071)
s22  1.832* (0.082) 1.634* (0.083) 1.645* (0.086) 1.643* (0.086) 1.648* (0.086) 1.647* (0.086)
Ri (γ
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. s11, s21, and s22 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,

\begin{align*} \log\frac{P(Y_{it1}\leq k)}{P(Y_{it1}>k)}&=\beta_{10k}+0.469\hbox{SEX}_i-0.026 \hbox{AGE}_{it} +0.071\hbox{DRINK1}_{it}+0.121\hbox{DRINK2}_{it}\\ &\quad -0.282\hbox{SMOKE1}_{it}-0.546\hbox{SMOKE2}_{it}, \end{align*}
where (β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,

\begin{align*} \log\frac{P(Y_{it2}\leq k)}{P(Y_{it2}>k)}&=\beta_{20k}-0.337\hbox{SEX}_i-0.028\hbox{AGE}_{it} +0.351\hbox{DRINK1}_{it}+0.111\hbox{DRINK2}_{it}\\ &\quad +0.184\hbox{SMOKE1}_{it}-0.223\hbox{SMOKE2}_{it}, \end{align*}
where (β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 bit1 and bit2 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.

## 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.

## References

Caflisch
R.
Monte Carlo and quasi-Monte Carlo methods.
Acta Numerica
,
1998
, vol.
7
(pg.
1
-
49
()
Chaganty
N. R.
Naik
D. N.
Analysis of multivariate longitudinal data using quasi-least squares.
Journal of Statistical Planning and Inference
,
2002
, vol.
103
(pg.
421
-
436
()
Daniels
M. J.
M.
Modeling covariance matrices via partial autocorrelations.
Journal of Multivariate Analysis
,
2009
, vol.
100
(pg.
2352
-
2363
()
Devaraj
S.
Singh
U.
Jialal
I.
Human C-reactive protein and the metabolic syndrome.
Current Opinion in Lipidology
,
2009
, vol.
20
(pg.
182
-
189
()
Fahrmeir
L.
Tutz
G.
Multivariate Statistical Modelling Based on Generalized Linear Models
,
2000
New York
Springer
Fitzmaurice
G. M.
Laird
N. M.
A likelihood-based method for analysing longitudinal binary responses.
Biometrika
,
1993
, vol.
80
(pg.
141
-
151
()
Galecki
A. T.
General class of covariance structures for two or more repeated factors in longitudinal data analysis.
Communication in Statistics-Theory Method
,
1994
, vol.
23
(pg.
3105
-
3119
()
Gonzalez
J.
Tuerlinckx
F.
Boeck
P. D.
Cools
R.
Numerical integration in logistic-normal models.
Computational Statistics and Data Analysis
,
2006
, vol.
51
(pg.
1535
-
1548
()
Gueorguieva
R.
A multivariate generalized linear mixed model for joint modelling of clustered outcomes in the exponential family.
Statistical Modelling
,
2001
, vol.
1
(pg.
177
-
193
()
Heagerty
P. J.
Marginally specified logistic-normal models for longitudinal binary data.
Biometrics
,
1999
, vol.
55
(pg.
688
-
698
()
Heagerty
P. J.
Marginalized transition models and likelihood inference for longitudinal categorical data.
Biometrics
,
2002
, vol.
58
(pg.
342
-
351
()
Ilk
O.
Daniels
M. J.
Marginalized transition random effects models for multivariate longitudinal binary data.
Canadian Journal of Statistics
,
2007
, vol.
35
(pg.
105
-
123
()
Joe
H.
Generating random correlation matrices based on partial correlations.
Journal of Multivariate Analysis
,
2006
, vol.
97
(pg.
2177
-
2189
()
Judd
K.
Numerical Methods in Economics
,
1998
Boston
MIT Press
Lee
K.
Daniels
M.
A class of Markov models for longitudinal ordinal data.
Biometrics
,
2007
, vol.
63
(pg.
1060
-
1067
()
Lee
K.
Daniels
M.
Marginalized models for longitudinal ordinal data with application to quality of life studies.
Statistics in Medicine
,
2008
, vol.
27
(pg.
4359
-
4380
()
Lee
K.
Joo
Y.
Yoo
J. K.
Lee
J.
Marginalized random effects models for multivariate longitudinal binary data.
Statistics in Medicine
,
2009
, vol.
28
(pg.
1284
-
1300
()
Lee
K.
Kang
S.
Liu
X.
Seo
D.
Likelihood-based approach for analysis of longitudinal nominal data using marginalized random effects models.
Journal of Applied Statistics
,
2011
, vol.
38
(pg.
1577
-
1590
()
Lee
K.
Mercante
D.
Longitudinal nominal data analysis using marginalized models.
Computational Statistics & Data Analysis
,
2010
, vol.
54
(pg.
208
-
218
()
Liu
L. C.
A model for incomplete longitudinal multivariate ordinal data.
Statistics in Medicine
,
2008
, vol.
27
(pg.
6299
-
6309
()
Liu
L. C.
Hedeker
D.
A mixed-effects regression model for longitudinal multivariate ordinal data.
Biometrics
,
2006
, vol.
62
(pg.
261
-
268
()
McCullagh
P.
Regression models for ordinal data.
Journal of the Royal Statistical Society, Series B
,
1980
, vol.
42
(pg.
109
-
142
)
Niederreiter
H.
Random Number Generation and Quasi-Monte Carlo Methods
,
1992
SIAM

CBMS-NSF Regional Conference Series in Applied Mathematics 63.
Post
J. M.
Beebe-Dimmer
J. L.
Morgenstern
H.
Neslund-Dudas
C.
Bock
C. H.
Nock
N.
Rundle
A.
Jankowski
M.
Rybicki
B. A.
The metabolic syndrome and biochemical recurrence following radical prostatectomy.
Prostate Cancer
,
2011
, vol.
2011
pg.
6

Advance Access published September 25, 2011, () ()
M.
Joint mean-covariance models with applications to longitudinal data: unconstrained parameterisation.
Biometrika
,
1999
, vol.
86
(pg.
677
-
690
()
Reilly
M. P.
Wolfe
M. L.
Localio
A. R.
D. J.
C-reactive protein and coronary artery calcification: The Study of Inherited Risk of Coronary Atherosclerosis (SIRCA).
Arteriosclerosis, Thrombosis, and Vascular Biology
,
2003
, vol.
23
(pg.
1851
-
1856
()
Ridker
P. M.
Buring
J. E.
Cook
N. R.
Rifai
N.
C-reactive protein, the metabolic syndrome, and risk of incident cardiovascular events.
Circulation
,
2003
, vol.
107
(pg.
391
-
397
()
Schildcrout
J.
Heagerty
P. J.
Marginalized models for moderate to long series of longitudinal binary response data.
Biometrics
,
2007
, vol.
63
(pg.
322
-
331
()
Todem
D.
Kim
K.
Lesaffre
E.
Latent-variable models for longitudinal data with bivariate ordinal outcomes.
Statistics in Medicine
,
2007
, vol.
26
(pg.
1034
-
1054
()
Wang
Y.
Daniels
M. J.
Bayesian modeling of the dependence in longitudinal data via partial autocorrelations and marginal variances.
Journal of Multivariate Analysis
,
2013
, vol.
116
(pg.
130
-
140
()
Ye
X.
Yu
Z.
Li
H.
Franco
O. H.
Liu
Y.
Lin
X.
Distribution of C-reactive protein and its association with metabolic syndrome in middle-aged and older Chinese people.
Journal of the American College of Cardiology.
,
2007
, vol.
49
(pg.
1798
-
1805
()