Robust estimation for partial functional linear regression models based on FPCA and weighted composite quantile regression

1 Introduction

In the literature of functional data analysis, functional linear regression models (FLRMs) provide an efficient method to analyze the relationship between a functional predictor and a scalar response. Therefore, much effort has been devoted to studying its estimation and other relevant inference problems, see [1,2, 3,4,5, 6,7,8, 9,10,11, 12,13]. However, in practice, we often see that a scalar response is related not only to functional covariates but also to scalar covariates. Thus, Shin considered partial functional linear regression models (PFLRMs) to balance the flexibility of FLRMs and interpretation of classical linear regression models [14]. Generally, PFLRM has the following form:

where Y is a real-valued response variable defined on a probability space ( Ω , ℬ , P ) , Z is a p -dimensional random vector and θ = ( θ 1 , … , θ p ) T is an unknown parameter vector, { X ( t ) ; t ∈ T } is a zero mean, second-order stochastic process defined on ( Ω , ℬ , P ) with sample paths in H = L 2 ( T ) , the Hilbert space containing square integrable functions on T with inner product ⟨ x , y ⟩ = ∫ T x ( t ) y ( t ) d t , ∀ x , y ∈ H and norm ‖ x ‖ = ⟨ x , x ⟩ 1 / 2 , β ( t ) is an unknown slope function belonging to H , ε is the random error independent of Z and X , E ( ε ∣ Z , X ) = 0 . T denotes the transport operation throughout this paper. Without loss of generality, we further assume T = [ 0 , 1 ] .

The PFLRMs have been studied by many authors. For example, Shin [14] proposed the estimation method based on functional principal component analysis (FPCA) and investigated the asymptotic properties of the estimators for model (1). Zhou et al. [15] considered the polynomial spline method to estimate model (1) and established the asymptotic normality for the parameter vector and the global convergence rate for the slope function. Furthermore, Kong et al. [16] studied variable selection in model (1) with multiple functional and ultrahigh-dimensional scalar predictors. Ling et al. [17] developed the k-nearest-neighbors estimation for PFLRMs. Kong et al. [18] proposed the partial functional linear cox regression model for censored outcomes. Cao et al. [19] developed a two-step estimation procedure via maximizing the quasi likelihood function in the framework of generalized PFLRMs. Yuan and Zhang [20] investigated the hypothesis test of the parametric components in PFLRMs based on B-spline. Zhu et al. [21] considered estimation and testing problems for PFLRMs when the covariates of the non-functional linear component are measured with additive error. Li et al. [22] proposed two test statistics for series correlation in PFLRMs and derived their asymptotic distributions under the null hypothesis.

Quantile regression (QR) is an effective approach for functional data analysis. This is due to the fact that QR models compared to the mean regression models [23] offer more robust merits especially the observations contain some outliers. For relevant literature, see [24,25,26]. Hence, QR method may result in an arbitrarily small relative efficiency compared with the mean regression. To conquer this drawback of QR method, Zou and Yuan proposed the composite quantile regression (CQR) method for linear models and showed that the relative efficiency of the CQR estimator is greater than 70% regardless of the error distribution compared with the least square (LS) estimator [27]. In many cases, CQR provides more efficient estimation compared with a single QR or the LS estimator for non-normal error distributions. This method has been deeply investigated in the literature and widely applied in many statistical models. For further references about the CQR method on nonparametric/semi-parametric models see [28,29,30, 31,32]. Even though CQR has been well developed, theory and methodology of PFLRMs are rare. Du et al. [33] used the penalized CQR method to study variable selection for parametric part in the PFLRMs, Yu et al. [34] considered composite quantile estimation for the PFLRMs with errors from short-range dependent and strictly stationary linear processes. However, the aforementioned CQR method for PFLRMs is a sum of different QRs with equal weights. Intuitively, using equal weights is not optimal in general, though CQR enjoys great advantages in terms of estimation efficiency. Therefore, Jiang et al. [35] addressed the issue of statistical inference for nonlinear models with large-dimensional covariates by using the weighted composite quantile regression (WCQR) method. Moreover, they pointed out that their method was more robust than some existing methods. Guo et al. [36] investigated group SCAD penalized WCQR estimation for varying coefficient models with diverging number of parameters. Jiang et al. [37] developed a data-driven WCQR for heteroscedastic partially linear varying coefficient models. Similar conclusions of WCQR method have been further confirmed in [6,38,39]. Until now, the WCQR method has not been used in functional regression analysis. Consequently, motivated by this fact, we devote to extending the WCQR method to PFLRMs in this paper. We show that the resulting estimators for both parametric components and slope function are more efficient in the case of asymmetrical error distribution, and as asymptotically efficient as the corresponding LS estimators when the error distribution is normal.

It is well known that many smoothing methods have been proposed for functional regressions, and two popular approaches for estimating the slope function β ( ⋅ ) include the FPCA approach based on spectral decompositions of both the covariance of X ( t ) and its estimator (see [3,7,40]) and the basis function expansion method with both β ( ⋅ ) and X ( t ) approximated by basis functions such as polynomial spline functions in [15]. In this paper, we propose a WCQR estimation method for PFLRMs through weighted composite quantile loss function with principal component basis function approximations. The main contributions of this paper are threefold:

1. This is the first attempt to provide robust estimate via WCQR and establish the asymptotic properties of the proposed estimators for model (1).

2. We develop a data-driven weighting strategy that maximizes the efficiency of the WCQR estimators. Practically, the proposed method is more efficient than the CQR method in [33] and rank regression method proposed under asymmetrical error distributed data set in [38].

3. Theoretically, the proposed WCQR method generalizes the results of the composite quantile estimators in [8]. Particularly, our method reduces to that of Du et al. when the weights are all equal.

The rest of this paper is organized as follows. In Section 2, we describe the implementation details of the proposed algorithm and establish the theoretical properties of the estimators. In Section 3, we discuss the selection of weights and tuning parameters. Simulation study and a real data analysis are conducted in Sections 4 and 5 to evaluate the performance of our proposed approach, and then followed by a short conclusion in Section 6. All technical proofs are collected in the Appendix.

2 Method and main results

In partial functional linear QR model, for a given quantile level τ ∈ ( 0 , 1 ) , ( Z i , X i ( ⋅ ) , Y i ) , i = 1 , … , n is an independent and identically distributed (i.i.d.) sample from model (1), that is

where Q τ ( Y ∣ Z , X ( t ) ) is the τ th conditional quantile of Y given ( Z , X ( t ) ) , Y is a real-valued random variable defined on a probability space ( Ω , ℬ , P ) , Z is a p -dimensional vector of random variables with finite second moments, θ τ is a p × 1 coefficient vector to be estimated, β τ ( t ) is an unknown square integrable function on [ 0 , 1 ], and ε τ is a random error whose τ th quantile conditional on ( Z , X ( t ) ) being zero. For the convenience of presentation, we will omit τ from θ τ , β τ ( t ) and ε τ in model (2) wherever clear from the context, but we should keep in mind that those quantities are τ -specific.

Denote the covariance function of the process X ( ⋅ ) and its empirical version, respectively, as

The covariance function K defines a linear operator which maps a function f to K f given by ( K f ) ( u ) = ∫ K ( u , v ) f ( v ) d v . We assume that the linear operator with kernel K is positive definite. By Mercer’s theorem, then

where λ 1 > λ 2 > ⋯ > 0 and λ ˆ 1 ≥ λ ˆ 2 ≥ ⋯ ≥ λ ˆ n + 1 = ⋯ = 0 are, respectively, the ordered eigenvalue sequences of the linear operators with kernels K and K ˆ , and { ϕ j } and { ϕ ˆ j } are the corresponding orthonormal eigenfunction sequences. Obviously, the vectors { ϕ j } and { ϕ ˆ j } each forms an orthonormal basis in H . According to the Karhunen-Loèvere presentation, then

where ξ i and γ i are defined by

where ξ i is referred to as the i th functional principal component score of X ( t ) . It follows that ξ i are uncorrelated random variables with mean 0 and variance E ( ξ i 2 ) = λ i . For more details see [41]. Substituting (3) into model (2), the model (2) can be written as

Therefore, the regression model in (4) can be well approximated by

where m ≤ n is the truncation level that trades off approximation error against variability and typically diverges with n . Replace ϕ j by ϕ ˆ j for j = 1 , … , m , model (5) can be rewritten as

where U = { ⟨ X , ϕ ˆ j ⟩ } j = 1 , … , m , γ = ( γ 1 , … , γ m ) T . Thus, the estimation problem of the parameters and the slope function are transformed into the estimation problem of the parameter vectors θ and γ .

Note that the CQR method used the same weight for different QR models. Intuitively, it will be more effective if different weights are used. According to Jiang et al. [35], the WCQR procedure estimates θ and γ via minimizing the following loss function:

where ρ τ k ( r ) = τ k r − r I ( r < 0 ) , k = 1 , … , q , q is the number of quantiles, τ k ∈ ( 0 , 1 ) . Typically, we use equally spaced quantile position: τ k = k / ( 1 + q ) , ω = ( ω 1 , … , ω q ) T is a vector of weight satisfying ω k ≥ 0 . c τ k is 100 τ k % quantile of ε , c = ( c τ 1 , … , c τ q ) T . Let the coefficients γ ^ = ( γ ^ 1 , … , γ ^ m ) T and θ ^ = ( θ ^ 1 , … , θ ^ p ) T be the solution of equation (6), then the resulting estimator of the slope function is given by

3 Assumptions and asymptotic results

To establish the asymptotic properties of the proposed estimators, first, we introduce some notations. Let θ 0 and β 0 ( t ) be the true values of θ and β ( t ) , respectively, the notation ‖ ⋅ ‖ is the L 2 norm for a function or the Euclidean norm for a vector. Let C denote a generic constant that might assume different values at different places, a n ∼ b n means that a n / b n is bounded away from 0 and infinity as n → ∞ . The following technical assumptions are imposed.

1. The random function X ( ⋅ ) satisfies E ∣ ∣ X ( ⋅ ) ∣ ∣ 4 < ∞ .

2. For each j , E [ U j 4 ] ≤ C λ j 2 . For the eigenvalues λ j and Fourier coefficients γ j , we require that λ j − λ j + 1 ≥ C − 1 j − a − 1 and ∣ γ j ∣ ≤ C j − b for j > 1 , a > 1 and b > a / 2 + 1 .

3. The tuning parameter m satisfies m ∼ n 1 / ( a + 2 b ) .

4. The random vector Z has bounded fourth moment, E ‖ Z ‖ 4 < ∞ .

5. Z = η + ⟨ g , X ⟩ , where η = ( η 1 , … , η p ) T is zero-mean random variable, g = ( g 1 , … , g p ) T with g j ∈ L 2 ( [ 0 , 1 ] ) , j = 1 , … , p .

6. E [ η ] = 0 and E [ η η T ] = Σ . Furthermore, we need that Σ is a positive definite matrix.

7. ε has the distribution function F ( ⋅ ) with density f ( ⋅ ) . The density function f ( ⋅ ) of ε is positive and continuous at the τ k th quantiles c τ k .

4 Computational issues

From Theorem 3.1, we find that the asymptotic variance of θ ^ depends on ω only through σ 2 ( ω ) . Thus, the optimal choice of weights for maximizing efficiency of the estimator θ ^ is ω opt = arg min ω σ 2 ( ω ) . Based on [36], the optimal choice of weights is to solve the following quadratic optimization problem:

where f = ( f ( c τ 1 ) , … , f ( c τ q ) ) T and Ω is a q × q matrix with the ( k , k ′ ) element Ω k , k ′ = min ( τ k , τ k ′ ) − τ k τ k ′ .

It can be seen from (8) that the optimal weight vector is rather complicated and involves the density of the errors c k = F − 1 ( τ k ) and f ( c k ) , k = 1 , … , q . In practice, the error density f ( ⋅ ) is generally unknown. Following [36], we propose an estimation procedure as follows.

Step 1. We give the initial estimators θ ˜ and β ˜ ( t ) by minimizing the LS objective function and estimate σ 2 ( ⋅ ) by

Step 2. Compute ε ˜ i = Y i − Z i T θ ˜ − ∫ 0 1 X i ( t ) β ˜ ( t ) d t σ ˜ and then we can use kernel density estimation to estimate f ( ⋅ ) as follows:

Similar to Silverman [42], the bandwidth h is taken as

K ( ⋅ ) is always chosen as the Gaussian kernel. Here std( ⋅ ) and IQR( ⋅ ) denote the sample standard deviation and sample interquartile, respectively.

Step 3. Estimate f { F − 1 ( τ k ) } by f ˜ { F ˜ − 1 ( τ k ) } , where F ˜ − 1 ( τ k ) is the sample τ k quantile of { ε ˜ i , i = 1 , … , n } .

The proposed estimators involve the tuning parameter m , which needs to be selected by minimizing some selection criteria such as cross validation (CV), generalized cross validation (GCV), AIC information criterion, BIC information criterion, etc. We use a BIC-type criterion to select m [36], that is

where p is the number of components in parameter vector. c ˆ τ k ( m ) , θ ˆ ( m ) , and γ ˆ ( m ) are obtained by minimizing equation (7) with the first m principal components.

5 Simulation study

In this section, we conduct a Monte Carlo simulation example to assess the finite sample performance of the proposed procedures. Moreover, we include four competitors in our comparisons: (1) the LS method [14], (2) the QR method (QR) with τ = 0.5 [24], (3) the rank regression (RR) method [43], (4) the composite quantile regression (CQR) method [33]. Following [30], the performance of CQR estimation does not depend sensitively on the choice of the number of quantiles q . For the sake of simplicity, we only consider q = 9 for CQR and WCQR method in our simulation. The data sets are generated from the following PFLRM:

where both Z i 1 and Z i 2 are standard normal with correlation coefficient 0.5, θ 1 = 0.8 and θ 2 = 2 . For the functional linear component, we take the same form as [40]. Specifically, the slope function β ( t ) = 2 sin ( π t / 2 ) + 3 2 sin ( 3 π t / 2 ) and X ( t ) = ∑ j = 1 100 ξ j v j ( t ) , where the ξ j are independently distributed as the normal with mean 0 and variances given by λ j = 10 ( ( j − 0.5 ) π ) − 2 , v j ( t ) = 2 sin ( ( j − 0.5 ) π t ) .

In order to show the robustness of our estimators, the following six different error distributions are considered: standard normal distribution N(0, 1), t (3) distribution that is used to produce heavy-tailed distribution, F distribution F(4, 6), exponential distribution Exp(1), Log normal distribution LN(0, 1), and Weibull distribution Wb(1, 1.5) are used to produce asymmetric error distributions.

To assess the performance of the estimated parameters, the following criteria are considered:

(1) the average value of the corresponding coefficients (MEAN);

(2) the sample standard deviation (SD);

(3) the median of absolute bias (ABISE). To assess the performance of the slope function, we consider the following square root of average square errors (RASE):

where { t k , k = 1 , … , n grid } are the grid points at which the slope function β ( ⋅ ) is evaluated. In our simulation experiments, we set n grid = 200 , 500 repetitions are carried out with sample size n = 100 , 200, 400. In addition, the range of the selected m is from 1 to 8. The simulation results are presented in Tables 1, 2, 3, 4, 5, 6. From Tables 1–6, we have the following findings: (1) As the sample size n increases, the ABISE, SD, and RASE values of all estimators decrease. As expected, increasing the sample size leads to a better performance of all estimators. (2) The standard errors become smaller and the estimators of MEANs for the parametric components are very close to the true values θ as n increases, due to the fact of root- n consistency of the parameter estimators. This result confirms the asymptotic property that the parameter estimators are asymptotically unbiased as given in Theorem 3.1. (3) When the error follows N(0, 1), LS is the best one among the five estimators and WCQR performs nearly as well as CQR and RR. Although the QR method with quantile being 0.5 have similar performances, the WCQR seems to be slightly better than QR in most scenarios. (4) When the error distribution is symmetrical, CQR, RR, and WCQR perform similarly and satisfactorily. For all asymmetric errors: F(4, 6), Exp(1), LN(0, 1), and Wb(1, 1.5), WCQR performs better than others.

6 Application to Tecator data

We illustrate the proposed approach by analyzing a real estate data set which is available at http://lib.stat.cmu.edu/datasets/tecator. The data have been widely used to predict fat content on samples of finely chopped meat [44,45,46]. For each food sample, the functional data consist of a 100-channel spectrum of absorbances recorded on a Tecator Infratec Food and Feed Analysis working in the wavelength range 850–1,050 nm by the near infrared transmission principle. More details on the data can be found in Ferraty and Vieu [47]. In this section, we construct the following PFLRM model:

We denote the scale of fat content as Y i , the protein content as Z i , the moisture content as U i , and the absorbances as X i ( t ) .

To evaluate the performance of the method, the sample is randomly divided into two subsamples: the training sample, I 1 = { ( X i , Y i , Z i , U i ) , ∣ I 1 ∣ = 165 } , where ∣ I 1 ∣ denotes the cardinality of I 1 and the remaining is test sample, I 2 = { ( X i , Y i , Z i , U i ) , ∣ I 2 ∣ = 50 } . To verify our robust estimation procedures, we reanalyzed this data set by including some outliers in the response variable. The case is considered through controlling the number of outliers and the values of outliers, we randomly generate κ outliers by increasing the value Y to Y + t ( κ = 20 , t = 5 ) . The training samples are used to estimate the parameters and the test samples are employed to verify the quality of predictions. For this, we compute the mean square error of prediction MSEP [44], which is defined by

where Y ˆ i is the predicted value based on the training sample and Var I 2 is the variance of response variables from test sample. We use different estimation methods to predict the fat content of a meat sample based on its protein and/or moisture contents and/or its NIT absorbance spectrum, and the corresponding average MSEPs based on 200 times randomly splits are shown in Table 7. Compared with other methods, the MSEP of our proposed method (WCQR) is relatively smaller. Meanwhile, based on random simulations, the average estimated weight is ω opt = ( 0.4287 , 0.07228 , 0.0119 , 0.0114 , 0.0624 , 0.0631 , 0.0244 , 0.0022 , 0 ) and the corresponding standard deviation is (0.1608, 0.0895, 0.0481, 0.0257, 0.0275, 0.0403, 0.0690, 0.1261, 0). From Figure 1, we can see that the estimated slope functions β ( t ) ˆ by our WCQR method and CQR method, respectively, and it demonstrates that the two estimation approaches perform similarly.

Figure 1

Curves of slope function β ( t ) ˆ with CQR (solid line) and WCQR (dotted line).

7 Conclusion

In this paper, a novel and robust procedure based on WCQR and principal component basis function approximations is developed for PFLRMs. Theoretical properties of the estimators of both slope function and linear parameters are derived under some mild assumptions. We show that the estimators do not require any specification for error distribution, and thus more robust than those based on the LS, QR and CQR in the case of asymmetric errors. Furthermore, we present an efficient algorithm for the selection of optimal weights and discuss the selection of tuning parameters. Numerical studies and the real data analysis illustrate that our proposed method performs very well for moderate sample size.

In addition, we could extend the proposed procedure to the functional regression models with diverse dimension of covariates. These deserve to be studied further.