2024-2 Contents

1.   Introduction

Longitudinal data arises frequently in the economics, social sciences, epidemiology and biological research. The observations of each subject are measured repeatedly over time and thus are intrinsically correlated^[1]. It is crucial to correctly address the within-subject covariation structure since ignoring this structure may lead to inefficient estimators of the mean parameter. Moreover, the covariation structure itself may be of scientific interest^[2], while covariance and precision matrices suffer from positive definiteness constraints and high-dimensional parameters, and the correlation matrix should satisfy the constraint that its diagonals must be 1’s in addition to the aforementioned two issues. Thus, there is a frequently-used strategy to estimate the correlation matrix, i.e., first estimating the covariance or precision matrix and then calculating the corresponding correlation matrix, which motivates us to develop the decomposition method proposed in this paper.

Now it is vital to decide to decompose and model whether the covariance matrix or the precision matrix. In fact, the covariance or precision matrix may possess typical structures, and the matrix inversion operation does not protect the original structure, such as sparsity and regression relationship, then decomposing and modeling the inverse of the covariance matrix or the inverse of the precision matrix may result in efficiency loss. For the precision matrix, Pourahmadi^[3] introduced the modified Cholesky decomposition (MCD), which leads to an unconstrained parameterization for the precision matrix that automatically guarantees its positive definiteness. The entries in MCD can be interpreted as generalized autoregressive parameters and innovation variances, which can be estimated by fitting a class of unconstrained joint mean-covariance models via maximum likelihood estimation. For the covariance matrix, Zhang and Leng^[4] analyzed within-subject covariation by decomposing the covariance matrix itself rather than its inverse by using the moving average Cholesky decomposition (MACD). The entries in this decomposition are moving average parameters and innovation variances. Based on MACD, the corresponding correlation expression depends on the innovation variances, and thus is not necessarily robust with respect to their model misspecification. In contrast, Chen and Dunson^[5] proposed an alternative Cholesky decomposition (ACD), which directly models the covariance matrix but in a way that the estimates of the correlation matrix do not depend on the quality of modeling and estimating the innovation variances; that is, estimation of the correlation matrix is robust with respect to model misspecification of the innovation variances, the components shared by both MACD and ACD. ACD is closely related to the moving average model of “standardized” measurements on a longitudinal subject^[6,7]. By directly considering the correlation matrix, Zhang et al.^[8] developed the hyperspherical parameterization of the Cholesky factor (HPC), which is very appealing since the resulting parameters are unconstrained on the support and directly interpretable in regard to the correlations. HPC can always lead to a more parsimonious model, however, at the cost of intensive computations. To the best of our knowledge, when the precision matrix has a typical structure, robust estimation of the correlation matrix against model misspecification of innovation variances has been much less investigated.

Moreover, longitudinal data may suffer from outliers, and several heavy distributed joint modeling approaches have been investigated. Lin and Wang^[9] proposed a joint mean and scale covariance model based on multivariate t distribution and MCD. Maadooliat et al.^[7] further investigated the joint model based on the t distribution and ACD. However, this kind of joint t-regression approach is computationally intensive since joint regression parameters together with the degree of freedom need to be estimated. As an alternative, Guney et al.^[10] adopted a joint model with multivariate Laplace distribution and MCD. The Laplace distribution is heavy-tailed and robust to outliers. It does not have the degree of freedom parameter; thus, its computation is more convenient than that of the t distribution. However, this joint modeling approach could not lead to the aforementioned robust estimation of the correlation matrix.

In this paper, motivated by ACD, we propose an alternative modified Cholesky decomposition (AMCD) of the precision matrix for longitudinal data. Specifically, the inverse of the diagonal matrix of innovation standard deviations is placed outside the two triangular matrices, which alone determines the correlations, and thus results in the aforementioned robust estimation of the correlation matrix. Then, we establish a joint mean-covariance model with multivariate normal distribution and AMCD, and the quasi-Fisher scoring algorithm, and show that the maximum likelihood estimators are consistent and asymptotically normally distributed. Furthermore, we present the double-robust joint modeling approach with multivariate Laplace distribution and AMCD, and the quasi-Newton algorithm for maximum likelihood estimation. We carry out the simulation studies, cattle data analysis, and sleep dose-response data analysis, which indicated that the proposed AMCD method performed well.

The rest of this paper is organized as follows. In Section 2, we review MACD, ACD and MCD, propose AMCD, and discuss their close relationships. In Section 3, we establish the joint normal mean covariance model with AMCD and give the iterative algorithm and asymptotic properties. In Section 4, we investigate the joint model with multivariate Laplace distribution and AMCD, and develop the quasi-Newton algorithm for maximum likelihood estimation. In Section 5, we carry out the simulations. In Section 6, real data analysis is conducted. In Section 7, the paper is summarized. All the calculations and proofs are in the Appendix.

2.   Existing and proposed Cholesky-type decompositions of the covariation structure

Assume that ${y}_{ij}$ is the measurement at the $j$th time point ${t}_{ij}$ for the $i$th subject ($i=\text{1},\dots ,n,  j=\text{1},\dots ,{m}_{i}$). Let ${\boldsymbol{y}}_{i}={\left({y}_{i1},\dots ,{y}_{i{m}_{i}}\right)}^{\text{T}}$, and ${\boldsymbol{t}}_{i}={\left({t}_{i1},\dots ,{t}_{i{m}_{i}}\right)}^{\text{T}}$. Moreover, suppose that ${\boldsymbol{y}}_{i}$ follows a normal distribution with $E\left({\boldsymbol{y}}_{i}\right)={\boldsymbol{\mu }}_{i}= {\left({\mu }_{i1},\dots ,{\mu }_{i{m}_{i}}\right)}^{\text{T}}$ and ${\rm var}\left({\boldsymbol{y}}_{i}\right)={\boldsymbol{\varSigma }}_{i}$. Then let ${\boldsymbol{y}}_{i}={\boldsymbol{\mu }}_{i}+{\boldsymbol{r}}_{i},$ where ${\boldsymbol{r}}_{i}={\left({r}_{i1},\dots ,{r}_{i{m}_{i}}\right)}^{\text{T}}$ is the normal random error, with $E\left({\boldsymbol{r}}_{i}\right)={\boldsymbol{0}}$ and ${\rm var}\left({\boldsymbol{r}}_{i}\right)={\boldsymbol{\varSigma }}_{i}$.

In some cases, the covariance matrix ${\boldsymbol{\varSigma }}_{i}$ (or its inverse, i.e., the precision matrix) is of scientific interest and possesses a typical structure. Then, a specific Cholesky-type decomposition can be established to adapt to the corresponding typical structure. Furthermore, the correlation matrix might be of particular interest; therefore, it is worthwhile to develop a robust estimation procedure based on the specific decomposition process, of course, with some suitable modifications. Along this line, several existing Cholesky-type decompositions are discussed, the proposed AMCD is introduced, and their close connections are investigated.

2.1   MACD of the covariance matrix

If the covariance matrix itself possesses a typical structure, it is reasonable to decompose the symmetric and positive definite covariance matrix as ${\boldsymbol{\varSigma }}_{i}={\boldsymbol{C}}_{i}{\boldsymbol{C}}_{i}^{\text{T}}$ based on the standard Cholesky decomposition. Here ${\boldsymbol{C}}_{i}$ is a unique lower triangular matrix with positive diagonal elements. Define $\text{diag}\left({\boldsymbol{C}}_{i}\right)={\left({c}_{i11},\dots ,{c}_{i{m}_{i}{m}_{i}}\right)}^{\text{T}}$, and ${\boldsymbol{\varLambda }}_{i}=\text{diag}\left({c}_{i11},\dots ,{c}_{i{m}_{i}{m}_{i}}\right)$. For simplification, let ${d}_{ij}={c}_{ijj},j=1,\dots ,{m}_{i}$; then, ${\boldsymbol{\varLambda }}_{i}=\text{diag}\left({d}_{i1},\dots ,{d}_{i{m}_{i}}\right)$. In Ref. [4], the standard Cholesky decomposition was transformed into MACD; that is, the matrix ${\boldsymbol{\varLambda }}_{i}$ is located inside two triangular matrices:

where ${\boldsymbol{L}}_{i}={\boldsymbol{C}}_{i}{\boldsymbol{\varLambda }}_{i}^{-1}$ results in a standardized ${\boldsymbol{C}}_{i}$, namely dividing each column of ${\boldsymbol{C}}_{i}$ by its diagonal elements. Let ${\boldsymbol{L}}_{i}={\left({l}_{ijk}\right)}_{{m}_{i}\times {m}_{i}}$ and ${\boldsymbol{D}}_{i}={\boldsymbol{\varLambda }}_{i}^{2}=\text{diag}\left({d}_{i1}^{2},\dots ,{d}_{i{m}_{i}}^{2}\right)$. Denote ${\boldsymbol{\varepsilon }}_{i}={\left({\varepsilon }_{i1},\dots ,{\varepsilon }_{i{m}_{i}}\right)}^{\text{T}}={\boldsymbol{L}}_{i}^{-1}{\boldsymbol{r}}_{i}$; then, $\text{var}\left({\boldsymbol{\varepsilon }}_{i}\right)={\boldsymbol{\varLambda }}_{i}^{2}$ and ${\boldsymbol{r}}_{i}={\boldsymbol{L}}_{i}{\boldsymbol{\varepsilon }}_{i}$. Note that ${\boldsymbol{L}}_{i}$ is a lower triangular matrix with diagonals being 1’s; then, a moving average representation for the residual ${r}_{ij}$ can be expressed as follows:

where $\displaystyle \sum _{t=1}^{0}{l}_{ijt}{\varepsilon }_{it}$ is denoted as 0. Then, from (2), for any $1\leqslant j,k\leqslant {m}_{i}$, with $j\wedge k=\mathrm{min}\{j,k\}$, it follows that

and thus, the correlation between ${r}_{ij}$ and ${r}_{ik}$ is given by

which is determined by the ${\boldsymbol{L}}_{i}$ and ${\boldsymbol{\varLambda }}_{i}$ matrices. Then it might not be robust against misspecification of the model for the innovation variances ${d}_{ij}^{2},j=1,\dots ,{m}_{i}$.

2.2   ACD of the covariance matrix

In contrast, the standard Cholesky decomposition can be transformed into ACD^[7]; that is, the matrix ${\boldsymbol{\varLambda }}_{i}$ is located outside the two triangular matrices:

where ${\boldsymbol{A}}_{i}={\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{C}}_{i}$ leads to a standardized ${\boldsymbol{C}}_{i}$, namely dividing each row of ${\boldsymbol{C}}_{i}$ by its diagonal elements. Denote ${\boldsymbol{A}}_{i}={\left({a}_{ijk}\right)}_{{m}_{i}\times {m}_{i}}$. It is obvious that ${\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{r}}_{i}$ has covariance matrix ${\boldsymbol{A}}_{i}{\boldsymbol{A}}_{i}^{\text{T}}$, within which ${\boldsymbol{\varLambda }}_{i}$ disappears. Specifically, let ${\boldsymbol{\varepsilon}}_{i}={\left({\varepsilon}_{i1},\dots ,{\varepsilon}_{i{m}_{i}}\right)}^{\text{T}}={\left({\boldsymbol{\varLambda }}_{i}{\boldsymbol{A}}_{i}\right)}^{-1}{\boldsymbol{r}}_{i}$; then, $\text{var}\left({\boldsymbol{\varepsilon}}_{i}\right)={\boldsymbol{I}}_{{m}_{i}}$ and ${\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{r}}_{i}={\boldsymbol{A}}_{i}{\boldsymbol{\varepsilon}}_{i}$. In fact, ${\boldsymbol{A}}_{i}$ is also a lower triangular matrix with diagonals being 1's; then, the “standardized” residual ${r}_{ij}/{d}_{ij}$ can be represented as

From (4), it follows that

and thus

which is determined by the ${\boldsymbol{A}}_{i}$ matrix alone. Modeling a correlation matrix via ACD is highly appealing since it is robust with respect to misspecification of the model for ${d}_{ij}^{2}$, $j=1,\dots ,{m}_{i}$.

2.3   MCD of the precision matrix

However, the precision matrix might possess a typical structure and should be decomposed appropriately. Based on the standard Cholesky decomposition, ${\boldsymbol{\varSigma }}_{i}^{-1}={\boldsymbol{C}}_{i}^{-\text{T}}{\boldsymbol{C}}_{i}^{-1}={\boldsymbol{B}}_{i}^{\text{T}}{\boldsymbol{B}}_{i}$, where ${\boldsymbol{B}}_{i}={\boldsymbol{C}}_{i}^{-1}$ is a unique lower triangular matrix with positive diagonals. Similar to (1),

where ${\boldsymbol{F}}_{i}^{\text{T}}={\boldsymbol{B}}_{i}^{\text{T}}{\boldsymbol{\varLambda }}_{i}$ leads to a standardized ${\boldsymbol{B}}_{i}^{\text{T}}$, namely dividing each column of ${\boldsymbol{B}}_{i}^{\text{T}}$ by its diagonal elements ${d}_{ij}^{-2}$. In fact, (5) can be transformed into ${\boldsymbol{F}}_{i}{\boldsymbol{\varSigma }}_{i}{\boldsymbol{F}}_{i}^{\text{T}}={\boldsymbol{\varLambda }}_{i}^{2}$, i.e., MCD^[6,9]. Notice that ${\boldsymbol{F}}_{i}$ can be expressed as a lower triangular matrix, with the diagonal entries being 1’s and the $\left(j,k\right)$th nonzero entries being $-{f}_{ijk},k < j$. Then let ${\boldsymbol{\varepsilon }}_{i}={\boldsymbol{F}}_{i}{r}_{i}$, and $\text{var}\left({\boldsymbol{\varepsilon }}_{i}\right)={\boldsymbol{\varLambda }}_{i}^{2}$, from which the autoregressive representation for the residuals ${r}_{ij}$ can be expressed as follows:

Then, by simple calculation, the correlation between ${r}_{ij}$ and ${r}_{ik}$ is determined by the ${\boldsymbol{F}}_{i}$ and ${\boldsymbol{\varLambda }}_{i}$ matrices, and thus is not robust against misspecification of the model for innovation variances.

2.4   AMCD of the precision matrix

Conversely, the proposed AMCD causes the matrix ${\boldsymbol{\varLambda }}_{i}^{-1}$ to be outside the two triangular matrices:

where ${\boldsymbol{T}}_{i}^{\text{T}}={\boldsymbol{\varLambda }}_{i}{\boldsymbol{B}}_{i}^{\text{T}}$ results in a standardized ${\boldsymbol{B}}_{i}^{\text{T}}$, namely dividing each row of ${\boldsymbol{B}}_{i}^{\text{T}}$ by its diagonals ${d}_{ij}^{-2}$. In fact, (6) is equivalent to ${\boldsymbol{T}}_{i}{\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{\varSigma }}_{i}{\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{T}}_{i}^{\text{T}}={\boldsymbol{I}}_{{m}_{i}}$. Apparently, the covariance matrix of ${\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{r}}_{i}$ is ${\boldsymbol{T}}_{i}^{-1}{\boldsymbol{T}}_{i}^{-\text{T}}$, where ${\boldsymbol{\varLambda }}_{i}$ disappears. Let ${\boldsymbol{\varepsilon}}_{i}={\boldsymbol{T}}_{i}{\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{r}}_{i}$; then, $\text{var}\left({\boldsymbol{\varepsilon}}_{i}\right)={\boldsymbol{I}}_{{m}_{i}}$. Note that

Then, the “standardized” residual ${r}_{ij}/{d}_{ij}$ can be modeled by

where the generalized autoregressive parameter ${\phi }_{ijt}$’s are unconstrained, the “standardized” factor ${d}_{ij}$’s are restricted to be positive, the innovation variances $\text{var}\left({\epsilon}_{ij}\right)=1$, and the innovation covariances $\text{cov}\left({\epsilon}_{ij},{\epsilon}_{ik}\right)=0,j\ne k$, implying their uncorrelation in general, more specifically, their independence when assuming normality. Note that ${d}_{ij}^{2}$’s are the innovation variances in MCD. For ease of comparison, they are also referred to as innovation variances rather than squared “standardized” factors, in AMCD.

Obviously, $\text{corr}\left({r}_{ij},{r}_{ik}\right)$ is determined by the ${\boldsymbol{T}}_{i}$ matrix alone, and thus is robust against misspecification of the model for ${d}_{ij}^{2}$, $j=1,\dots ,{m}_{i}$. This advantage makes the proposed AMCD attractive for modeling the correlation matrix. Moreover, it is more suitable to adopt AMCD than ACD when the precision matrix, rather than the covariance matrix, possesses a typical structure.

3.   Normal joint modeling approach

3.1   Joint mean-covariance model

Based on the proposed AMCD in (6), the estimation of ${\boldsymbol{\varSigma }}_{i}^{-1}$ is equivalent to that of ${\boldsymbol{T}}_{i}$ and ${\boldsymbol{\varLambda }}_{i}$. Define ${\text{log}}\left({\boldsymbol{\varLambda }}_{i}^{2}\right)= {\text{diag}}\left\{{\text{log}}\left({d}_{i1}^{2}\right),\dots ,{\text{log}}\left({d}_{i{m}_{i}}^{2}\right)\right\}$. Following the general approach in Ref. [11], the unconstrained nonredundant entries of $\text{log}\left({\boldsymbol{\varLambda }}_{i}^{2}\right)$ and ${\boldsymbol{T}}_{i}$ can be modeled together with ${\boldsymbol{\mu }}_{i}$ by generalized linear regression models

where $g\left(\cdot \right)$ is assumed to be a monotone and differentiable known link function; ${\boldsymbol{x}}_{ij}$, ${\boldsymbol{z}}_{ij}$ and ${\boldsymbol{w}}_{ijk}$ are the vectors of covariates; and ${\boldsymbol{\beta}} ,{\boldsymbol{\lambda}}$ and ${\boldsymbol{\gamma}}$ are the $p\times 1$, $s\times 1$ and $q\times 1$ vectors, respectively, of the corresponding associated parameters. The covariates ${\boldsymbol{x}}_{ij}$, ${\boldsymbol{z}}_{ij}$ and ${\boldsymbol{w}}_{ijk}$ may consist of the baseline covariates, polynomials in time (for ${\boldsymbol{x}}_{ij}$ and ${\boldsymbol{z}}_{ij}$) or time lag (for ${\boldsymbol{w}}_{ijk}$), and even their interactions. For instance, when the entries of ${\boldsymbol{\mu }}_{i}$, $\text{log}\left({\boldsymbol{\varLambda }}_{i}^{2}\right)$ and ${\boldsymbol{L}}_{i}$ are modeled by polynomials in time and time lag, the covariates may be constructed as

assuming that the correlation between ${r}_{ij}$ and ${r}_{ik}$ relies only on the time lag between ${t}_{ij}$ and ${t}_{ik}$ $\left(1\leqslant k < j\leqslant {m}_{i}\right)$.

3.2   Maximum likelihood estimation

According to the proposed AMCD, inversing (6), it can be obtained that

Let ${\boldsymbol{\theta}} ={\left({\boldsymbol{\beta }}^{\text{T}},{\boldsymbol{\lambda }}^{\text{T}},{\boldsymbol{\gamma }}^{\text{T}}\right)}^{\text{T}}$; then, from (6) and (8), twice the negative log-likelihood function of the multivariate normal distribution, except for a constant, can be written as

where ${\boldsymbol{r}}_{i}={\boldsymbol{y}}_{i}-{\boldsymbol{\mu }}_{i}$. In this case, the joint mean-covariance model can be referred to as the normal joint modeling model (NJMM). By taking partial derivatives of $l\left({\boldsymbol{\theta}} \right)$ with respect to ${\boldsymbol{\beta}}$, ${\boldsymbol{\lambda}}$ and ${\boldsymbol{\gamma}}$, respectively, the maximum likelihood estimating equations for these parameters become

Here $\dfrac{\partial {\boldsymbol{\mu }}_{i}^{\text{T}}}{\partial {\boldsymbol{\beta}} }$ is the $p\times {m}_{i}$ matrix with $j$th column $\dfrac{\partial {\mu }_{ij}}{\partial {\boldsymbol{\beta}} }={\dot{g}}^{-1}\left({\boldsymbol{x}}_{ij}^{\text{T}}{\boldsymbol{\beta}} \right){\boldsymbol{x}}_{ij}$, ${\dot{g}}^{-1}\left(\cdot \right)$ is the derivative of the inverse function ${g}^{-1}\left(\cdot \right)$, and we denote $\mu \left(\cdot \right)={g}^{-1}\left(\cdot \right)$; ${\boldsymbol{Z}}_{i}={\left({\boldsymbol{z}}_{i1}^{\text{T}},\dots ,{\boldsymbol{z}}_{i{m}_{i}}^{\text{T}}\right)}^{\text{T}}$, ${\boldsymbol{h}}_{i}=\text{diag}\left({\boldsymbol{T}}_{i}^{\text{T}}{\boldsymbol{T}}_{i}{\boldsymbol{\varLambda }}_{i}^{-1}{\boldsymbol{r}}_{i}{\boldsymbol{r}}_{i}^{\text{T}}{\boldsymbol{\varLambda }}_{i}^{-1}\right)$, and ${\boldsymbol{1}}_{{m}_{i}}$ is a ${m}_{i}\times 1$ vector of 1’s; $\dfrac{\partial {\boldsymbol{T}}_{i}^{-1}}{\partial {\gamma }_{j}}=-{\boldsymbol{T}}_{i}^{-1}\dfrac{\partial {\boldsymbol{T}}_{i}}{\partial {\gamma }_{j}}{\boldsymbol{T}}_{i}^{-1}$ with

for $j=1,\dots ,q$, and ${\boldsymbol{w}}_{ijk}={\left({w}_{ijk,1},\dots ,{w}_{ijk,q}\right)}^{\text{T}}$.

Since the solutions satisfy equations (9), the parameters ${\boldsymbol{\beta}}$, ${\boldsymbol{\lambda}}$ and ${\boldsymbol{\gamma}}$ can be solved iteratively with the others kept fixed. More specifically, the numerical solutions for these parameters can be calculated by the quasi-Fisher scoring algorithm.

First, the expectations of the Hessian matrices are listed below and discussed explicitly in Appendix A.1. That is,

where

with ${\boldsymbol{A}}\circ {\boldsymbol{B}}$ denoting the Hadamard product of matrix ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$, ${\boldsymbol{v}}_{ijk}=-\displaystyle \sum _{t=k+1}^{j-1}\dfrac{\partial {a}_{itk}}{\partial {\boldsymbol{\gamma}} }{\phi }_{ijt}$, and ${a}_{itk}$ being the $\left(t,k\right)$th element of ${\boldsymbol{T}}_{i}^{-1}$.

Then the iterative algorithm is as follows.

Step 1. Initialize the starting values ${\boldsymbol{\beta }}^{\left(0\right)}$, ${\boldsymbol{\lambda }}^{\left(0\right)}$ and ${\boldsymbol{\gamma }}^{\left(0\right)}$. Set $k=0$.

Step 2. Compute ${\boldsymbol{\varSigma }}_{i}$ with given ${\boldsymbol{\lambda }}^{\left(k\right)}$ and ${\boldsymbol{\gamma }}^{\left(k\right)}$. Update ${\boldsymbol{\beta}}$ via

Step 3. Given ${\boldsymbol{\beta}} ={\boldsymbol{\beta }}^{\left(k+1\right)}$, update ${\boldsymbol{\lambda}}$ and ${\boldsymbol{\gamma}}$ as

Step 4. Set $k\leftarrow k+1$. Repeat steps 2 and 3 until a prespecified convergence criterion is satisfied.

Note that ${\boldsymbol{\lambda}}$ and ${\boldsymbol\gamma}$ are updated together, which is caused by their asymptotic dependence, as seen from Theorem 1 in Section 3.3. It can only be ensured that this algorithm converges to a local optimum which relies critically on the starting values. It is natural to choose the starting value of ${\boldsymbol{\beta}}$ as the least square estimator in equation (10) with ${\boldsymbol{\varSigma }}_{i}$’s set as identity matrices. Then, the starting value of ${\boldsymbol{\gamma}}$ can be determined assuming ${\boldsymbol{T}}_{i}={\boldsymbol{I}}_{{m}_{i}}$, and the starting value of ${\boldsymbol{\lambda}}$ can be given as the least square estimator based on the residuals. The starting estimates of ${\boldsymbol{\beta}}$ and ${\boldsymbol{\lambda}}$ are clearly $\sqrt{n}$-consistent. Moreover, the negative log-likelihood function is asymptotically convex around a small neighbourhood of true parameters, according to the theoretical results in Theorem 1 in Section 3.3 and the proofs in Appendix A.2. Then the final estimates via the iterative algorithm are guaranteed to be the global optima and more efficient than the starting values in terms of the asymptotical viewpoint. For simulation studies and real data analysis, convergence was usually achieved within several iterations.

3.3   Asymptotic properties

In this section, we establish the strong consistency and asymptotic normality of the maximum likelihood estimators. The following regularity conditions are imposed for the theoretical analysis.

(C1) The dimensions $p$, $s$ and $q$ of covariates ${\boldsymbol{x}}_{ij}$, ${\boldsymbol{z}}_{ij}$ and ${\boldsymbol{\varOmega }}_{ijk}$ are fixed, $n\to \text{∞}$, and $\underset{1\leqslant i\leqslant n}{\mathrm{max}}\;{m}_{i}$ is bounded.

(C2) The parameter space $\bf{\varTheta }$ of ${\boldsymbol{\theta}}$ is a compact subset of ${\boldsymbol{R}}^{p+s+q}$, and the true parameter ${\boldsymbol{\theta }}_{0}={\left({\boldsymbol{\beta }}_{0}^{\text{T}},{\boldsymbol{\lambda }}_{0}^{\text{T}},{\boldsymbol{\gamma }}_{0}^{\text{T}}\right)}^{\text{T}}$ lies in the interior of $\bf{\varTheta }$.

(C3) When $n\to \mathrm{\infty }$, ${\boldsymbol{I}}\left({\boldsymbol{\theta }}_{0}\right)/n$ converges to a positive definite matrix ${\boldsymbol{\mathcal{I}}}\left({\boldsymbol{\theta }}_{0}\right)$.

Condition (C1) is standard for practical longitudinal data analysis. Condition (C2) is natural in the theoretical study of the maximum likelihood estimation. Condition (C3) is conventional in regression analysis for modeling unbalanced longitudinal data.

Theorem 1. If $n\to \mathrm{\infty }$ and regularity conditions (C1)–(C3) hold, then (a) the maximum likelihood estimator $\widehat{\boldsymbol{\theta }}={\left({\widehat{\boldsymbol{\beta }}}^{T},{\widehat{\boldsymbol{\lambda }}}^{T},{\widehat{\boldsymbol{\gamma }}}^{T}\right)}^{\text{T}}$ is strongly consistent for the true value ${\boldsymbol{\theta }}_{0}={\left({\boldsymbol{\beta }}_{0}^{\text{T}},{\boldsymbol{\lambda }}_{0}^{\text{T}},{\boldsymbol{\gamma }}_{0}^{\text{T}}\right)}^{\text{T}}$, and (b) $\widehat{\boldsymbol{\theta }}={\left({\widehat{\boldsymbol{\beta }}}^{T},{\widehat{\boldsymbol{\lambda }}}^{T},{\widehat{\boldsymbol{\gamma }}}^{T}\right)}^{\text{T}}$ is asymptotically normal, that is, $\sqrt{n}\left(\widehat{\boldsymbol{\theta }}-{\boldsymbol{\theta }}_{0}\right)\xrightarrow{\text{d}} N\left\{{\boldsymbol{0}},{\boldsymbol{{\mathcal{I}}}}^{-1}\left({\boldsymbol{\theta }}_{0}\right)\right\}$.

It can be easily seen that the Fisher information matrix ${\boldsymbol{\mathcal{I}}}\left({\boldsymbol{\theta }}_{0}\right)$ is a block matrix, more specifically,

from which $\widehat{\boldsymbol{\beta }}$ is asymptotically independent of $\widehat{\boldsymbol{\lambda} }$ and $\widehat{\boldsymbol{\gamma }}$, whereas $\widehat{\boldsymbol{\lambda}}$ and $\widehat{\boldsymbol{\gamma }}$ are not asymptotically independent. Since $\widehat{\boldsymbol{\theta }}={\left({\widehat{\boldsymbol{\beta }}}^{T},{\widehat{\boldsymbol{\lambda }}}^{T},{\widehat{\boldsymbol{\gamma }}}^{T}\right)}^{\text{T}}$ is a consistent estimator for ${\boldsymbol{\theta }}_{0}$, the asymptotic covariance matrix ${\boldsymbol{\mathcal{I}}}^{-1}$ can be consistently estimated by the inverse of a matrix with block components

4.   Laplace joint modeling approach

4.1   Laplace joint modeling model

In this section, we investigate the joint modeling approach based on the multivariate Laplace distribution, which is useful when the response variable has heavier tails. However, there are several kinds of multivariate Laplace distributions, which can be respectively regarded as a particular multivariate Linnik distribution^[12], a special case of the multivariate power exponential (PE) distribution^[13,14], a Gaussian scale mixture^[15], and so on^[16]. Among them, we adopt the special case of the multivariate PE distribution.

If ${\boldsymbol{y}}_{i}$ follows a ${m}_{i}$-dimensional PE distribution, $i=1,\dots ,n,$ then its density function is

where

with the location vector ${\boldsymbol{\mu }}_{i}\in {\boldsymbol{R}}$, the positive definite dispersion matrix ${\boldsymbol{\varSigma }}_{i}$ and $\upsilon$. When $\upsilon =1$, (12) corresponds to a multivariate normal distribution. When $\upsilon =0.5$, (12) is the density function of a multivariate Laplace distribution, i.e.,

then E(${\boldsymbol{y}}_{i}$)= ${\boldsymbol{\mu }}_{i}$ and var(${\boldsymbol{y}}_{i}$)=4(${m}_{i}$+1)${\boldsymbol{\varSigma }}_{i}$. Based on the multivariate Laplace distribution (13), the proposed AMCD of ${\boldsymbol{\varSigma }}_{i}^{-1}$ (6) and the joint model in Section 3.1, we can obtain the Laplace joint modeling model (LJMM).

4.2   The quasi-Newton algorithm for maximum likelihood estimation

Twice the negative log-likelihood function of the multivariate Laplace distribution, except for a constant, can be written as

Let ${{{\boldsymbol{\Delta }}_{i}={\boldsymbol{r}}}_{i}^{\text{T}}{\boldsymbol{\varLambda}} }_{i}^{-1}{{\boldsymbol{T}}}_{i}^{\text{T}}{{\boldsymbol{T}}}_{i}{{\boldsymbol{\varLambda}} }_{i}^{-1}{{\boldsymbol{r}}}_{i}$. By taking partial derivatives of ${l}_{L}\left({\boldsymbol{\theta}} \right)$ with respect to ${\boldsymbol{\beta}}$, ${\boldsymbol{\lambda}}$ and ${\boldsymbol{\gamma}}$, respectively, the maximum likelihood estimating equations are

where ${\bf{\Delta }}_{i}^{-\frac{1}{2}}$ can be regarded as a weight, providing robustness against outliers in the data.

We then estimate ${\boldsymbol{\theta}}$ by minimizing expression (14) via the quasi-Newton algorithm. More specifically, the parameters in ${\boldsymbol{\theta}}$ can be split into ${\boldsymbol{\theta }}_{1}={\boldsymbol{\beta}}$ and ${\boldsymbol{\theta }}_{23}={\left({\boldsymbol{\lambda }}^{\text{T}},{\boldsymbol{\gamma }}^{\text{T}}\right)}^{\text{T}}$, and solved sequentially with the other parameters held fixed. The detailed algorithm is as follows.

Step 1. Initialize the parameters ${\boldsymbol{\beta }}^{\left(0\right)}$, ${\boldsymbol{\lambda }}^{\left(0\right)}$and ${\boldsymbol{\gamma }}^{\left(0\right)}$ as in Section 3.2, the inverse Hessian matrix ${\boldsymbol{H}}_{1}^{\left(0\right)}$ for ${\boldsymbol{\theta }}_{1}$ and ${\boldsymbol{H}}_{23}^{\left(0\right)}$ for ${\boldsymbol{\theta }}_{23}$ as identity matrix. Set $k=0$.

Step 2. For ${\boldsymbol{\theta }}_{1}={\boldsymbol{\beta}}$, compute the score function

Step 3. Update the search direction

and compute the step size by an approximate line minimization

Step 4. Update ${\boldsymbol{\theta }}_{1}$ and the score function as

Step 5. Update the inverse Hessian matrix via the BFGS formula

where

Step 6. For ${\boldsymbol{\theta }}_{23}={\left({\boldsymbol{\lambda }}^{\text{T}},{\boldsymbol{\gamma }}^{\text{T}}\right)}^{\text{T}}$, compute the score function

Step 7. Update the search direction

and compute the step size by an approximate line minimization:

Step 8. Update ${{\boldsymbol{\theta}} }_{23}$ and the score function as

Step 9. Update the inverse Hessian matrix via the BFGS formula as in Step 5, with the subscript “1” replaced by “23”.

Step 10. Set $k=k+1$ and repeat Steps 2 to 9 until a prespecified convergence criterion is met.

Note that the BFGS algorithm is implemented in this section since it is demonstrated to be one of the best quasi-Newton algorithms for solving unconstrained smooth optimization problems and performs very well here. Alternative quasi-Newton algorithms are also deserved to be investigated in the future. For more details, see Refs. [17, 18].

5.   Simulation studies

In this section, we investigate the finite-sample performance of the proposed AMCD, compare the robustness of AMCD, MCD and ACD for modeling correlation matrices based on various typical covariation structures, and compare the modeling capacities of NJMM and LJMM with AMCD based on different distributions.

Study 1. The 1000 replications of the datasets were generated from the following NJMM:

where ${r}_{ij}$ follows the proposed AMCD, i.e., expression (7), with

Here, $n=50, 100$, or 200, ${m}_{i}-1\sim$binomial(6,0.8), and ${t}_{ij}$’s were generated from the uniform distribution in the unit interval. The covariates ${\boldsymbol{x}}_{ij}={\left({x}_{ij1},{x}_{ij2},{x}_{ij3}\right)}^{\text{T}}$ were generated from the multivariate normal distribution with mean zero, marginal variance 1 and correlation 0.5. In addition, we take ${\boldsymbol{z}}_{ij}={\boldsymbol{x}}_{ij}$ and ${\boldsymbol{w}}_{ijk}=\left\{1,{t}_{ij}-{t}_{ik},{\left({t}_{ij}-{t}_{ik}\right)}^{2}\right\}$. In Table 1, “Bias” denotes the average bias of the parameter estimates, and “SD” represents the sample standard deviation of the estimates for the parameters, which can be regarded as the true standard deviation of the resulting estimates. “SE” denotes the sample average of the estimated standard errors by the formula, and “Std” represents the standard deviation of these standard errors. Table 1 indicates that the proposed AMCD approach results in unbiased parameter estimates, and the standard error formula works well, especially when $n$ is large.

Table  1.  Simulation results for Study 1. All the results are multiplied by a factor ${10}^{2}$.

		$\mathit{n}=50$				$\mathit{n}=100$				$\mathit{n}=200$
	True	Bias	SD	SE (Std)		Bias	SD	SE (Std)		Bias	SD	SE (Std)
${\mathit{\beta }}_{0}$	1	−0.11	2.00	1.88 (0.19)		0.02	1.38	1.35 (0.09)		−0.01	0.96	0.96 (0.04)
${\mathit{\beta }}_{1}$	−0.6	−0.17	3.22	3.07 (0.38)		0.02	2.37	2.19 (0.19)		−0.02	1.53	1.54 (0.10)
${\mathit{\beta }}_{2}$	0.6	0.08	3.42	3.12 (0.38)		0.01	2.32	2.22 (0.18)		0.04	1.55	1.57 (0.09)
${\mathit{\beta }}_{3}$	0.4	0.09	3.44	3.13 (0.35)		−0.03	2.43	2.25 (0.18)		0.06	1.53	1.58 (0.09)
${\mathit{\lambda }}_{0}$	−0.8	−3.88	8.64	7.73 (0.22)		−2.21	6.42	5.50 (0.09)		−1.07	4.22	3.89 (0.05)
${\mathit{\lambda }}_{1}$	0.8	0.50	8.91	7.98 (0.51)		0.27	6.06	5.71 (0.26)		0.41	4.23	4.05 (0.13)
${\mathit{\lambda }}_{2}$	−0.5	−0.59	8.66	7.99 (0.52)		−0.10	5.85	5.71 (0.26)		−0.09	4.22	4.04 (0.12)
${\mathit{\lambda }}_{3}$	0.25	0.15	8.81	8.00 (0.53)		−0.10	6.27	5.72 (0.26)		−0.14	4.27	4.05 (0.12)
${\mathit{\gamma }}_{0}$	−0.5	−2.55	9.41	8.96 (0.32)		−0.83	6.63	6.35 (0.16)		−0.50	4.63	4.48 (0.08)
${\mathit{\gamma }}_{1}$	0.3	2.61	58.08	52.66 (2.37)		0.32	39.65	37.27 (1.18)		−0.72	26.29	26.27 (0.56)
${\mathit{\gamma }}_{2}$	−0.5	−6.35	72.12	64.96 (3.87)		−1.71	49.57	45.92 (1.90)		−0.22	32.53	32.30 (0.93)

 | Show Table

DownLoad: CSV

Study 2. We compare the performances of the proposed AMCD with those of MCD and ACD based on different datasets generated from various types of covariation structures in terms of the estimation robustness of the correlation matrices against model misspecification of the innovation variances. To assess the estimation accuracy of the correlation matrices, we define the entropy loss function ${\Delta }_{\text{1}}\left({\boldsymbol{R}}_{i},{\widehat{\boldsymbol{R}}}_{i}\right)=\text{trace} \left({\boldsymbol{R}}_{i}^{-1}{\widehat{\boldsymbol{R}}}_{i}\right)-\text{log}\left|{\boldsymbol{R}}_{i}^{-1}\right|-{m}_{i}$ and the quadratic loss function ${\Delta}_{\text{2}}\left({\boldsymbol{R}}_{i},{\widehat{\boldsymbol{R}}}_{i}\right)= \text{trace} {\left({\boldsymbol{R}}_{i}^{-1}{\widehat{\boldsymbol{R}}}_{i}-{\boldsymbol{I}}_{{m}_{i}}\right)}^{2}$, where ${R}_{i}$ is the true correlation matrix and ${\widehat{\boldsymbol{R}}}_{i}$ is an estimated correlation matrix. Both of the loss functions are 0 when ${\boldsymbol{R}}_{i}={\widehat{\boldsymbol{R}}}_{i}$ and positive when ${\boldsymbol{R}}_{i}\ne {\widehat{\boldsymbol{R}}}_{i}$. It is obvious that a smaller loss implies better estimates of ${\boldsymbol{R}}_{i}$. Furthermore, the loss functions for all subjects are defined by ${L}_{l}=\dfrac{1}{n}\displaystyle \sum _{i=1}^{n}{\Delta}_{l}\left({\boldsymbol{R}}_{i},{\widehat{\boldsymbol{R}}}_{i}\right),i=\mathrm{1,2}$.

In the following, we generated 1000 replications of the datasets with $n=100$ from the joint mean-covariance model similar to that in Study 1, except that ${r}_{ij}$ follows respectively the proposed AMCD, MCD and ACD. Then the estimates and losses of the correlation matrices were computed respectively based on the true model $\mathrm{log}\left({d}_{ij}^{2}\right)={\mathrm{\lambda }}_{0}+{z}_{ij1}{\mathrm{\lambda }}_{1}+{z}_{ij2}{\mathrm{\lambda }}_{2}+{z}_{ij3}{\mathrm{\lambda }}_{3}$ and the misspecified model $\mathrm{log}\left({d}_{ij}^{2}\right)={\mathrm{\lambda }}_{0}+{z}_{ij1}{\mathrm{\lambda }}_{1}$ for the innovation variances. For comparison, the losses obtained from fitting the true and misspecified models for the innovation variances are exhibited in Tables 2 and 3, respectively, where ENL and QUL represent respectively the average of the entropy loss ${L}_{1}$ and quadratic loss ${L}_{2}$, with empirical standard errors in parentheses.

Table  2.  Simulation results for Study 2. Fitting the true model for log(${d}_{ij}^{2}$).

True structure	AMCD			MCD			ACD
	ENL	QUL		ENL	QUL		ENL	QUL
AMCD	0.0342	0.1239		0.3999	1.6446		0.1611	0.4953
	(0.0287)	(0.1415)		(0.0949)	(0.7985)		(0.0355)	(0.3158)
MCD	0.5481	2.7191		0.0432	0.1413		0.7785	4.9234
	(0.0868)	(1.3806)		(0.0307)	(0.1545)		(0.1191)	(2.2544)
ACD	0.0807	0.1274		0.2208	0.4502		0.0359	0.0701
	(0.0300)	(0.1133)		(0.0444)	(0.2350)		(0.0299)	(0.0852)

 | Show Table

DownLoad: CSV

Table  3.  Simulation results for Study 2. Fitting ${\lambda }_{\text{0}}+{z}_{ij\text{1}}{\lambda }_{\text{1}}$for log(${d}_{ij}^{2}$).

True structure	AMCD			MCD			ACD
	ENL	QUL		ENL	QUL		ENL	QUL
AMCD	0.0364	0.1348		0.4551	1.9058		0.1638	0.5323
	(0.0308)	(0.1584)		(0.0974)	(0.8514)		(0.0375)	(0.3441)
MCD	0.5551	2.8642		0.1051	0.3242		0.7857	5.1367
	(0.0908)	(1.4454)		(0.0350)	(0.2610)		(0.1215)	(2.3220)
ACD	0.0819	0.1362		0.2411	0.4912		0.0364	0.0728
	(0.0310)	(0.1220)		(0.0448)	(0.2397)		(0.0301)	(0.0898)

 | Show Table

DownLoad: CSV

In terms of the average losses in Table 2, it is vital to know the true covariation structure; that is, when the true covariation structure follows AMCD, the losses of estimating the correlation matrices increase when the covariation structure is incorrectly decomposed based on MCD or ACD, and vice versa. By comparing Table 2 with Table 3, we can analyze the influence of model misspecification for the innovation variances on estimating the correlation matrices. Specifically, by comparing the left two AMCD columns of Table 2 with those of Table 3, ENL and QUL both vary little, no matter whether the true covariation structure is AMCD, MCD or ACD; the losses in the right two ACD columns of Tables 2 and 3 also exhibit a similar pattern; however, the losses in the middle two MCD columns of Tables 2 and 3 vary significantly. Thus, the estimation of the correlation matrices is robust with respect to the model misspecification for the innovation variances when fitting the AMCD or ACD structure, but is nevertheless not robust when fitting the MCD structure.

Study 3. We generated 1000 replications of the datasets with $n=100$ respectively from several distributions, where the joint model was the same as that in Study 1. Specifically, ${y}_{i}$ follows respectively the multivariate normal distribution (Scenario 1), Laplace distribution (Scenario 2) and PE distribution with $\upsilon =0.7$ (Scenario 3), and we fitted NJMM and LJMM under each scenario to compare thier modeling capacities.

To assess the estimation accuracy of the parameters, Bias, SD and mean squared errors (MSE) were calculated, and the simulation results are displayed in Tables 4–6. From Table 4, NJMM yields smaller MSE than does LJMM for all parameters, when the data follows the multivariate normal distribution (Scenario 1). Moreover, Bias of NJMM and LJMM are generally small, except that Bias for ${\boldsymbol{\lambda }}_{0}$ when fitting LJMM are rather large under Scenario 1, which might be due to the influence of distribution misspecification and the difficulty of estimating the constant term in innovation variance models. From Table 5, we can carry out a similar discussion for Scenario 2, and find out that LJMM outperforms NJMM when the data follows the multivariate Laplace distribution. Table 6 shows that LJMM has slightly smaller MSE than NJMM, when the data follows the multivariate PE distribution with $\upsilon =0.7$ (Scenario 3), which is between the multivariate normal distribution and the multivariate Laplace distribution.

Table  4.  Simulation results for Study 3. Fitting NJMM and LJMM based on the data generated from the multivariate normal distribution (Scenario 1).

		NJMM				LJMM
	True	Bias	SD	MSE		Bias	SD	MSE
${\mathit{\beta }}_{0}$	1	0.000293	0.013820	0.000191		0.000389	0.014523	0.000211
${\mathit{\beta }}_{1}$	−0.6	0.000956	0.023084	0.000534		0.001186	0.023680	0.000562
${\mathit{\beta }}_{2}$	0.6	−0.001044	0.022847	0.000523		−0.001132	0.023518	0.000554
${\mathit{\beta }}_{3}$	0.4	−0.000942	0.024125	0.000583		−0.001063	0.025185	0.000635
${\mathit{\lambda }}_{0}$	−0.8	−0.018044	0.057580	0.003641		−3.284723	0.059723	10.792970
${\mathit{\lambda }}_{1}$	0.8	0.002058	0.061634	0.003803		0.002269	0.064317	0.004142
${\mathit{\lambda }}_{2}$	−0.5	−0.001617	0.060490	0.003662		−0.002075	0.064486	0.004163
${\mathit{\lambda }}_{3}$	0.25	0.002351	0.057719	0.003337		0.002540	0.061003	0.003728
${\mathit{\gamma }}_{0}$	−0.5	−0.010487	0.063410	0.004131		−0.010345	0.064888	0.004318
${\mathit{\gamma }}_{1}$	0.3	−0.000795	0.379933	0.144350		0.000034	0.387449	0.150117
${\mathit{\gamma }}_{2}$	−0.5	−0.015246	0.473920	0.224833		−0.018321	0.486507	0.237025

 | Show Table

DownLoad: CSV

Table  6.  Simulation results for Study 3. Fitting NJMM and LJMM based on the data generated from the multivariate PE distribution with $\upsilon$=0.7 (Scenario 3).

		NJMM				LJMM
	True	Bias	SD	MSE		Bias	SD	MSE
${\mathit{\beta }}_{0}$	1	0.000279	0.031374	0.000984		−0.000075	0.031226	0.000975
${\mathit{\beta }}_{1}$	−0.6	−0.000363	0.053292	0.002840		−0.001130	0.052419	0.002749
${\mathit{\beta }}_{2}$	0.6	0.000324	0.053262	0.002837		0.001559	0.053214	0.002834
${\mathit{\beta }}_{3}$	0.4	−0.000064	0.054821	0.003005		−0.000298	0.053428	0.002855
${\mathit{\lambda }}_{0}$	−0.8	1.697678	0.076425	2.887951		−1.610802	0.075379	2.600366
${\mathit{\lambda }}_{1}$	0.8	0.005810	0.065800	0.004363		0.006267	0.063930	0.004126
${\mathit{\lambda }}_{2}$	−0.5	−0.001172	0.067431	0.004548		−0.001043	0.066281	0.004394
${\mathit{\lambda }}_{3}$	0.25	0.000540	0.064320	0.004137		0.000625	0.063059	0.003977
${\mathit{\gamma }}_{0}$	−0.5	−0.007865	0.068495	0.004753		−0.007884	0.068537	0.004759
${\mathit{\gamma }}_{1}$	0.3	−0.000655	0.422247	0.178293		0.000933	0.421153	0.177371
${\mathit{\gamma }}_{2}$	−0.5	−0.019594	0.521888	0.272751		−0.022436	0.515721	0.266472

 | Show Table

DownLoad: CSV

Table  5.  Simulation results for Study 3. Fitting NJMM and LJMM based on the data generated from the multivariate Laplace distribution (Scenario 2).

		NJMM				LJMM
	True	Bias	SD	MSE		Bias	SD	MSE
${\mathit{\beta }}_{0}$	1	0.000443	0.074431	0.005540		−0.001020	0.070935	0.005033
${\mathit{\beta }}_{1}$	−0.6	0.001874	0.125683	0.015800		−0.000377	0.115993	0.013454
${\mathit{\beta }}_{2}$	0.6	−0.006187	0.120606	0.014584		−0.004163	0.112968	0.012779
${\mathit{\beta }}_{3}$	0.4	0.003073	0.124357	0.015474		0.004704	0.116733	0.013649
${\mathit{\lambda }}_{0}$	−0.8	3.325543	0.090872	11.067496		−0.019765	0.085474	0.007696
${\mathit{\lambda }}_{1}$	0.8	0.000266	0.068353	0.004672		0.000632	0.063271	0.004004
${\mathit{\lambda }}_{2}$	−0.5	−0.000353	0.068978	0.004758		−0.000868	0.065576	0.004301
${\mathit{\lambda }}_{3}$	0.25	0.003889	0.067912	0.004627		0.003261	0.063900	0.004094
${\mathit{\gamma }}_{0}$	−0.5	−0.009415	0.072748	0.005381		−0.007799	0.069382	0.004875
${\mathit{\gamma }}_{1}$	0.3	−0.012116	0.428129	0.183441		−0.016835	0.406756	0.165733
${\mathit{\gamma }}_{2}$	−0.5	0.003418	0.527299	0.278056		0.011848	0.507051	0.257241

 | Show Table

DownLoad: CSV

6.   Real data analysis

6.1   The cattle data

We applied the proposed AMCD approach to the balanced cattle data, which was initially introduced in Ref. [19]. This dataset consists of two treatment groups, A and B, and the weights of each cattle were measured 11 times over a period of 133 d. Notably, the 30 animals in group A were analyzed in Refs. [11, 20] and [8]. Pan and Mackenzie^[20] found out that it is reasonable to adopt three polynomials for modeling jointly the mean, the log-innovation variances and the autoregressive coefficients based on MCD, and the optimal triplet of the polynomial orders is determined as (8,3,4) in terms of the Bayesian information criterion (BIC).

From expression (8) based on AMCD, we can calculate the corresponding innovation variances ${d}_{ij}^{2}$ and generalized autoregressive parameters ${\phi }_{ijk}\left(j > k\right)$ of the sample covariance matrix. Fig. 1 displays $\text{log}\left({d}_{ij}^{2}\right)$ versus time and ${\phi }_{ijk}$ versus time lag, both indicate trends which might be properly characterized by polynomials. Then we adopt three polynomials

Figure  1.  Cattle data: sample regressograms and fitted curves for (a) log-innovation variances and (b) autoregressive coefficients. Solid lines, curves fitted by the proposed AMCD method; dashed lines, 95% pointwise confidence intervals using the bootstrapping method.

DownLoad: Full-Size Img PowerPoint

for modelling jointly the mean and the covariation structure based on the proposed AMCD, and select the optimal model in terms of

where $p$, $s$ and $q$ are the three polynomial orders, and ${\widehat{l}}_{\text{max}}$ is the corresponding maximum log-likelihood. Then some model with the smallest BIC is generally judged to be the optimal model, which can capture the dynamics much more precisely and parsimoniously. Due to the asymptotic orthogonality of the mean parameter and the covariation parameters, these two kinds of parameters can be searched separately, which is more computationally efficient. Note that $p$ within the optimal model in Refs. [20] and [8] are both 8, which motivates us to first fix $p=\text{8}$ and search for $s$ and $q$. The optimal $\left(s,q\right)$ is determined to be (3,3) by comparing several candidate models, and the major searching process is displayed in Table 7. Then, we fix $\left(s,q\right)=\left(\text{3},\text{3}\right)$, search for $p$, and determine the optimal order for AMCD as $\left(p,s,q\right)=\left(\text{8,3,3}\right)$. Fig. 1 shows the optimal AMCD-based fitted curves for $\text{log}\left({d}_{ij}^{2}\right)$’s and ${\phi }_{ijk}$’s and their corresponding $\text{95} \%$ pointwise confidence intervals obtained using the bootstrapping method; these curves capture the pattern well for the log-innovation variances and rather well for the generalized autoregressive parameters.

Table  7.  Cattle data: comparison of various models based on the proposed AMCD approach.

$\text{Poly}\left(\text{8}\text{,}\text{}{s}\text{,}\text{}{q}\right)$	Number of parameters	${\widehat{l}}_{\text{max}}$	BIC
(8,1,1)	13	−808.07	55.35
(8,1,2)	14	−764.56	52.56
(8,1,3)	15	−750.13	51.71
(8,1,4)	16	−749.86	51.80
(8,2,1)	14	−806.79	55.37
(8,2,2)	15	−763.07	52.57
(8,2,3)	16	−749.16	51.76
(8,2,4)	17	−749.00	51.86
(8,3,1)	15	−806.52	55.47
(8,3,2)	16	−760.96	52.54
(8,3,3)†	17	−746.07	51.67
(8,3,4)	18	−745.82	51.76
(8,4,1)	16	−806.45	55.58
(8,4,2)	17	−759.88	52.59
(8,4,3)	18	−746.06	51.78
(8,4,4)	19	−745.80	51.87

 | Show Table

DownLoad: CSV

The optimal order based on AMCD is slightly less than that based on MCD, but identical to that based on ACD, which is also determined by using a similar searching process. To compare the performances of the aforementioned three decompositions, we regard the sample correlation matrix as the true matrix, and still assess the estimation accuracy of the correlation matrix by the entropy loss ${L}_{1}$ and quadratic loss ${L}_{2}$. The $\left({L}_{1},{L}_{2}\right)$ values for the optimal models based on AMCD, ACD and MCD are (3.0866,35.6300), (2.8660,30.8312), and (1.9237,12.9392), respectively. Obviously, the MCD-based model yields more accurate estimates of the correlation matrix, and the AMCD and ACD-based models have similar performances. This indicates that the MCD structure might be more suitable for this dataset than AMCD and ACD, in terms of the estimation accuracy of the correlation matrix. However, when we keep $p$ and $q$ fixed, and take the polynomial order for the innovation variances as $s=\text{1}$ rather than the optimal order $s=\text{3}$, the $\left({L}_{1},{L}_{2}\right)$’s are respectively (3.0743, 31.8634), (2.9381, 30.0376), and (2.8819, 29.2350). The estimation accuracy of the correlation matrix for AMCD and ACD-based models varies slightly, again implying robustness with respect to the model misspecification of the innovation variances; while the losses for the MCD-based model increase significantly which might lead to the opposite conclusion.

6.2   The sleep dose-response data

We apply NJMM and LJMM based on the proposed AMCD to the sleep dose-response data, which was initially investigated by Ref. [21] and subsequently analyzed in Ref. [22] using Bayesian inference for the joint model based on the t distribution and MCD. The days of sleep deprivation and the corresponding average reaction times were recorded for each of the 18 participants in the 3-h group. Fig. 2a shows the trajectories of the average reaction times over an equally spaced 10-d period, together with the mean profile plot and the corresponding $\pm 1$ standard deviations across the period. There seem to be sudden jumps and drops in the trajectories of the first and sixth subjects, which may be obvious outliers and thus not able to be properly handled using NJMM.

Figure  2.  Sleep dose-response data: (a) trajectories of average reaction time; (b) sample regressograms for log-innovation variances; (c) sample regressograms for autoregressive coefficients.

DownLoad: Full-Size Img PowerPoint

The mean response varies linearly over time, and thus can be modeled by a 1st-degree polynomial in time. Fig. 2b and c show the sample regressograms for log-innovation variances and autoregressive coefficients, implying higher-order polynomial functions in time or time lag. Hence, it is suitable to adopt three polynomials expressed as (16) with $p=1$ and $s,q$ to be determined. We fitted the sleep dose-response data using AMCD-based NJMM and LJMM for various degrees of the Poly($1,s,q$) models. The numbers of parameters, ${\widehat{l}}_{\text{max}}$ and BIC($1,s,q$) values for the searched Poly($1,s,q$)’s are listed in Table 8, where the BIC values are calculated via formula (17). According to the BIC values, Poly(1,3,4) and Poly(1,3,2) are the best for NJMM and LJMM, respectively, and this data is fitted significantly better when using LJMM. Table 9 displays the maximum likelihood estimates and standard errors when fitting respectively NJMM and LJMM with the above two best polynomials. Note that NJMM and LJMM have similar parameter estimates. However, the standard errors of $\widehat{\boldsymbol{\lambda } }$ and $\widehat{\boldsymbol{\gamma } }$ for LJMM are smaller than those of NJMM with Poly(1,3,2); the standard errors of $\widehat{\boldsymbol{\lambda } }$ for LJMM are slightly larger than those of NJMM with Poly(1,3,4), and the standard errors of $\widehat{\boldsymbol{\gamma }}$ for LJMM are still smaller in this case. This indicates convincingly that parameter estimates of LJMM always possess lower variability for this data.

Table  8.  Sleep dose-response data: comparison of various Poly($1,{s,q}$) choices between NJMM and LJMM.

$\text{Poly}\left({1},{s,q}\right)$	Number of parameters	${\widehat{l}}_{\text{max}}$			BIC
		NJMM	LJMM		NJMM	LJMM
(1,1,1)	6	−710.32	−454.85		79.89	51.50
(1,1,2)	7	−698.82	−441.85		78.77	50.22
(1,1,3)	8	−698.07	−441.25		78.85	50.31
(1,1,4)	9	−696.49	−440.01		78.83	50.34
(1,1,5)	10	−695.85	−439.72		78.92	50.46
(1,2,1)	7	−710.13	−454.47		80.03	51.62
(1,2,2)	8	−698.66	−441.54		78.91	50.34
(1,2,3)	9	−697.86	−440.87		78.99	50.43
(1,2,4)	10	−696.40	−439.79		78.98	50.47
(1,2,5)	11	−695.78	−439.52		79.08	50.60
(1,3,1)	8	−709.45	−454.21		80.11	51.75
(1,3,2)	9	−695.00	−438.03		78.67	50.11
(1,3,3)	10	−694.00	−437.36		78.72	50.20
(1,3,4)	11	−691.02	−435.40		78.55	50.14
(1,3,5)	12	−690.22	−435.04		78.62	50.27
(1,4,1)	9	−708.21	−453.47		80.13	51.83
(1,4,2)	10	−694.34	−437.76		78.76	50.25
(1,4,3)	11	−693.20	−436.96		78.79	50.32
(1,4,4)	12	−690.52	−435.11		78.65	50.27
(1,4,5)	13	−689.79	−434.79		78.73	50.40
(1,5,1)	10	−707.14	−452.99		80.18	51.94
(1,5,2)	11	−693.42	−437.21		78.81	50.35
(1,5,3)	12	−692.06	−436.24		78.82	50.40
(1,5,4)	13	−689.65	−434.64		78.72	50.38
(1,5,5)	14	−688.83	−434.26		78.78	50.50

 | Show Table

DownLoad: CSV

Table  9.  Sleep dose-response data: parameter estimates for NJMM and LJMM with the two best polynomials.

	Poly(1,3,2)						Poly(1,3,4)
	NJMM			LJMM			NJMM			LJMM
	MLE	SE		MLE	SE		MLE	SE		MLE	SE
${\beta }_{0}$	240.9311	6.7036		241.1218	7.7784		241.6758	6.3496		241.1834	7.7073
${\beta }_{1}$	10.1691	1.6517		9.7037	1.6478		10.2647	1.6540		9.8365	1.6582
${\lambda }_{0}$	7.5886	0.6562		4.1020	0.6333		7.8139	0.5422		4.1793	0.5746
${\lambda }_{1}$	−1.0052	0.4517		−1.0521	0.4055		−1.2803	0.3616		−1.1990	0.3734
${\lambda }_{2}$	0.2199	0.0875		0.2107	0.0783		0.2827	0.0753		0.2470	0.0767
${\lambda }_{3}$	−0.0126	0.0050		−0.0118	0.0045		−0.0165	0.0044		−0.0141	0.0045
${\gamma }_{0}$	0.8835	0.1715		0.9978	0.1456		2.1703	0.5384		2.0808	0.4754
${\gamma }_{1}$	−0.3654	0.0978		−0.4222	0.0795		−2.0456	0.7714		−1.8478	0.6405
${\gamma }_{2}$	0.0349	0.0112		0.0407	0.0090		0.6838	0.3205		0.5931	0.2577
${\gamma }_{3}$	−	−		−	−		−0.0956	0.0495		−0.0814	0.0393
${\gamma }_{4}$	−	−		−	−		0.0047	0.0025		0.0040	0.0020

 | Show Table

DownLoad: CSV

7.   Conclusions

In longitudinal data analysis, we always suffer from various kinds of covariation structure and need to choose the most suitable decomposition among the candidates. If the covariance matrix possesses a typical structure, we might apply the MACD method; otherwise, we might adopt the MCD method for the precision matrix. Furthermore, when we are interested in robust estimation of the correlation matrix against model misspecification of the innovation variances, we can apply the ACD method to the covariance matrix. However, the decomposition of the precision matrix targeting the above robustness has been rather less investigated. In this paper, we have proposed AMCD of the precision matrix and established its role in providing robust estimator for the correlation matrix. In addition, we have investigated the AMCD-based LJMM which may achieve both the aforementioned robust estimation and robustness to outliers in the data.

Recently, an autoregressive moving average Cholesky decomposition (ARMACD) has been proposed in Ref. [23] ; this decomposition is more general than both MACD and MCD. Therefore, it is worthwhile to develop a robust estimation of the correlation matrix on the basis of ARMACD in the future. Furthermore, these decompositions correspond to various time series models, then various robust time series approaches could also be generalized for robust estimation of the correlation, covariance or precision matrix. In addition, longitudinal data may suffer from missingness, such as informative dropouts^[24]. Then we need to specify a suitable dropout mechanism, establish the complete data log-likelihood, and implement an EM algorithm to calculate the maximum likelihood estimates. Moreover, heterogeneity is also deserved to be investigated, not only for the mean parameters, but also for the covariation structures. The finite mixture model^[25] based on NJMM or LJMM with AMCD may perform well, which is left for future research.