A new GEE method to account for heteroscedasticity, using asymmetric least-square regressions
WWeighted asymmetric least squares regression forlongitudinal data using GEE
Amadou Barry ∗1 , Karim Oualkacha , Arthur Charpentier Department of Mathematics, Université du Québec à Montréal, Canada, Faculty of Economics, Université de Rennes 1, France
Abstract
The well-known generalized estimating equations (GEE) is widely used to estimate the effect of thecovariates on the mean of the response variable. We apply the GEE method using the asymmetricleast-square regression (expectile) to analyze the longitudinal data. Expectile regression naturally extendsthe classical least squares method and has properties similar to quantile regression. Expectile regressionallows the study of the heterogeneity of the effects of the covariates over the entire distribution of theresponse variable, while also accounting for unobserved heterogeneity. In this paper, we present thegeneralized expectile estimating equations estimators, derive their asymptotic properties and propose arobust estimator of their variance-covariance matrix for inference. The performance of the new estimatorsis evaluated through exhaustive simulation studies, and their advantages in relation to existing methodsare highlighted. Finally, the labor pain dataset is analyzed to illustrate the usefulness of the proposedmodel. keywords:
Expectile regression, quantile regression, GEE, working correlation, cluster data, longitu-dinal data.
Longitudinal and clustered data arise in many application fields such as epidemiology (Smith et al. 2015),genetics (Furlotte, Eskin, and Eyheramendy 2012), economics (Hsiao 2007), and other fields of biological andsocial sciences. They are characterized by the presence of a within-subject dependence which gives themdesirable properties over the cross-Sectional data. However, that dependency makes the statistical analysischallenging and needs to be addressed in order to generate unbiased and high efficient estimators. Generalizedestimating equations (GEE) (Liang and Zeger 1986) is a commonly used method for the analysis of such datawithin a marginal model framework.A marginal model estimates the expectation of the marginal distribution of the outcome without specifyingthe within-subject dependence, like cross-Sectional models (Fitzmaurice et al. 2008, Diggle et al. (2013)). TheGEE model completes the marginal model by introducing a “working” covariance matrix in the estimationprocess to account for the within-subject dependence. As a result, the GEE yields a consistent estimatorwith high efficiency even with misspecification of the true covariance structure (Liang and Zeger 1986). TheGEE model estimates only the effects of the covariates on the expectation of the marginal distribution of theoutcome. The expectation is a very important summary statistic, but limiting the study of the effects of thecovariates to this is not enough unless the covariates uniformly affect the whole distribution of the responsevariable. With its favorable properties, the GEE can be extended beyond the mean using the expectileregression (ER).ER models the relationship between the covariates and the response variable by estimating the effect ofthe predictors at different points of the conditional cumulative distribution function (c.d.f.) of the responsevariable. These points are generally the mean (expectile of level 0.5), and other expectiles above and belowthe mean. ER estimates the impact of the covariates on the location, scale, and shape of the responsedistribution. By doing so, the ER captures the heterogeneity of the effects of the covariates on the responsevariable; for example, when covariates affect the mean and the tail of the distribution in different ways. ∗ Corresponding author: barryhafi[email protected] a r X i v : . [ s t a t . M E ] O c t R was introduced by Aigner, Amemiya, and Poirier (1976) under a likelihood-based approach. A decadelater, Newey and Powell (1987) presented a detailed study of this new class of estimator. They presented itsfavorable properties, like its location and scale equivariance property, and derived its asymptotic properties.After the Efron paper (Efron 1991), expectile lived quietly in the shadows for decades, as mentioned in (Eilers2013). But recently it has comeback into the spotlight. After its re-discovery, early contributions to ERfocused on the application of the ER method to spline and smoothing model (Schnabel and Eilers 2009, Rossiand Harvey (2009), Sobotka and Kneib (2012), Sobotka et al. (2013)). Other works focused on contrastingER and QR, on showing how to get quantiles from a fine grid of expectiles, on the crossing curves problemand on promoting application of ER (Kneib 2013a, Schnabel and Eilers (2013), Waltrup et al. (2015)). TodayER is extended to many classes of models, such as Bayesian (Majumdar and Paul 2016,Waldmann, Sobotka,and Kneib (2016), Xing and Qian (2017)), nonparametric (Righi, Yang, and Ceretta 2014, Yang and Zou(2015)), nonlinear (Kim and Lee 2016), neural network (Xu et al. 2016, Jiang et al. (2017)), and supportvector machine (Farooq and Steinwart 2017). Recently, Waltrup and Kauermann (Schulze Waltrup andKauermann 2015) combined smoothing and random intercept to fit clustered data with penalized splines.ER generalizes mean regression in the same way that quantile regression (QR) (Koenker and Bassett 1978)generalizes median regression. The QR method was adapted to longitudinal data using GEE approach. Themain idea consists of smoothing the QR estimating functions in order to make them differentiable withrespect to regression parameters. Jung (1996) proposed the quasi-likelihood approach to analyze the medianregression model for dependent data. Chen, Wei, and Parzen (2004) derived a QR estimator for correlateddata using GEE approach based on independence working correlation. Along the same lines, Fu and Wang(2012) combined the between- and within-subject estimating functions to account for the correlations betweenrepeated measurements in the estimation of QR model. Lu and Fan (2015) proposed a general stationaryauto-correlation matrix for the working correlation to enhance the efficiency of the QR inference.Both QR and ER provide an overview of the effects of the covariates on the distribution of the responsevariable. Their resemblance and usefulness have already been discussed in the literature, see for examples(Efron 1991, Kneib (2013b), Waltrup et al. (2015)).This paper makes its contribution by introducing a new class of estimators for the analysis of dependentdata. Section 2 defines the expectile statistic, and introduces the ER method for cross-Sectional data and thegeneralized expectile estimating equation (GEEE) method for longitudinal data. In Section 3, the asymptoticproperties of the GEEE estimator of the model parameters and an estimator of its variance-covariance matrixare presented. The evaluation of the small sample performance of the estimators, through extensive simulationstudies, is presented in Section 4. The GEEE estimator is applied to a real data set and the results arepresented in Section 5; the conclusion is presented in Section 6. All proofs are in the appendix.
This Section introduces the univariate expectile and the ER model.
Expectile of a random variable Y is defined as the solution µ τ ( Y ) which minimize the loss function E { ρ τ ( Y − θ ) } (1)over θ ∈ R for a fixed value of τ ∈ (0 , . The function ρ τ ( · ) , of the form ρ τ ( t ) = | τ − ( t ≤ | · t
2s the asymmetric square loss function that assigns weights τ and 1 − τ to positive and negative deviations,respectively.By equating the first derivative of (1) to zero, the expectile can also be defined as solution of µ τ ( Y ) = µ τ = µ − − τ − τ E (cid:2) { Y − µ τ ( Y ) } { Y > µ τ ( Y ) } (cid:3) , (2)where µ = µ . ( Y ) = E ( Y ) . This definition, presented by Newey and Powell (1987), shows that expectile isdetermined by the tail expectations of the distribution of Y. Interestingly, we found that expectile can bedefined as µ τ = E " ψ τ ( Y − µ τ ) E (cid:2) ψ τ ( Y − µ τ ) (cid:3) Y , where ψ τ ( t ) = | τ − ( t ≤ | is the check function. This latter definition, which is much more meaningful inthe context of regression, reveals that expectiles, like the mean, are weighted averages.Given a random sample, { ( y i ) } ni =1 , the τ -th empirical expectile b µ τ = n X i =1 ψ τ ( y i − b µ τ ) P ni =1 ψ τ ( y i − b µ τ ) y i is the solution which minimizes the empirical loss function1 n n X i =1 ρ τ ( y i − θ ) . (3)Newey and Powell (1987) have shown that expectile function has attractive properties. Expectile summarizesthe c.d.f. of Y in the same way that quantile does. Moreover, expectile is location and scale equivariant, thatis for s > t ∈ R , µ τ ( sY + t ) = sµ τ ( Y ) + t. More details about the properties of expectile and resultson ER are given by Efron (1991).To introduce the ER method, consider the classical linear regression y i = x i T β + ε i , (4)where y i is the scalar response, ε i is the random error, x i ∈ R p is the vector of covariates and β ∈ R p is theunknown parameter that needs to be estimated. Under this framework, Newey and Powell (1987) introducedthe ER model for a fixed τ ∈ (0 ,
1) as µ τ ( y i | x i ) = x i T β τ , with µ τ ( ε i ) = 0 . (5)The assumption µ τ ( ε i ) = 0 ensures that the random error is centered on the τ -th expectile. The correspondingER estimator is defined as the unique solution which minimizes the objective function1 n n X i =1 ρ τ ( y i − x i T β τ ) (6)over β τ ∈ R p . The asymmetric loss function associated with the expectile function is continuously differentiable,and solving equation (6) gives 3 β τ = (cid:16) n X i =1 x i T ψ τ ( b ε it ) x i (cid:17) − (cid:16) n X i =1 x i ψ τ ( b ε it ) y i (cid:17) , (7)where b ε it = y i − x i T b β τ . The ER estimator is easily computed with the iterated weighted least squaresalgorithm. In addition to deriving the asymptotic properties of the above ER estimator, Newey and Powell(1987) proposed a robust estimator of the variance-covariance matrix of b β τ . Note that the ER estimator was presented previously by Aigner, Amemiya, and Poirier (1976) through alikelihood-based approach. The likelihood is derived by assuming an asymmetric normal distribution (AND)for the disturbance, u ∼ AND( u ; µ τ , σ , τ ) , with density f ( u ; µ τ , σ , τ ) = 2 √ πσ p τ (1 − τ ) √ τ + √ − τ exp ( − ρ τ u − µ τ σ !) , (8)where µ τ , σ, and τ are respectively the location, scale and asymmetric parameters. Now substitute µ iτ = x i T β τ and assume the n observations are independents, then the ER estimator is equivalent to the maximumof the likelihood function L ( β ; σ, τ, y ) ∝ σ − n exp ( − n X i =1 ρ τ y i − x i T β τ σ !) , where y = ( y , . . . , y n ) T . The AND distribution is not to be confused with the class of density functionsrelated to the standard density function and proposed by Azzalini (1985).
This Section presents the model and method of the GEEE for longitudinal data. Consider the data { y it , x it } ≤ i ≤ n, ≤ t ≤ m i generated by the following model y it = x it T β + ε it , (9)where y it is the t -th observation of the response variable for the i -th individual, x it = ( x it , . . . , x pit ) is the p × ε it the random error and β the p × y i = X i β + ε i , (10)where y i is the dependent observations of the response variable of the individual i, X i the corresponding m i × p matrix covariates and ε i the vector error. Individual observations can also be stacked and presentedin matrix form as y = Xβ + ε , (11)where y and ε are N × X is N × p matrix and N = P ni =1 m i . Using the location-scale equivariance property of the expectile function, the corresponding conditionalexpectile of level τ of the model equation (9) is defined as µ τ ( y it | x it ) = x it T β τ , µ τ ( ε it ) = 0 . (12)4he assumption, µ τ ( ε it ) = 0 , is introduced to guarantee that the random error is centered on the τ -thexpectile. The parameter β τ measures the effect of the covariates x it on the location, scale and shape of theconditional distribution of the response.A practical estimator of the parameter can be obtained by looking for the solution of the following expectileestimating equations S I ( β τ ) = n X i =1 X i T Ψ τ ( y i − X i β τ ) h y i − X i β τ i = , (13)where Ψ τ ( y i − X i β τ ) = diag (cid:16) ψ τ ( y i − x i T β τ ) , . . . , ψ τ ( y im i − x im i T β τ ) (cid:17) . The resulting estimator b β Iτ canalso be derived as the minimizer of the following objective function1 N n X i =1 m i X t =1 ρ τ (cid:16) y it − x it T β τ (cid:17) (14)over β τ ∈ R p . The explicit form of the resulting estimator b β Iτ is similar to (7).When τ = 0 . , the estimator b β Iτ corresponds to the GEE estimator introduced by Liang and Zeger (1986)with an independent working correlation between observations from the same subject. This fact is exploitedto extend the GEE to the generalized expectile estimating equation (GEEE).The GEEE method models the underlying correlation structure from the same subject by formally includinga hypothesized structure with the within-subject correlation. For a fixed τ, the GEEE estimator b β τ is derivedby solving the following GEEE equations S ( β τ ) = n X i =1 X i T V − iτ Ψ τ ( y i − X i β τ ) h y i − X i β τ i = , (15)where V iτ is a working covariance matrix represented as V iτ = σ τ A iτ R i ( α τ ) A iτ , (16)and σ τ is the nuisance parameter. A iτ is the m i × m i diagonal matrix with the variance function ν ( µ i ) asdiagonal elements and R i ( α τ ) as the working correlation matrix.The working correlation matrix R i ( α τ ) describes the correlation pattern of within-subject observationswith the K × α τ . Liang and Zeger (1986) proposed several types of working correlationstructures (independent, exchangeable, autoregressive, unstructured, etc.) for the case τ = 0 . . These workingcorrelations are adapted and extended to the GEEE approach. The extension of some of the most commonand popular ones are presented below.The GEEE independent working correlation structure is the simplest form of working correlation withidentity matrix and is the structure assumed by the expectile estimating equations model (13). The GEEEexchangeable structure is a simple extension of the independence working correlation. It assumes a commoncorrelation, ρ tsτ = α τ , ∀ t = s, between any pair of measurements. The GEEE AR1 structure correlationdefines the correlation of a pair of observations as a decreasing function of their distance in time, ρ tsτ = α | t − s | τ . This structure assigns the highest correlation to adjacent pairs of observation and the lowest correlationto distant pairs. The GEEE unstructured, as its name suggests, imposes no structure to the correlationmatrix and defines the correlations of all pairs of measurements differently, without any explicit pattern, ρ tsτ = α tsτ , ∀ t = s. All these types of working correlation are usually unknown and must be estimated. They are estimated inthe iterative fitting process using the current value of the parameter vector. Indeed, the estimators can be5omputed as iterated weighted least squares estimators. The estimation algorithm for the GEEE exchangeableworking correlation can be summarized as the following stepwise procedure.
Algorithm: The GEEE algorithmStep
1. Let e β (0) τ = b β Iτ , the estimator obtained from (13). Step
2. Given e β ( r − τ from the r − b σ r ) τ ← N − p n X i =1 m i X t =1 ψ τ ( b ε itτ ) b ε itτ , b α ( r ) τ ← N − p ) b σ r ) τ n X i =1 m i X t
1) and b ε itτ = y it − x it T b β ( r − τ . Step
3. Update b β ( r ) τ by b β ( r ) τ ← b β ( r − τ + h n X i =1 X i T V − iτ ( b α ( r − τ ) Ψ τ ( b β ( r − τ ) X i i − S ( b α ( r − τ , b β ( r − τ ) , where Ψ τ ( b β ( r − τ ) = Ψ τ ( y i − X i b β ( r − τ ) . Step
4. Repeat the above iteration,
Steps α which is either a scalar or a vector, depending on the type of correlation. For example, for aGEEE autoregressive AR1 working correlation structure, the scalar parameter α τ is estimated by b α τ = 1( N − p ) σ τ n X i =1 X t Assume that b β I τ is the solution of the estimating function (13) and suppose the data aregenerated by model (9), and that conditions A1 - A3 are satisfied. If E | ψ τ k ( ε itτ k ) | ν < ∆ and E | ε itτ k | ν < ∆ for some ν > and ∆ > , then for every fixed sequence of expectiles τ = ( τ , . . . , τ q ) √ N (cid:0)b β I τ − β τ (cid:1) d −→ N (cid:18) , D − I ( τ ) D I ( τ ) D − I ( τ ) (cid:19) . In order to use this new estimator b β I τ to make inference, an estimator of its VC-matrix is presented in Theorem 2 . This will make it possible to construct large sample confidence intervals or hypotheses tests.This estimator is a generalization of the robust VC estimator proposed by White (1980) and used in, amongother things, multilevel analysis (Liang and Zeger 1986). This estimator inherits the same property namely inthat it takes into account the within-subject-correlation and the heteroscedasticity between subjects. In sum,the proposed VC-matrix estimator is a commonly advocated covariance matrix estimator for longitudinaldata. Let, b D I ( τ ) = N − n X i =1 ( W ⊗ X i ) T Ψ τ ( b ε i τ )( I q ⊗ X i ) , b D I ( τ ) = N − n X i =1 ( W ⊗ X i ) T b Σ i τ ( W ⊗ X i )where b Σ i τ = Ψ τ ( b ε i τ ) b ε i τ b ε i τ T Ψ τ ( b ε i τ ) and b ε i τ is obtained by replacing β τ with b β I τ in the expression of ε i τ . Then, we have Theorem 2 . Theorem 2. Suppose the data are generated by model (9) and that conditions A1 - A3 are satisfied. If E | ψ τ k ( b ε itτ k ) | ν < ∆ and E | ε itτ k | ν < ∆ for some ν > and ∆ > , then for every fixed sequence ofexpectiles τ = ( τ , . . . , τ q ) b D − I ( τ ) b D I ( τ ) b D − I ( τ ) p −→ D − I ( τ ) D I ( τ ) D − I ( τ ) . .2 Asymptotic properties for the general GEEE estimator After presenting the asymptotic properties of the GEEE-independent working correlation estimator, thissubSection presents the asymptotic properties of the GEEE-estimator for a general working correlation.Assume that B1 . The data { ( y i , X i ) } ni =1 are independent across i and Var h Ψ τ ( ε i τ ) ε i τ i = Σ i τ . B2 . The limiting forms of the following matrices are positive definite D ( τ ) = lim n →∞ N − n X i =1 ( W ⊗ X i ) T V − i τ E [ Ψ τ ( ε i τ )]( I q ⊗ X i ) , D ( τ ) = lim n →∞ N − n X i =1 ( W ⊗ X i ) T V − i τ Σ i τ V − i τ ( W ⊗ X i ) . B3 . max ≤ i ≤ n ≤ t ≤ m i k x it k < M. The following Theorem derives the asymptotic properties of the GEEE estimator with a general workingcorrelation, under the above conditions. Theorem 3. Suppose the data are generated by model (9) and that conditions B1 - B3 are satisfied. If E | ψ τ k ( ε itτ k ) | ν < ∆ and E | ε itτ k | ν < ∆ for some ν > and ∆ > , then for every fixed sequence ofexpectiles τ = ( τ , . . . , τ q ) √ N (cid:0)b β τ − β τ (cid:1) d −→ N (cid:18) , D − ( τ ) D ( τ ) D − ( τ ) (cid:19) . In the same way as with the GEEE-independent working correlation estimator, the following Theorem 4 proposes an estimator of the VC-matrix of estimator b β τ . Consider b V i τ to be a consistent estimator of V i τ . Then, under the above conditions, Theorem 4 is stated as follows Theorem 4. Suppose the data are generated by model (9) and that conditions B1 - B3 are satisfied. Assume E | ψ τ k ( b ε itτ k ) | ν < ∆ and E | ε itτ k | ν < ∆ for some ν > and ∆ > . Then for every fixed sequence ofexpectiles τ = ( τ , . . . , τ q ) b D − ( τ ) b D ( τ ) b D − ( τ ) p −→ D − ( τ ) D ( τ ) D − ( τ ) where b D ( τ ) = N − n X i =1 ( W ⊗ X i ) T b V − i τ Ψ τ ( b ε i τ )( I q ⊗ X i ) , b D ( τ ) = N − n X i =1 ( W ⊗ X i ) T b V − i τ b Σ i τ b V − i τ ( W ⊗ X i ) , and b Σ i τ = Ψ τ ( b ε i τ ) b ε i τ b ε i τ T Ψ τ ( b ε i τ ) . In this Section, the small sample performance of the estimators is evaluated through extensive simulationstudies. The random samples are generated from the following linear model ( M γ ) :9 it = β + x it β + (1 + γx it ) ε it , i = 1 , . . . , n and t = 1 , . . . , m i . (18)Two versions of model (18) are considered with respect to the parameter γ ∈ { , / } : a location-shift model( M ) corresponding to γ = 0 , which helps assess the performance of the estimators for an homoscedasticscenario, and a location-scale-shift model ( M / ) corresponding to γ = 1 / , serving to assess the performanceof the estimators in the presence of heteroscedasticity.The corresponding GEEE model of ( M ) is µ τ ( y it ) = β τ + x it β where β τ = β + µ τ ( ε it ) , so that only theintercept term varies with τ and the expectiles functions are parallel lines. The GEEE model related to the( M / ) model is µ τ ( y it ) = β τ + x it β τ where β τ = β + µ τ ( ε it ) and β τ = β + γµ τ ( ε it ) . Therefore, in thepresence of heteroscedasticity both the intercept and the slope vary with τ. We generate the regressor x it from a Gaussian distribution in ( M ) and from a Chi-square distribution in( M / ) and set the parameters β and β to 0. In order to allow for simulation of dependent errors withdifferent marginal distributions, we first simulate dependent uniform margins from a Gaussian copula with anAR1 correlation structure. We then generate the dependent random errors as quantiles of the uniform marginsfrom three distinct marginal distributions: Normal, Student with three degrees of freedom and Chi-squarewith three degrees of freedom. We also centered the random errors on the τ -th expectile. Specifically, wegenerate the data as follows1. Generate x it from a Gaussian distribution in ( M ) and from a Chi-square in ( M / );2. Generate a uniform sample: ( u , . . . , u m i ) from a Gaussian Copula with AR1 correlation structure;3. For t = 1 , . . . , m i , generate the dependent random error ε it = F − ( u it ) , where F ( . ) is one of the threemarginal distributions: Gaussian, Student or Chi-square distribution;4. Center the random error: ε it = ε it − µ τ ( ε it );5. Generate the final sample: y it = β + x it β + (1 + γx it ) ε it . We used three different values for the AR1 correlation structure: low ρ = 0 . , medium ρ = 0 . , and high ρ = 0 . n ∈ { , } . Finally,for the number of repeated measurements m i , a balanced design with m i = 4 and an unbalanced design arestudied.In the unbalanced design, m i is an integer number randomly generated between 3 and 8 with equal probability.The extensive simulation is carried out with 400 replications for each parameter-combination scenario. In eachscenario, the focus is on the effect of the regressor, x it , at the expectiles τ ∈ { . , . , . } . All computationsare implemented in R code language (R Core Team 2018).The results of the GEEE regression are analyzed using four different and popular working correlation matrices:independence, exchangeable, AR1 and unstructured correlation. The average bias (Bias) and relative efficiency(EFF) of the estimates are reported for the measurement of the quality of the different related estimators.The standard deviation (SD) of the 400 parameter estimates is used as a benchmark to evaluate the averageasymptotic standard errors (SE).We use the quasi-likelihood criterion (QIC) as a model-selection criteria to choose among the different workingcorrelation structures. The QIC is a criteria developed by Pan (2001) for model selection and selection ofworking correlation structures. The QIC is a modification of the Akaike Information Criterion (AIC) for theGEE model. In our case, the statistic is defined as QIC ( R ) = 12 n X i =1 m i X t =1 b ε it b σ + 2 Trace (cid:16) b Ω I b V R ( b β ) (cid:17) , (19)where b V R ( b β ) is the robust covariance estimate and b Ω I is the inverse of the covariance estimate under theindependent working correlation evaluated at β τ ( R ) , the parameter estimate under the working correlationof interest. 10o compare with the quantile regression approach, the simulation results of the linear quantile mixed model(lqmm) (Geraci and Bottai 2007,Geraci and Bottai (2014)) were reported. The lqmm is a conditional quantileregression model with random effects parameters included to account for the within-subject dependence. Thechoice of the lqmm is motivated by the fact that the linear mixed model (lmm) estimate is equivalent to theexchangeable correlation estimate in the linear Gaussian setting, when τ = 0 . τ = (0 . , . , . 67) correspond to theGaussian expectiles of τ = (0 . , . , . . For the sake of brevity, we present in this paper only the simulation results for the normal distribution. Thesimulation results for the other distributions (Student and Chi-Square) are in the supplementary materialI . We also published these results and the codes on our GitHub repository (https://github.com/AmBarry/expectgee). Table 1 and Table 2 report the Bias and EFF results, respectively, for the ( M ) and ( M / ) models whenthe error follows a multivariate normal distribution with an AR1 correlation structure and ρ ∈ { . , . , . } . Overall, the estimation biases are all very close to 0 for the location-shift and location-scale-shift scenarios.The bias of the three estimators is comparable for the three values 0 . , . . ρ. The Un and AR1estimators do slightly better than the Ind estimator in term of efficiency, particularly when the correlation ishigher ρ ∈ (0 . , . . The Un estimator do much better in general than the other three estimators (Ind, Exc,AR1).To evaluate the asymptotic standard error (SE), we use the standard deviation (SD) as a benchmark. Theresults are presented in Table 3 and Table 4 when the error follows a normal distribution respectively forthe ( M ) and ( M / ) models. We observe that the values of SD and SE decrease as n becomes large. Ingeneral, the values of SD and SE are identical for each of the estimators. This identity is more pronouncedin the case of the independent and unstructured working correlation. Similar performances are observedwhen the error is generated by a Student or a Chi-Square distribution. These results can be found in the supplementary material I .Overall, the different estimators are efficient and have small biases regardless of the correlation structure.Hence our results confirm that the GEEE method yields a consistent and highly efficient estimator even witha misspecification of the true covariance structure (AR1). Table 5 presents the QIC results of the different correlation structures with respect to the balanced/unbalanceddata, the M /M / model and the sample sizes n ∈ { , } . The QIC is most likely to correctly select theAR1 structure from the four given correlation structures in the M scenario, particularly for the unbalanceddata. In the M / scenario, the QIC is most likely to select either the AR1 or the Un structure. This lastresult is unexpected but is not surprising. Similar results have been reported in (Jang 2011).As in Pan’s paper (Jang 2011), many did not include the unrestricted structure correlation as a candidate inthe evaluation of the QIC or other criteria (Jang 2011). But when included, the results showed that the QICwas strongly biased towards selecting the unrestricted structure. Please, see (Jang 2011) and the referencetherein for further details.The last Tables 6-7 report the simulation results (Bias and RMSE) of the lqmm estimator and the GEEEestimator with exchangeable working correlation. The results show that both methods are competitive interm of Bias and RMSE. 11 able 1: Bias and relative efficiency of GEEE estimator with different correlation structures at 3 percentileswith ρ ∈ (0 . , . , . , and ε ∼ N (0 , 1) under a location-shift scenario. m = 4 m ∼ U (3 , τ ρ GEEE Bias EFF Bias EFF Bias EFF Bias EFF0.25 0.1 Ind -0.0006 1.000 -0.0003 1.000 0.0006 1.000 -0.0007 1.000AR1 -0.0006 1.000 -0.0003 1.010 0.0007 1.000 -0.0006 1.000Un -0.0003 1.170 -0.0001 1.099 0.0006 1.254 -0.0006 1.1670.5 Ind -0.0006 1.000 0.0000 1.000 0.0004 1.000 -0.0009 1.000AR1 -0.0006 1.007 0.0000 1.000 0.0005 1.008 -0.0008 1.000Un -0.0006 1.142 0.0001 1.083 0.0002 1.233 -0.0008 1.1400.9 Ind -0.0006 1.000 0.0001 1.000 0.0001 1.000 -0.0010 1.000AR1 -0.0006 1.007 0.0001 1.000 0.0001 1.000 -0.0010 1.000Un -0.0007 1.164 0.0002 1.089 -0.0002 1.272 -0.0010 1.1670.50 0.1 Ind 0.0001 1.000 -0.0007 1.000 -0.0002 1.000 0.0002 1.000AR1 0.0001 1.079 -0.0007 1.070 -0.0002 1.081 0.0002 1.089Un -0.0003 1.664 -0.0002 1.600 -0.0003 1.790 0.0004 1.7220.5 Ind 0.0000 1.000 -0.0005 1.000 -0.0003 1.000 0.0002 1.000AR1 0.0001 1.104 -0.0005 1.105 -0.0003 1.118 0.0002 1.116Un -0.0002 1.687 -0.0001 1.653 -0.0004 1.807 0.0003 1.7330.9 Ind 0.0000 1.000 -0.0003 1.000 -0.0003 1.000 0.0001 1.000AR1 -0.0001 1.080 -0.0003 1.070 -0.0004 1.080 0.0001 1.067Un -0.0003 1.681 0.0000 1.590 -0.0007 1.784 0.0001 1.6780.75 0.1 Ind -0.0009 1.000 0.0003 1.000 0.0003 1.000 0.0001 1.000AR1 -0.0005 1.277 0.0002 1.265 0.0004 1.294 0.0000 1.284Un -0.0007 3.328 -0.0001 3.316 0.0001 3.698 0.0005 3.6820.5 Ind -0.0009 1.000 0.0004 1.000 0.0005 1.000 0.0001 1.000AR1 -0.0005 1.376 0.0003 1.362 0.0005 1.408 0.0000 1.405Un -0.0001 3.391 0.0001 3.351 0.0001 3.783 0.0002 3.7620.9 Ind -0.0008 1.000 0.0006 1.000 0.0006 1.000 0.0001 1.000AR1 -0.0006 1.273 0.0005 1.265 0.0005 1.298 0.0000 1.284Un -0.0007 3.367 0.0000 3.306 0.0005 3.823 -0.0001 3.67012 able 2: Bias and relative efficiency of GEEE estimator with different correlation structures at 3 percentileswith ρ ∈ (0 . , . , . , and ε ∼ N (0 , 1) under a location-scale-shift scenario. m = 4 m ∼ U (3 , τ ρ GEEE Bias EFF Bias EFF Bias EFF Bias EFF0.25 0.1 Ind 0.0013 1.000 0.0005 1.000 0.0007 1.000 0.0006 1.000AR1 -0.0014 1.761 -0.0078 1.611 0.0074 2.272 -0.0101 1.717Un 0.0035 2.608 0.0006 2.365 0.0048 2.693 -0.0002 2.9740.5 Ind -0.0006 1.000 -0.0001 1.000 0.0001 1.000 0.0003 1.000AR1 0.0040 1.853 -0.0003 1.484 0.0107 2.806 0.0023 1.671Un 0.0021 2.338 0.0002 2.081 0.0017 2.352 -0.0008 2.4700.9 Ind -0.0025 1.000 -0.0009 1.000 -0.0008 1.000 0.0002 1.000AR1 0.0120 1.830 0.0083 1.739 0.0128 3.492 0.0080 17.936Un 0.0009 2.752 0.0007 2.636 -0.0009 3.020 -0.0017 3.1670.50 0.1 Ind -0.0022 1.000 0.0012 1.000 0.0015 1.000 0.0003 1.000AR1 -0.0039 1.758 -0.0023 1.780 0.0045 3.080 -0.0005 2.646Un 0.0042 2.390 0.0021 2.304 -0.0011 3.219 -0.0008 2.9940.5 Ind -0.0031 1.000 0.0001 1.000 -0.0004 1.000 -0.0001 1.000AR1 0.0044 1.678 0.0028 1.699 0.0094 3.172 0.0054 2.766Un 0.0034 2.075 0.0017 2.030 0.0006 2.778 0.0002 2.4810.9 Ind -0.0040 1.000 -0.0009 1.000 -0.0020 1.000 -0.0006 1.000AR1 0.0128 1.813 0.0084 1.836 0.0195 3.217 0.0131 2.726Un 0.0026 2.461 0.0025 2.591 0.0029 3.429 0.0008 3.1400.75 0.1 Ind 0.0000 1.000 0.0024 1.000 -0.0010 1.000 0.0013 1.000AR1 0.0002 2.312 0.0032 2.207 -0.0165 3.860 -0.0049 3.845Un -0.0009 2.192 0.0017 2.151 0.0000 2.634 -0.0019 2.8050.5 Ind -0.0015 1.000 0.0009 1.000 -0.0023 1.000 0.0006 1.000AR1 -0.0019 2.150 0.0009 2.114 -0.0197 3.792 -0.0038 3.700Un -0.0029 1.850 -0.0006 1.892 -0.0032 2.248 -0.0022 2.1590.9 Ind -0.0024 1.000 -0.0006 1.000 -0.0039 1.000 -0.0001 1.000AR1 -0.0042 2.210 -0.0001 2.251 -0.0235 3.785 -0.0027 3.701Un -0.0079 1.931 -0.0013 2.328 -0.0069 2.544 -0.0044 2.49413 able 3: Standard deviation and asymptotic standard errors of the GEEE estimator with different correlationstructures at 3 percentiles with ρ ∈ (0 . , . , . , and ε ∼ N (0 , 1) under a location-shift scenario. m = 4 m ∼ U (3 , τ ρ GEEE SD SE SD SE SD SE SD SE0.25 0.1 Ind 0.014 0.014 0.010 0.010 0.014 0.013 0.009 0.009AR1 0.014 0.014 0.010 0.010 0.014 0.013 0.009 0.009Un 0.014 0.016 0.010 0.011 0.014 0.016 0.009 0.0100.5 Ind 0.013 0.013 0.010 0.010 0.013 0.012 0.009 0.009AR1 0.013 0.014 0.010 0.010 0.013 0.012 0.009 0.009Un 0.014 0.015 0.010 0.010 0.013 0.015 0.009 0.0100.9 Ind 0.014 0.014 0.010 0.010 0.014 0.012 0.009 0.009AR1 0.014 0.014 0.010 0.010 0.014 0.012 0.009 0.009Un 0.015 0.016 0.010 0.011 0.014 0.016 0.009 0.0100.50 0.1 Ind 0.014 0.014 0.010 0.010 0.012 0.012 0.009 0.009AR1 0.013 0.015 0.009 0.011 0.012 0.013 0.008 0.010Un 0.013 0.023 0.008 0.016 0.011 0.022 0.008 0.0160.5 Ind 0.014 0.013 0.010 0.010 0.012 0.012 0.009 0.009AR1 0.012 0.015 0.009 0.010 0.011 0.013 0.008 0.010Un 0.012 0.023 0.008 0.016 0.010 0.022 0.007 0.0150.9 Ind 0.014 0.014 0.010 0.010 0.012 0.012 0.009 0.009AR1 0.013 0.015 0.010 0.011 0.012 0.014 0.008 0.010Un 0.013 0.023 0.009 0.016 0.011 0.022 0.008 0.0150.75 0.1 Ind 0.014 0.014 0.010 0.010 0.014 0.013 0.009 0.009AR1 0.011 0.018 0.008 0.012 0.011 0.016 0.007 0.011Un 0.011 0.046 0.007 0.032 0.009 0.047 0.006 0.0320.5 Ind 0.013 0.013 0.009 0.009 0.013 0.012 0.009 0.008AR1 0.010 0.018 0.007 0.013 0.010 0.017 0.006 0.012Un 0.008 0.045 0.005 0.032 0.006 0.045 0.004 0.0320.9 Ind 0.014 0.014 0.010 0.010 0.014 0.012 0.009 0.009AR1 0.011 0.018 0.008 0.012 0.011 0.016 0.007 0.011Un 0.011 0.047 0.006 0.032 0.009 0.047 0.006 0.03214 able 4: Standard deviation and asymptotic standard errors of the GEEE estimator with different correlationstructures at 3 percentiles with ρ ∈ (0 . , . , . , and ε ∼ N (0 , 1) under a location-scale-shift scenario. m = 4 m ∼ U (3 , τ ρ GEEE SD SE SD SE SD SE SD SE0.25 0.1 Ind 0.026 0.021 0.019 0.017 0.022 0.020 0.017 0.015AR1 0.062 0.037 0.041 0.027 0.253 0.046 0.052 0.026Un 0.048 0.054 0.033 0.040 0.043 0.054 0.032 0.0450.5 Ind 0.024 0.020 0.018 0.016 0.022 0.020 0.016 0.015AR1 0.060 0.038 0.038 0.024 0.241 0.055 0.048 0.025Un 0.042 0.048 0.030 0.034 0.036 0.046 0.027 0.0370.9 Ind 0.024 0.021 0.019 0.016 0.024 0.020 0.018 0.016AR1 0.085 0.038 0.050 0.029 0.240 0.069 0.075 0.280Un 0.055 0.057 0.039 0.044 0.043 0.060 0.035 0.0490.50 0.1 Ind 0.027 0.023 0.019 0.017 0.024 0.020 0.017 0.016AR1 0.070 0.041 0.045 0.030 0.237 0.062 0.068 0.043Un 0.053 0.055 0.036 0.039 0.059 0.065 0.038 0.0480.5 Ind 0.025 0.023 0.019 0.017 0.023 0.020 0.017 0.015AR1 0.064 0.038 0.040 0.028 0.239 0.063 0.063 0.043Un 0.045 0.047 0.032 0.034 0.050 0.055 0.028 0.0380.9 Ind 0.025 0.023 0.020 0.017 0.024 0.020 0.018 0.016AR1 0.095 0.042 0.050 0.031 0.278 0.065 0.070 0.043Un 0.060 0.057 0.042 0.044 0.082 0.070 0.036 0.0490.75 0.1 Ind 0.028 0.024 0.020 0.018 0.028 0.024 0.020 0.017AR1 0.068 0.056 0.045 0.040 0.295 0.091 0.075 0.067Un 0.062 0.053 0.038 0.038 0.059 0.062 0.042 0.0490.5 Ind 0.028 0.024 0.020 0.018 0.026 0.023 0.019 0.017AR1 0.057 0.052 0.039 0.037 0.296 0.086 0.065 0.063Un 0.048 0.044 0.033 0.033 0.043 0.051 0.030 0.0370.9 Ind 0.030 0.025 0.020 0.018 0.027 0.023 0.019 0.017AR1 0.065 0.055 0.049 0.041 0.309 0.086 0.072 0.064Un 0.053 0.048 0.042 0.043 0.048 0.058 0.033 0.04315 able 5: Total of the frequency of the working correlation matrix selected by QIC for the different correlations ρ ∈ { . , . , . } from 1200 independent replications. The true correlation matrix is AR1.Balanced Data Unbalanced Data n = 50 n = 100 n = 50 n = 100Ind Exc AR1 Un Ind Exc AR1 Un Ind Exc AR1 Un Ind Exc AR1 UnLocation-shift model N 43 82 550 525 41 61 581 517 46 99 759 296 32 49 665 454 T 98 168 560 374 101 142 587 370 96 169 710 225 70 135 702 293 χ 69 138 548 445 60 113 577 450 48 130 743 279 49 79 704 368Location-scale-shift model N 201 200 245 554 232 185 279 504 164 219 325 492 158 194 284 564 T 99 164 519 418 129 177 466 428 125 143 474 458 124 172 442 462 χ 51 115 455 579 41 82 458 619 42 112 468 578 53 84 485 57816 able 6: Bias and RMSE of the GEEE estimator with exchangeable working correlation and the lqmm estimatorat 3 percentiles when ρ ∈ (0 . , . , . 9) for a balanced panel with ε ∼ N (0 , .n = 50 n = 100Bias RMSE Bias RMSE ρ τ Exc lqmm Exc lqmm Exc lqmm Exc lqmmLocation-shift model0.1 τ -0.0005 -0.0008 0.0142 0.0154 -0.0002 -0.0001 0.0104 0.0121 τ -0.0006 -0.0008 0.0132 0.0147 0.0000 0.0000 0.0099 0.0116 τ -0.0006 -0.0008 0.0140 0.0161 0.0001 0.0001 0.0102 0.01150.5 τ τ τ τ τ τ τ -0.0008 -0.0003 0.0025 0.0243 -0.0007 -0.0016 0.0020 0.0168 τ -0.0001 0.0004 0.0023 0.0189 0.0000 -0.0010 0.0017 0.0150 τ τ -0.0019 0.0027 0.0033 0.0263 -0.0018 0.0023 0.0026 0.0189 τ -0.0002 0.0017 0.0026 0.0228 -0.0001 0.0011 0.0018 0.0162 τ τ -0.0018 0.0155 0.0030 0.0322 -0.0017 0.0149 0.0023 0.0249 τ -0.0001 -0.0001 0.0023 0.0242 0.0000 0.0001 0.0016 0.0174 τ able 7: Bias and RMSE of the GEEE estimator with exchangeable working correlation and the lqmm estimatorat 3 percentiles when ρ ∈ (0 . , . , . 9) for an unbalanced panel with ε ∼ N (0 , .n = 50 n = 100Bias RMSE Bias RMSE ρ τ Exc lqmm Exc lqmm Exc lqmm Exc lqmmLocation-shift model τ τ τ τ -0.0002 -0.0003 0.0114 0.0135 0.0002 0.0001 0.0084 0.0098 τ -0.0003 -0.0006 0.0107 0.0127 0.0002 0.0001 0.0079 0.0093 τ -0.0004 0.0004 0.0115 0.0136 0.0002 0.0000 0.0082 0.01000.9 τ τ τ τ -0.0007 0.0007 0.0023 0.0218 -0.0008 -0.0008 0.0018 0.0152 τ τ τ -0.0015 0.0019 0.0027 0.0243 -0.0015 0.0001 0.0022 0.0171 τ -0.0001 0.0005 0.0021 0.0215 0.0000 -0.0005 0.0015 0.0146 τ τ -0.0017 0.0134 0.0028 0.0303 -0.0017 0.0161 0.0023 0.0251 τ τ Application In this Section, the proposed method is applied to the repeated measurements labor pain dataset previouslyreported by Davis (1991). It is a commonly used dataset in biostatistics, and used several times in thequantile regression framework (Jung 1996, Geraci and Bottai (2007), Lu and Fan (2015)). The dataset comesfrom a clinical trial on the effectiveness of a medication for labor pain relief for 83 women in labor. Thetreatment group (43 women) and the placebo group (40 women) were randomly assigned. The responsevariable is a self-reported score measured every 30 min on a 100-mm line, where 0 means no pain and 100means extreme pain. A nearly monotone pattern of missing data was found for the response variable, andthe maximum number of measurements per woman was six. Figure 1 shows the box-plot of the responsevariable for all women by treatment group. At first glance, we can determine that the response variable isnon-normal. We also observe the evolution of the mean and the median over time. The objective is to studythe effect of medication on the self-reported pain score from the following equation y it = β + β R i + β T it + β R i T it + ε it , (20)where y it is the t -th measure of the pain for the i -th subject. R i is the treatment variable with value 0 forplacebo and 1 for treatment, and T it is the measurement time divided by 30min. The corresponding GEEEequation, for a fixed τ µ τ ( y it ) = β τ + β R i + β T it + β R i T it , (21)with fourth working correlation (Ind, Exc, AR1, Un), was estimated for three asymmetric points,(0 . , . , . . Table 8 presents the results of the estimated parameters, their standard errors, as wellas their 95% confidence intervals. It is observed that the different GEEE methods produce comparableestimates. The estimated parameter b β is not significant in any of the models and the estimated parameter b β is not significant except for the percentile τ = 0 . . This indicates that the baseline pain does not differsignificantly between the two groups. The estimated parameters b β and b β are significant at 5% level for theGEEE methods, except for b β of the GEEE with an Exc working correlation. This means that time and itsinteraction with treatment affect the amount of pain. To investigate the effect of the drug on pain over time,we study the evolution of this difference µ τ ( y it | R i = 1) − µ τ ( y it | R i = 0) = β + β T it , for which the result is presented in Figure 2 . We see that medication helps women relieve their labor pain.Indeed, the pain of women in the placebo group grows faster with time than that of the treated group, andthis is for all the GEEE methods and at different percentiles (0 . , . , . . Using the QIC measure to choose among the 4 working correlation structures leads us to select the Un( QIC = 2416 . QIC = 2416 . QIC = 2418 . QIC = 2419 . otal Medication Placebo30 60 90 120 150 180 30 60 90 120 150 180 30 60 90 120 150 180020406080100 Time Groups Labo r P a i n Figure 1: Box-plot of measured labor pain for all women in placebo and medication groups. The solid linesconnect the medians and the dashed lines connect the means.20 able 8: Parameters estimates (Est) with their standard errors (SE) and 95% confidence interval (CI) obtained from the GEEEindependent, exchangeable, AR1 and unstructured working correlation at three percentiles, τ = 0 . , . , . .τ = 0.25 0.5 0.75GEEE Est SE CI Est SE CI Est SE CIInd β β β β -9.65 2.12 (-13.80, -5.49) -9.58 2.03 (-13.56, -5.59) -7.32 2.22 (-11.67, -2.97)Exc β β β β -8.71 4.94 (-18.39, 0.96) -9.15 4.96 (-18.87, 0.56) -7.03 4.34 (-15.53, 1.47)AR1 β β β β -7.62 1.64 (-10.84, -4.40) -8.14 2.21 (-12.46, -3.81) -6.52 2.91 (-12.23, -0.82)Un β β β β -9.13 2.97 (-14.94, -3.32) -9.21 2.97 (-15.03, -3.39) -6.95 2.94 (-12.70, -1.19)21 Conclusion We combined weighted asymmetric least squares regression and generalized estimating equations to derive anew class of estimators: generalized expectile estimating equations estimators. This new GEEE class modelsthe underlying correlation structure from one subject by formally including a hypothesized structure with thewithin-subject correlation. In addition, this new model captures the heterogeneity of covariate effects andaccounts for unobserved heterogeneity. We also showed how to extend and adapt some of the most commonand popular GEE working correlations.We derived the asymptotic properties of this new estimator and proposed a robust estimator of its variance-covariance matrix. The results of the exhaustive simulations displayed its favorable qualities under variousscenarios and its advantages in relation to existing methods. The QIC is most likely to select the correctworking correlation (AR1) among the four working correlation structures used in the simulation. Finally, wefit the GEEE model to the labor pain data. The results revealed a strong association of treatment and timeon the labor pain of the two groups of women. This result is consistent with the results obtained in otherstudies (Lu and Fan 2015, Leng and Zhang (2014)). In addition, the results show that the heterogeneity ofthe evolution of pain according to time depends on whether one is in the center or on the left/right of the tailof the distribution response. The application of the QIC criterion to choose between the four correlationstructures leads to the selection of either the AR1 or the Un working correlation structures.The proposed model opens the door to other avenues of research. Unlike the quantile regression model, theexpectile regression and the GEEE method will naturally generalize to the dichotomous or count data, or toother longitudinal models already used to estimate the effect of covariates on the average of the response.22 nd UnAr1 Exc50 100 150 50 100 150255075100255075100 time pa i n group PlaceboTreatment tau Figure 2: Representation of the estimated labor pain obtained from the GEEE method with the differentworking correlations (Ind, Exc, AR1, Un) at three percentiles, (0 . , . , . . APPENDIX Sketches of the proofs are provided below. A detailed version is available in the supplementary materialII . Proof of Theorem 1 . The proof of Theorem 1 is based on the following result from Hjort (2011). Corollary 1. (Basic Corollary of Hjort (2011) ) Consider A n ( s ) is convex and can be represented as s T V s + U n T s + C n + r n ( s ) , where V is symmetric and positive definite, U n is stochastically bounded, C n is arbitrary, and r n ( s ) goes to zero in probability for each s. Then α n the argmin of A n , is only o p (1) awayfrom β n = − V − U n , the argmin of s T V s + U n T s + C n . If also U n d −→ U then α n d −→ − V − U n . With this result, the proof of Theorem 1 is obtained by approximating the objective function R Nq using aconvex quadratic function with a unique minimum. This approximation is possible through the following Lemma : Lemma 1. Let r ( t ) = ρ τ ( u − t ) − ρ τ ( t ) + 2 tψ τ ( u ) u then r ( t ) = O ( t ) . Proof of Lemma 1 . Replacing the functions ρ τ ( · ) , and ψ τ ( · ) by their expression, the function r ( · ) is r ( t ) = | τ − ( u < t ) | ( u − t ) − | τ − ( u < | u + 2 | τ − ( u < | ut. According to the sign of t, we have: r ( t ) = (1 − τ ) t if u < < t (1 − τ )( u − t ) + τ t if 0 < u < tτ t if 0 < t < u if t > , and r ( t ) = (1 − τ ) t if u < t < τ − u − t ) + (1 − τ ) t if t < u < τ t if t < < u if t < . Hence we have r ( t ) = O ( t ) . The objective function R Nq is defined as R Nq ( δ ) = q X k =1 n X i =1 m i X t =1 w k (cid:26) ρ τ (cid:16) ε itτ k − x it T δ τ k / √ N (cid:17) − ρ τ (cid:0) ε itτ k (cid:1)(cid:27) (22)where δ = ( δ τ T , . . . , δ τ q T ) T is a pq × δ τ k T a p × ε itτ k = y it − µ itτ k and µ itτ k = x it T β τ k .R Nq is a convex function of δ and is minimized by b δ δ = δ τ ... δ τ q = √ N (cid:0)b β τ − β τ (cid:1) ... √ N (cid:0)b β τ q − β τ q (cid:1) . Now using Lemma 1 , we are able to split the objective function in two functions R (1) Nq and R (2) Nq R Nq ( δ ) ’ − √ N q X k =1 w k δ τ k T X T Ψ τ k ( ε τ k ) ε τ k + 1 N q X k =1 w k δ τ k T X T E [ Ψ τ k ( ε τ k )] Xδ τ k = R (1) Nq ( δ ) + R (2) Nq ( δ ) . Conditions A2 and A3 imply a Lindeberg condition and we have R (1) Nq ( δ ) = − √ N δ T ( W ⊗ X ) T Ψ τ ( ε τ ) ε τ d −→ − δ T B where B is a zero mean Gaussian vector with covariance matrix D I ( τ ) . For the second term to the right of R Nq we have, by condition A2 , R (2) Nq ( δ ) = N − δ T n X i =1 ( W ⊗ X i ) T E h Ψ τ ( ε i τ ) i ( I q ⊗ X i ) δ → δ T D I ( τ ) δ . Thus the limiting form of the objective function is R q ( δ ) = − δ T B + δ T D I ( τ ) δ where B is a zero mean Gaussian vector with covariance matrix D I ( τ ) . Application of Corollary 1 givesthe result of Theorem 1 . Proof of Theorem 2 . For the proof of Theorem 2 , we must show the convergence of b D I ( τ ) and b D I ( τ ) . The matrix b D I ( τ ) is a diagonal matrix of general term N − w k P ni =1 X i T Ψ τ k ( b ε iτ k ) X i , < k < q. Itsuffices to show the convergence of N − P ni =1 X i T Ψ τ k ( b ε iτ k ) X i to obtain that of b D I ( τ ) . Using the resultplim b β Iτ = β τ , the convergence of b D I ( τ ) p −→ D I ( τ ) follows by application of the Markov’s inequality andthe Markov law of large numbers (LLN).The matrix b D I ( τ ) is a block matrix of dimension pq × pq and of general term w k w j N − P ni =1 X i T b Σ iτ k τ j X i , Theorem 2 , the proof of Theorem 4 consists in showing the convergence of b D ( τ ) and b D ( τ ) . Thisis done following the same steps as Theorem 2 . 26 eferences Aigner, D.J., T. Amemiya, and D.J. Poirier. 1976. “ON the Estimation of Production Frontiers: MAXIMUMLikelihood Estimation of the Parameters of a Discontinuous Density Function.” International EconomicReview 17 (2): 377.Azzalini, A. 1985. “A Class of Distributions Which Includes the Normal Ones.” Scandinavian Journalof Statistics Proceedings of the Second Seattle Symposium in Biostatistics: Analysis of Correlated Data , edited by D. Y.Lin and P. J. Heagerty, 51–69. New York, NY: Springer New York. doi:10.1007/978-1-4419-9076-1_4.Davis, Charles S. 1991. “Semi-Parametric and Non-Parametric Methods for the Analysis of RepeatedMeasurements with Applications to Clinical Trials.” Statistics in Medicine 10 (12). Wiley SubscriptionServices, Inc., A Wiley Company: 1959–80. doi:10.1002/sim.4780101210.Diggle, P., P. Heagerty, K.Y. Liang, and S. Zeger. 2013. Analysis of Longitudinal Data . Oxford StatisticalScience Series. OUP Oxford. https://books.google.ca/books?id=zAiK-gWUqDUC.Efron, B. 1991. “REGRESSION Percentiles Using Asymmetric Squared Error Loss.” Statistica Sinica Statistical Modelling 13 (4): 317–22.Farooq, M., and I. Steinwart. 2017. “An Svm-Like Approach for Expectile Regression.” ComputationalStatistics and Data Analysis Longitudinal Data Analysis . Chapman& Hall/Crc Handbooks of Modern Statistical Methods. CRC Press. https://books.google.ca/books?id=zVBjCvQCoGQC.Fu, Liya, and You-Gan Wang. 2012. “Quantile Regression for Longitudinal Data with a Working CorrelationModel.” Computational Statistics & Data Analysis 56 (8): 2526–38. doi:https://doi.org/10.1016/j.csda.2012.02.005.Furlotte, Nicholas A., Eleazar Eskin, and Susana Eyheramendy. 2012. “Genome-Wide Association Mappingwith Longitudinal Data.” Genetic Epidemiology 36 (5): 463–71.Geraci, M., and M. Bottai. 2007. “Quantile Regression for Longitudinal Data Using the Asymmetric LaplaceDistribution.” Biostatistics Statistics and Computing 24 (3):461–79. doi:10.1007/s11222-013-9381-9.Hsiao, Cheng. 2007. “Panel Data Analysis—advantages and Challenges.” TEST 16 (1): 1–22.Jang, Mi. 2011. “Working Correlation Selection in Generalized Estimating Equations.” ProQuest DissertationsPublishing. http://search.proquest.com/docview/920555336/.Jiang, C., M. Jiang, Q. Xu, and X. Huang. 2017. “Expectile Regression Neural Network Model withApplications.” Neurocomputing Journal of the American StatisticalAssociation 91 (433). Taylor & Francis: 251–57.Kim, M., and S. Lee. 2016. “Nonlinear Expectile Regression with Application to Value-at-Riskand Expected Shortfall Estimation.” Computational Statistics and Data Analysis 94. Elsevier: 1–19.doi:10.1016/j.csda.2015.07.011.Kneib, Thomas. 2013a. “Beyond Mean Regression.” Statistical Modelling 13 (4): 275–303.27oi:10.1177/1471082X13494159.———. 2013b. “Rejoinder.” Statistical Modelling 13 (4): 373–85.Koenker, Roger, and Gilbert Bassett Jr. 1978. “Regression Quantiles.” Econometrica. Journal of theEconometric Society 46 (1): 33–50.Leng, Chenlei, and Weiping Zhang. 2014. “Smoothing Combined Estimating Equations in Quantile Regressionfor Longitudinal Data.” Statistics and Computing 24 (1). Boston: Springer US: 123–36.Liang, Kung-Yee, and Scott L. Zeger. 1986. “Longitudinal Data Analysis Using Generalized Linear Models.” Biometrika 73 (1): 13–22.Liu, Yufeng, and Yichao Wu. 2011. “Simultaneous Multiple Non-Crossing Quantile Regression EstimationUsing Kernel Constraints.” Journal of Nonparametric Statistics 23 (2): 415–37.Lu, Xiaoming, and Zhaozhi Fan. 2015. “Weighted Quantile Regression for Longitudinal Data.” ComputationalStatistics 30 (2): 569–92. doi:10.1007/s00180-014-0550-x.Majumdar, A., and D. Paul. 2016. “Zero Expectile Processes and Bayesian Spatial Regression.” Journal of Computational and Graphical Statistics 25 (3). American Statistical Association: 727–47.doi:10.1080/10618600.2015.1062014.Newey, Whitney K., and James L. Powell. 1987. “Asymmetric Least Squares Estimation and Testing.” Econometrica 55 (4): 819–47.Pan, Wei. 2001. “Akaike’s Information Criterion in Generalized Estimating Equations.” Biometrics 57 (1).Oxford, UK: Blackwell Publishing Ltd: 120–25.R Core Team. 2018. R: A Language and Environment for Statistical Computing . Vienna, Austria: RFoundation for Statistical Computing.Righi, M.B., Y. Yang, and P.S. Ceretta. 2014. “Nonparametric Expectile Regression for ConditionalAutoregressive Expected Shortfall Estimation.” Contemporary Studies in Economic and Financial Analysis 96. Emerald Group Publishing Ltd.: 83–95. doi:10.1108/S1569-375920140000096003.Rossi, Giuliano De, and Andrew Harvey. 2009. “Quantiles, Expectiles and Splines.” Journal of Econometrics 152 (2): 179–85. doi:http://dx.doi.org/10.1016/j.jeconom.2009.01.001.Schnabel, S.K., and P.H.C. Eilers. 2009. “Optimal Expectile Smoothing.” Computational Statistics and DataAnalysis 53 (12): 4168–77. doi:10.1016/j.csda.2009.05.002.———. 2013. “A Location-Scale Model for Non-Crossing Expectile Curves.” Stat Statistics and Computing 27 (1): 271–82. doi:10.1007/s11222-015-9621-2.Smith, Luke B., Montserrat Fuentes, Penny Gordon-Larsen, and Brian J. Reich. 2015. “Quantile Regressionfor Mixed Models with an Application to Examine Blood Pressure Trends in China.” Annals of AppliedStatistics Computational Statistics and DataAnalysis 56 (4): 755–67. doi:10.1016/j.csda.2010.11.015.Sobotka, Fabian, Göran Kauermann, Linda Schulze Waltrup, and Thomas Kneib. 2013. “On Confidence Inter-vals for Semiparametric Expectile Regression.” Statistics and Computing 23 (2): 135–48. doi:10.1007/s11222-011-9297-1.Waldmann, E., F. Sobotka, and T. Kneib. 2016. “Bayesian Regularisation in Geoadditive Expectile Regression.” Statistics and Computing . Springer New York LLC, 1–15. doi:10.1007/s11222-016-9703-9.Waltrup, L.S., F. Sobotka, T. Kneib, and G. Kauermann. 2015. “Expectile and Quantile Regression—David28nd Goliath?” Statistical Modelling 15 (5): 433–56. doi:10.1177/1471082X14561155.White, Halbert. 1980. “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test forHeteroskedasticity.” Econometrica 48 (4): 817–38. https://EconPapers.repec.org/RePEc:ecm:emetrp:v:48:y:1980:i:4:p:817-38.Xing, J.-J., and X.-Y. Qian. 2017. “Bayesian Expectile Regression with Asymmetric Normal Distri-bution.” Communications in Statistics - Theory and Methods 46 (9). Taylor; Francis Inc.: 4545–55.doi:10.1080/03610926.2015.1088030.Xu, Q., X. Liu, C. Jiang, and K. Yu. 2016. “Nonparametric Conditional Autoregressive Expectile Model viaNeural Network with Applications to Estimating Financial Risk.” Applied Stochastic Models in Business andIndustry 32 (6). John Wiley; Sons Ltd: 882–908. doi:10.1002/asmb.2212.Yang, Y., and H. Zou. 2015. “Nonparametric Multiple Expectile Regression via Er-Boost.”