ESL-Based Robust Estimation for Mean-Covariance Regression with Longitudinal Data

Fei Lu^*, Liugen Xue, Xiong Cai

College of Applied Sciences, Beijing University of Technology, Beijing, China.

DOI: 10.4236/ojs.2020.101002                

Abstract

When longitudinal data contains outliers, the classical least-squares approach is known to be not robust. To solve this issue, the exponential squared loss (ESL) function with a tuning parameter has been investigated for longitudinal data. However, to our knowledge, there is no paper to investigate the robust estimation procedure against outliers within the framework of mean-covariance regression analysis for longitudinal data using the ESL function. In this paper, we propose a robust estimation approach for the model parameters of the mean and generalized autoregressive parameters with longitudinal data based on the ESL function. The proposed estimators can be shown to be asymptotically normal under certain conditions. Moreover, we develop an iteratively reweighted least squares (IRLS) algorithm to calculate the parameter estimates, and the balance between the robustness and efficiency can be achieved by choosing appropriate data adaptive tuning parameters. Simulation studies and real data analysis are carried out to illustrate the finite sample performance of the proposed approach.

Keywords

Exponential Squared Loss Function, Within-Subject Correlation, Longitudinal Data, Modified Cholesky Decomposition, Robustness

1. Introduction

Longitudinal data arises frequently in many fields, such as biological research, social science and other fields. The observations on the same subject are measured repeatedly over time, and thus intrinsically correlated [1]. The covariance matrix of such data is important since ignoring the correlation structure may lead to inefficient estimators of the mean parameter. Furthermore, the covariance matrix itself may be of scientific interest [2]. However, it is challenging to model the covariance matrix which suffers from the positive-definiteness constraint and includes many unknown parameters. To avoid this challenge, a common strategy is to specify the working correlation structure [3], which does not permit more general structures and cannot flexibly incorporate covariates that may be helpful to explain the covariations. To overcome this limitation, joint modelling for the mean and covariance of longitudinal data has received increasing interest recently; see, for example, [4] - [12]. Among these joint modeling approaches, the modified Cholesky decomposition (MCD) for the covariance matrix proposed in [4] [5] is attractive owing to the fact that it leads automatically to positive definite covariance matrix, and the parameters in it can be interpreted by suitable statistical concepts. As a result, the regression techniques and model based inference can be adopted to the parameters in this decomposition, see [7] for more details.

However, the aforementioned approaches are very sensitive to outliers or heavy-tailed distributions. [13] proposed a robust procedure for modeling the correlation matrix of longitudinal data based on an alternative Cholesky decomposition and heavy-tailed multivariate t-distributions with unknown degrees of freedom. It should be pointed out that the use of the multivariate t-distribution alone does not necessarily guarantee robustness. In addition, [14] [15] developed robust generalized estimating equations (GEE) for regression parameters in joint mean-covariance model by employing the Huber’s function and leverage-based weights. [16] developed an efficient parameter estimation via MCD for quantile regression with longitudinal data. [17] further proposed a moving average Cholesky factor model, which is transformed from MCD, in covariance modeling for composite quantile regression with longitudinal data. Then, [18] carried out smoothed empirical likelihood inference via MCD for quantile varying coefficient models with longitudinal data. Later, [19] developed quantile estimation via MCD for longitudinal single-index models.

Although M-type regression and quantile regression procedures can overcome outliers and heavy-tail errors, they may lose efficiency under the normal distribution. To overcome this difficulty, [20] recently proposed a robust variable selection approach by adopting the exponential squared loss (ESL) function with a tuning parameter. They have showed that, with properly selected tuning parameter, the proposed approach not only achieves good robustness with respect to outliers in the dataset, but also is as asymptotically efficient as the least squares estimation without outliers under normal error. Later, some authors employed the ESL funtion for longitudinal data. For example, [21] propsed an efficient and robust variable selection method for longitudinal generalized linear models based on GEE. [22] proposed a robust and efficient estimation procedure for simultaneous model structure identification and variable selection in generalized partial linear varying coefficient models for longitudinal data. [23] developed GEE-based robust estimation and empirical likelihood inference approach with ESL for panel data models. With a similar loss function, [24] proposed a robust variable selection method in modal varying-coefficient models with longitudinal data. [25] proposed modal regression statistical inference for longitudinal semivarying coefficient models, including GEE, empirical likelihood and variable selection. However, to our knowledge, there is no paper to investigate the robust estimation procedure against outliers within the framework of mean-covariance regression analysis for longitudinal data employing the ESL function.

In this paper, we propose a robust estimation approach for the model parameters of the mean and generalized autoregressive parameters in the within subject covariance matrices for longitudinal data based on the ESL function. We begin with the ESL-based estimation for the mean parameters pretending that the repeated measurements within a subject are independent. Then based on the roughly estimated mean parameters, the simultaneous estimation for the mean and generalized autoregressive parameters is carried out using the ESL function. The proposed estimators can be shown to be asymptotically normal under certain conditions. Moreover, we develop an iteratively reweighted least squares (IRLS) algorithm [26] to calculate the parameter estimates. Numerical studies are carried out to illustrate the finite sample performance of our approach.

The outline of this paper is organized as follows. Section 2 develops the robust estimation procedure. Section 3 establishes the asymptotic properties of the proposed ESL estimators. The IRLS algorithm and a data-driven method for the selection of tuning parameters are presented in Section 4. Simulations are carried out in Section 5. Section 6 analyses a real data set. A discussion is given in Section 7. The technical proofs are provided in the Appendix.

2. Robust Estimation Procedure

2.1. Initial Estimate for the Mean Parameters

Consider n subjects, where each subject is measured repeatedly over time. For the ith subject, suppose that $Y_{ij}$ is the observed scale response variable at time $t_{ij}$, and $X_{ij}$ is the corresponding $p\times 1$ covariate vector, $i=1,\cdots ,n,j=1,\cdots ,m_i$. Denote $N={\displaystyle \underset{i=1}{\overset{n}{\sum }}{m}_i}$. Furthermore, let $Y_i={\left(Y_{i1}\mathrm{,}\cdots \mathrm{,}{Y}_{i{m}_i}\right)}^{{\rm T}}$, $X_i={\left(X_{i1}\mathrm{,}\cdots \mathrm{,}{X}_{i{m}_i}\right)}^{{\rm T}}$, $t_i={\left(t_{i1}\mathrm{,}\cdots \mathrm{,}{t}_{i{m}_i}\right)}^{{\rm T}}$. A longitudinal linear regression model has the form

$Y_{ij}=X_{ij}^{{\rm T}}\beta +{\epsilon }_{ij}\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}i=\mathrm{1,}\cdots \mathrm{,}n,j=\mathrm{1,}\cdots \mathrm{,}{m}_i\mathrm{,}$ (1)

where $\beta$ is the $p\times 1$ vector of associated parameters, ${\epsilon }_i={\left({\epsilon }_{i1}\mathrm{,}\cdots \mathrm{,}{\epsilon }_{i{m}_i}\right)}^{\text{T}}$ is the random error satisfying that $E\left({\epsilon }_{ij}\mathrm{<mo>\; |\; </mo>}{X}_{ij}\right)=0$ and $\mathrm{var}\left({\epsilon }_i|X_i\right)={\Sigma }_i$.

Define the ESL function $1-{\varphi }_{\tau }\left(t\right)=1-\exp \left(-t^2/\tau \right)$. We first estimate $\beta$ pretending that the random error ${\epsilon }_{ij}$ ’s are independent. More specifically, we estimate $\beta$ by maximizing the objective function

$L_1\left(\beta \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta \right)\mathrm{,}$ (2)

where ${\tau }_1$ is a tuning parameter. The resulting initial estimate of $\beta$ is denoted by $\tilde{\beta }$.

2.2. Simultaneous Estimate for the Mean and Generalized Autoregressive Parameters

Based on the Cholesky decomposition, there exists a lower triangle matrix ${\Phi }_i$ with diagonal ones such that

$\mathrm{var}\left({\Phi }_i{\epsilon }_i\mathrm{<mo>\; |\; </mo>}{X}_i\right)={\Phi }_i{\Sigma }_i{\Phi }_i^{{\rm T}}=D_i\mathrm{,}$

where $D_i=\text{diag}\left(d_{i1}^2\mathrm{,}\cdots \mathrm{,}{d}_{i{m}_i}^2\right)$. In other words, let ${ϵ}_i={\left({ϵ}_{i1}\mathrm{,}\cdots \mathrm{,}{ϵ}_{i{m}_i}\right)}^{\text{T}}={\Phi }_i{\epsilon }_i$, we obtain that

$\begin{array}{l}{\epsilon }_{i1}={ϵ}_{i1}\mathrm{,}\\ {\epsilon }_{ij}={\varphi }_{j\mathrm{,1}}^{\left(i\right)}{\epsilon }_{i1}+\cdots +{\varphi }_{j\mathrm{,}j-1}^{\left(i\right)}{\epsilon }_{i\mathrm{,}j-1}+{ϵ}_{ij}\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n,j=\mathrm{2,}\cdots \mathrm{,}{m}_i\mathrm{,}\end{array}$ (3)

where ${\varphi }_{j\mathrm{,}k}^{\left(i\right)}$ is the negative of the $\left(j\mathrm{,}k\right)\text{th}$ element of ${\Phi }_i$. It’s obvious that ${ϵ}_{ij}$ ’s are uncorrelated with $E\left({ϵ}_{ij}\right)=0$ and $\mathrm{var}\left({ϵ}_{ij}\right)=d_{ij}^2$, $j=1,\cdots ,m_i$. If ${\epsilon }_i$ ’s were available, then (1) could be transformed as the following linear model with uncorrelated error ${ϵ}_{ij}$ ’s:

$\begin{array}{l}{Y}_{i1}=X_{i1}^{{\rm T}}\beta +{ϵ}_{i1}\mathrm{,}\\ Y_{ij}=X_{ij}^{{\rm T}}\beta +{\varphi }_{j\mathrm{,1}}^{\left(i\right)}{\epsilon }_{i1}+\cdots +{\varphi }_{j\mathrm{,}j-1}^{\left(i\right)}{\epsilon }_{i\mathrm{,}j-1}+{ϵ}_{ij}\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n\mathrm{,\; <mi>\; </mi>}j=\mathrm{2,}\cdots \mathrm{,}{m}_i\mathrm{.}\end{array}$ (4)

[4] pointed out that the MCD has a well-founded statistical interpretation, and it has the advantage that the generalized autoregressive parameter ${\varphi }_{j\mathrm{,}k}^{\left(i\right)}$ ’s and log-innovation variance $\log \left(d_{ij}^2\right)$ ’s are unconstrained. For simplicity, we assume that $d_{ij}^2=d^2$ for all $i=1,\cdots ,n,\mathrm{<mi>\; </mi>}j=2,\cdots ,m_i$. Since ${\Sigma }_i=\mathrm{var}\left({\epsilon }_i|X_i\right)$ may depend on $X_i$, we adopt a more parsimonious structure,

${\varphi }_{j\mathrm{,}k}^{\left(i\right)}=Z_{j\mathrm{,}k}^{\left(i\right){\rm T}}\gamma \mathrm{,}$ (5)

where $\gamma ={\left({\gamma }_1\mathrm{,}\cdots \mathrm{,}{\gamma }_q\right)}^{{\rm T}}$ is the regression coefficient, and the covariates $Z_{j\mathrm{,}k}^{\left(i\right)}={\left(Z_{j\mathrm{,}k\mathrm{,1}}^{\left(i\right)}\mathrm{,}\cdots \mathrm{,}{Z}_{j\mathrm{,}k\mathrm{,}q}^{\left(i\right)}\right)}^{{\rm T}}$ may contain the time, the baseline covariates $X_i$, the interactions and so on.

Based on $\tilde{\beta }$, we can obtain the estimated residuals

${\tilde{\epsilon }}_{ij}=Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n\mathrm{,\; <mi>\; </mi>}j=\mathrm{1,}\cdots \mathrm{,}{m}_i-1.$ (6)

From (4), (5) and (6), we can obtain the simultaneous estimate $\hat{\theta }={\left({\hat{\beta }}^{{\rm T}}\mathrm{,}{\hat{\gamma }}^{{\rm T}}\right)}^{{\rm T}}$ of $\theta ={\left({\beta }^{{\rm T}}\mathrm{,}{\gamma }^{{\rm T}}\right)}^{{\rm T}}$ by maximizing the following objective function:

$\begin{array}{c}{L}_2\left(\theta \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta -{\tilde{\epsilon }}_{i1}{Z}_{j\mathrm{,1}}^{\left(i\right){\rm T}}\gamma -\cdots -{\tilde{\epsilon }}_{i\mathrm{,}j-1}{Z}_{j\mathrm{,}j-1}^{\left(i\right){\rm T}}\gamma \right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta -{\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\tilde{\epsilon }}_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,1}}^{\left(i\right)}{\gamma }_1-\cdots -{\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\tilde{\epsilon }}_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,}q}^{\left(i\right)}{\gamma }_q\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta -{\tilde{\zeta }}_{ij}^{{\rm T}}\gamma \right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-{\tilde{\delta }}_{ij}^{{\rm T}}\theta \right)\mathrm{,}\end{array}$ (7)

where ${\tilde{\delta }}_{ij}={\left(X_{ij}^{{\rm T}}\mathrm{,}{\tilde{\zeta }}_{ij}^{{\rm T}}\right)}^{{\rm T}}$ with ${\tilde{\zeta }}_{ij}={\left({\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\tilde{\epsilon }}_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,1}}^{\left(i\right)}\mathrm{,}\cdots \mathrm{,}{\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\tilde{\epsilon }}_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,}q}^{\left(i\right)}\right)}^{{\rm T}}$. Then we can obtain the estimates of ${\varphi }_{j\mathrm{,}k}^{\left(i\right)}$ ’s in model (3) by combining (5) and $\hat{\gamma }$.

3. Asymptotic Properties

Define $F_1\left(x,{\tau }_1\right)=E\left({{\varphi }^{″}}_{{\tau }_1}\left({\epsilon }_{ij}\right)|X_{ij}=x\right)$, $G_1\left(x,{\tau }_1\right)=E\left({{\varphi }^{\prime }}_{{\tau }_1}{\left({\epsilon }_{ij}\right)}^2|X_{ij}=x\right)$, $F_2\left(x,{\tau }_2\right)=E\left({{\varphi }^{″}}_{{\tau }_2}\left({ϵ}_{ij}\right)|X_{ij}=x\right)$, $G_2\left(x,{\tau }_2\right)=E\left({{\varphi }^{\prime }}_{{\tau }_2}{\left({ϵ}_{ij}\right)}^2|X_{ij}=x\right)$, $H\left({\bar{x}}_{ij}\right)=E\left({\zeta }_{ij}{\zeta }_{ij}^{{\rm T}}\mathrm{<mo>\; |\; </mo>}{\bar{X}}_{ij}={\bar{x}}_{ij}\right)$, ${\bar{X}}_{ij}={\left(X_{i1}\mathrm{,}\cdots \mathrm{,}{X}_{ij}\right)}^{{\rm T}}$, ${\bar{x}}_{ij}={\left(x_{i1}\mathrm{,}\cdots \mathrm{,}{x}_{ij}\right)}^{{\rm T}}$, ${\delta }_{ij}={\left(X_{ij}^{{\rm T}}\mathrm{,}{\zeta }_{ij}^{{\rm T}}\right)}^{{\rm T}}$, ${\zeta }_{ij}={\left({\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\epsilon }_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,1}}^{\left(i\right)}\mathrm{,}\cdots \mathrm{,}{\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\epsilon }_{i\mathrm{,}k}{Z}_{j\mathrm{,}k\mathrm{,}q}^{\left(i\right)}\right)}^{{\rm T}}$. To establish the asymptotic properties of the proposed ESL estimator, assume that the following regularity conditions hold:

(C1) There exists a positive integer M such that $\underset{1\le i\le n}{\max }{m}_i\le M<\infty$. This means that n and N have the same order.

(C2) There exists a positive constant C such that $\left|X_{ij}^{\left(k\right)}\right|\le C$ for $1\le i\le n$, $1\le j\le m_i$ and $0\le k\le p$. In addition, $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{X}_{ij}^{{\rm T}}$ converges to a finite positive definite matrix in probability.

(C3) $F_1\left(x\mathrm{,}{\tau }_1\right)$ and $G_1\left(x\mathrm{,}{\tau }_1\right)$ are continuous with respect to x. Moreover, for any ${\tau }_1>0$, $F_1\left(x,{\tau }_1\right)<0$.

(C4) $E\left({{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)|X_{ij}=x\right)=0$, $E\left({\left|{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)\right|}^3|X_{ij}=x\right)$, $E\left({{\varphi }^{″}}_{{\tau }_1}{\left({\epsilon }_{ij}\right)}^2|X_{ij}=x\right)$, $E\left({{\varphi }^{‴}}_{{\tau }_1}\left({\epsilon }_{ij}\right)|X_{ij}=x\right)$ and $E\left({{\varphi }^{‴}}_{{\tau }_1}{\left({\epsilon }_{ij}\right)}^2\mathrm{<mo>\; |\; </mo>}{X}_{ij}=x\right)$ are continuous with respect to x.

(C5) $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E\left(X_i^{{\rm T}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_i\right){{\varphi }^{\prime }}_{{\tau }_1}{\left({\epsilon }_i\right)}^{{\rm T}}{X}_i\right)$ converges to a finite positive definite matrix.

(C6) For $1\le i\le n\mathrm{,\; <mi>\; </mi>1}\le j_1\mathrm{,}\cdots \mathrm{,}{j}_6\le m_i-1$, $E\left\{\left|{\epsilon }_{i{j}_1}{\epsilon }_{i{j}_2}{\epsilon }_{i{j}_3}\right|\mathrm{<mo>\; |\; </mo>}\left(X_{i{j}_1}\mathrm{,}{X}_{i{j}_2}\mathrm{,}{X}_{i{j}_3}\right)=x\right\}$ and $E\left\{{\epsilon }_{i{j}_1}\cdots {\epsilon }_{i{j}_k}\mathrm{<mo>\; |\; </mo>}\left(X_{i{j}_1}\mathrm{,}\cdots \mathrm{,}{X}_{i{j}_k}\right)=x\right\}\left(k=\mathrm{2,3,4,6}\right)$ are continuous with respect to x.

(C7) $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\zeta }_{ij}{\zeta }_{ij}^{{\rm T}}$ converges to a finite positive definite matrix in probability.

(C8) $F_2\left(x\mathrm{,}{\tau }_2\right)$ and $G_2\left(x\mathrm{,}{\tau }_2\right)$ are continuous with respect to x. Moreover, for any ${\tau }_2>0$, $F_2\left(x,{\tau }_2\right)<0$.

(C9) $E\left({{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)|X_{ij}=x\right)=0$, $E\left({\left|{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)\right|}^3|X_{ij}=x\right)$, $E\left({{\varphi }^{″}}_{{\tau }_2}{\left({ϵ}_{ij}\right)}^2|X_{ij}=x\right)$, $E\left({{\varphi }^{‴}}_{{\tau }_2}\left({ϵ}_{ij}\right)|X_{ij}=x\right)$ and $E\left({{\varphi }^{‴}}_{{\tau }_2}{\left({ϵ}_{ij}\right)}^2|X_{ij}=x\right)$ are continuous with respect to x.

Theorem 1 If regularity conditions (C1)-(C5) hold, then

$\sqrt{n}\left(\tilde{\beta }-{\beta }_0\right)\to N\left(\mathrm{0,}{\Sigma }_1^{-1}{\Sigma }_2{\Sigma }_1^{-1}\right)$

in distribution as $n\to \infty$, where

${\Sigma }_1=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{F_1\left(X_{ij}\mathrm{,}{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}\mathrm{,}$

${\Sigma }_2=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E\left(X_i^{{\rm T}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_i\right){{\varphi }^{\prime }}_{{\tau }_1}{\left({\epsilon }_i\right)}^{{\rm T}}{X}_i\right)\mathrm{.}$

Then if the random error ${\epsilon }_{ij}$ ’s are independent, we have the following corollary, which is similar to Corollary 2 and useful for the choice of tuning parameters in Section 4.2.

Corollary 1 If regularity conditions (C1)-(C4) hold and ${\epsilon }_{ij}$ ’s are independent, then

$\sqrt{n}\left(\tilde{\beta }-{\beta }_0\right)\to N\left(\mathrm{0,}{\Sigma }_1^{-1}{\Sigma }_2^{\ast }{\Sigma }_1^{-1}\right)$

in distribution as $n\to \infty$, where

${\Sigma }_1=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{F_1\left(X_{ij}\mathrm{,}{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}\mathrm{,}$

${\Sigma }_2^{\ast }=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{G_1\left(X_{ij}\mathrm{,}{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}\mathrm{.}$

Theorem 2 If regularity conditions (C1)-(C9) hold, then

$\sqrt{n}\left(\hat{\theta }-{\theta }_0\right)\to N\left(\mathrm{0,}{\Sigma }_3^{-1}{\Sigma }_4{\Sigma }_3^{-1}\right)$

in distribution as $n\to \infty$, where

${\Sigma }_3=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{F_2\left(X_{ij}\mathrm{,}{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}\mathrm{<mo>\; |\; </mo>}{\bar{X}}_{ij}\right)\right\}\mathrm{,}$

${\Sigma }_4=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{G_2\left(X_{ij}\mathrm{,}{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}\mathrm{<mo>\; |\; </mo>}{\bar{X}}_{ij}\right)\right\}\mathrm{.}$

In fact, we have

$E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}\mathrm{<mo>\; |\; </mo>}{\bar{X}}_{ij}\right)=\left(\begin{array}{cc}{X}_{ij}{X}_{ij}^{{\rm T}}& 0\\ 0& H\left({\bar{X}}_{ij}\right)\end{array}\right)\mathrm{.}$ (8)

Then it can be deduced that $\hat{\beta }$ and $\hat{\gamma }$ are asymptotically independent, that is, the following corollaries hold.

Corollary 2 If regularity conditions (C1)-(C9) hold, then

$\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)\to N\left(\mathrm{0,}{\Sigma }_5^{-1}{\Sigma }_6{\Sigma }_5^{-1}\right)$

in distribution as $n\to \infty$, where

${\Sigma }_5=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{F_2\left(X_{ij}\mathrm{,}{\tau }_2\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}\mathrm{,}$

${\Sigma }_6=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{G_2\left(X_{ij}\mathrm{,}{\tau }_2\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}\mathrm{.}$

Corollary 3 If regularity conditions (C1)-(C9) hold, then

$\sqrt{n}\left(\hat{\gamma }-{\gamma }_0\right)\to N\left(\mathrm{0,}{\Sigma }_7^{-1}{\Sigma }_8{\Sigma }_7^{-1}\right)$

in distribution as $n\to \infty$, where

${\Sigma }_7=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>2}}{\overset{m_i}{\sum }}}E\left\{F_2\left(X_{ij}\mathrm{,}{\tau }_2\right)H\left({\bar{X}}_{ij}\right)\right\}\mathrm{,}$

${\Sigma }_8=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>2}}{\overset{m_i}{\sum }}}E\left\{G_2\left(X_{ij}\mathrm{,}{\tau }_2\right)H\left({\bar{X}}_{ij}\right)\right\}\mathrm{.}$

4. Implementation of the ESL Estimator

4.1. IRLS Algorithm

In this subsection, we develop an IRLS algorithm to calculate the parameter estimates. The IRLS algorithm has been commonly adopted for general M-estimators. Since the maximizers of (2) and (7) can be regarded as special M-estimators, the IRLS algorithm can be carried out to find $\tilde{\beta }$ and $\hat{\theta }$. In the following, we first develop the IRLS algorithm to find the maximizer of (2), and then we can calculate the maximizer of (7) in a similar way. Later, we summarize the algorithm in detail.

Because $\tilde{\beta }$ maximizes (2), we have the following normal equation:

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{{\varphi }^{\prime }}_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\right)=\mathrm{0,}$

or

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\right)\left(Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\right)=0.$ (9)

Let $W_{ij}={\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\right)$, then (9) can be transformed as

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{W}_{ij}\left(Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }\right)=\mathrm{0,}$

or $X^{{\rm T}}WX\tilde{\beta }=X^{{\rm T}}WY$, where $Y={\left(Y_1\mathrm{,}\cdots \mathrm{,}{Y}_n\right)}^{{\rm T}}$, $X={\left(X_1\mathrm{,}\cdots \mathrm{,}{X}_n\right)}^{{\rm T}}$ and $W=\text{diag}\left(W_{11}\mathrm{,}\cdots \mathrm{,}{W}_{n\mathrm{,}{m}_n}\right)$. Given the k-th approximation ${\tilde{\beta }}^{\left(k\right)}$, we can compute the corresponding weight matrix $W^{\left(k\right)}$ with $W_{ij}^{\left(k\right)}={\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}{\tilde{\beta }}^{\left(k\right)}\right)$. Then we have

${\beta }^{\left(k+1\right)}={\left(X^{{\rm T}}{W}^{\left(k\right)}X\right)}^{-1}{X}^{{\rm T}}{W}^{\left(k\right)}Y\mathrm{.}$ (10)

This iteration of (10) will monotonically non-decrease the objective function (2), that is, $L_1\left({\tilde{\beta }}^{\left(k+1\right)}\right)-L_1\left({\tilde{\beta }}^{\left(k\right)}\right)\ge 0$. In fact,

$\begin{array}{l}\log \left\{L_1\left({\tilde{\beta }}^{\left(k+1\right)}\right)\right\}-\log \left\{L_1\left({\tilde{\beta }}^{\left(k\right)}\right)\right\}\\ =\log \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{\left(k+1\right)}\right)-\log \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{\left(k\right)}\right)=\log \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}\frac{W_{ij}^{\left(k+1\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{(k)}}\right)\end{array}$

$=\log \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}\frac{W_{ij}^{\left(k\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{\left(k\right)}}\frac{W_{ij}^{\left(k+1\right)}}{W_{ij}^{\left(k\right)}}\right)=\log \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\pi }_{ij}^{\left(k\right)}\frac{W_{ij}^{\left(k+1\right)}}{W_{ij}^{\left(k\right)}}\right)\mathrm{,}$

where ${\pi }_{ij}^{\left(k\right)}=\frac{W_{ij}^{\left(k\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{\left(k\right)}}$. Based on the Jensen’s inequality, we have

$\begin{array}{l}\log \left\{L_1\left({\tilde{\beta }}^{\left(k+1\right)}\right)\right\}-\log \left\{L_1\left({\tilde{\beta }}^{\left(k\right)}\right)\right\}\\ \ge {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\pi }_{ij}^{\left(k\right)}\log \left(\frac{W_{ij}^{\left(k+1\right)}}{W_{ij}^{\left(k\right)}}\right)\\ =\frac{1}{{\tau }_1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\pi }_{ij}^{\left(k\right)}\left\{{\left(Y_{ij}-X_{ij}^{{\rm T}}{\tilde{\beta }}^{\left(k\right)}\right)}^2-{\left(Y_{ij}-X_{ij}^{{\rm T}}{\tilde{\beta }}^{\left(k+1\right)}\right)}^2\right\}.\end{array}$

From expression (10), we can find out that ${\tilde{\beta }}^{\left(k+1\right)}$ minimizes ${\left(Y-X\beta \right)}^{{\rm T}}{W}^{\left(k\right)}\left(Y-X\beta \right)$, or ${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{W}_{ij}^{\left(k\right)}{\left(Y_{ij}-X_{ij}^{{\rm T}}\beta \right)}^2$. Then we have $\log \left\{L_1\left({\tilde{\beta }}^{\left(k+1\right)}\right)\right\}-\log \left\{L_1\left({\tilde{\beta }}^{\left(k\right)}\right)\right\}\ge 0$.

The IRLS algorithm is summarized as follows:

Step 1. Computation of $\tilde{\beta }$ by maximizing (2). Take the initial value ${\tilde{\beta }}^{\left(0\right)}$ as the ordinary least squares (OLS) estimator. Given the k-th approximation ${\tilde{\beta }}^{\left(k\right)}$, the IRLS iteration updates ${\tilde{\beta }}^{\left(k+1\right)}$ through (10). Repeat this iteration until the convergence occurs. The resulting estimator is denoted as $\tilde{\beta }$.

Step 2. Computation of $\hat{\theta }={\left({\hat{\beta }}^{{\rm T}}\mathrm{,}{\hat{\gamma }}^{{\rm T}}\right)}^{{\rm T}}$ by maximizing (7). Take ${\hat{\theta }}^{\left(0\right)}$ as the OLS estimator by minimizing ${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\left(Y_{ij}-{\tilde{\delta }}_{ij}^{{\rm T}}\theta \right)}^2$. Similar to (10), given the k-th approximation ${\hat{\theta }}^{\left(k\right)}$, the IRLS iteration updates

${\hat{\theta }}^{\left(k+1\right)}={\left({\tilde{\Delta }}^{{\rm T}}{V}^{\left(k\right)}\tilde{\Delta }\right)}^{-1}{\tilde{\Delta }}^{{\rm T}}{V}^{\left(k\right)}Y\mathrm{,}$

where $\tilde{\Delta }={\left({\tilde{\delta }}_1\mathrm{,}\cdots \mathrm{,}{\tilde{\delta }}_n\right)}^{{\rm T}}$, ${\tilde{\delta }}_i={\left({\tilde{\delta }}_{i2}\mathrm{,}\cdots \mathrm{,}{\tilde{\delta }}_{i{m}_i}\right)}^{{\rm T}}$, $V^{\left(k\right)}=\text{diag}\left(V_{12}^{\left(k\right)}\mathrm{,}\cdots \mathrm{,}{V}_{n\mathrm{,}{m}_n}^{\left(k\right)}\right)$, and $V_{ij}^{\left(k\right)}={\varphi }_{{\tau }_2}\left(Y_{ij}-{\tilde{\delta }}_{ij}^{{\rm T}}{\hat{\theta }}^{\left(k\right)}\right)$. Repeat this iteration until the convergence occurs. The resulting estimator is denoted as $\hat{\theta }$.

4.2. The Choice of Tuning Parameters

In this subsection, we give a data driving method to determine the tuning parameters $\left\{{\tau }_1\mathrm{,}{\tau }_2\right\}$. In order to simplify the calculation with respect to ${\tau }_1$, we assume that the random error ${\epsilon }_{ij}$ ’s are independent of each other and $X_{ij}$ ’s. Then from Corollary 1, we can obtain that the ratio between the asymptotic variance of the initial ESL estimator and that of the OLS estimator for $\beta$ is

$r\left({\tau }_1\right)\triangleq \frac{G_1\left({\tau }_1\right)F_1^{-2}\left({\tau }_1\right)}{{\sigma }^2},$

where ${\sigma }^2=\mathrm{var}\left({\epsilon }_{ij}\right)$. Therefore, the ideal choice of ${\tau }_1$ is

${\tau }_{\mathrm{1,}\text{opt}}=\arg \underset{{\tau }_1}{\min }r\left({\tau }_1\right)=\arg \underset{{\tau }_1}{\min }{G}_1\left({\tau }_1\right)F_1^{-2}\left({\tau }_1\right)/{\sigma }^2\mathrm{.}$

Then $r\left({\tau }_1\right)$ can be estimated by ${\tilde{G}}_1\left({\tau }_1\right){\tilde{F}}_1^{-2}\left({\tau }_1\right)/{\tilde{\sigma }}^2$, where

${\tilde{F}}_1\left({\tau }_1\right)=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{″}}_{{\tau }_1}\left({\tilde{\epsilon }}_{ij}\right)\text{and}{\tilde{G}}_1\left({\tau }_1\right)=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{\prime }}_{{\tau }_1}{\left({\tilde{\epsilon }}_{ij}\right)}^2$

with ${\tilde{\epsilon }}_{ij}=Y_{ij}-X_{ij}^{{\rm T}}\tilde{\beta }$, and $\tilde{\sigma }$ is the standard deviation of ${\tilde{\epsilon }}_{ij}$. Then ${\tau }_{\mathrm{1,}\text{opt}}$ can be easily obtained using the grid search approach. ${\tau }_2$ can also be chosen in a similar way.

5. Simulation Studies

In this section, we conduct some simulation studies to investigate the finite sample performance of the proposed approach. We generate 200 datasets and consider sample sizes $n=50$, 100 and 200. In particular, the datasets are generated from the following model:

$Y_{ij}=X_{ij}^{{\rm T}}\beta +{\epsilon }_{ij},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}i=1,\cdots ,n,\mathrm{<mi>\; </mi>}j=1,\cdots ,m_i,$

where $m_i-1\sim \text{binomial}\left(\mathrm{6,0.8}\right)$, $X_{ij1}\mathrm{,}{X}_{ij2}\sim N\left(\mathrm{0,1}\right)$, $\beta ={\left(\mathrm{1,2}\right)}^{{\rm T}}$,

$t_{ij}\sim U\left(\mathrm{0,1}\right)$, $Z_{j\mathrm{,}k}^{\left(i\right)}={\left(\mathrm{1,}{t}_{ij}-t_{ik}\right)}^{{\rm T}}$, $\gamma ={\left(\mathrm{0.5,0.2}\right)}^{{\rm T}}$, $D_i=I_{m_i}$, and ${\epsilon }_i\sim N\left(\mathrm{0,}{\Sigma }_i\right)$ with ${\Phi }_i{\Sigma }_i{\Phi }_i^{{\rm T}}=D_i$.

To investigate robustness, we denote the above datasets as no contamination (NC) situation and consider the following four contaminations:

1) ${\epsilon }_i$ follows the multivariate t-distribution with 2 degrees of freedom and covariance matrix ${\Sigma }_i$.

2) ${\epsilon }_i$ follows the multivariate t-distribution with 2 degrees of freedom and covariance matrix ${\Sigma }_i$, and randomly choose 2% of $X_{ij}$ to be $X_{ij}+{\left(\mathrm{6,6}\right)}^{{\rm T}}$.

3) ${\epsilon }_i\sim N\left(\mathrm{0,}{\Sigma }_i\right)$, and randomly choose 2% of ${\epsilon }_{ij}$ to be ${\epsilon }_{ij}+15$.

4) ${\epsilon }_i\sim N\left(\mathrm{0,}{\Sigma }_i\right)$, randomly choose 2% of ${\epsilon }_{ij}$ to be ${\epsilon }_{ij}+15$ and 2% of $X_{ij}$ to be $X_{ij}+{\left(\mathrm{6,6}\right)}^{{\rm T}}$.

We compare the proposed ESL method with the OLS method, the M-estimator (M) in [26] and the quantile regression (QR) method in [27]. Note that the OLS, M and QR methods follow the estimation procedure similar to that of ESL, while the main difference is that the objective function is respectively replaced by their counterparts. To assess the finite sample performance, we calculate the mean and standard deviation (SD) for the estimators of $\beta$ and $\gamma$. The corresponding simulation results are displayed in Tables 1-5.

From Table 1, it can be observed that the performance of the M, QR and ESL methods is comparable to that of the OLS method, when the error follows a normal distribution and there are no outliers in the data. From Tables 2-5, it can be found out that the M, QR and ESL methods outperform significantly the OLS method, particularly in terms of SD, in several contamination cases; moreover, the ESL method always perform best in these cases. More specifically, Table 2 and Table 3 indicate that the M, QR and ESL methods outperform significantly the OLS method with respect to both $\beta$ and $\gamma$, when the error follows a heavy-tailed error distribution. Table 4 and Table 5 indicate that the M, QR and ESL methods outperform significantly the OLS method with respect to $\beta$, when there are outliers in responses; the OLS, M and QR methods perform rather poorly with respect to $\gamma$, however, the ESL method sill performs well in this case.

6. Real Data Analysis

In this section, we analyse the CD4 cell study, which was previously analysed by [7] [28]. This dataset comprises CD4 cell counts of 369 HIV-infected men, and there are totally 2376 values collected at different times for each individual, over a period of approximately eight and a half years. The number of measurements for each individual varies from 1 to 12 and the time points are not equally spaced. We use square roots of the CD4 counts [28] to make the response variable closer to the normal distribution, and the related six covariates are respectively time since seroconversion $t_{ij}$, age relative to arbitrary origin $X_{ij1}$, packs of cigarettes smoked per day $X_{ij2}$, recreational drug use $X_{ij3}$, number of sexual partners $X_{ij4}$ and mental illness score $X_{ij5}$. Note that Figure 1 displays the sample regressogram and local linear fitted curve for the square root of CD4 count over time, which reflects polynomial trend with respect to time. In the following, we use the mean model [1].

Table 1. Parameter estimates (with SD in parentheses) for the NC situation.

Table 2. Parameter estimates (with SD in parentheses) for the first contamination situation.

Table 3. Parameter estimates (with SD in parentheses) for the second contamination situation.

Table 4. Parameter estimates (with SD in parentheses) for the third contamination situation.

Table 5. Parameter estimates (with SD in parentheses) for the fourth contamination situation.

Figure 1. Sample regressogram and local linear fitted curve for the CD4 data: square root of CD4+ number versus time.

$Y_{ij}={\beta }_1{X}_{ij1}+{\beta }_2{X}_{ij2}+{\beta }_3{X}_{ij3}+{\beta }_4{X}_{ij4}+{\beta }_5{X}_{ij5}+f\left(t_{ij}\right)+{\epsilon }_{ij},$

where $f\left(t\right)={\beta }_0+{\beta }_6{t}_{+}+{\beta }_7{t}_{+}^2$ with $t_{+}=tI\left(t>0\right)$. We use cubic polynomial to model the generalized autoregressive parameters, that is,

$Z_{j\mathrm{,}k}^{\left(i\right)}={\left(\mathrm{1,}{t}_{ij}-t_{ik}\mathrm{,}{\left(t_{ij}-t_{ik}\right)}^2\mathrm{,}{\left(t_{ij}-t_{ik}\right)}^3\right)}^{{\rm T}}$.

We apply the OLS, M, QR methods and the proposed ESL method to the CD4 cell study. To assess the prediction performance, we randomly split the data into three parts, each with 369/3 = 123 subjects. We use the first two parts as the training dataset to fit the model, and then assess the out of sample performance on the testing dataset (defined as TD) which is left out. This process is repeated 200 times. We define the median absolute prediction error (MAPE) as median of $\left\{\left|Y_{ij}-{\hat{Y}}_{ij}\right|\mathrm{,}i\in \text{TD}\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}{m}_i\right\}$. To illustrate the estimation robustness of the proposed ESL method compared with the other methods, we re-analyse the dataset by including 5% outliers in the dataset, which are randomly generated by replacing $Y_{ij}$ with $Y_{ij}+300$. Moreover, we also re-analyse the dataset by including 5% outliers only in the training data set, in order to assess the robustness in terms of prediction performance. The results are displayed in Table 6. It can be observed that the estimates for $\gamma$ are very similar based on different methods, but the estimates for $\beta$ based on the ESL method are somewhat distinct with those of the OLS, M and QR methods in the no outlier case. Moreover, the MAPE of the proposed ESL method is slightly larger than the others, indicating that these methods possess comparative prediction performance in the case of no outlier. In the 5% outliers case, the MAPEs indicate that the ESL, M and QR methods perform similarly but much better than the OLS method in the prediction performance. Comparing the no outlier with 5% outliers case, we can see that the estimates based on the ESL method varies more slightly than the other methods, especially in terms of the estimates for $\gamma$ ; the ESL, M and QR methods are robust to the outliers, while the OLS method is adversely affected by the outliers.

7. Conclusions

In this paper, based on the ESL function, we proposed a robust estimation approach for the mean and generalized autoregressive parameters with longitudinal

Table 6. Parameter estimates and prediction results for the CD4 data.

data. The generalized autoregressive parameters that resulted from the MCD of the covariance matrix are unconstrained and can be well interpreted by statistical concepts in the framework of time series. Then the mean and generalized autoregressive parameters can be estimated via linear regression models using the ESL function. Moreover, the balance between the robustness and efficiency can be achieved by choosing appropriate data adaptive tuning parameters. Under certain conditions, we established the theoretical properties. Simulation studies and real data analysis were also carried out to illustrate the finite sample performance of our approach.

Several further problems need to be investigated. First, the dimension of the covariates in regression models is assumed to be fixed, thus it is interesting to extend our approach to the high-dimensional settings. Second, the models can be extended to nonparametric and semiparametric models. For more discussion along this line, references including [10] [29] may be helpful. Moreover, this paper targets the conditional mean of the response given covariates, which suffers from difficulties when the conditional distribution of the response is asymmetric. In this case, the conditional mode may be a more useful summary than the conditional mean, and thus modal linear regression [30] may be an interesting problem.

Acknowledgements

We thank the editor and the referee for their comments. This work was supported by the National Natural Science Foundation of China (No. 11971001) and the Beijing Natural Science Foundation (No. 1182002).

Appendix

Lemma 1 If regularity conditions (C1)-(C5) hold, then with probability approaching to 1, there exists a local maximizer of (2), denoted as $\tilde{\beta }$, such that

$\Vert \tilde{\beta }-{\beta }_0\Vert =O_P\left(\frac{1}{\sqrt{n}}\right)\mathrm{.}$

Proof of Lemma 1. Let $R_1\left(\beta \right)=\frac{1}{n}{L}_1\left(\beta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta \right)$. It suffices to show that for any given $\rho >0$, there exists a large constant $C>0$ such that

$P\left\{\underset{\Vert v\Vert =C}{\sup }{R}_1\left({\beta }_0+v/\sqrt{n}\right)<R_1\left({\beta }_0\right)\right\}\ge 1-\rho \mathrm{,}$ (1)

for any p-dimensional vector v satisfying that $\Vert v\Vert =C$. Based on the Taylor expansion, we have

$\begin{array}{l}{R}_1\left({\beta }_0+v/\sqrt{n}\right)-R_1\left({\beta }_0\right)\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}\left\{{\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}\left({\beta }_0+v/\sqrt{n}\right)\right)-{\varphi }_{{\tau }_1}\left(Y_{ij}-X_{ij}^{{\rm T}}{\beta }_0\right)\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}\left\{{\varphi }_{{\tau }_1}\left({\epsilon }_{ij}-v^{{\rm T}}{X}_{ij}/\sqrt{n}\right)-{\varphi }_{{\tau }_1}\left({\epsilon }_{ij}\right)\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}\left\{-{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}/\sqrt{n}+\frac{1}{2n}{{\varphi }^{″}}_{{\tau }_1}\left({\epsilon }_{ij}\right){\left(v^{{\rm T}}{X}_{ij}\right)}^2-\frac{1}{6n^{\frac{3}{2}}}{{\varphi }^{‴}}_{{\tau }_1}\left({\epsilon }_{ij}^{\ast }\right){\left(v^{{\rm T}}{X}_{ij}\right)}^3\right\}\\ \triangleq I_1+I_2+I_3,\end{array}$

where ${\epsilon }_{ij}^{\ast }$ lies between ${\epsilon }_{ij}-v^{{\rm T}}{X}_{ij}/\sqrt{n}$ and ${\epsilon }_{ij}$. Then, for $I_1$, we have

$\begin{array}{c}E\left(I_1\right)=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}/\sqrt{n}\right\}\\ =-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{E\left({{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)|X_{ij}\right)v^{{\rm T}}{X}_{ij}/\sqrt{n}\right\}=0,\end{array}$

and

$\mathrm{var}\left(I_1\right)=E{\left\{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}/\sqrt{n}\right\}}^2=\frac{1}{n^3}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E{\left\{{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}\right\}}^2.$

Using the C_r inequality and condition (C1), we can obtain that

$E{\left\{{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}\right\}}^2\le m_i{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E{\left\{{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}\right\}}^2\le M{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E{\left\{{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)v^{{\rm T}}{X}_{ij}\right\}}^2$

for $i=1,\cdots ,n$, then we have

$\begin{array}{c}\mathrm{var}\left(I_1\right)\le \frac{M}{n^3}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{v}^{{\rm T}}E\left\{{{\varphi }^{\prime }}_{{\tau }_1}{\left({\epsilon }_{ij}\right)}^2{X}_{ij}{X}_{ij}^{{\rm T}}\right\}v\\ =\frac{M}{n^3}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{v}^{{\rm T}}E\left\{G_1\left(X_{ij},{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right\}v\\ =O\left({\Vert v\Vert }^2{n}^{-2}\right)=O\left(C^2{n}^{-2}\right).\end{array}$

Therefore, we can deduce that $I_1=E\left(I_1\right)+O_P\left(\sqrt{\mathrm{var}\left(I_1\right)}\right)=O_P\left(C{n}^{-1}\right)$. Similarly, we get $I_3=O_P\left(n^{-3/2}\right)$. For $I_2$, we have

$\begin{array}{c}{I}_2=\frac{1}{2n^2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{v}^{{\rm T}}E\left(F_1\left(X_{ij},{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right)v\left(1+o_P\left(1\right)\right)\\ =O\left({\Vert v\Vert }^2{n}^{-1}\right)=O\left(C^2{n}^{-1}\right).\end{array}$

Note that $\Vert v\Vert =C$, we can choose a sufficiently large C such that $I_2$ dominates both $I_1$ and $I_3$ with a probability of at least $1-\rho$. From condition (C3), we get $F_1\left(x,{\tau }_1\right)<0$, then (1) holds. The proof of Lemma 1 is finished.

Proof of Theorem 1. Let ${\tilde{\phi }}_{ij}=X_{ij}^{{\rm T}}\left(\tilde{\beta }-{\beta }_0\right)$, then $\tilde{\beta }$ satisfies the following equation:

$\begin{array}{c}0=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}-{\tilde{\phi }}_{ij}\right)\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}\left\{{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)-{{\varphi }^{″}}_{{\tau }_1}\left({\epsilon }_{ij}\right){\tilde{\phi }}_{ij}+\frac{1}{2}{{\varphi }^{‴}}_{{\tau }_1}\left({\epsilon }_{ij}^{\ast }\right){\tilde{\phi }}_{ij}^2\right\}\\ \triangleq I_4+I_5+I_6,\end{array}$

where ${\epsilon }_{ij}^{\ast }$ locates between ${\epsilon }_{ij}-{\tilde{\phi }}_{ij}$ and ${\epsilon }_{ij}$. It can be easily shown that

$\begin{array}{c}{I}_5=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{{\varphi }^{″}}_{{\tau }_1}\left({\epsilon }_{ij}\right)X_{ij}{X}_{ij}^{{\rm T}}\left(\tilde{\beta }-{\beta }_0\right)\\ =-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left(F_1\left(X_{ij},{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right)\left(\tilde{\beta }-{\beta }_0\right)+o_P\left(1\right).\end{array}$

Since $\left|{\tilde{\phi }}_{ij}\right|=\left|X_{ij}^{{\rm T}}\left(\tilde{\beta }-{\beta }_0\right)\right|=O_P\left(\Vert \tilde{\beta }-{\beta }_0\Vert \right)$, then ${\left|{\tilde{\phi }}_{ij}\right|}^2=O_P\left({\Vert \tilde{\beta }-{\beta }_0\Vert }^2\right)$. Thus, $I_6=o_P\left(I_5\right)$. Then,

$\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left(F_1\left(X_{ij},{\tau }_1\right)X_{ij}{X}_{ij}^{{\rm T}}\right)\sqrt{n}\left(\tilde{\beta }-{\beta }_0\right)=\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)+o_P\left(1\right).$

Notice that $E\left(X_{ij}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)\right)=0$. For any $p\times 1$ vector $\kappa$, let $e_i={\kappa }^{{\rm T}}{U}_i$, where $U_i={\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}{X}_{ij}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)=X_i^{{\rm T}}{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_i\right)$. By C_r inequality and condition (C1), we have

$\begin{array}{c}E{\left|e_i\right|}^3\le M{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E{\left|\left({\kappa }^{{\rm T}}{X}_{ij}\right){{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)\right|}^3\\ =M{\displaystyle \underset{j=1}{\overset{m_i}{\sum }}}E\left\{E\left({\left|{{\varphi }^{\prime }}_{{\tau }_1}\left({\epsilon }_{ij}\right)\right|}^3|X_{ij}\right){\left|{\kappa }^{{\rm T}}{X}_{ij}\right|}^3\right\}\\ \le K_1,\end{array}$

where $K_1$ is a constant independent of i. Then by condition (C5), $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E\left(U_i{U}_i^{{\rm T}}\right)$ converges to a finite positive definite matrix, which is denoted by ${\Sigma }_2$ here. Thus, we obtain that

$\begin{array}{l}\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\mathrm{var}\left(e_i\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E\left(e_i^2\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\kappa }^{{\rm T}}E\left(U_i{U}_i^{{\rm T}}\right)\kappa \to {\kappa }^{{\rm T}}{\Sigma }_2\kappa >0.\hfill \end{array}$

So the proof of theorem 1 is completed by the Lyapunov central limit theorem and Slutsky’s theorem.

Lemma 2 If regularity conditions (C1)-(C9) hold, then with probability approaching to 1, there exists a local maximizer of (7), denoted as $\hat{\theta }$, such that

$\Vert \hat{\theta }-{\theta }_0\Vert =O_P\left(\frac{1}{\sqrt{n}}\right)\mathrm{.}$

Proof of Lemma 2. Since $\Vert \tilde{\beta }-{\beta }_0\Vert =O_P\left(\frac{1}{\sqrt{n}}\right)$, we can replace $\tilde{\beta }$ with ${\beta }_0$ during the proof from now on. Then ${\tilde{\delta }}_{ij}={\left(X_{ij}^{{\rm T}}\mathrm{,}{\tilde{\zeta }}_{ij}^{{\rm T}}\right)}^{{\rm T}}$ can be substituted by ${\delta }_{ij}={\left(X_{ij}^{{\rm T}}\mathrm{,}{\zeta }_{ij}^{{\rm T}}\right)}^{{\rm T}}$, where ${\zeta }_{ij}={\left({\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\epsilon }_{i,k}{Z}_{j,k,1}^{\left(i\right)},\cdots ,{\displaystyle \underset{k=1}{\overset{j-1}{\sum }}}{\epsilon }_{i,k}{Z}_{j,k,q}^{\left(i\right)}\right)}^{{\rm T}}$. Let

$R_2\left(\theta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-X_{ij}^{{\rm T}}\beta -{\zeta }_{ij}^{{\rm T}}\gamma \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\varphi }_{{\tau }_2}\left(Y_{ij}-{\delta }_{ij}^{{\rm T}}\theta \right).$

It suffices to show that for any given $\rho >0$, there exists a large constant $C>0$ such that

$P\left\{\underset{\Vert \nu \Vert =C}{\lim }{R}_2\left({\theta }_0+\nu /\sqrt{n}\right)<R_2\left({\theta }_0\right)\right\}\ge 1-\rho \mathrm{,}$ (2)

for any $\left(p+q\right)$ -dimensional vector $\nu$ satisfying that $\Vert \nu \Vert =C$. Based on the Taylor expansion, we have

$\begin{array}{l}{R}_2\left({\theta }_0+\nu /\sqrt{n}\right)-R_2\left({\theta }_0\right)\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}\left\{{\varphi }_{{\tau }_2}\left(Y_{ij}-{\delta }_{ij}^{{\rm T}}\left({\theta }_0+\nu /\sqrt{n}\right)\right)-{\varphi }_{{\tau }_2}\left(Y_{ij}-{\delta }_{ij}^{{\rm T}}{\theta }_0\right)\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}\left\{{\varphi }_{{\tau }_2}\left({ϵ}_{ij}-{\nu }^{{\rm T}}{\delta }_{ij}/\sqrt{n}\right)-{\varphi }_{{\tau }_2}\left({ϵ}_{ij}\right)\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}\left\{-{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right){\nu }^{{\rm T}}{\delta }_{ij}/\sqrt{n}+\frac{1}{2n}{{\varphi }^{″}}_{{\tau }_2}\left({ϵ}_{ij}\right){\left({\nu }^{{\rm T}}{\delta }_{ij}\right)}^2-\frac{1}{6n^{\frac{3}{2}}}{{\varphi }^{‴}}_{{\tau }_2}\left({ϵ}_{ij}^{\ast }\right){\left({\nu }^{{\rm T}}{\delta }_{ij}\right)}^3\right\}\\ \triangleq I_7+I_8+I_9,\end{array}$

where ${ϵ}_{ij}^{\ast }$ lies between ${ϵ}_{ij}-{\nu }^{{\rm T}}{\delta }_{ij}/\sqrt{n}$ and ${ϵ}_{ij}$. Recall that ${\bar{X}}_{ij}={\left(X_{i1}\mathrm{,}\cdots \mathrm{,}{X}_{ij}\right)}^{{\rm T}}$. It should be emphasized that ${\delta }_{ij}$ consists of both covariates and ${\epsilon }_{i1}\mathrm{,}\cdots \mathrm{,}{\epsilon }_{i\mathrm{,}j-1}$, and thus ${ϵ}_{ij}$ is independent of ${\delta }_{ij}$ conditional on ${\bar{X}}_{ij}$. In addition, from (8), we have

$E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}\right)=\left(\begin{array}{cc}E\left(X_{ij}{X}_{ij}^{{\rm T}}\right)& 0\\ 0& E\left({\zeta }_{ij}{\zeta }_{ij}^{{\rm T}}\right)\end{array}\right)\mathrm{.}$

Then, for $I_7$, we have

$\begin{array}{c}E\left(I_7\right)=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right){\nu }^{{\rm T}}{\delta }_{ij}/\sqrt{n}\right\}\\ =-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{E\left({{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)|X_{ij}\right)E\left({\nu }^{{\rm T}}{\delta }_{ij}|{\bar{X}}_{ij}\right)/\sqrt{n}\right\}\\ =0,\end{array}$

and

$\begin{array}{c}\mathrm{var}\left(I_7\right)=E{\left\{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right){\nu }^{{\rm T}}{\delta }_{ij}/\sqrt{n}\right\}}^2\\ =\frac{1}{n^3}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\nu }^{{\rm T}}E\left\{{{\varphi }^{\prime }}_{{\tau }_2}{\left({ϵ}_{ij}\right)}^2{\delta }_{ij}{\delta }_{ij}^{{\rm T}}\right\}\nu \\ =\frac{1}{n^3}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\nu }^{{\rm T}}E\left\{G_2\left(X_{ij},{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}|{\bar{X}}_{ij}\right)\right\}\nu \\ =O\left({\Vert \nu \Vert }^2{n}^{-2}\right)=O\left(C^2{n}^{-2}\right).\end{array}$

Therefore, we can deduce that $I_7=E\left(I_7\right)+O_P\left(\sqrt{\mathrm{var}\left(I_7\right)}\right)=O_P\left(C{n}^{-1}\right)$. Similarly, we get $I_9=O_P\left(n^{-3/2}\right)$. For $I_8$, we have

$\begin{array}{c}{I}_8=\frac{1}{2n^2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\nu }^{{\rm T}}E\left\{F_2\left(X_{ij},{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}|{\bar{X}}_{ij}\right)\right\}\nu \left(1+o_P\left(1\right)\right)\\ =O\left({\Vert \nu \Vert }^2{n}^{-1}\right)=O\left(C^2{n}^{-1}\right).\end{array}$

Note that $\Vert \nu \Vert =C$, we can choose a sufficiently large C such that $I_8$ dominates both $I_7$ and $I_9$ with a probability of at least $1-\rho$. From condition (C8), we get $F_2\left(x,{\tau }_2\right)<0$, then (2) holds. The proof of Lemma 2 is finished.

Proof of Theorem 2. Let ${\hat{\alpha }}_{ij}={\delta }_{ij}^{{\rm T}}\left(\hat{\theta }-{\theta }_0\right)$, then $\hat{\theta }$ satisfies the following equation:

$\begin{array}{c}0=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\delta }_{ij}{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}-{\hat{\alpha }}_{ij}\right)\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\delta }_{ij}\left\{{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)-{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right){\hat{\alpha }}_{ij}+\frac{1}{2}{{\varphi }^{‴}}_{{\tau }_2}\left({ϵ}_{ij}^{\ast }\right){\hat{\alpha }}_{ij}^2\right\}\\ \triangleq I_{10}+I_{11}+I_{12},\end{array}$

where ${ϵ}_{ij}^{\ast }$ locates between ${ϵ}_{ij}-{\hat{\alpha }}_{ij}$ and ${ϵ}_{ij}$. It can be easily obtained that

$\begin{array}{c}{I}_{11}=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{{\varphi }^{″}}_{{\tau }_2}\left({ϵ}_{ij}\right){\delta }_{ij}{\delta }_{ij}^{{\rm T}}\left(\hat{\theta }-{\theta }_0\right)\\ =-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{F_2\left(X_{ij},{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}|{\bar{X}}_{ij}\right)\right\}\left(\hat{\theta }-{\theta }_0\right)+o_P\left(1\right),\end{array}$

and

$I_{12}=\frac{1}{2n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{{\varphi }^{‴}}_{{\tau }_2}\left({ϵ}_{ij}^{\ast }\right){\delta }_{ij}{\left\{{\delta }_{ij}^{{\rm T}}\left(\hat{\theta }-{\theta }_0\right)\right\}}^2=O_P\left({\Vert \hat{\theta }-{\theta }_0\Vert }^2\right)=o_P\left(I_{11}\right).$

Then

$\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{F_2\left(X_{ij},{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}|{\bar{X}}_{ij}\right)\right\}\sqrt{n}\left(\hat{\theta }-{\theta }_0\right)=\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\delta }_{ij}{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)+o_P\left(1\right).$

Notice that $E\left({\delta }_{ij}{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)\right)=0$. For any $\left(p+q\right)\times 1$ vector $\eta$, let ${\psi }_i={\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}\left({\eta }^{{\rm T}}{\delta }_{ij}\right){{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)$. By C_r inequality and condition (C1), we have

$\begin{array}{c}E{\left|{\psi }_i\right|}^3\le M{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E{\left|\left({\eta }^{{\rm T}}{\delta }_{ij}\right){{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)\right|}^3\\ =M{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{E\left({\left|{{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)\right|}^3|X_{ij}\right)E\left({\left|{\eta }^{{\rm T}}{\delta }_{ij}\right|}^3|{\bar{X}}_{ij}\right)\right\}\\ \le K_2,\end{array}$

where $K_2$ is a constant independent of i. Moreover,

$\begin{array}{c}\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\mathrm{var}\left({\psi }_i\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E\left({\psi }_i^2\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}E{\left\{{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}\left({\eta }^{{\rm T}}{\delta }_{ij}\right){{\varphi }^{\prime }}_{{\tau }_2}\left({ϵ}_{ij}\right)\right\}}^2\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}E\left\{{\left({\eta }^{{\rm T}}{\delta }_{ij}\right)}^2{{\varphi }^{\prime }}_{{\tau }_2}{\left({ϵ}_{ij}\right)}^2\right\}\\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=2}{\overset{m_i}{\sum }}}{\eta }^{{\rm T}}E\left\{G_2\left(X_{ij},{\tau }_2\right)E\left({\delta }_{ij}{\delta }_{ij}^{{\rm T}}|{\bar{X}}_{ij}\right)\right\}\eta \\ \to K_3>0.\end{array}$

So the proof is completed by the Lyapunov central limit theorem and Slutsky’s theorem.

List of Main Symbols

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Diggle, P.J., Heagerty, P.J., Liang, K.-Y. and Zeger, S. (2002) Analysis of Longitudinal Data. Oxford University Press, Oxford.
[2]	Diggle, P.J. and Verbyla, A.P. (1998) Nonparametric Estimation of Covariance Structure in Longitudinal Data. Biometrics, 91, 403-415. https://doi.org/10.2307/3109751
[3]	Liang, K.-Y. and Zeger, S.L. (1986) Longitudinal Data Analysis Using Generalized Linear Models. Biometrika, 73, 13-22. https://doi.org/10.1093/biomet/73.1.13
[4]	Pourahmadi, M. (1999) Joint Mean-Covariance Models with Applications to Longitudinal Data: Unconstrained Parameterisation. Biometrika, 86, 677-690. https://doi.org/10.1093/biomet/86.3.677
[5]	Pourahmadi, M. (2000) Maximum Likelihood Estimation of Generalised Linear Models for Multivariate Normal Covariance Matrix. Biometrika, 87, 425-435. https://doi.org/10.1093/biomet/87.2.425
[6]	Pan, J. and Mackenzie, G. (2003) On Modelling Mean-Covariance Structures in Longitudinal Studies. Biometrika, 90, 239-244. https://doi.org/10.1093/biomet/90.1.239
[7]	Ye, H. and Pan, J. (2006) Modelling of Covariance Structures in Generalised Estimating Equations for Longitudinal Data. Biometrika, 93, 927-941. https://doi.org/10.1093/biomet/93.4.927
[8]	Pourahmadi, M. (2007) Cholesky Decompositions and Estimation of a Covariance Matrix: Orthogonality of Variance-Correlation Parameters. Biometrika, 94, 1006-1013. https://doi.org/10.1093/biomet/asm073
[9]	Daniels, M.J. and Pourahmadi, M. (2009) Modeling Covariance Matrices via Partial Autocorrelations. Journal of Multivariate Analysis, 100, 2352-2363. https://doi.org/10.1016/j.jmva.2009.04.015
[10]	Leng, C., Zhang, W. and Pan, J. (2010) Semiparametric Mean-Covariance Regression Analysis for Longitudinal Data. Journal of the American Statistical Association, 105, 181-193. https://doi.org/10.1198/jasa.2009.tm08485
[11]	Zhang, W. and Leng, C. (2012) A Moving Average Cholesky Factor Model in Covariance Modelling for Longitudinal Data. Biometrika, 99, 141-150. https://doi.org/10.1093/biomet/asr068
[12]	Zhang, W., Leng, C. and Tang, C.Y. (2015) A Joint Modelling Approach for Longitudinal Studies. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77, 219-238. https://doi.org/10.1111/rssb.12065
[13]	Maadooliat, M., Pourahmadi, M. and Huang, J.Z. (2013) Robust Estimation of the Correlation Matrix of Longitudinal Data. Statistics and Computing, 23, 17-28. https://doi.org/10.1007/s11222-011-9284-6
[14]	Zheng, X., Fung, W.K. and Zhu, Z. (2013) Robust Estimation in Joint Mean-Covariance Regression Model for Longitudinal Data. Annals of the Institute of Statistical Mathematics, 65, 617-638. https://doi.org/10.1007/s10463-012-0383-8
[15]	Zheng, X., Fung, W.K. and Zhu, Z. (2014) Variable Selection in Robust Joint Mean and Covariance Model for Longitudinal Data Analysis. Statistica Sinica, 24, 515-531.
[16]	Lv, J. and Guo, C. (2017) Efficient Parameter Estimation via Modified Cholesky Decomposition for Quantile Regression with Longitudinal Data. Computational Statistics, 32, 947-975. https://doi.org/10.1007/s00180-017-0714-6
[17]	Lv, J., Guo, C., Yang, H. and Li, Y. (2017) A Moving Average Cholesky Factor Model in Covariance Modeling for Composite Quantile Regression with Longitudinal Data. Computational Statistics & Data Analysis, 112, 129-144. https://doi.org/10.1016/j.csda.2017.02.015
[18]	Lv, J., Guo, C. and Wu, J. (2019) Smoothed Empirical Likelihood Inference via the Modified Cholesky Decomposition for Quantile Varying Coefficient Models with Longitudinal Data. TEST, 28, 999-1032. https://doi.org/10.1007/s11749-018-0616-0
[19]	Lv, J. and Guo, C. (2019) Quantile Estimations via Modified Cholesky Decomposition for Longitudinal Single-Index Models. Annals of the Institute of Statistical Mathematics, 71, 1163-1199. https://doi.org/10.1007/s10463-018-0673-x
[20]	Wang, X., Jiang, Y., Huang, M. and Zhang, H. (2013) Robust Variable Selection with Exponential Squared Loss. Journal of the American Statistical Association, 108, 632-643. https://doi.org/10.1080/01621459.2013.766613
[21]	Lv, J., Yang, H. and Guo, C. (2015) An Efficient and Robust Variable Selection Method for Longitudinal Generalized Linear Models. Computational Statistics & Data Analysis, 82, 74-88. https://doi.org/10.1016/j.csda.2014.08.006
[22]	Wang, K. and Lin, L. (2016) Robust and Efficient Estimator for Simultaneous Model Structure Identification and Variable Selection in Generalized Partial Linear Varying Coefficient Models with Longitudinal Data. Statistical Papers, 1-28.
[23]	Li, S., Wang, K. and Ren, Y. (2018) Robust Estimation and Empirical Likelihood Inference with Exponential Squared Loss for Panel Data Models. Economics Letters, 164, 19-23. https://doi.org/10.1016/j.econlet.2017.12.029
[24]	Yang, H., Lv, J. and Guo, C. (2015) Robust Variable Selection in Modal Varying-Coefficient Models with Longitudinal. Journal of Statistical Computation and Simulation, 85, 3064-3079. https://doi.org/10.1080/00949655.2014.949716
[25]	Wang, K., Li, S., Sun, X. and Lin, L. (2019) Modal Regression Statistical Inference for Longitudinal Data Semivarying Coefficient Models: Generalized Estimating Equations, Empirical Likelihood and Variable Selection. Computational Statistics & Data Analysis, 133, 257-276. https://doi.org/10.1016/j.csda.2018.10.010
[26]	Huber, P.J. (2011) Robust Statistics. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-04898-2_594
[27]	Koenker, R. (2005) Quantile Regression (Econometric Society Monographs; No. 38). Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511754098
[28]	Zeger, S.L. and Diggle, P.J. (1994) Semiparametric Models for Longitudinal Data with Application to CD4 Cell Numbers in HIV Seroconverters. Biometrics, 50, 689-699. https://doi.org/10.2307/2532783
[29]	Yao, W. and Li, R. (2013) New Local Estimation Procedure for a Non-Parametric Regression Function for Longitudinal Data. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75, 123-138. https://doi.org/10.1111/j.1467-9868.2012.01038.x
[30]	Yao, W. and Li, L. (2014) A New Regression Model: Modal Linear Regression. Scandinavian Journal of Statistics, 41, 656-671. https://doi.org/10.1111/sjos.12054