Robust estimation for varying coefficient partially functional linear regression models based on exponential squared loss function

1 Introduction

Varying coefficient models have become an important research topic in modern statistics and have been frequently applied for addressing problems in various scientific areas, such as economics, biomedical science, finance, medicine, and engineering. One superior characteristic that makes the varying coefficient models attractive is that they have the ability to explore the dynamic pattern by allowing coefficients to vary over the so-called index variable. Moreover, varying coefficient models retain the virtues of both parametric and nonparametric models; therefore, the possible modeling biases can be reduced and the curse of dimensionality can be avoided. Recently, much effort has been devoted to studying the methodological, theoretical, and applied sides around varying coefficient models, and we refer, for instance, to [1,2,3, 4,5,6].

Nowadays, technological innovations allow us to create and store large-scale data sets. In many cases, these data sets can be viewed as functions of time or spatial points or some other continua. Such type of data is termed functional data. In recent years, functional data analysis has attracted much attention among statisticians and practitioners. As the fundamental tool and mathematical framework for functional data analysis, the functional linear model (FLM), which characterizes the relationship between a scalar response and a functional predictor with a linear operator, is seeing a rise in popularity in practice. For a systematic coverage of statistical methods and inferences of FLM, see [7,8, 9,10,11, 12,13] and others. It is well known that the FLM model is useful. But in many practical situations, the simple relationship between the response variable and functional covariate may not be sufficient. To improve the power of prediction and interpretation of FLM, several extensions of FLM have been proposed for the analysis of mixed data with complex structures. Recently, Peng et al. [14] incorporated some real-valued covariates in FLM and proposed a mixed data model, namely, a varying coefficient partially functional linear regression model (VCPFLRM). The VCPFLRM takes the form

where Y is a real-valued random variable defined on a probability space ( Ω , ℬ , P ) , Z = ( Z 1 , … , Z p ) T and U are observed explanatory variables and so-called index variables defined on the same probability space, respectively. β ( U ) = ( β 1 ( U ) , … , β p ( U ) ) T consists of p unknown varying-coefficient functions, { X ( t ) : t ∈ T } is a zero mean, the second-order stochastic process valued in H = L 2 ( T ) , the set of all square-integrable functions on T with the inner product ⟨ x , y ⟩ = ∫ T x ( t ) y ( t ) d t , ∀ x , y ∈ H and norm ‖ x ‖ = ⟨ x , x ⟩ 1 ∕ 2 , α ( t ) ∈ H is unknown slope function, the error term ε is assumed to be independent of ( Z , U , X ) with E ( ε ∣ Z , U , X ) = 0 . Without loss of generality, we assume that U takes values from the unit interval [ 0 , 1 ] and suppose throughout that T = [ 0 , 1 ] . The VCPFLRM is flexible and includes some existing statistical models. For example, when α ( t ) = 0 , it reduces to the varying coefficient model; when β ( ⋅ ) are all set as constant, it reduces to the partially FLM [15]; when p ≡ 1 and Z 1 = 1 , it reduces to the semi-FLM [16].

To estimate the slope function α ( ⋅ ) and the coefficient functions β j ( ⋅ ) , j = 1 , … , p in VCPFLRM, Peng et al. [14] proposed the least squares estimation procedure based on B-spline basis function approximations. Feng and Xue [17] proposed a profile least squares estimation by using functional principal component (FPC) analysis and local linear smoothing technique. However, all the above existing articles are built on least square method and the assumptions that the error is normally distributed, which are sensitive to outliers and heavy-tailed error distributions. Hence, robustness against outliers is a very important issue in the VCPFLRM. Recently, Wang et al. [18] investigated a new estimation method for linear regression models by minimizing exponential squared loss (ESL) ρ ( ⋅ ) , i.e., ρ ( t ) = 1 − exp ( − t 2 ∕ h ) , where the tuning parameter h ∈ ( 0 , ∞ ) controls the degree of robustness and efficiency for the estimators. For large h , 1 − exp ( − t 2 ∕ h ) ≈ t 2 ∕ h , and therefore the proposed estimators are similar to the least squares estimators in the extreme case. For a small h , the outliers in the observations t will receive a small impact on the estimators. Therefore, this method is more robust to outliers than the least squares method. Due to its nice property, this method has been quickly extended to the statistical models for the sake of robust inference, the recent literature include but are not restricted to Jiang et al. [19] for the varying coefficient partially nonlinear model, Song et al. [20] for the single index varying coefficient model, Yu et al. [21] for the semi-functional linear regression model, Lv et al. [22] for the generalized varying coefficient partially linear model. Even though exponential-squared-based estimation has been well developed, the theory and methodology for VCPFLRM are none. It is therefore our impetus for applying this method to VCPFLRM.

The main objective of this article is to propose a new robust estimation procedure based on the ESL function for the VCPFLRM. To estimate α ( ⋅ ) in VCPFLRM, two commonly used approaches have been proposed, including the FPC basis method and the basis function expansions (such as B-spline functions approximation) method. Since FPC has some desirable properties in approximating a smooth function, for example, it is a data-driven technique where the data can be represented as linear combinations of functions estimated from the data. Besides, by requiring the order of generalized Fourier coefficients based on the FPC basis, we can control the smoothness of the slope function. Due to these advantages mentioned above, we first use B splines and FPC basis to approximate the varying coefficient functions and slope function in model (1), respectively. With the help of these approximations, the estimation problem based on the ESL function for VCPFLRM becomes a problem of estimating the coefficients in the linear combinations. Under some mild regularity conditions, we show that the estimators of coefficient functions and slope functions can achieve the optimal convergence rate. Moreover, our numerical studies and a real data application in the sequel demonstrate that our proposed method is at least comparable to some existing works. It is worth pointing out that there exist some essential differences between our proposed method and that of Wang et al. [18] and Yu et al. [21]. First, our studied VCPFLRM is flexible enough to cover the classical linear regression model and the semi-functional linear regression model as special cases. Accordingly, a different approach is required to deal with varying coefficient functions β ( ⋅ ) , this difference leads to different theoretical properties for proposed estimators. Second, compared with the works of [18] and [21], though the estimate of β ( ⋅ ) and α ( ⋅ ) can be easily obtained by transforming the original VCPFLRM to a standard linear regression model, we need to choose extra tuning parameters by the data-driven procedure for spline estimates of β ( ⋅ ) , such as the number of interior knots, which controls the dimensions of the spline spaces used to approximate the true varying coefficient functions.

The rest of the article is organized as follows. In Section 2, we describe the proposed estimation procedure for VCPFLRM. In Section 3, we establish the theoretical properties of the resulted estimators. In Section 4, we develop an estimation algorithm and discuss the selection of tuning parameters. Section 5 illustrates the numerical performance of the proposed method through simulation experiments. A real data analysis is presented in Section 6, and some concluding remarks are followed in Section 7. All proofs of theorems are provided in the Appendix.

2 Estimation procedure

Let { Y i , Z i , U i , X i } i = 1 n be an independent and identically distributed sample that is generated from model (1). Motivated by [18], we estimate the slope function α ( t ) and nonparametric coefficients β ( U ) in VCPFLRMs by maximizing the following ESL function

In what follows, we first approximate the nonparametric components β ( U ) by B splines. Let B ( u ) = ( B 1 ( u ) , … , B K ( u ) ) T be the B-spline basis functions, where K = k n + q + 1 is the number of basis functions, q is the degree of the spline, and k n is the number of internal knots that control the smoothness of nonparametric functions. Thus, for any function β j ( u ) , j = 1 , … , p , we have

Note that it is possible to specify different K for each component, but we assume that they are the same for simplicity. Next, we construct the covariance function and the empirical covariance function for the random process X ( ⋅ ) as follows:

Obviously, the covariance function K ( s , t ) defines a linear operator that maps a function f to K f given by ( K f ) ( s ) = ∫ K ( t , s ) f ( t ) d t . We assume that the linear operator with kernel K ( s , t ) is positive definite. Let λ 1 > λ 2 > ⋯ > 0 and λ ˆ 1 ≥ λ ˆ 2 ≥ ⋯ ≥ λ ˆ n + 1 = ⋯ = 0 be the ordered eigenvalue sequences of the linear operators with kernels K ( s , t ) and K ˆ ( s , t ) . { ϕ j } and { ϕ ˆ j } are an orthonormal sequence of eigenfunctions, and it is clear that the vectors ϕ 1 , ϕ 2 , … form an orthonormal basis in H . Then, consider the spectral decomposition of the covariance function K ( s , t ) and K ˆ ( s , t ) as follows:

According to the Karhunen-Loèvere presentation, two L 2 -valued functions X ( t ) and α ( t ) have the following expansions:

where ξ i = ∫ 0 1 X ( t ) ϕ i ( t ) d t , γ j = ∫ 0 1 α ( t ) ϕ j ( t ) d t . It follows that ξ i s are uncorrelated random variables with mean 0 and variance E ( ξ i 2 ) = λ i . For more details, see [23].

Let Π i = ( Z i 1 B 1 ( U i ) , … , Z i 1 B K ( U i ) , … , Z i p B K ( U i ) ) T , η = ( η 1 T , … , η p T ) T with η j = ( η j 1 , … , η j K ) T . Based on equations (3)–(6), model (1) can be rewritten as follows:

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 (7) can be rewritten as follows:

where V i = { ⟨ X i , ϕ ˆ j ⟩ } j = 1 , … , m , γ = ( γ 1 , … , γ m ) T . Thus, the objective function (2) becomes

Denote the maximizer of (9) by η ˆ and γ ˆ . Consequently, the estimated component functions are given by β ˆ j ( u ) = B T ( u ) η ˆ j , α ˆ ( t ) = ∑ j = 1 m γ ˆ j ϕ ˆ j ( t ) .

3 Theoretical results

In this section, we aim to establish the theoretical properties of the resulting estimators. First, we need to make some definitions and assumptions. We assume that the true functions in model (1) are α 0 ( ⋅ ) and β 0 j ( ⋅ ) , j = 1 , … , p . We use C to denote a generic positive constant that can vary from line to line. a n ∼ b n means that a n ∕ b n is bounded away from zero and infinity as n → ∞ . The notation ‖ ⋅ ‖ is the L 2 norm for a function or the Euclidean norm for a vector. Let S = ⟨ X ( t ) , α 0 ( t ) ⟩ , ψ h ( ς ) = exp ( − ς 2 ∕ h ) , F ( z , u , s , h ) = E ( ψ h ″ ( ε ) ∣ Z = z , U = u , S = s ) , and G ( z , u , s , h ) = E ( ψ h ′ ( ε ) 2 ∣ Z = z , U = u , S = s ) . The following regularity conditions will be required, which are common and some can be found in [14], [24].

1. E ‖ X ( ⋅ ) ‖ 4 < C < ∞ , E [ ξ i 4 ] ≤ C λ i 2 , for i ≥ 1 , for ι > 2 , E ( ∣ Z j ∣ ι ) < ∞ for j = 1 , … , p .

2. The eigenvalues λ i and score coefficients γ j satisfy the following conditions, respectively.

   1. There exist some constant C and a > 1 such that

      C − 1 i − a ≤ λ i ≤ C i − a , λ i − λ i + 1 ≥ C i − a − 1 , i ≥ 1 .

   2. There exist some constant C and b > a ∕ 2 + 1 such that

      ∣ γ j ∣ ≤ C j − b , j ≥ 1 .

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

4. The density function of U , f U ( u ) is bounded away from 0 and infinity on [ a , b ] . Furthermore, we assume that f U ( u ) is continuously differentiable on ( a , b ) .

5. β j ( t ) ∈ ℱ r , j = 1 , … , p for some r > 2 , where ℱ r is the collection of all functions on [ 0 , 1 ] whose d th order derivative satisfies the Hölder condition of the order ς with r ≡ d + ς and 0 < ς ≤ 1 .

6. F ( z , u , s , h ) and G ( z , u , s , h ) are continuous with respect to ( z , u , s ) and F ( z , u , s , h ) < 0 for any h > 0 .

7. E ( ψ h ′ ( ε ) ∣ z , u , s ) = 0 , E ( ψ h ″ 2 ( ε ) ∣ z , u , s ) , E ( ψ h ′ 3 ( ε ) ∣ z , u , s ) , E ( ψ h ‴ ( ε ) ∣ z , u , s ) are continuous with respect to ( z , u , s ) .

4 Computational issue

5 Simulation experiments

In this section, we conduct simulations to compare the finite sample performance of our proposed estimation method (ESL) with the other two methods: the quantile regression method [24] with quantile τ = 0.5 (QR) and the least square (LS) method [14]. We generated data from the model as follows:

where the coefficient functions are β 1 ( u ) = 2 sin ( 2 π u ) , β 2 ( u ) = 8 u ( 1 − u ) , and β 3 ( u ) = 5.5 + 10 exp ( 2 u − 1 ) , the covariate Z = ( Z 1 , Z 2 , Z 3 ) T , with each Z i , i = 1 , 2 , 3 , being standard normal distribution and the correlation between Z i and Z i being 0 . 5 ∣ i − j ∣ , U ∼ U ( 0 , 1 ) . For the functional linear component, we take the same form as [21], that is, α ( t ) = 2 sin ( π t ∕ 2 ) + 3 2 sin ( 3 π t ∕ 2 ) and X ( t ) = ∑ j = 1 50 ξ j v j ( t ) , where the ξ j are independently distributed as the normal with mean 0 and variance λ j = 10 ( ( j − 0.5 ) π ) − 2 , v j ( t ) = 2 sin ( ( j − 0.5 ) π t ) . In order to examine the robustness and efficiency of our proposed method, the following three different error distributions are considered:

Case 1. ε follows the standard normal distribution, abbreviated as N ( 0 , 1 ) .

Case 2. ε follows the standard t distribution with three degrees of freedom, abbreviated as t ( 3 ) , which is used to produce heavy-tailed distribution.

Case 3. ε follows the contaminated normal distribution with k % data from N ( 0 , ℓ ) and remaining 1 − k % data from N ( 0 , 1 ) , abbreviated as 0.9 N ( 0 , 1 ) + 0.1 N ( 0 , 10 ) . Obviously, the k % data can be regarded as outliers, and we set k = 10 and ℓ = 10 in this case.

To evaluate the performance of our estimation for slope function α ( ⋅ ) and coefficient functions β j ( ⋅ ) , we define the following square root of the average square errors (RASE):

where { t i , i = 1 , … , n 1 , grid } and { u i , i = 1 , … , n 2 , grid } are the grid points at which the functions α ˆ ( ⋅ ) and β ˆ j ( ⋅ ) are evaluated, respectively. In our simulation experiments, we simply set n 1 , grid = 201 and n 2 , grid equals to the sample size n in each scenario. A total of 200 replications with n = 100 , 200 , 600 are performed. The mean, median, and standard error (SE) of the RASE values for three error distribution cases are reported in Tables 1, 2, 3. Several observations can be found in Tables 1–3. First, when the error follows a standard normal distribution, the LS estimate slightly outperforms the QR and ESL estimates, and ESL performs slightly better than QR. Second, ESL is significantly superior to LS in the other two cases, and this is expected because the performance of LS is sensitive to heavy-tailed distributions or outliers, ESL seems to be slightly better than QR in most cases. Third, for the given error distribution, all the considered estimators become better and better as the sample size n increase gradually. Figure 1 demonstrates the performances of average curve estimates of β j ( ⋅ ) and α ( ⋅ ) when t ∼ 0.9 N ( 0 , 1 ) + 0.1 N ( 0 , 10 ) and n = 100 . The plots for other cases are similar and are omitted. All these results indicate that our proposed method works well under the considered settings. To further demonstrate the performance of three different methods, we also consider the Cauchy error distribution. The results are given in Figure 2. We find LS has poor performance since the LS estimates will break down and suffer from consistency when the error is Cauchy.

Figure 1

The estimated curves of functions β j ( u ) , j = 1 , 2 , 3 and α ( t ) when ε ∼ 0.9 N ( 0 , 1 ) + 0.1 N ( 0 , 10 ) and n = 100 .

Figure 2

The estimated curves of functions β j ( u ) , j = 1 , 2 , 3 and α ( t ) when ε ∼ Cauchy and n = 100 .

To further compare the efficiency of the LS, QR, and ESL estimation of the nonparametric functions g ˆ = ( α ˆ ( ⋅ ) , β ˆ j ( ⋅ ) ) , we consider the ratio of average squared error (RT) defined as

where g ˆ LS denotes the LS estimate of g , and g ˆ could be taken as the QR or ESL estimate. The corresponding results are summarized in Table 4.

From Table 4, we can see that the RTs of the ESL method are slightly smaller than 1 compared with the RTs of the QR method when the error is normally distributed, for the other two error distributions, the RTs of the ESL method are greater than 1, which indicates that our proposed estimation method can significantly improve the estimation efficiency.

6 A real data analysis

For illustration purposes, we apply the proposed methodology to analyze the spectral data. The data set contains 215 chopped meat samples, each data sample contains fat, protein, moisture contents, and spectral curve. The three contents are measured in percent and determined by analytic chemistry. The spectral curve consists of 100 wavelengths of absorbance spectrum records. Our main interest is to explore the relationship between fat content and other factors. In the subsequent analysis, we treat the fat content as the response variable Y , the protein content as Z , the moisture content as U , and the spectral curve as X ( t ) . Beforehand, we standardize the response variable Y and two other contents Z and U . As done in [17], we use the following VCPFLRM to fit the data:

For a fair evaluation, the randomly selected 129 samples are used as training data to build the model, and the remaining 86 observations are left as the validation data for prediction purposes. The prediction performance is measured by the mean squared error proposed by [28], which is defined as MSE = 1 86 ∑ i ∈ G ( Y i − Y ˆ i ) 2 , where G is the validation data set. To verify the proposed robust estimation procedure, we re-analyzed this data set by including some outliers in the training response variable; in other words, the response Y of training data with ι high leverage outliers being replaced by Y + κ (Case 1: κ = 5 , ι = 10 , Case 2: κ = 10 , ι = 10 and Case 3: κ = 5 , ι = 20 ). Then, we apply the three different estimation procedures (LS, QR with quantile τ = 0.5 , and ESL) to analyze the data set by VCPFLRM stated as (10). The mean values of MSE of the proposed ESL procedure with those of QR and LS methods based on 200-time repetitions are shown to be 0.022, 0.023, and 0.706 for Case 1, 0.019, 0.032, and 2.378 for Case 2, 0.035, 0.042, and 1.323 for Case 3, respectively. Consequently, our proposed ESL method is better than the others in terms of performance of prediction when data set contains outliers. For comparison, we also use the functional linear regression model to fit the data. We only consider case 2, that is, the response Y of training data has 10 high leverage outliers, with 200-time repetitions, the mean values of MSE based on ESL, QR, and LS are shown to be 0.067, 0.102, and 3.187. This result reveals that our proposed model and estimation method can improve the accuracy of the prediction when data sets are subjected to high outliers in the response.

7 Concluding remarks

In this article, we provide a robust and efficient estimation approach for VCPFLRMs based on the ESL function. The varying coefficients are approximated by B-splines, and the slope function is estimated by the FPC basis. The theoretical properties of proposed estimators are established under some conditions. Furthermore, we develop a computation algorithm to solve the optimization problems and then introduce a data-driven procedure to select the tuning parameters. The outstanding feature of the newly proposed method is that it can achieve good performance for a wide class of errors by selecting an appropriate tuning parameter h based on observed real data. The good properties of this method have also been demonstrated through both numeric examples and real data analysis. As discussed in [28], the FPC-based methods select data-driven basis functions relying on the decomposition of covariance function of the functional predictor X ( t ) . To compensate for this defect, we could use B-spline basis functions to approximate the slope function for its efficiency in function approximation and numerical computation. In addition, the method proposed here may be extended to other functional regression models, which could be considered in future research.