Robustness of Variable Selection Based on Functional Response Regression Model

Ziqi Yan, Wanzhou Ye

Department of Mathematics, College of Science, Shanghai University, Shanghai, China.

DOI: 10.4236/apm.2025.153010                

Abstract

Function-on-scalar regression is a type of function response regression used to analyze the relationship between function response and a set of scalar predictor factors. The variable selection methods of FOSR models mostly focus on the linear effects of scalar predictor factors. Therefore, in this paper, we perform robust variable selection for nonlinear FOSR models with the presence of multiple continuous covariates to further explain the behavior of function response over time. We project functional data into a low dimensional principal component space via the principal component score of the function response, in order to use the principal component score for variable selection and regression modeling. In this paper, we develop a regularized iterative algorithm based on exponential squared loss and group smoothly clipped absolute for predicting estimates of scalar factors and function coefficients, and use tuning parameters chosen by data-driven methods to achieve robustness in variable selection. The robustness of the proposed method is verified through simulation studies and demonstrated on real datasets.

Keywords

Functional Data Analysis, Functional Response Regression, B-Splines, Robust Estimation, Variable Selection

1. Introduction

With the advancement of modern technology, more and more data are recorded continuously at intervals of time, or intermittently at several discrete points in time, often stored in the form of curves and images. Functional data are essentially infinite dimensional its basic unit is the function. A common type of functional data are curves recorded during a time interval. Functional Data Analysis (FDA) includes statistical methods for such data. FDA deals with the analysis of data which is in the form of functions. The functions could be curves, images, or other types of forms, depending on the set on which the functions are recorded. Functional data analysis is mainly divided into the following three regression models: 1) scalar-on-function regression: a continuous response variable regressed on functional covariates; 2) function-on-function regression: a functional response regressed on functional predictors; 3) Function-on-scalar regression: a functional response regressed on scalar predictors. The aim of this paper is to provide robust variable selection for function-on-scalar regression models containing nonlinear scalar predictors.

Function-on-scalar regression (FOSR) is an increasingly popular area of research in the analysis of functional data, e.g. Barber et al. [1] for genome-wide association studies (GWAS), Chen et al. [2] for the effect of stroke severity on motor control, Goldsmith et al. [3] for physical activity (PA), and so on. In FOSR, a continuous functional response is typically modeled by the additive effects of a set of scalar predictors of interest, which are captured using smooth univariate coefficient functions. This has led to several methods of variable selection in the FOSR model. Chen et al. [2] used a pre-whitening approach to account for temporal correlation within curves and proposed a variable selection procedure via group minima and maxima concave penalties (MCP) regularization with a penalized least squares technique. Barber et al. [1] used a penalized least squares objective function with the LASSO penalty on the basis coefficients and proposed a function-on-scalar LASSO (FSL). Fan and Reimherr [4] further developed an adaptive FSL in high dimensions. Parodi and Reimherr [5] proposed a functional linear adaptive mixed estimation method to achieve variable selection and smoothing for function-on-scalar regression.

The vast majority of FOSR methods focus on linear effects of scalar predictors, but this may not necessarily hold true in many applications. Ghosal and Maity [6] proposed a nonlinear function-on-scalar regression approach considering the additional nonparametric effects of the scalar predictors into the FOSR model. This paper is based on this nonlinear model for robust variable selection. It is well known that the least squares method is sensitive to outliers in the data, however, the variable selection methods discussed above are based on the least squares criterion. Therefore, it is necessary to propose a more robust variable selection method to replace the least squares method in the presence of outliers. Many loss functions are resistant to outliers, such as minimum absolute deviation loss, quantile loss, or Huber loss, which can be used for robust estimation and variable selection. In this paper, we introduce the exponential squared loss (ESL) introduced by Wang et al. [7] in nonlinear FOSR to reduce the effect of outliers. ESL is a bounded continuous function, and we can control the robustness of the estimator by selecting the appropriate adjustment parameter $h$ , that is, reducing the influence of outliers that lead to large errors. The ESL function is defined as ${\varphi }_h\left(t\right)=1-\exp \left(-t^2/h\right)$ . When $h$ is very large, ${\varphi }_h\left(x\right)\approx x^2/h$ , and therefore the proposed estimator is similar to the least squares estimator in the extreme case. When $h$ is small, observations with larger absolute values will have a small impact on the loss function ${\varphi }_h\left(x\right)$ and therefore have a smaller impact on the estimator.

In this paper, we investigate the use of functional principal component analysis (FPCA) to project the functional response into a low-dimensional principal component space in a nonlinear FOSR model. The principal component scores are used as target variables in a regression model to analyze the effect of scalar covariates on the functional response and perform robust variable selection for the nonlinear covariates. The nonlinear covariates are approximated by a B-spline basis function, and the estimation of the coefficient function combines the FPCA eigenfunctions and the B-spline basis function. We integrate the exponential squared loss function into the nonlinear FOSR model and propose a robust variable selection process using the group SCAD regularization method. By locally linearly approximating the SCAD penalty function, it is replaced by a set of weighted LASSO penalty functions. We utilize a data-driven approach to select the smoothing parameters, which achieves robustness and establishes consistency of the proposed estimators. The estimated coefficients are projected back using univariate basis functions and eigenfunctions of the functional response to recover the effects of the scalar predictors. Finally, we show through simulations and real data applications that the proposed method performs well. When outliers are present in the dataset, our method is robust to outliers. Meanwhile, in the absence of outliers, the proposed estimation method is comparable to the least squares estimation method.

The rest of this paper is organized as follows. In Section 2, we describe the estimation methods for nonlinear function-on-scalar regression and discuss the theoretical properties of the proposed estimators. In Section 3, we present the implementation algorithm for computing the estimator based on nonlinear FOSR and the criterion used to select the smoothing parameter. To demonstrate the superiority of the proposed method, in Section 4, we report and compare the simulation results and illustrate the proposed method with a real diffusion tensor imaging data example. Finally, in Section 5, the summary and conclusion of the paper are given in Section 5.

2. Estimation Method and Theoretical Properties

2.1. Estimation Method

Suppose that the observed data for the $i\text{th}$ subject is $\left\{Y_i\left(t\right),M_{i1},M_{i2},\cdots ,M_{iq}\right\}$ $\left(i=1,2,\cdots ,n\right)$ where $Y_i\left(t\right)$ represents the functional response and $M_{i1},M_{i2},\cdots ,M_{iq}$ are the corresponding scalar predictors of interest. Sometimes there might be additional control or confounding variables $X_{i1},X_{i2},\cdots ,X_{ip}$ which we want to adjust for. In this paper, we assume the functions are observed on a dense and regular grid of points $S=\left\{t_1,t_2,\cdots ,t_m\right\}\subset T=\left[a,b\right]$ for some $a,b\in R$ where the functional response is observed on an irregular and sparse domain. Now, the linear FOSR model is extended to the nonlinear dynamic effect of scalar predictors, and we use the following generalized function on scalar regression (GFOSR) model proposed by Ghosal and Maity [6]:

$Y_i\left(t\right)=\mu \left(t\right)+{\displaystyle \sum _{j=1}^p{X}_{ij}{\beta }_j\left(t\right)}+{\displaystyle \sum _{l=1}^q{\theta }_l\left(M_{il},t\right)}+{\epsilon }_i\left(t\right),$ (1)

Here, ${\beta }_j\left(t\right)$ represents the functional effect of the scalar predictor $X_{ij}$ , and ${\theta }_l(,)$ is an unknown smooth function on $R\times T$ (twice differentiable with respect to both parameters), used to capture the functional effect of the predictor $M_{il}$ . The coefficient functions ${\beta }_j(\cdot )$ , ${\theta }_l\left(\cdot ,\cdot \right)$ , functional response $Y_i(\cdot )$ , and the error process ${\epsilon }_i(\cdot )$ are assumed to lie in a real separable Hilbert space [1]. In this paper, we focus our attention to $L^2\left(T\right)$ . We assume that the covariates are centered and $\mu (\cdot )$ is a smooth function capturing marginal mean of the functional response $Y_i(\cdot )$ , where the error functions ${\epsilon }_i(\cdot )$ are i.i.d. of ${\epsilon }_i(\cdot )$ which is a mean zero stochastic process with unknown nontrivial covariance structure.

Let $Y_1\left(t\right),\cdots ,Y_n\left(t\right)$ be independent realizations of a smooth random function $Y\left(t\right)$ with mean function $\mu \left(t\right)$ and covariance function $\mathrm{cov}\left\{Y\left(t\right),Y\left(s\right)\right\}=$ $C\left(t,s\right)$ . By the Karhunen-Loève theorem by Ash and Gardner [8], the random

function $Y_i\left(t\right)=\mu \left(t\right)+{\displaystyle \sum _{k=1}^{\infty }{\xi }_{ik}{\varphi }_k\left(t\right)}$ where $\mu \left(t\right)=E\left\{Y_i\left(t\right)\right\}$ is the overall mean

function; ${\varphi }_k\left(t\right)$ is the $k\text{th}$ orthonormal eigenfunction of the covariance function, satisfying ${\displaystyle \int {\varphi }_k\left(t\right){\varphi }_j\left(t\right)\text{d}t}=1$ if $j=k$ and 0 otherwise; ${\xi }_{ik}$ denotes FPC scores with $E\left({\xi }_{ik}\right)=0$ , $\mathrm{var}\left({\xi }_{ik}\right)={\rho }_k$ and $\mathrm{cov}\left({\xi }_{ij},{\xi }_{ik}\right)=0$ if $j\ne k$ , and ${\rho }_k$ is the eigenvalue corresponding to the eigenfunction ${\varphi }_k\left(t\right)$ . Since

$\underset{t\in \left[0,1\right]}{\sup }E{\left\{{\displaystyle \sum _{k=1}^{\infty }{\xi }_{ik}{\varphi }_k\left(t\right)}-{\displaystyle \sum _{k=1}^{K_n}{\xi }_{ik}{\varphi }_k\left(t\right)}\right\}}^2\to 0$ , one can reasonably suppose that

$Y_i\left(t\right)\approx \mu \left(t\right)+{\displaystyle \sum _{k=1}^{K_n}{\xi }_{ik}{\varphi }_k\left(t\right)}$ (2)

as $K_n\to \infty$ . The approximation of (2) with a fixed $K_n$ is commonly adopted for longitudinal and functional data analysis (see Yao et al. [9]). To be more flexible, Hall and Mohammad [10] considered model (2) with $K_n\to \infty$ as $n\to \infty$ . In this paper, we adopt the FPCA method of Yao et al. [9], which is used to estimate the FPCA score ${\xi }_{ik}$ , the mean function $\mu \left(t\right)$ , and the eigenfunction ${\varphi }_k\left(t\right)$ . Specifically, a spectral decomposition of the covariance function yields $\mu \left(t\right)$ and ${\varphi }_k\left(t\right)$ , and thus ${\xi }_{ik}$ is estimated by numerical integration. Thus, the model in (1) becomes:

$\begin{array}{c}{\xi }_{ik}={\displaystyle \underset{T}{\int }\left(Y_i\left(t\right)-\mu \left(t\right)\right){\varphi }_k\left(t\right)\text{d}t}\\ ={\displaystyle \sum _{j=1}^p{X}_{ij}{\displaystyle \underset{T}{\int }{\beta }_j\left(t\right){\varphi }_k\left(t\right)\text{d}t}}+{\displaystyle \sum _{l=1}^q{\displaystyle \underset{T}{\int }{\theta }_l\left(M_{il},t\right){\varphi }_k\left(t\right)\text{d}t}}+{\displaystyle \underset{T}{\int }{\epsilon }_i\left(t\right){\varphi }_k\left(t\right)\text{d}t}\\ ={\displaystyle \sum _{j=1}^p{X}_{ij}{\beta }_{jk}}+{\displaystyle \sum _{l=1}^q{\theta }_{lk}\left(M_{il}\right)}+{\epsilon }_{ik}\end{array}$ (3)

for $k=1,\cdots ,K$ , where we define ${\beta }_{jk}={\displaystyle \underset{s}{\int }{\beta }_j\left(t\right){\varphi }_k\left(t\right)\text{d}t}$ and

${\theta }_{lk}\left(u\right)={\displaystyle \underset{T}{\int }{\theta }_l\left(u,t\right){\varphi }_k\left(t\right)\text{d}t}$ .

We project the nonlinear variable $M_{il}$ onto the B-spline basis function space, and expand the nonparametric function ${\theta }_l(\cdot )$ with B-spline basis functions as

${\theta }_{lk}\left(u\right)={\displaystyle \underset{T}{\int }{\theta }_l\left(u,t\right){\varphi }_k\left(t\right)\text{d}t}$ . Let $B_j\left(t\right)={\left(B_j\left(t\right):1\le j\le J_n\right)}^{\text{T}}$ denote the qth-order

B-spline basis functions, where $J_n=N_n+q$ , and $N_n$ is the number of interior knots for a knot sequence ${\alpha }_1=\cdots =0={\alpha }_q<{\alpha }_{q+1}<\text{}\cdots <{\alpha }_{N_n+q}<1={\alpha }_{N_n+q+1}=\cdots ={\alpha }_{N_n+2q}$ . So, model (3) becomes as follows:

${\xi }_{ik}=X_i^{\text{T}}{\beta }_K+{\displaystyle \sum _{l=1}^q{\displaystyle \sum _{j=1}^{J_n}{B}_j\left(M_{il}\right){\delta }_{l,k,j}+{\epsilon }_{ik}}}$ (4)

where we define $X_i^{\text{T}}=\left(X_{i1},\cdots ,X_{ip}\right)$ , ${\beta }_k^{\text{T}}=\left({\beta }_{1k},\cdots ,{\beta }_{pk}\right)$ , ${\delta }_{\cdot k\cdot }^{\text{T}}=\left({\delta }_{1,k,1},\cdots ,{\delta }_{1,k,J_n},\cdots ,{\delta }_{q,k,J_n}\right)$ , ${\beta }_k^{\text{T}}=\left({\beta }_{1k},\cdots ,{\beta }_{pk}\right)$ , and

${\eta }_i^{\text{T}}=\left(B_1\left(M_{i1}\right),\cdots ,B_1\left(M_{iq}\right),B_2\left(M_{i1}\right),\cdots ,B_{J_n}\left(M_{iq}\right)\right)$ . Let $Z_i\left[U_i,V_i\right]$ dimension

of $K\times \left(Kp+Kq{J}_n\right)$ , $\gamma =\left(\beta ,\delta \right)$ , where $U_i=I_{K\times K}\otimes X_i^{\text{T}}$ dimension of $K\times \left(Kp\right)$ , $V_i=I_{K\times K}\otimes {\eta }_i^{\text{T}}$ dimension of $K\times \left(Kq{J}_n\right)$ , $\delta ={\left({\delta }_{\cdot 1}^{\text{T}},\cdots ,{\delta }_{\cdot k}^{\text{T}}\right)}^{\text{T}}$ , ${\xi }_1=\left({\xi }_{i1},{\xi }_{i2},\cdots ,{\xi }_{ik}\right)$ , and ${\epsilon }_1=\left({\epsilon }_{i1},{\epsilon }_{i2},\cdots ,{\epsilon }_{ik}\right)$ . Thus, we can write the FOSR model as,

${\xi }_i=Z_i\gamma +{\epsilon }_i.$ (5)

Then, the robust estimator of $\gamma =\left(\beta ,\delta \right)$ can be obtained by minimizing

$Q\left(\gamma \right)={\displaystyle \sum _{i=1}^n{\varphi }_h}\left({\xi }_i-Z_i\gamma \right),$ (6)

where ${\varphi }_h(\cdot )$ is the ESL function. We select variables for all scalar predictors, including linear and nonlinear variables, and group them. We group the variables $\gamma$ to be selected into, $\gamma =\left({\gamma }_1,\cdots ,{\gamma }_{l+1}\right)$ where each ${\gamma }_i$ represents a group of variables. Specifically, when $l=0$ , ${\gamma }_1=\left({\beta }_1,\cdots ,{\beta }_p\right)$ , the effect of all control or confounding variables $X_j,j=1,\cdots ,p$ , is divided into a group that is either simultaneously selected or simultaneously 0. When $l=1,\cdots ,q$ , ${\gamma }_{l+1}=\left({\delta }_{l,1,1},\cdots ,{\delta }_{l,k,J_n}\right)$ , $l=0,1,\cdots ,q$ each nonlinear predictor $M_l$ is divided into a group, with the projection coefficients under the B-spline basis forming a group for variable selection. We use the group penalty function $p_{\lambda }\left({\Vert {\gamma }_l\Vert }_2\right)$ for group variable selection, where

${\Vert {\gamma }_l\Vert }_2=\sqrt{{\displaystyle \sum _{j=1}^d{\gamma }_{l,j}^2}}$ , where d is the number of variables within theth group. In this

paper, we use the group smoothly clipped absolute deviation (SCAD) regularization method, thus proposing to minimizing:

$L\left(\gamma \right)={\displaystyle \sum _{i=1}^n{\varphi }_h\left({\xi }_i-Z_i\gamma \right)+n{\displaystyle \sum _{l=0}^q{p}_{\lambda }\left(\Vert {\gamma }_l\Vert \right)}}$ (7)

The SCAD penalty of Fan and Li [11] is defined in the following way,

$p_{\lambda }\left(\theta \right)=\{\begin{array}{l}\lambda \left|\theta \right|,\left|\theta \right|\le \lambda ,\\ -\frac{{\theta }^2-2a\lambda \left|\theta \right|+{\lambda }^2}{2\left(a-1\right)},\lambda <\left|\theta \right|\le a\lambda ,\\ \frac{\left(a+1\right){\lambda }^2}{2},\left|\theta \right|>a\lambda .\end{array}$

for some $a>2$ . Where $\Vert {\gamma }_j\Vert$ is the Euclidean norm of ${\gamma }_j$ , and $p_{\lambda }\left(\Vert {\gamma }_j\Vert \right)$ is the penalty function with $\lambda$ as a regularization parameter. When the loss function ${\varphi }_h\left(x\right)=x^2$ is the least square loss, it is the least squares (LS) method proposed by Wang, et al. [12]. When the loss function ${\varphi }_h\left(x\right)=\left|x\right|$ is the least absolute deviation loss, there is a LAD estimation method.

2.2. Theoretical Properties

In this subsection, we discuss the theoretical properties of the proposed estimators. We apply local linear smoothing to the mean function $\mu$ . Let $\widehat{C}\left(s,t\right)$ be the estimated smooth surface of $C\left(t,s\right)=\mathrm{cov}\left\{Y\left(t\right),Y\left(s\right)\right\}$ obtained through local quadratic fitting, and $h_{\mu }$ and $h_G$ be the bandwidths of the estimated values $\widehat{\mu }$ and $\widehat{G}$ , respectively. We use the symbol $\Vert \cdot \Vert$ to denote the Euclidean norm for vectors and the $L^2\left(T\right)$ norm for functions defined on $T$ . Let $f^{\prime }$ and $f^{″}$ present the first and second derivatives of $f$ , respectively. We impose the following regularity conditions:

(C1) The mean function $\mu \left(t\right)$ and the covariance function $G\left(s,t\right)$ are

smooth, and the estimation error satisfies: $\underset{t\in T}{\sup }\left|\widehat{\mu }-\mu \left(t\right)\right|=O_p\left(1/\sqrt{n{h}_{\mu }}\right)$ , $\underset{t\in T}{\sup }\left|\widehat{G}\left(s,t\right)-G\left(s,t\right)\right|=O_p\left(1/\sqrt{n{h}_G^2}\right)$ , thus FPCA estimation characteristic function

${\widehat{\varphi }}_k\left(t\right)$ converge to the real characteristic function:

$\sup \left|{\widehat{\varphi }}_k\left(t\right)-{\varphi }_k\left(t\right)\right|=O_p\left(1/\sqrt{n{h}_{\mu }}\right)$ .

(C2) Error function $\epsilon \left(t\right)$ is a finite variance with a mean of zero, satisfy the$\mathrm{cov}\left({\epsilon }_i\left(t\right),{\epsilon }_j\left(t\right)\right)=0$ , $\forall j\ne i$ , so that the FPCA score ${\xi }_{ik}$ satisfies consistency:${\widehat{\xi }}_{ik}\stackrel{p}{\to }{\xi }_{ik}$ .

(C3) The matrix $E\left(Z{Z}^{\text{T}}\right)$ is positively definite and the eigenvalues are bounded, where the dimension of $Z$ is fixed. There exists a positive constant $M$ such that $\left|Z_j\right|\le M$ , for all $1\le j\le p$ .

(C4) The functional coefficient ${\delta }_l(\cdot )\in C^r\left[0,1\right]$ for some integer $r\ge 2$ , and the spline order $q$ satisfies $q\ge r$ . $C^r\left[0,1\right]=\left\{\phi |{\phi }^{\left(r\right)}\in C\left[0,1\right]\right\}$ means the space of $r$ -th order smooth function.

(C5) The distance between neighboring knots $H_i={\tau }_{i+1}-{\tau }_i$ , and $H={\max }_{1\le i\le 1+N_n}{H}_i$ . For some constant $0<C<\infty$ such that

$H/{\min }_{1\le i\le 1+N_n}{H}_i\le C$ and ${\max }_{1\le i\le N_n}\left\{H_{i+1}-H_i\right\}=o\left(N_n^{-1}\right)$ .

(C6) $E\left[{{\varphi }^{\prime }}_h\left(\epsilon \left(t\right)\right)\right]=0$ , and $E\left[{{\varphi }^{″}}_h\left(\epsilon \left(t\right)\right)\right]>0$ , for all $t$ and any $h>0$ .

Remark 1 The conditions (C1) and (C2) are the consistency and asymptotic results derived by Yao et al. [9] for estimating the FPC score ${\widehat{\xi }}_{ik}$ . Therefore, in this paper, the function FPCA scores ${\xi }_{ik}$ can be used as the target variable, and the estimated coefficients can be projected back using the univariate basis functions and the eigenfunctions of the functional response to recover the influence of the scalar predictor. The conditions (C3)-(C5) are standard in spline estimation literature, where the number of nodes $J_n$ for B-splines should satisfy

$\underset{n\to \infty }{\lim }{J}_n/n^{1/\left(2n+1\right)}=C$ , where $C$ is some nonzero constant. Condition (C6) is introduced by Yao et al. [13] as an additional tuning parameter to automatically

select based on observed data, ensuring that the objective function has a local minimum.

Theorem 1 Suppose that the conditions (C1)-(C6) are satisfied, if ${\lambda }_n\to \infty$ and $n^{r/\left(2r+1\right)}{\lambda }_n\to \infty$ as $n\to \infty$ , then we have

$\Vert {\widehat{\gamma }}_i\left(t\right)-{\gamma }_i\left(t\right)\Vert =O_p\left(n^{-r/\left(2r+1\right)}\right),j=1,\cdots ,p+q.$

Remark 2 Based on the convergence rate in Theorem 1, we can see that our proposed method provides the optimal rate of convergence $O_p\left(n^{-r/\left(2r+1\right)}\right)$ established by Stone [14] for estimating the functional coefficients. When $r=2$ , the optimal convergence rate becomes $O_p\left(n^{-2/5}\right)$ . In the studies of Fan and Li [11] and Wang et al. [7], similar variable selection methods have been proven to be consistent and asymptotically normal. Therefore, the asymptotic property of our method still holds if the above conditions are satisfied.

3. Computational Algorithm and Tuning Parameter Selection

3.1. Computational Algorithm

In this paper, we develop a new iterative algorithm using the local linear approximation (LLA) [15] to the group SCAD method and combining it with the exponential square loss function. The minimization problem of the objective function is solved through iterative reweighted least squares. Since the SCAD penalty function itself is non-convex, the LLA approximation is used to transform it into a weighted LASSO, which penalizes different groups of variables with different weights, and the weights are updated at each iteration, thereby enhancing the robustness of variable selection. The ESL function is applied to the residuals to reduce the effect of outliers on the objective function.

In the LLA method, we perform a local linear approximation on the SCAD

penalty function: $p_{\lambda }\left({\Vert {\gamma }_l\Vert }_2\right)\approx {p^{\prime }}_{\lambda }\left({\Vert {\gamma }_l^{\left(m\right)}\Vert }_2\right){\Vert {\gamma }_l\Vert }_2$ ,where:

$w_l={p^{\prime }}_{\lambda }\left({\Vert {\gamma }_l^{\left(m\right)}\Vert }_2\right)=\{\begin{array}{l}\lambda ,\text{if}{\Vert {\gamma }_l^{\left(m\right)}\Vert }_2\le \lambda ,\\ \frac{a\lambda -\left|\gamma \right|}{a-1},\text{if}\lambda <{\Vert {\gamma }_l^{\left(m\right)}\Vert }_2\le a\lambda ,\\ 0,\text{if}{\Vert {\gamma }_l^{\left(m\right)}\Vert }_2>a\lambda ,\end{array}$

Thus, the SCAD penalty function is approximated as a weighted group LASSO:

${\displaystyle \sum _{l=0}^q{p}_{\lambda }\left({\Vert {\gamma }_l\Vert }_2\right)\approx {\displaystyle \sum _{l=0}^q{w}_l{\Vert {\gamma }_l\Vert }_2}}$ where $w_l$ is the weight updated at each iteration, let$W_{SCAD}=diag\left(\frac{w_1}{{\Vert {\gamma }_l\Vert }_2},\cdots ,\frac{w_l}{{\Vert {\gamma }_l\Vert }_2}\right)$ . Then our objective Function (7) becomes:

$L\left(\gamma \right)={\displaystyle \sum _{i=1}^n{\varphi }_h\left({\xi }_i-Z_i\gamma \right)+n{\displaystyle \sum _{l=0}^q{w}_l{\Vert {\gamma }_l\Vert }_2}}$ (8)

Derive this equation with respect to $\gamma$ and set it to 0, resulting in the following equation:

$-{\displaystyle \sum _{i=1}^n{Z}_i^{\text{T}}{{\varphi }^{\prime }}_h\left({\xi }_i-Z_i\gamma \right)+n{W}_{SCAD}\gamma =0.}$ (9)

Under ESL loss, let the weight matrix $W_{ESL}$ be calculated from the derivatives

of ESL: $w_i=\frac{{{\varphi }^{\prime }}_h\left(r_i\right)}{r_i}=\frac{2}{h}\exp \left(-\frac{r_i^2}{h}\right)$ , where $r_i={\xi }_i-Z_i^{\text{T}}\gamma$ , $i=1,\cdots ,n$ . With

$W_{ESL}=diag\left(w_1,\cdots ,w_n\right)$ . Then Equation (9) becomes:

$\left(Z^{\text{T}}{W}_{ESL}Z+n{W}_{SCAD}\right)\gamma =Z^{\text{T}}{W}_{ESL}\xi$ (10)

Update $\gamma$ as in each iteration step:

$\widehat{\gamma }={\left(Z^{\text{T}}{W}_{ESL}Z+n{W}_{SCAD}\right)}^{-1}{Z}^{\text{T}}{W}_{ESL}\xi .$

Therefore, the iterative calculation process for $\widehat{\gamma }={\left({\widehat{\gamma }}_0^{\text{T}},\cdots ,{\widehat{\gamma }}_q^{\text{T}}\right)}^{\text{T}}$ is as follows:

Step 1: Obtain initial values using least squares estimation ${\widehat{\gamma }}^{\left(0\right)}={\left(Z^{\text{T}}Z\right)}^{-1}{Z}^{\text{T}}\xi$.

Step 2: Given ${\widehat{\gamma }}^{\left(m\right)}$ : Calculate the residual weight matrix ${\widehat{W}}^{\left(m\right)}$ , where the residual is ${\widehat{\gamma }}_i^{\left(m\right)}={\xi }_i-Z_i^{\text{T}}{\widehat{\gamma }}^{\left(r\right)}$.

Step 3: Given ${\widehat{\gamma }}^{\left(m\right)}$ , calculate the LLA approximation weights for SCAD ${\widehat{w}}^{\left(m\right)}$ .

Step 4: Update ${\widehat{\gamma }}^{\left(m+1\right)}$ through iterative reweighted least squares: ${\widehat{\gamma }}^{\left(m+1\right)}={\left(Z^{\text{T}}{W}_{ESL}Z+n{W}_{SCAD}\right)}^{-1}{Z}^{\text{T}}{W}_{ESL}\xi$.

Step 5: Repeat steps 2 to 4 until convergence. If $\Vert {\stackrel{⌢}{\gamma }}^{\left(m+1\right)}-{\stackrel{⌢}{\gamma }}^{\left(m\right)}\Vert$ is less than a predetermined threshold, the final estimate of γ is obtained.

3.2. The Choice of Tuning Parameters

The parameter $K$ is the number of eigenbases in the function principal component analysis, and the parameter $J$ is the number of B-spline bases used to approximate the nonlinear parameters as a linear model. To select the number of these two basis functions more robustly, we use the following weighted generalized cross-validation (GCV) statistic:

$GCV\left(K,J\right)=\frac{{\left(\xi -Z\eta \right)}^{\text{T}}\stackrel{⌢}{W}\left(Y-Z\eta \right)/nk}{{\left(1-t\left(K,J\right)/nk\right)}^2},$

where $t\left(K,J\right)=trace\left\{Z{\left(Z^{\text{T}}\widehat{W}Z\right)}^{-1}{Z}^{\text{T}}\stackrel{⌢}{W}\right\}$. Thus, the optimal number of $K$ and

$J$ is obtained by $\left(\widehat{K},\widehat{J}\right)=\min GCV\left(K,J\right)$ . In order to achieve consistency in variable selection, we use the extended Bayesian information criterion (EBIC) [16] of a corresponding Gaussian likelihood for choosing the penalty parameter $\lambda$ . For the parameter $a$ in the SCAD penalty function, following the suggestion of Fan and Li [11], we adopt $a=3.7$ for implementation. The tuning parameter $h$ impacts the efficiency of the proposed estimator. To attain high efficiency, we adopt the approach of Yao et al. [13] to select $h$ .

4. Simulation and Application

4.1. Simulation Studies

In this subsection, we implement several simulation studies to investigate the finite sample performance of the proposed method. We take the sample sizes of $n=50,100,200$ and generate the data from the nonlinear FOSR model:

$Y_i\left(t\right)=\mu \left(t\right)+{\displaystyle \sum _{j=1}^2{X}_{ij}{\beta }_j\left(t\right)+{\displaystyle \sum _{l=1}^{20}{\theta }_l\left(M_{il},t\right)+{\epsilon }_i\left(t\right),}}$

Here, we have the coefficient functions given by $\mu \left(t\right)=8\sin \left(\pi t/50\right)$ , ${\beta }_1\left(t\right)=3+5t/100$ , ${\beta }_2\left(t\right)=4\sin \left(\pi t/50\right)+4\cos \left(\pi t/50\right)$ , the nonparametric functions are given by ${\theta }_1\left(x,t\right)=2x^3\exp \left(t/100\right)$ , ${\theta }_2\left(x,t\right)=5\left(x+x^3\right)\sin \left(\pi t/100\right)$ ,${\theta }_3\left(x,t\right)=16\sin \left(xt/50\right)$ . The rest of the nonlinear functions ${\theta }_l\left(x,t\right)$ are set to zero. The exogenous covariates $X_{ij}~Unif\left(-2j,-2i\right)$ , and the predictors of interest$M_{ij}~N\left(0,1^2\right)$ . The error process ${\epsilon }_i\left(t\right)$ is generated as ${\epsilon }_i\left(t\right)={\xi }_{i1}\cos \left(t\right)+{\xi }_{i2}\sin \left(t\right)+N\left(0,1^2{I}_m\right)$ , where ${\xi }_{i1}~N\left(0,{0.5}^2\right)$ and ${\xi }_{i2}~N\left(0,{0.75}^2\right)$ . The response $Y_i\left(t\right)$ is observed on a grid of $m=80$ . The values of $K$ and $J$ (number of B-spline basis) were chosen from a two-dimensional grid with $K\in \left\{3,4,5,6,7\right\}$ and $J\in \left\{5,7,9\right\}$ .

In order to examine the accuracy of the proposed method, we use the standard deviation of the integrated squared error (ISE), which is defined as

$\text{ISE}={\displaystyle \sum _{j=1}^p{\displaystyle \underset{T}{\int }{\left\{{\widehat{\eta }}_j\left(t\right)-{\eta }_j\left(t\right)\right\}}^2\text{d}t}}$. Besides, we use the mean square prediction error

(MSPE) to assess the accuracy of prediction. MSPE is given by

$\text{MSPE}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\displaystyle \underset{T}{\int }{\left\{{\widehat{Y}}_i^{*}\left(t\right)-{\displaystyle \sum _{j=1}^p{\eta }_j\left(t\right)Z_{ij}^{*}}\right\}}^2\text{d}t}}$, where the independent samples

$Z_i=\left[X_i,M_i\right]$ are generated for the test set, ${\widehat{Y}}_i^{*}\left(t\right)={\displaystyle \sum _{j=1}^2{X}_{ij}{\beta }_j\left(t\right)+{\displaystyle \sum _{l=1}^{20}{\theta }_l\left(M_{il},t\right)}}$

are estimated from the training data.

To demonstrate the effectiveness of the algorithm described in Section 4, we generated data with a sample size of $n=100$ and used cubic splines of order $q=4$ to approximate the nonlinear predictor variables. The tuning parameters $K,J$ and $h$ selected by the proposed process are approximately 5, 7, and 3.5, respectively.

In order to show the robustness of the proposed approach, we simulated three types of outlier contamination: outliers in the function response, outliers in the predictor variables, and outliers in both the function response and scalar predictor. To introduce outliers in the function response, we assumed that 10% of the response curves of the original sample were contaminated by outlier curves, denoted as ${\widehat{Y}}_i^{*}\left(t\right)$ , so that the contaminated function response data $Y_i^0\left(t\right)=\left(1-R_i\right)Y_i\left(t\right)+R_i{Y}_i^{*}\left(t\right)$ is treated as a numerical simulation, where $R_i$ is a Bernoulli random variable such that $P\left(R_i=1\right)=0.1$ . To introduce outliers in the predictor variables, we randomly selected subsets of the predictors from the original sample and contaminated them with peaks values. The detail mechanisms for producing the outlier response $Y_i^{*}\left(t\right)$ and outliers in scalar predictors are described as follows:

Scenario 1. No outliers. In this setting, $Y_i^{*}\left(t\right)=Y_i\left(t\right)$

Scenario 2. When the outlier is in the function response. The outlier curve $Y_i^{*}\left(t\right)$ is generated as $Y_i^{*}\left(t_k\right)=Y_i\left(t_k\right)+r_{ik}{I}_{\left[U_i,U_i+v\right]}\left(t_k\right),k=1,\cdots ,m$ , where $r_{ik}$ are random values taken from a uniform distribution on $\left[-r_u,-r_l\right]\cup \left[r_l,r_u\right]$ and $I(\cdot )$ is the indicator function, $\left[U_i,U_i+v\right]$ be an interval in $\left[0,1\right]$ for the $i\text{th}$ subject, where $U_i$ is randomly chosen from $\left[0,1-v\right]$ , and $v$ is a constant. The $r_{ik}$ control the strength of outliers while the constant $v$ determines the contamination length of each outlier curve. We set $r_l=3$ , $r_u=10$ and $v=0.5$ in the scenario.

Scenario 3. When the outliers are in the scalar predictors. We randomly select 5% of the predictors that are contaminated by the high leverage outliers $\left(X_{i1}^{*},\cdots ,X_{ip}^{*}\right)=\left(X_{i1},\cdots ,X_{ip}\right)+5$ .

Scenario 4. When the outliers are both in the function response and in the scalar predictor. We consider the contaminated error distribution $0.8N\left(0,{0.5}^2\right)+0.2Cauchy\left(0,1\right)$ , which is used to produce outliers with a heavy-tailed error distribution for data set.

For each scenario, our method is compared with the other two methods. Since the method proposed in this paper is based on the exponential square loss function ${\varphi }_h\left(t\right)=1-\exp \left(-t^2/h\right)$ , our model is called ESL. When ${\varphi }_h\left(x\right)$ is the square loss function ${\varphi }_h\left(x\right)=x^2$ , it is the least squares (LS) method proposed by Wang et al. [12]. When ${\varphi }_h\left(x\right)=\left|x\right|$ , it leads to the least absolute deviation method (LAD). For each scenario, 100 repetitions are used, and the numerical results are recorded in the following four tables. First, as shown in Table 1, when there are no outliers in the dataset and the error distribution is normal, the ISE and MSPE values of the ESL estimator are close to those of the LS estimator. This means that for datasets without outliers, the performance of these two methods is comparable. Moreover, the LS estimator performs best because its sample mean and the standard deviations of ISE and MSPE are smaller than those of the proposed ESL and LAD estimators. Second, from Tables 2-4, we can see that the ESL estimator produces smaller ISE and MSPE values than the LS and LAD estimators. This indicates that for datasets with outliers, the proposed ESL method outperforms the LS and LAD methods. This is mainly because when there are some very large outliers in the dataset, the proposed ESL method places more weight on the “most likely” data around the true value, thereby achieving a robust estimator. Finally, the ISE and MSPE values of the ESL estimator decrease with increasing sample size $n$ , which confirms the consistency result established in Theorem 1. In summary, the proposed ESL method has certain advantages in variable selection and estimation.

4.2. Application

In this subsection, we apply the proposed method to analyze the real DTI dataset collected from the NIH Alzheimer’s Disease Neuroimaging Initiative (ADNI) study, which is available at http://adni.loni.usc.edu/. The DTI data includes 213 subjects. For each subject, we calculate the fractional anisotropy (FA) curve along the midsagittal skeleton of the corpus callosum at m = 83 positions, as shown in

Table 1. Simulation results of the averages (×10⁻²) and the standard deviations (×10⁻²) of integrated squared error (ISE) and mean square prediction error (MSPE) under Scenario 1.

n	Methods	ISE		MSPE
		Aver.	std	Aver.	std
50	ESL	0.532	0.214	1.763	0.473
	LS	0.438	0.186	1.516	0.432
	LAD	0.732	0.278	1.928	0.675
100	ESL	0.273	0.106	0.316	0.127
	LS	0.258	0.103	0.308	0.119
	LAD	0.376	0.145	0.531	0.193
200	ESL	0.124	0.032	0.078	0.026
	LS	0.121	0.032	0.076	0.025
	LAD	0.183	0.054	1.117	0.041

Table 2. Simulation results of the averages (×10⁻²) and the standard deviations (×10⁻²) of integrated squared error (ISE) and mean square prediction error (MSPE) under Scenario 2.

n	Methods	ISE		MSPE
		Aver.	std	Aver.	std
50	ESL	0.574	0.203	1.837	0.413
	LS	3.486	2.693	13.723	7.516
	LAD	0.904	0.523	3.841	1.967
100	ESL	0.241	0.098	0.564	0.139
	LS	2.773	1.392	5.103	2.877
	LAD	0.339	0.124	1.072	0.438
200	ESL	0.134	0.045	0..341	0.113
	LS	1.747	1.219	2.914	1.485
	LAD	0.202	0.074	0.517	0.149

Table 3. Simulation results of the averages (×10⁻²) and the standard deviations (×10⁻²) of integrated squared error (ISE) and mean square prediction error (MSPE) under Scenario 3.

n	Methods	ISE		MSPE
		Aver.	std	Aver.	std
50	ESL	0.728	0.165	2.137	0.413
	LS	87.72	31.46	486.7	103.2
	LAD	14.47	16.82	63.59	89.74
100	ESL	0.341	0.105	0.536	0.153
	LS	83.45	16.82	253.6	41.67
	LAD	12.97	8.183	39.26	23.92
200	ESL	0.179	0.063	0.218	0.072
	LS	76.81	11.74	165.3	17.23
	LAD	11.34	3.769	14.35	6.138

Table 4. Simulation results of the averages (×10⁻²) and the standard deviations (×10⁻²) of integrated squared error (ISE) and mean square prediction error (MSPE) under Scenario 4.

n	Methods	ISE		MSPE
		Aver.	std	Aver.	std
50	ESL	0.931	0.273	2.928	0.773
	LS	673.5	703.4	2156	2013
	LAD	32.18	132.6	107.3	388.5
100	ESL	0.452	0.156	0.763	0.248
	LS	756.1	514.9	1235	891.4
	LAD	15.14	14.42	42.37	28.36
200	ESL	0.237	0.083	0.153	0.076
	LS	548.6	413.7	409.6	314.9
	LAD	12.16	6.297	13.54	9.433

Figure 1. We refer to Smith et al. [17] for a detailed description of the application of this data.

Figure 1. For 213 subjects, fractional anisotropy (FA) curves were calculated for $m=83$ locations along the sagittal skeleton in the corpus callosum.

We applied the model proposed in this paper,

$Y_i\left(t\right)=\mu \left(t\right)+{\displaystyle \sum _{j=1}^p{X}_{ij}{\beta }_j\left(t\right)}+{\displaystyle \sum _{l=1}^q{\theta }_l\left(M_{il},t\right)}+{\epsilon }_i\left(t\right)$ , to this dataset. Here, $X_{ij}$ are the

categorical predictors (dummy encoded) with linear effects ${\beta }_j\left(t\right)$ , and $M_{il}$ are the continuous predictors. In this study, we regarded gender (coded by a dummy variable indicating males), handedness (coded by a dummy variable indicating left-handedness), education level, Alzheimer’s disease (AD) status, and mild cognitive impairment (MCI) status as categorical predictors, and age as a continuous variable. We centered all the functional responses and scalar predictors at zero mean. Then, we use the proposed ESL method to analyze the ADNI data, which approximates the nonlinear scalar function using cubic splines of order $q=4$ . The first five panels of Figure 2 show the estimated function coefficients of these five linear variables. We note that the function coefficient corresponding to gender takes positive values for the most part, while the functional coefficients of the other four predictors take negative values. This indicates that males tend to have higher FA values than females, and that patients with AD and MCI have lower FA values. Left-handedness or higher education level may lead to smaller FA values. The last graph shows the effect of continuous predicted age on FA values, which decrease with age and decline more rapidly than FA values during normal aging. These results are consistent with the previous analysis by Cai et al. [18].

Figure 2. Estimated functional coefficients (×10⁻³) for gender, handedness, education, AD, age and MCI for the ADNI data.

To evaluate the predictive performance of the model, we randomly divided 213 samples into a training set containing 149 samples and a testing set containing 64 samples. We use the training set to estimate the function coefficients and then predict the response in the test set. We use mean squared prediction error (MSPE) to verify the quality of the prediction, and the smaller the MSPE, the better the prediction. Based on 100 random splits, we obtained LS estimation average MSPEs and their standard deviations of $0.74\times {10}^{-2}$ and $0.21\times {10}^{-2}$ , respectively. We obtained the average MSPEs and standard deviation of our proposed ESL estimator, which are $0.69\times {10}^{-2}$ and $0.18\times {10}^{-2}$ , respectively. These results indicate that our proposed method is preferred when considering the robustness of predictive performance.

5. Conclusions

In order to better explain the behavior of the function response over time, this paper develops a robust variable selection procedure for the nonlinear FOSP model containing continuous type. The function response is projected into a low-dimensional principal component space by using functional principal component analysis, and the principal component scores are used for variable selection and regression modeling. A robust variable selection for all scalar predictors is achieved by combining the ESL function with the group SCAD regularization method. To implement the computational procedure of the proposed method, we propose an iterative algorithm based on the LLA method. The tuning parameters are selected through a data driven program to achieve the robustness of variable selection. Simulation studies and real data analysis show that when there are outliers in the function response or scalar predictors, the proposed method can select the relevant predictors.

There are still some problems worth solving in the future research. In this paper, all linear covariates are divided into a group, which can either be selected or rejected at the same time. This condition is too strict and has limitations in practical application. Therefore, the variable selection of paper for scalar predictors with nonlinear covariates needs further study. When nonlinear covariates have high dimensional characteristics, robust variable selection for functional response regression is a significant topic for continued research.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Barber, R.F., Reimherr, M. and Schill, T. (2017) The Function-on-Scalar LASSO with Applications to Longitudinal GWAS. Electronic Journal of Statistics, 11, 1351-1389. https://doi.org/10.1214/17-ejs1260
[2]	Chen, Y., Goldsmith, J. and Ogden, R.T. (2016) Variable Selection in Function-on-Scalar Regression. Stat, 5, 88-101. https://doi.org/10.1002/sta4.106
[3]	Goldsmith, J., Liu, X., Jacobson, J.S. And Rundle, A. (2016) New Insights into Activity Patterns in Children, Found Using Functional Data Analyses. Medicine & Science in Sports & Exercise, 48, 1723-1729. https://doi.org/10.1249/mss.0000000000000968
[4]	Fan, Z. and Reimherr, M. (2017) High-Dimensional Adaptive Function-on-Scalar Regression. Econometrics and Statistics, 1, 167-183. https://doi.org/10.1016/j.ecosta.2016.08.001
[5]	Parodi, A. and Reimherr, M. (2018) Simultaneous Variable Selection and Smoothing for High-Dimensional Function-On-Scalar Regression. Electronic Journal of Statistics, 12, 4602-4639. https://doi.org/10.1214/18-ejs1509
[6]	Ghosal, R. and Maity, A. (2021) Variable Selection in Nonlinear Function-on-Scalar Regression. Biometrics, 79, 292-303. https://doi.org/10.1111/biom.13564
[7]	Wang, X., Jiang, Y., Huang, M. and Zhang, H. (2013) Robust Variable Selection with Exponential Squared Loss. Journal of the American Statistical Association, 108, 632-643. https://doi.org/10.1080/01621459.2013.766613
[8]	Ash, R.B. and Gardner, M.F. (1978) Topics in Stochastic Processes. Academic Press.
[9]	Yao, F., Müller, H. and Wang, J. (2005) Functional Data Analysis for Sparse Longitudinal Data. Journal of the American Statistical Association, 100, 577-590. https://doi.org/10.1198/016214504000001745
[10]	Hall, P. and Hosseini-Nasab, M. (2005) On Properties of Functional Principal Components Analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology, 68, 109-126. https://doi.org/10.1111/j.1467-9868.2005.00535.x
[11]	Fan, J. and Li, R. (2001) Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties. Journal of the American Statistical Association, 96, 1348-1360. https://doi.org/10.1198/016214501753382273
[12]	Wang, L., Chen, G. and Li, H. (2007) Group SCAD Regression Analysis for Microarray Time Course Gene Expression Data. Bioinformatics, 23, 1486-1494. https://doi.org/10.1093/bioinformatics/btm125
[13]	Yao, W., Lindsay, B.G. and Li, R. (2012) Local Modal Regression. Journal of Nonparametric Statistics, 24, 647-663. https://doi.org/10.1080/10485252.2012.678848
[14]	Stone, C.J. (1982) Optimal Global Rates of Convergence for Nonparametric Regression. The Annals of Statistics, 10, 1040-1053. https://doi.org/10.1214/aos/1176345969
[15]	Zou, H. and Li, R. (2008) One-Step Sparse Estimates in Nonconcave Penalized Likelihood Models. The Annals of Statistics, 36, 1509-1533. https://doi.org/10.1214/009053607000000802
[16]	Chen, J. and Chen, Z. (2008) Extended Bayesian Information Criteria for Model Selection with Large Model Spaces. Biometrika, 95, 759-771. https://doi.org/10.1093/biomet/asn034
[17]	Smith, S.M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T.E., Mackay, C.E., et al. (2006) Tract-Based Spatial Statistics: Voxelwise Analysis of Multi-Subject Diffusion Data. NeuroImage, 31, 1487-1505. https://doi.org/10.1016/j.neuroimage.2006.02.024
[18]	Cai, X., Xue, L., Pu, X. and Yan, X. (2020) Efficient Estimation for Varying-Coefficient Mixed Effects Models with Functional Response Data. Metrika, 84, 467-495. https://doi.org/10.1007/s00184-020-00776-0