Remodeling and Estimation for Sparse Partially Linear Regression Models

Abstract

When the dimension of covariates in the regression model is high, one usually uses a submodel as a working model that contains significant variables. But it may be highly biased and the resulting estimator of the parameter of interest may be very poor when the coefficients of removed variables are not exactly zero. In this paper, based on the selected submodel, we introduce a two-stage remodeling method to get the consistent estimator for the parameter of interest. More precisely, in the first stage, by a multistep adjustment, we reconstruct an unbiased model based on the correlation information between the covariates; in the second stage, we further reduce the adjusted model by a semiparametric variable selection method and get a new estimator of the parameter of interest simultaneously. Its convergence rate and asymptotic normality are also obtained. The simulation results further illustrate that the new estimator outperforms those obtained by the submodel and the full model in the sense of mean square errors of point estimation and mean square prediction errors of model prediction.

1. Introduction

Consider the following partially linear regression model: 𝑌=𝛽^𝑇𝑋+𝛾^𝑇𝑍+𝑓(𝑈)+𝜀, (1) where 𝑌 is a scalar response variable, 𝑋 and 𝑍 are, respectively, 𝑝-dimensional and 𝑞-dimensional continuous-valued covariates with 𝑝 being finite and 𝑝≪𝑞, 𝛽 is the parameter vector of interest and 𝛾 is the nuisance parameter vector which is supposed to be sparse in the sense that ‖𝛾‖₂ is small, 𝑓(⋅) is an unknown function satisfying 𝐸𝑓(𝑈)=0 for identification, 𝜀 is the random error satisfying 𝐸(𝜀∣𝑋,𝑍,𝑈)=0. For simplicity, we assume that 𝑈 is univariate. Let (𝑌_𝑖,𝑋_𝑖,𝑍_𝑖,𝑈_𝑖), 𝑖=1,…,𝑛, be i.i.d. observations of (𝑌,𝑋,𝑍,𝑈) obtained from the above model.

A feature of the model is that the parametric part contains both the parameter vector of interest and nuisance parameter vector. The reason for this coefficient separation is as follows. In practice we often use such a model to distinguish the main treatment variables of interest from the state variables. For instance, in a clinical trial, 𝑋 consists of treatment variables and can be easily controlled, 𝑍 is a vector of many clinical variables, such as patient ages and body weights. The variables in 𝑍 may have an impact on 𝑌 but are not of primary interest and the effects may be small. In order to make up for potentially nonnegligible effects on the response 𝑌, the nuisance covariate 𝑍 are introduced into model (1); see Shen et al. [1]. Model (1) contains all relevant covariates and in this paper we call it full model.

The purpose of this paper is to estimate 𝛽, the parameter of interest, when 𝛾^𝑇𝑍 is removed from the model. The main idea is remodeling based on the following working model: 𝑌=𝛽^𝑇𝑋+𝑓(𝑈)+𝜂. (2) As is known, 𝐸(𝜂∣𝑋=𝑥,𝑈=𝑢) is a nonzero function if 𝛾^𝑇𝐸(𝑍∣𝑋,𝑈)≠0, which relies on two elements, one is 𝐸(𝑍∣𝑋,𝑈), related with the correlation between the covariates of 𝑍 and (𝑋,𝑈), the other is 𝛾, determined by the nuisance parameter in the removed part. Thus the least squares estimator based on model (2) may be inconsistent. In the following, we will make use of the above two elements. Specifically, in the first stage, we shall construct a remodeled model by a multistep-adjustment to correct the submodel bias based on the correlation information between the covariates. This adjustment is motivated by Gai et al. [2]. In the paper, they proposed a nonparametric adjustment by adding a univariate nonparametric estimation to the working model (2), and it can dramatically reduce the bias of the working model. But this only holds in a subset of the covariates, although the subset may be fairly large. In order to obtain a globally unbiased working model for linear regression model, Zeng et al. [3] adjusted the working model by multiple steps. Because only those variables in 𝑍 correlated with (𝑋,𝑈) may have impact on estimation of 𝛽, in each step a univariate nonparametric part was added to the working model and consequently a globally unbiased working model was obtained.

However, when many components of 𝑍 are correlated with (𝑋,𝑈), the number of nonparametric functions added in the above working model is large. Such a model is improper in practice. Thus, in the second stage, we further simplify the above adjusted model by a semiparametric variable selection procedure proposed by Zhao and Xue [4]. Their method can select significant parametric and nonparametric components simultaneously under sparsity condition for semiparametric varying coefficient partially linear models. The relevant papers include Fan and Li [5], Wang et al. [6, 7], among others. After two-stage remodeling, the final model is conditionally unbiased. Based on this model, the estimation and model prediction are significantly improved.

The rest of this paper is organized as follows. In Section 2, a multistep adjustment and remodeled models are firstly proposed, then the models are further simplified via the semiparametric SCAD variable selection procedure. A new estimator of the parameter of interest based on the simplified model is derived, its convergence rate and asymptotic normality are also obtained. Simulations are given in Section 3. A short conclusion and some remarks are contained in Section 4. Some regular conditions and theoretical proofs are presented in the appendix.

2. New Estimator for the Parameter of Interest

In this paper, we suppose that covariate 𝑍 has zero mean, 𝑝 is finite and 𝑝≪𝑞, 𝐸(𝜀∣𝑋,𝑍,𝑈)=0 and Var(𝜀∣𝑋,𝑍,𝑈)=𝜎². We also assume that covariates 𝑋 and 𝑈 and parameter 𝛽 are prespecified, so that the submodel (2) is a fixed model.

2.1. Multistep-Adjustment by Correlation

In this subsection, we first adjust the submodel to be conditionally unbiased by a multistep-adjustment.

When 𝑍 is normally distributed, the principal component analysis (PCA) method will be used. Let Σ_𝑍 be the covariance matrix of 𝑍, then there exists an orthogonal 𝑞×𝑞 matrix 𝑄 such that 𝑄Σ_𝑍𝑄^𝑇=Λ, where Λ is the diagonal matrix diag(𝜆₁,𝜆₂,…,𝜆_𝑞) with 𝜆₁≥𝜆₂≥⋯≥𝜆_𝑞≥0 being eigenvalues of Σ_𝑍. Denote 𝑄^𝑇=(𝜏₁,𝜏₂,…,𝜏_𝑞) and 𝑍^(𝑗)=𝜏^𝑇_𝑗𝑍.

When 𝑍 is centered but nonnormally distributed, we shall apply independent component analysis (ICA) method. Assume that 𝑍 is generated by a nonlinear combination of independent components 𝑍^(𝑖), that is 𝑍=𝐹(𝑍), where 𝐹(⋅) is an unknown nonlinear mapping from 𝑅^𝑞 to 𝑅^𝑞, 𝑍 is an unknown random vector with independent components. By imposing some constraints on the nonlinear mixing mapping 𝐹 or the independent components 𝑍^(𝑖), the independent components 𝑍^(𝑖) can be properly estimated. See Simas Filho and Seixas [8] for an overview of the main statistical principles and some algorithms for estimating the independent components. For simplicity, in this paper we suppose that 𝑍=(𝑍⁽¹⁾,…,𝑍^(𝑞))^𝑇 with 𝑍^(𝑙)=Σ^𝑞_𝑗=1𝐹_𝑙𝑗(𝑍^(𝑗)), 𝑙=1,…,𝑞, and 𝐹_𝑙𝑗(⋅) are scalar functions.

In the above two cases, 𝑍^(𝑗)'s are independent of each other. Set 𝐾₀ to be the size of set 𝑀₀={𝑗∶𝐸(𝑍^(𝑗)∣𝑋,𝑈)≠0,1≤𝑗≤𝑞}. Without loss of generality, let 𝑀₀={1,…,𝐾₀}.

We construct the following adjusted model: 𝑌=𝛽^𝑇𝑋+^𝐾0_𝑗=1𝑔_𝑗𝑍^(𝑗)+𝑓(𝑈)+𝜁_𝐾0, (3) where 𝑔_𝑗(𝑍^(𝑗))=𝐸(𝑌−𝛽^𝑇𝑋−𝑓(𝑈)∣𝑍^(𝑗))=𝛾^𝑇𝐸(𝑍∣𝑍^(𝑗)), 𝑗=1,…,𝐾₀ and 𝜁_𝐾0=𝑌−𝛽^𝑇𝑋−𝑔₁(𝑍⁽¹⁾)−⋯−𝑔_𝐾0(𝑍^(𝐾0))−𝑓(𝑈). The model (3) is based on 𝑍's population and depends on the distributions of 𝑋, 𝑈 and 𝑍. It is easy to see that model (3) is conditionally unbiased, that is, 𝐸(𝜁_𝐾0∣𝑋,𝑈,𝑍^(𝑗),1≤𝑗≤𝐾₀)=0.

The adjusted model (3) is an additive partially linear model, in which 𝛽^𝑇𝑋 is the parametric part, 𝑓(𝑈) and 𝑔_𝑗(𝑍^(𝑗)), 𝑗=1,…,𝐾₀, are the nonparametric parts and 𝜁_𝐾0 is the random error. Compared with the submodel (2), the nonparametric parts 𝑔_𝑗(𝑍^(𝑗)), 𝑗=1,…,𝐾₀, may be regarded as bias-corrected terms for the random error 𝜂. For centered 𝑍, 𝐸(𝑔_𝑗(𝑍^(𝑗)))=0, 𝑗=1,…,𝐾₀, the nonparametric components 𝑔₁(𝑍⁽¹⁾),…,𝑔_𝐾0(𝑍^(𝐾0)) can be properly identified. In fact, centered 𝑍 can be relaxed to any 𝑍 such that satisfies 𝛾^𝑇𝐸(𝑍)=0.

When 𝑍 is centered and normally distributed, the nonparametric parts 𝑔_𝑗(𝑍^(𝑗))=𝛼_𝑗𝜏^𝑇_𝑗𝑍=𝛼_𝑗𝑍^(𝑗), 𝑗=1,…,𝐾₀. So the multistep adjusted model (3) is really a partially linear model 𝑌=𝛽^𝑇𝑋+𝛼^𝑇𝑍_𝐾0+𝑓(𝑈)+𝜁_𝐾0 (4) with 𝛼=(𝛼₁,…,𝛼_𝐾0)^𝑇 and 𝑍_𝐾0=(𝑍⁽¹⁾,…,𝑍^(𝐾0))^𝑇. Specially, when 𝑓(𝑈)≡0, the full model is a linear model, the multistep adjusted model is also a linear model 𝑌=𝛽^𝑇𝑋+𝛼^𝑇𝑍_𝐾0+𝜁_𝐾0. (5) But when the variables in 𝑍 are not jointly normal, the nonparametric parts 𝑔_𝑗 can be highly nonlinear, which are similar to the results of marginal regression; see Fan et al. [9].

2.2. Model Simplification

When the most of the features in the full model are correlated, then 𝐾₀ is very large and even is close to 𝑞. In this case, the adjusted model (3) is improper in practice, so we shall use the group SCAD regression procedure, proposed by Wang et al. [6], and the semiparametric variable selection procedure, proposed by Zhao and Xue [4], to further simplify the model.

Let 𝑠=|ℳ_∗| with ℳ_∗={1≤𝑗≤𝐾₀∶𝐸(𝑔_𝑗(𝑍^(𝑗)))²>0}, and assume that the model (3) is sparse, that is, 𝑠 is small. We define the semiparametric penalized least squares as 𝐹𝛽,𝑔(⋅),𝑓(⋅)=^𝑛_𝑖=1𝑌_𝑖−𝛽^𝑇𝑋_𝑖−^𝐾0_𝑗=1𝑔_𝑗𝑍_𝑖^(𝑗)−𝑓𝑈_𝑖²+𝑛^𝐾0_𝑗=1𝑝_𝜆𝑗𝑔_𝑗𝑍^(𝑗), (6) where ‖𝑔_𝑗(𝑍^(𝑗))‖=(𝐸(𝑔_𝑗(𝑍^(𝑗)))²)^1/2, and 𝑝_𝜆(⋅) is the SCAD penalty function with 𝜆 being a tuning parameter defined as 𝑝^′_𝜆(𝑤)𝐼()+()=𝜆𝑤≤𝜆𝑎𝜆−𝑤₊()𝜆𝐼(,𝑎−1𝑤>𝜆) (7) with 𝑎>2, 𝑤>0 and 𝑝_𝜆(0)=0. In (6), 𝑔(⋅) denotes the set {𝑔_𝑗(𝑍^(𝑗)),𝑗=1,…,𝐾₀}. Because 𝑔_𝑗 are nonparametric functions, thus they cannot be directly applied for minimization. Here we will replace 𝑓(⋅) and 𝑔(⋅) by basis function approximations. For 1≤𝑗≤𝐾₀, let {Ψ_𝑗𝑘,𝑘=1,…,𝐿} be orthogonal basis functions satisfying 𝐸Ψ_𝑗𝑘Ψ_𝑗𝑙≡_suppΨ_𝑗𝑘𝑍^(𝑗)Ψ_𝑗𝑙𝑍^(𝑗)𝑟_𝑗𝑍^(𝑗)𝑑𝑍=𝛿_𝑘𝑙=0,𝑘≠𝑙;1,𝑘=𝑙, (8) where 𝑟_𝑗(𝑍^(𝑗)) is the density function of 𝑍^(𝑗). Similarly, let {Ψ_0𝑘,𝑘=1,…,𝐿} be orthogonal basis functions satisfying the above condition which is only replaced by the support and density function of 𝑈. Denote Ψ_𝑗(𝑍^(𝑗))=(Ψ_𝑗1(𝑍^(𝑗)),…,Ψ_𝑗𝐿(𝑍^(𝑗)))^𝑇, Ψ₀(𝑈)=(Ψ₀₁(𝑈),…,Ψ_0𝐿(𝑈))^𝑇. Then 𝑔_𝑗(𝑍^(𝑗)) and 𝑓(𝑈) can be approximated by 𝑔_𝑗𝑍^(𝑗)≈𝜽^𝑇_𝑗Ψ_𝑗𝑍^(𝑗),𝑓(𝑈)≈𝜈^𝑇Ψ₀(𝑈). (9) Denote ‖𝜃_𝑗‖₂=(𝜃^𝑇_𝑗𝜃_𝑗)^1/2, invoking that 𝐸(Ψ_𝑗(𝑍^(𝑗))⋅Ψ^𝑇_𝑗(𝑍^(𝑗)))=𝐈_𝐿 the identity matrix, we get 𝐹𝛽,𝜃,𝜈=^𝑛_𝑖=1𝑌_𝑖−𝛽^𝑇𝑋_𝑖−𝜃^𝑇Ψ_𝐢−𝜈^𝑇Ψ_𝟎𝐢²+𝑛^𝐾0_𝑗=1𝑝_𝜆𝑗𝜽_𝑗2, (10) where 𝜃=(𝜃^𝑇₁,…,𝜃^𝑇_𝐾0)^𝑇, Ψ_𝐢≡Ψ(𝑍_𝑖)=Vec(Ψ₁(𝑍_𝑖⁽¹⁾),…,Ψ_𝐾0(𝑍^(𝐾0)_𝑖)), Ψ_𝟎𝐢≡Ψ₀(𝑈_𝑖).

Denote by 𝛽, 𝜽=(𝜃^𝑇₁,…,𝜃^𝑇_𝐾0)^𝑇 and 𝜈 the least squares estimators based on the penalized function (10), that is (𝛽,𝜽,𝜈)=argmin_𝛽∈𝑅^𝑝_{,𝜽𝑗∈𝑅}^𝐿_,𝜈∈𝑅^𝐿𝐹(𝛽,𝜃,𝜈). Let 𝑔_𝑗≡𝑔_𝑗(𝑍^(𝑗))=𝜽^𝑇_𝑗Ψ_𝑗(𝑍^(𝑗)) and 𝑓≡𝑓(𝑈)=𝜈^𝑇Ψ₀(𝑈), then 𝑔_𝑗 is an estimator of 𝑔_𝑗(𝑍^(𝑗)), 𝑓 is an estimator of 𝑓(𝑈).

Let ℳ_𝑛={1≤𝑗≤𝐾₀∶𝜃_𝑗≠0} and 𝐾_𝑛=|ℳ_𝑛|. For simplicity, we assume that ℳ_∗={1,2,…,𝑠} and ℳ_𝑛={1,2,…,𝐾_𝑛}. So we get the following simplified working model 𝑌=𝛽^𝑇𝑋+^𝐾𝑛_𝑗=1𝑔_𝑗𝑍^(𝑗)+𝑓(𝑈)+𝜁_𝐾𝑛, (11) where 𝑔_𝑗(𝑍^(𝑗))=𝐸(𝛾^𝑇𝑍∣𝑍^(𝑗)), 𝑗=1,…,𝐾_𝑛 and 𝜁_𝐾𝑛=𝑌−𝛽^𝑇𝑋−𝑔₁(𝑍⁽¹⁾)−⋯−𝑔_𝐾𝑛(𝑍^(𝐾𝑛))−𝑓(𝑈). Under the assumption of sparsity, the model (11) contains all of significant nonparametric functions and fully utilizes both the correlation information of covariates and the model sparsity on nuisance covariate.

If 𝑍 is centered and normally distributed with covariance matrix Σ_𝑍=𝐼_𝑞 the identity matrix, then 𝜏_𝑗=𝑒_𝑗, 𝑗=1,…,𝑞, where 𝑒_𝑗 denotes the unit vector with 1 at position 𝑗, and 𝛼 is sparse with 𝛼_𝑗=𝛾^𝑇𝜏_𝑗=𝛾_𝑗. So the model (4) is sparse. For model (5), the special case of model (4), we can apply the SCAD penalty method proposed by Fan and Li [5] to select variables in 𝑍_𝐾0 and estimate parameters 𝛼 and 𝛽 simultaneously. The selected covariate and the corresponding parameter are denoted by 𝑍_𝐾𝑛 and 𝛼_𝐾𝑛, the resulting parameter estimators are denoted by 𝛼_𝐾𝑛 and 𝛽, respectively. Finally, we can use the simplified model 𝑌=𝛽^𝑇𝑋+𝛼^𝑇_𝐾𝑛𝑍_𝐾𝑛+𝜁_𝐾𝑛 (12) for model prediction. Under the condition of sparsity, its model size is much smaller than those of the multistep adjusted model (5) and the full model (1).

2.3. Asymptotic Property of Point Estimator

Let 𝛽₀, 𝜃₀, 𝜈₀, and 𝑔_𝑗0(⋅), 𝑓₀(⋅) be the true values of 𝛽, 𝜃, 𝜈, and 𝑔_𝑗(⋅), 𝑓(⋅), respectively, in model (3). Without loss of generality, we assume that 𝑔_𝑗0(𝑍^(𝑗))=0, 𝑗=𝑠+1,…,𝐾₀, and 𝑔_𝑗0(𝑍^(𝑗)), 𝑗=1,…,𝑠, are all nonzero components.

We suppose that 𝑔_𝑗(𝑍^(𝑗)), 𝑗=1,…,𝐾₀ can be expressed as ∑^∞_𝑘=1𝜃_𝑗𝑘Ψ_𝑗𝑘(𝑍^(𝑗)) and 𝑓(𝑈) can be expressed as ∑^∞_𝑘=1𝜈_𝑘Ψ_0𝑘(𝑈), 𝜃_𝑗 and 𝜈 belong to the Sobolev ellipsoid 𝑆(𝑟,𝑀)={𝜃∶∑^∞_𝑘=1𝜃²_𝑘𝑘^2𝑟≤𝑀,𝑀>0,𝑟>0}.

The following theorem gives the consistency of the penalized SCAD estimators.

Theorem 1. Suppose that the regularity conditions (C1)–(C5) in the appendix hold and the number of terms 𝐿=𝑂_𝑝(𝑛^1/(2𝑟+1)). Then,(i)‖𝛽−𝛽₀‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛),(ii)‖𝑔_𝑗(⋅)−𝑔_𝑗0(⋅)‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛),𝑗=1,…,𝐾₀,(iii)‖𝑓(⋅)−𝑓₀(⋅)‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛),
where 𝑎_𝑛=max_𝑗{|𝑝^′_𝜆𝑗(‖𝜃_𝑗0‖₂)|∶𝜃_𝑗0≠0}.

From the last paragraph of Section 2.2 we know that, for linear regression model and normally distributed 𝑍, the multistep adjusted model (5) is a linear model. By orthogonal basis functions, such as power series, we have 𝑟=∞, then ‖𝛽−𝛽₀‖=𝑂_𝑝(𝑛^−1/2+𝑎_𝑛), implying the estimator 𝛽 has the same convergence rate as that of the SCAD estimator in Fan and Li [5].

Theorem 2. Suppose that the regularity conditions (C1)–(C6) in the appendix hold and the number of terms 𝐿=𝑂_𝑝(𝑛^1/(2𝑟+1)). Let 𝜆_max=max_𝑗{𝜆_𝑗} and 𝜆_min=min_𝑗{𝜆_𝑗}. If 𝜆_max→0 and 𝑛^{𝑟/(2𝑟+1)}𝜆_min→∞ as 𝑛→∞, then, with probability tending to 1, 𝑔_𝑗(⋅)≡0, 𝑗=𝑠+1,…,𝐾₀.

Remark 3. By Remark  1 of Fan and Li [5], we have that, if 𝜆_max→0 as 𝑛→∞, then 𝑎_𝑛→0. Hence from Theorems 1 and 2, by choosing proper tuning parameters, the variable selection method is consistent and the estimators of nonparametric components achieve the optimal convergence rate as if the subset of true zero coefficients was already known; see Stone [10].

Let 𝜃^∗=(𝜃^𝑇₁,…,𝜃^𝑇_𝑠)^𝑇 be the nonzero components of 𝜃, corresponding covariates are denoted by Ψ^∗_𝑖,𝑖=1,…,𝑛. In addition, let 1Σ=𝜎²_𝐾0𝐸𝑋𝑋^𝑇−𝐸𝑋Ψ^∗𝑇𝐸⁻¹Ψ^∗Ψ^∗𝑇𝐸Ψ^∗𝑋^𝑇−𝐸𝑋Ψ̆^𝑇₀𝐸⁻¹Ψ̆₀Ψ̆^𝑇₀𝐸Ψ̆₀𝑋^𝑇, (13) where 𝜎²_𝐾0=Var(𝜁_𝐾0𝑖) for homoscedastic case, Ψ̆₀=Ψ₀−𝐸(Ψ₀Ψ^∗𝑇)𝐸⁻¹(Ψ^∗Ψ^∗𝑇)Ψ^∗.

Theorem 4. Suppose that the regularity conditions (C1)–(C6) in the appendix hold and the number of terms 𝐿=𝑂_𝑝(𝑛^1/(2𝑟+1)). If Σ is invertible, then √𝑛𝛽−𝛽−−→₀^𝒟𝑁0,Σ⁻¹, (14) where “−−→^𝒟” denotes the convergence in distribution.

Remark 5. From Theorems 1 and 4, it can be found that the penalized estimators have the oracle property. Furthermore, the estimator of the parameter of interest has the same asymptotic distribution as that based on the correct submodel.

2.4. Some Issues on Implementation

In the adjusted model (4), 𝜏_𝑗,𝑗=1,…,𝐾₀ are used. When the population distribution is not available, they need to be approximated by estimators. When 𝑍 is normally distributed and eigenvalues 𝑟_𝑗,𝑗=1,…,𝑞 of the covariance matrix Σ_𝑍 are different from each other, then √𝑛(𝑢_𝑗−𝜏_𝑗) is asymptotically 𝑁(𝟎,𝑉_𝑗) with 𝑉_𝑗=∑^𝑞_𝑙≠𝑗(𝑟_𝑗𝑟_𝑙/(𝑟_𝑗−𝑟_𝑙)²)𝜏_𝑙𝜏^𝑇_𝑙, where 𝑢_𝑗 is the 𝑗th eigenvector of 𝑆=(1/(𝑛−1))∑^𝑛_𝑖=1(𝑍_𝑖−𝑍)(𝑍_𝑖−𝑍)^𝑇 with 𝑍=(1/𝑛)∑^𝑛_𝑖=1𝑍_𝑖; see Anderson [11]. For the case when the population size is large and comparable with the sample size, if the covariance matrix is sparse, we can use the method in Rütimann and Bühlmann [12] or Cai and Liu [13] to estimate the covariance matrix. So we can use 𝑢_𝑗 to approximate 𝜏_𝑗. When 𝜏_𝑗 in model (4) are replaced by these consistent estimators, one can see that the approximation error can be neglected without changing the asymptotic property.

The nonparametric parts 𝑔_𝑙(𝑍^(𝑙)) in the adjusted model depend on the univariate variable 𝑍^(𝑙), for 𝑙=1,…,𝐾₀. So it needs to choose the steps 𝐾₀ firstly. In real implementation, we compute all the 𝑞 multiple correlation coefficients of 𝑍^(𝑙) (𝑙=1,…,𝑞) with 𝑋 and 𝑈. Then we choose the components 𝑅={𝑍^(𝑙):|mcorr(𝑍^(𝑙),(𝑋,𝑈))|≥𝛿,𝑙=1,…,𝑞} for given small number 𝛿>0, where mcorr(𝑢,𝑉) denotes the multicorrelation coefficient between 𝑢 and 𝑉 and can be approximated by its sample form; see Anderson [11].

There are some tuning parameters needing to choose in order to implement the two-stage remodeling procedure. Fan and Li [5] showed that the SCAD penalty with 𝑎=3.7 performs well in a variety of situations. Hence, we use their suggestion throughout this paper. We still need to choose the positive integer 𝐿 for basis functions and the tuning parameter 𝜆_𝑗 of the penalty functions. Similar to the adaptive lasso of Zou [14], we suggest taking 𝜆_𝑗=𝜆/‖𝜽_𝑗⁽⁰⁾‖₂, where 𝜽_𝑗⁽⁰⁾ is initial estimator of 𝜽_𝑗 by using ordinary least squares method based on the first term in (10). So the two remaining parameters 𝐿 and 𝜆 can be selected simultaneously using the leave-one-out CV or GCV method; see Zhao and Xue [4] for more details.

3. Simulation Studies

In this section, we investigate the behavior of the newly proposed method by simulation studies.

3.1. Linear Model with Normally Distributed Covariates

The dimensions of the full model (1) and the submodel (2) are chosen to be 100 and 5, respectively. We set 𝛽=(0.5,3.5,2.5,1.5,4.0)^𝑇 and 𝛾=(𝛾₁,𝛾₂,𝟎^𝑇₅₅)^𝑇, where 𝛾₂∼Unif[−0.5,0.5]³⁰, a 30-dimensional uniform distribution on [−0.5,0.5]³⁰, and 𝛾₁ is chosen in the following ways:

Case (I). 𝛾₁∼Unif[0.5,1.0]¹⁰.

Case (II). 𝛾₁=(1.0,1.0,1.0,1.5,1.5,1.5,2.0,2.0,2.0,2.0).

We assume that (𝑋^𝑇,𝑍^𝑇)^𝑇∼𝑁((𝟏^𝑇₅,𝟎^𝑇₄₀,𝟏^𝑇₅₅)^𝑇,ΣΣ^𝑇), where Σ=𝜎_𝑖𝑗,𝜎_𝑖𝑗=𝜎_𝑗𝑖=1.0,𝑗=𝑖,𝑖=1,…,𝑝+𝑞;𝑐,𝑗=𝑖+𝑝,𝑖=1,3,…,𝑞;0,otherwise, (15) with 𝑐=0.5 or 𝑐=0.8. The error term 𝜀 is assumed to be normally distributed as 𝑁(0,0.3²).

Here we denote the submodel (2) as model (I), the multistep adjusted linear model (5) as model (II), the two-stage model (12) as model (III), and the full model (1) as model (IV). We compare mean square errors (MSEs) of the new two-stage estimator 𝛽_TS based on model (III) with the estimator 𝛽_𝑆 based on model (I), the multistep estimator 𝛽_𝑀 based on model (II), the SCAD estimator 𝛽_SCAD and the least squares estimator 𝛽_𝐹 based on model (IV). We also compare mean square prediction errors (MSPEs) of the above mentioned models with corresponding estimators.

The data are simulated from the full model (1) with sample size 𝑛=100 and simulation times 𝑚=1000. We use the sample-based PCA approximations to substitute 𝜏_𝑗's. The parameter 𝑎 in the SCAD penalty function is set to be 3.7 and 𝜆 is selected by leave-one-out CV method.

Table 1 reports the MSEs of point estimators on the parameter 𝛽 and the MSPEs of model predictions. From the table, we have the following findings: (1)  𝛽_𝐹 has the largest MSEs and 𝛽_𝑆 takes the second place, nearly all the new estimator 𝛽_TS has the smallest MSEs. (2) When 𝑐=0.5, the MSEs of 𝛽_SCAD are smaller than those of 𝛽_𝑀, while when 𝑐=0.8 they are larger than those of 𝛽_𝑀. These show that if the correlation between the covariates is strong, the MSEs of 𝛽_SCAD are larger than those of 𝛽_𝑀, the multistep-adjustment is necessary, so the estimations and model predictions based on two-stage model are significantly improved. (3) In case (I) and (II) the simulation results have the similar performance. (4) Similar to the trend of the MSEs of the five estimators, the MSPE of the two-stage adjusted model is the smallest among the mentioned five models.

Table 1 

MSEs on the parameter 𝛽 and MSPEs of the two-stage adjusted linear model (12) compared with the submodel, the SCAD-penalized model, the multistep adjusted model and the full model.

In summary, Table 1 indicates that the two-stage adjusted linear model (12) performs much better than the full model, and better than the submodel, the SCAD-penalized model and the multistep adjusted model.

3.2. Partially Linear Model with Nonnormally Distributed Covariates

The dimensions of the linear part in the full model (1) and the submodel (2) are chosen to be 50 and 5, respectively. We set 𝛽=(0.5,3.5,2.5,1.5,4.0)^𝑇, 𝛾=(𝛾₁,𝛾₂,𝟎^𝑇₂₅)^𝑇, 𝑓(𝑢)=𝑢²∗sin(3𝑢), where  𝛾₁=(0.5,0.1,0.8,0.2,0.5,0.2,0.6,0.5,0.1,0.9),  𝛾₂∼Unif[−0.3,0.3]¹⁰, a 10-dimensional uniform distribution on [−0.3,0.3].

We assume that the covariates are distributed in the following two ways.

Case (I). (𝑋^𝑇,𝑍^𝑇,𝑈)^𝑇∼𝑡(𝟎^𝑇₅₁,ΣΣ^𝑇), a 51-dimensional student distribution with degree of freedom df=5, where Σ=𝜎_𝑖𝑗,𝜎_𝑖𝑗=𝜎_𝑗𝑖=1.0,𝑗=𝑖,𝑖=1,…,𝑝+𝑞+1;0.95,𝑗=𝑖+𝑝,𝑖=1,2,…,𝑞+1;0.9,𝑗=𝑖+𝑝−2,𝑖=1,2,…,𝑞+3;0,otherwise. (16)

Case (II). 𝑋=(1/(1+𝑐))(𝑊₁+𝑐𝑉), 𝑍=(𝑍^𝑇₁,𝑍^𝑇₂,𝑍^𝑇₃,𝑍^𝑇₄)^𝑇 with 𝑍₁=(1/(1+𝑐))(𝑊₂+𝑐𝑉), 𝑍₂=𝑊₃, 𝑍₃=(1/(1+𝑐))(𝑊₄+𝑐𝑉), 𝑍₄=𝑊₅, 𝑈=𝑊₅⁽¹⁾, where 𝑊₁,𝑊₂,𝑊₃,𝑊₄∼Unif[−1.0,1.0]⁵, 𝑊₅∼Unif[−1.0,1.0]³⁰, 𝑉∼Unif[−1.0,1.0]⁵, uniform distributions on [−1.0,1.0] and constant 𝑐=0.1. All 𝑊₁, 𝑊₂, 𝑊₃, 𝑊₄, 𝑊₅, and 𝑉 are independent.

The error term 𝜀 is assumed to be normally distributed as 𝑁(0,0.3²).

Here we denote the submodel (2) as model (I)′, the multistep adjusted additive partially linear model (3) as model (II)′, the two-stage model (11) as model (III)′ and the full model (1) as model (IV)′. We compare mean square errors (MSEs) of the new two-stage estimator 𝛽_TS based on model (III)′ with the estimator 𝛽_𝑆 based on model (I)′, the estimator 𝛽_𝑀 based on model (II)′ and the least squares estimator 𝛽_𝐹 based on model (IV)′. We also compare the mean average square errors (MASEs) of the nonparametric estimators of 𝑓(⋅) and the mean square prediction errors (MSPEs) of different models with corresponding estimators.

The data are simulated from the full model (1) with sample size 𝑛=100 and simulation times 𝑚=500. We use the sample-based approximations of ICA, see Hyvärinen and Oja [15]. The parameter 𝑎 in the SCAD penalty function is set to be 3.7, the number 𝐿 and the parameter 𝜆 is selected by GCV method. We use the standard Fourier orthogonal basis as the basis functions.

Table 2 reports the MSEs of point estimators on the parameter 𝛽, the MASEs of 𝑓(⋅) and the MSPEs of model predictions. From the table, we have the following results: (1)  𝛽_𝐹 has the largest MSEs, its MSEs are much larger than the MSEs of the other estimators, and the new estimator 𝛽_TS always has the smallest MSEs. (2) The MASEs of 𝑓(⋅) have similar trend to the MSEs of the four estimators, while the differences are not very noticeable. (3) Similar to the MSEs of the estimators, the MSPEs of the two-stage adjusted model are the smallest among the four models. (4) In Case (II), the simulation results of models (I)′, (II)′ and (III)′, perform a little better than those in Case (I) because of the correlation structure among the covariates.

Table 2 

MSEs on the parameter 𝛽, MASEs of 𝑓(⋅) and MSPEs of the two-stage adjusted model (11) compared with the submodel, multistep adjusted model and the full model.

In summary, Table 2 indicates that the two-stage adjusted model (11) performs much better than the full model and the multistep adjusted model, and better than the submodel.

4. Some Remarks

In this paper, the main objective is to consistently estimate the parameter of interest 𝛽. When estimating the parameter of interest, its bias is mainly determined by the relevant variables, and its variance may be impacted by other variables. Because variable selection much relies on the sparsity of the parameter, when we directly consider the partially linear model, some irrelevant variables with nonzero coefficients may be selected in the final model. This may affect the estimation of the parameter 𝛽 on its efficiency and stability. Thus based on the prespecified submodel, a two-stage remodeling method is proposed. In the new remodeling procedure, the correlation among the covariates (𝑋,𝑍) and the sparsity of the regression structure are fully used. So the final model is sufficiently simplified and conditionally unbiased. Based on the simplified model, the estimation and model prediction are significantly improved. Generally, after the first stage the adjusted model is an additive partially linear model. Therefore, the remodeling method can be applied to partially linear regression model with linear regression model as a special case.

From the remodeling procedure, we can see that it can be directly applied to additive partially linear model, in which the nonparametric function 𝑓(𝑈) has component-wise additive form. As for general partially linear model with multivariate nonparametric function, we should resort to multivariate nonparametric estimation method. If the dimension of covariate 𝑈 is high, it may be faced with “the curse of dimensionality”.

In the procedure of model simplification, orthogonal series estimation method is used. This is only for technical convenience, because the semiparametric penalized least squares (6) can be easily transformed into parametric penalized least squares (10) and then the theoretic results are obtained. Although other nonparametric methods such as kernel and spline can be used without any essential difficulty, they can not directly achieve this goal. Compared with kernel method, it is somewhat difficult for series method to establish the asymptotic normality result for the nonparametric component 𝑓(𝑈) under primitive conditions.

Appendix

A. Some Conditions and Proofs

A.1. Regularity Conditions (C1)–(C6)

(C1)(𝑍,𝑈) has finite nondegenerate compact support, denoted as supp(𝑍,𝑈).(C2) The density function 𝑟_𝑗(𝑡) of 𝑍(𝑗) and 𝑟₀(𝑡) of 𝑈 satisfies 0<𝐿₁≤𝑟_𝑗(𝑡)≤𝐿₂<∞ on its support for 0≤𝑗≤𝐾₀ for some constants 𝐿₁ and 𝐿₂, and it is continuously differentiable.(C3)𝐺(𝑍,𝑈)=𝐸(𝑋𝑋^𝑇∣𝑍,𝑈) and 𝐸(𝜁²_𝐾0∣𝑍,𝑈) are continuous. For given 𝑧 and 𝑢, 𝐺(𝑧,𝑢) is positive definite, and its eigenvalues are bounded.(C4)sup₍𝑧_{,𝑢)∈supp(}𝑍_,𝑈)𝐸(‖𝑋‖³∣𝑍=𝑧,𝑈=𝑢)<∞, 𝐸𝑓(𝑈)=0, the first two derivatives of 𝑓(⋅) are Lipschitz continuous of order one.(C5)𝑏_𝑛=max_𝑗{𝑝_𝜆^′′_𝑗(‖𝜃_𝑗0‖₂)∶𝜃_𝑗0≠0}→0 as 𝑛→∞.(C6)liminf_𝑛→∞liminf_{‖𝜃𝑗0‖2→0}𝜆_𝑗⁻¹𝑝^′_𝜆𝑗(‖𝜃_𝑗0‖₂)>0 for 𝑗=𝑠+1,…,𝐾₀ where 𝑠 satisfies 𝛾^𝑇𝐸(𝐸(𝑍∣𝑍^(𝑗))𝐸(𝑍^𝑇∣𝑍^(𝑗)))𝛾>0 for 1≤𝑗≤𝑠; 𝛾^𝑇𝐸(𝐸(𝑍∣𝑍^(𝑗))𝐸(𝑍^𝑇∣𝑍^(𝑗)))𝛾=0 for 𝑠<𝑗≤𝐾₀.

Conditions (C1)–(C3) are some regular constraints on the covariates and condition (C4) is some constraints on the regression structure as those in Härdle et al. [16]. Conditions (C5)-(C6) are assumptions on the penalty function which are similar to those used in Fan and Li [5] and Wang et al. [7].

A.2. Proof for Theorem 1

Let 𝛿=𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛, 𝛽=𝛽₀+𝛿𝑇₁, 𝜃=𝜃₀+𝛿𝑇₂, 𝜈=𝜈₀+𝛿𝑇₃ and 𝑇=(𝑇^𝑇₁,𝑇^𝑇₂,𝑇^𝑇₃)^𝑇. Firstly, we shall prove that, ∀𝜖>0, ∃𝐶>0, 𝑃{inf_{‖𝑇‖=𝐶}𝐹(𝛽,𝜃,𝜈)>𝐹(𝛽₀,𝜃₀,𝜈₀)}≥1−𝜖.

Denote 𝐷(𝛽,𝜃,𝜈)=𝐿⁻¹{𝐹(𝛽,𝜃,𝜈)−𝐹(𝛽₀,𝜃₀,𝜈₀)}, then we have 𝐷=1𝛽,𝜃,𝜈𝐿^𝑛_𝑖=1𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖×−2𝛿𝑌_𝑖+2𝛿𝛽^𝑇₀𝑋_𝑖+𝜃^𝑇₀Ψ𝑍_𝑖+𝜈^𝑇₀Ψ₀𝑈_𝑖×𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖+𝛿²𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖²+𝑛𝐿^𝐾0_𝑗=1𝑝_𝜆𝑗𝜃_𝑗2−𝑝_𝜆𝑗𝜃_𝑗02≥−2𝛿𝐿^𝑛_𝑖=1𝜁_𝐾0𝑖+𝑅𝑍_𝑖,𝑈_𝑖×𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖+𝛿²𝐿^𝑛_𝑖=1𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖²+𝑛𝐿^𝑠_𝑗=1𝑝_𝜆𝑗𝜃_𝑗2−𝑝_𝜆𝑗𝜃_𝑗02≡𝐼₁+𝐼₂+𝐼₃, (A.1) where 𝑅(𝑍_𝑖,𝑈_𝑖)=∑^𝐾0_𝑗=1𝑅_𝑗(𝑍_𝑖)+𝑅₀(𝑈_𝑖) with 𝑅_𝑗(𝑍_𝑖)=𝑔_𝑗(𝑍_𝑖^(𝑗))−𝜃^𝑇_𝑗Ψ_𝑗(𝑍_𝑖^(𝑗)), 𝑗=1,…,𝐾₀ and 𝑅₀(𝑈_𝑖)=𝑓(𝑈_𝑖)−𝜈^𝑇Ψ₀(𝑈_𝑖).

By the conditions (C1) and (C2), the maximal squared bias of 𝑔_𝑗(𝑍^(𝑗)) is equal to 𝑅_𝑗𝑍^(𝑗)2=^∞_𝑘=𝐿+1𝜃²_𝑗𝑘≤^∞_𝑘=𝐿+1𝜃²_𝑗𝑘𝑘𝐿^2𝑟≤𝑀𝐿^−2𝑟, (A.2) so ‖𝑅_𝑗(𝑍^(𝑗))‖=𝑂(𝐿^−𝑟). Similarly, ‖𝑅₀(𝑈)‖=𝑂(𝐿^−𝑟). Then, ^𝑛_𝑖=1𝑅𝑍_𝑖,𝑈_𝑖𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖=𝑂_𝑝𝑛𝐾₀𝐿^−𝑟‖𝑇‖. (A.3) Noticing that 𝐸(𝜁_𝐾0∣𝑋,𝑍,𝑈)=0, by Zhao and Xue [4], we have 1√𝑛^𝑛_𝑖=1𝜁_𝐾0𝑖𝑇^𝑇₁𝑋_𝑖+𝑇^𝑇₂Ψ𝑍_𝑖+𝑇^𝑇₃Ψ₀𝑈_𝑖=𝑂_𝑝(‖𝑇‖). (A.4) So 𝐼₁=−2𝛿𝐿𝑂_𝑝𝑛𝐾₀𝐿^−𝑟‖𝑇‖+𝑂_𝑝√𝑛‖𝑇‖=𝑂_𝑝1+𝑛^{𝑟/(2𝑟+1)}𝑎_𝑛‖𝑇‖. (A.5)

Similarly, we have 0<𝐼₂=𝑂_𝑝𝑛𝐿⁻¹𝛿²‖𝑇‖²=𝑂_𝑝1+2𝑛^{𝑟/(2𝑟+1)}𝑎_𝑛+𝑛^{2𝑟/(2𝑟+1)}𝑎²_𝑛‖𝑇‖². (A.6) By properly choosing a sufficiently large 𝐶, 𝐼₂ dominates 𝐼₁ uniformly in ‖𝑇‖=𝐶.

Using Taylor expansion, 𝐼₃=^𝑠_𝑗=1𝑛𝐿𝑝^′_𝜆𝑗𝜃_𝑗02𝜃_𝑗02^′𝛿𝑇_2𝑗+𝑛𝑝2𝐿^{𝜆′′𝑗}𝜃_𝑗02𝜃_𝑗02^′𝛿𝑇_2𝑗²×(1+𝑜(1))≡𝐼₃₁+𝐼₃₂. (A.7) By simple calculations, we have that 𝐼₃₁≤𝑛𝐿𝛿𝑎_𝑛𝑙₁^𝑠_𝑗=1𝑇_2𝑗≤√𝑠𝑛𝐿𝛿𝑎_𝑛𝑙₁‖𝑇‖=𝑂_𝑝𝑛^{𝑟/(2𝑟+1)}𝑎_𝑛+𝑛^{2𝑟/(2𝑟+1)}𝑎²_𝑛‖𝑇‖,𝐼₃₂≤𝑛𝛿2𝐿²𝑏_𝑛𝑙₂^𝑠_𝑗=1𝑇_2𝑗²≤𝑛𝛿2𝐿²𝑏_𝑛𝑙₂‖𝑇‖², (A.8) where 𝑙₁ and 𝑙₂ are some positive constants. We can find that 𝐼₃₁ is also dominated by 𝐼₂ uniformly in ‖𝑇‖=𝐶, and under the condition (C5), we have 0<𝐼₃₂≤𝑜_𝑝1+2𝑛^{𝑟/(2𝑟+1)}𝑎_𝑛+𝑛^{2𝑟/(2𝑟+1)}𝑎²_𝑛‖𝑇‖². (A.9) Hence, by choosing a sufficiently large 𝐶, 𝑃{inf_{‖𝑇‖=𝐶}𝐹(𝛽,𝜃,𝜈)>𝐹(𝛽₀,𝜃₀,𝜈₀)}≥1−𝜖, which implies that with probability at least 1−𝜖 there exists a local minimum of 𝐹(𝛽,𝜃,𝜈) in the ball {𝛽₀+𝛿𝑇₁∶‖𝑇₁‖≤𝐶}. Denote the local minimizer as 𝛽, then 𝛽−𝛽₀=𝑂_𝑝(𝛿)=𝑂_𝑝𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛. (A.10) With the same argument as above, there exists a local minimum in the ball {𝜃₀+𝛿𝑇₂∶‖𝑇₂‖≤𝐶}, and the local minimizer 𝜃 satisfies that 𝜃−𝜃₀=𝑂_𝑝𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛. (A.11) For the nonparametric component 𝑔(⋅), noticing that 𝑔_𝑗−𝑔_𝑗0²=𝐸𝑔_𝑗𝑍^(𝑗)−𝑔_𝑗0𝑍^(𝑗)2=𝐸Ψ_𝑗𝑍^(𝑗)𝜃_𝑗−Ψ_𝑗𝑍^(𝑗)𝜃_𝑗0+𝑅_𝑗0𝑍^(𝑗)2≤2𝐸Ψ_𝑗𝑍^(𝑗)𝜃_𝑗−Ψ_𝑗𝑍^(𝑗)𝜃_𝑗0²+2𝐸𝑅_𝑗0𝑍^(𝑗)2=2𝜃_𝑗−𝜃_𝑗0^𝑇𝜃_𝑗−𝜃_𝑗0+2𝐸𝑅_𝑗0𝑍^(𝑗)2, (A.12) it is known that ‖𝑅_𝑗(𝑍^(𝑗))‖=𝑂(𝐿^−𝑟), so 𝐸𝑅_𝑗0𝑍^(𝑗)2=𝑂_𝑝𝑛^{−2𝑟/(2𝑟+1)}. (A.13) Thus, we get 𝑔_𝑗−𝑔_𝑗0=𝑂_𝑝𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛. (A.14)

Similarly, there exists a local minimizer 𝜈 satisfies that ‖𝜈−𝜈₀‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛). Then we can get ‖𝑓−𝑓₀‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}+𝑎_𝑛).

A.3. Proof for Theorem 2

When 𝜆_max→0, 𝑎_𝑛=0 for large 𝑛 by the form of 𝑝^′_𝜆(𝑤). Then by Theorem 1, it is sufficient to show that: with probability tending to 1 as 𝑛→∞, for any 𝛽, it satisfies ‖𝛽−𝛽₀‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}), 𝜃_𝑗 satisfies ‖𝜃_𝑗−𝜃_𝑗0‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}) with 𝑗=1,…,𝑠, and 𝜈 satisfies ‖𝜈−𝜈₀‖=𝑂_𝑝(𝑛^{−𝑟/(2𝑟+1)}), for some small 𝜄_𝑛=𝐶𝑛^{−𝑟/(2𝑟+1)}, 𝜕𝐹𝛽,𝜃,𝜈𝜕𝜃_𝑗>0,for0<𝜃_𝑗<𝜄_𝑛,𝑗=𝑠+1,…,𝐾₀,𝜕𝐹𝛽,𝜃,𝜈𝜕𝜃_𝑗<0,for−𝜄_𝑛<𝜃_𝑗<0,𝑗=𝑠+1,…,𝐾₀. (A.15) So the minimizer of 𝐹(𝛽,𝜃,𝜈) is obtained at 𝜃_𝑗=0,𝑗=𝑠+1,…,𝐾₀.

In fact, 𝜕𝐹𝛽,𝜃,𝜈𝜕𝜃_𝑗=−2^𝑛_𝑖=1Ψ_𝑗𝑍_𝑖^(𝑗)𝑌_𝑖−𝛽^𝑇𝑋_𝑖−𝜃^𝑇Ψ𝑍_𝑖−𝜈^𝑇Ψ₀𝑈_𝑖+𝑛𝑝^′_𝜆𝑗𝜃_𝑗2𝜃_𝑗2^′=−2^𝑛_𝑖=1Ψ_𝑗𝑍_𝑖^(𝑗)𝜁_𝐾0𝑖+𝑅𝑍_𝑖,𝑈_𝑖−2^𝑛_𝑖=1Ψ_𝑗𝑍_𝑖^(𝑗)𝑋^𝑇_𝑖𝛽₀−𝛽−2^𝑛_𝑖=1Ψ_𝑗𝑍_𝑖^(𝑗)Ψ^𝑇𝑍_𝑖𝜃₀−𝜃−2^𝑛_𝑖=1Ψ_𝑗𝑍_𝑖^(𝑗)Ψ^𝑇₀𝑈_𝑖𝜈₀−𝜈+𝑛𝑝^′_𝜆𝑗𝜃_𝑗2𝜃_𝑗2^′=𝑛𝜆_𝑗𝑂_𝑝𝜆_𝑗⁻¹𝑛^{−𝑟/(2𝑟+1)}+𝜆_𝑗⁻¹𝑝^′_𝜆𝜃_𝑗2𝜃_𝑗𝜃_𝑗2. (A.16)

Under the conditions liminf_𝑛→∞liminf_{‖𝜃𝑗0‖2→0}𝜆_𝑗⁻¹𝑝^′_𝜆𝑗⋅(‖𝜃_𝑗0‖₂)=𝐶>0 and 𝜆_𝑗𝑛^{𝑟/(2𝑟+1)}>𝜆_min𝑛^{𝑟/(2𝑟+1)}→∞, then 𝜕𝐹(𝛽,𝜃,𝜈)/𝜕𝜃_𝑗=𝑂_𝑝(𝑛𝜆_𝑗(𝜃_𝑗/‖𝜃_𝑗‖₂)). So the sign of the derivative is determined by 𝜃_𝑗.

So with probability tending to 1, 𝜃_𝑗=0, 𝑗=𝑠+1,…,𝐾₀. Then under sup_𝑍‖Ψ_𝑗(𝑍^(𝑗))‖=𝑂(1), 𝑔_𝑗(𝑍^(𝑗))=𝜃^𝑇_𝑗Ψ_𝑗(𝑍^(𝑗))≡0, 𝑗=𝑠+1,…,𝐾₀.

A.4. Proof for Theorem 4

By Theorems 1 and 2, we know that, as 𝑛→∞, with probability tending to 1, 𝐹(𝛽,𝜃,𝜈) attains the local minimum value at 𝛽  and (𝜃^∗𝑇,0)^𝑇 and 𝜈. Let 𝐹_1𝑛(𝛽,𝜃,𝜈)=𝜕𝐹(𝛽,𝜃,𝜈)/𝜕𝛽, 𝐹_2𝑛(𝛽,𝜃,𝜈)=𝜕𝐹(𝛽,𝜃,𝜈)/𝜕𝜃^∗ and 𝐹_3𝑛(𝛽,𝜃,𝜈)=𝜕𝐹(𝛽,𝜃,𝜈)/𝜕𝜈, then 1𝑛𝐹_1𝑛𝛽,𝜃^∗𝑇,0^𝑇,𝜈=1𝑛^𝑛_𝑖=1𝑋_𝑖𝑌_𝑖−𝛽^𝑇𝑋_𝑖−𝜃^∗𝑇Ψ^∗_𝑖−𝜈^𝑇Ψ_0𝑖1=0,𝑛𝐹_2𝑛𝛽,𝜃^∗𝑇,0^𝑇,𝜈=1𝑛^𝑛_𝑖=1Ψ^∗_𝑖𝑌_𝑖−𝛽^𝑇𝑋_𝑖−𝜃^∗𝑇Ψ^∗_𝑖−𝜈^𝑇Ψ_0𝑖+^𝑠_𝑗=1𝑝^′_𝜆𝑗𝜃_𝑗2𝜃_𝑗𝜃_𝑗21=0,𝑛𝐹_3𝑛𝛽,𝜃^∗𝑇,0^𝑇,𝜈=1𝑛^𝑛_𝑖=1Ψ_0𝑖𝑌_𝑖−𝛽^𝑇𝑋_𝑖−𝜃^∗𝑇Ψ^∗_𝑖−𝜈^𝑇Ψ_0𝑖=0. (A.17)(A.18)(A.19) From (A.17), it yields that 1𝑛^𝑛_𝑖=1𝑋_𝑖𝛽₀−𝛽^𝑇𝑋_𝑖+𝜃^∗₀−𝜃^∗𝑇Ψ^∗_𝑖+𝜈₀−𝜈^𝑇Ψ_0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖+𝜁_𝐾0𝑖=0, (A.20) where 𝑅^∗(𝑍_𝑖,𝑈_𝑖)=∑^𝑠_𝑗=1𝑅^∗_𝑗(𝑍_𝑖)+𝑅₀(𝑈_𝑖). Applying the Taylor expansion, we get 𝑝^′_𝜆𝑗𝜃_𝑗2=𝑝^′_𝜆𝑗𝜃_𝑗02+𝑝_𝜆^′′_𝑗𝜃_𝑗02𝜃_𝑗𝜃_𝑗2+𝑜_𝑝(1)𝜃^∗−𝜃^∗₀. (A.21) Furthermore, condition (C5) implies that 𝑝_𝜆^′′_𝑗(‖𝜃_𝑗0‖₂)=𝑜_𝑝(1), and noting that 𝑝^′_𝜆𝑗(‖𝜃_𝑗0‖₂)=0 as 𝜆_max→0, then 𝑝^′_𝜆𝑗(‖𝜃_𝑗‖₂)=𝑜_𝑝(𝜃^∗−𝜃^∗₀). So from (A.18), it yields 1𝑛^𝑛_𝑖=1Ψ^∗_𝑖𝛽₀−𝛽^𝑇𝑋_𝑖+𝜃^∗₀−𝜃^∗𝑇Ψ^∗_𝑖+𝜈₀−𝜈^𝑇Ψ_0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖+𝜁_𝐾0𝑖+𝑜_𝑝𝜃^∗₀−𝜃^∗=0. (A.22)

Let Φ_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ^∗_𝑖Ψ_𝑖^∗𝑇, Γ_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ^∗_𝑖𝑋^𝑇_𝑖 and Π_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ^∗_𝑖Ψ^𝑇_0𝑖, then we have 𝜃^∗−𝜃^∗₀=Φ_𝑛+𝑜_𝑝(1)⁻¹×Γ_𝑛𝛽₀−𝛽+Π_𝑛𝜈₀−𝜈×1𝑛^𝑛_𝑖=1Ψ^∗_𝑖𝑅^∗𝑍_𝑖+𝜁_𝐾0𝑖. (A.23) Substituting (A.23) into (A.20), it yields 1𝑛^𝑛_𝑖=1𝑋_𝑖Ψ_0𝑖−Π^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖^𝑇𝛽−𝛽₀+𝑋_𝑖−Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖^𝑇𝜈−𝜈₀+𝑜_𝑝𝛽−𝛽₀+𝑜_𝑝𝜈−𝜈₀=1𝑛^𝑛_𝑖=1𝑋_𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹+𝑜_𝑝(1)Λ_𝑛, (A.24) where Λ_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ^∗_𝑖(𝑅^∗(𝑍_𝑖,𝑈_𝑖)+𝜁_𝐾0𝑖).

From (A.19), it yields that 1𝑛^𝑛_𝑖=1Ψ_0𝑖𝛽₀−𝛽^𝑇𝑋_𝑖+𝜃^∗₀−𝜃^∗𝑇Ψ^∗_𝑖+𝜈₀−𝜈^𝑇Ψ_0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖+𝜁_𝐾0𝑖=0. (A.25) Substituting (A.23) into (A.25), it yields 1𝑛^𝑛_𝑖=1Ψ_0𝑖𝑋_𝑖−Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖^𝑇𝛽−𝛽₀+Ψ_0𝑖−Π^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖^𝑇𝜈−𝜈₀+𝑜_𝑝𝛽−𝛽₀+𝑜_𝑝𝜈−𝜈₀=1𝑛^𝑛_𝑖=1Ψ_0𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹+𝑜_𝑝(1)Λ_𝑛. (A.26)

Noting that 𝑛^𝑛−1_𝑖=1Π^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖Ψ^𝑇_0𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹Π_𝑛𝑛=0,_{𝑛−1𝑖=1}Π^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹Λ_𝑛=0, (A.27) Equation (A.26) can be rewritten as 1𝑛^𝑛_𝑖=1Ψ_0𝑖𝑋̆^𝑇_𝑖𝛽−𝛽₀+𝑜_𝑝𝛽−𝛽₀+1𝑛^𝑛_𝑖=1Ψ̆_0𝑖Ψ̆^𝑇_0𝑖𝜈−𝜈₀+𝑜_𝑝𝜈−𝜈₀=1𝑛^𝑛_𝑖=1Ψ̆_0𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹+𝑜_𝑝(1)Λ_𝑛, (A.28) where 𝑋̆_𝑖=𝑋_𝑖−Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖, Ψ̆_0𝑖=Ψ_0𝑖−Π^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖. Let Ξ_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ̆_0𝑖Ψ̆^𝑇_0𝑖, then we have 𝜈−𝜈₀=Ξ_𝑛⁻¹1𝑛^𝑛_𝑖=1Ψ_0𝑖𝑋̆^𝑇_𝑖𝛽₀−𝛽+Ξ_𝑛⁻¹1𝑛^𝑛_𝑖=1Ψ̆_0𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ^∗_𝑖Φ_𝑛⁻¹Λ_𝑛+𝑜_𝑝𝛽−𝛽₀. (A.29) Substituting (A.29) into (A.24), and noting that 𝑛^𝑛−1_𝑖=1Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖𝑋_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹Γ_𝑛𝑛=0,_{𝑛−1𝑖=1}Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹Λ_𝑛𝑛=0,_{𝑛−1𝑖=1}Ψ̆_0𝑖𝑋^𝑇_𝑖=𝑛^𝑛−1_𝑖=1Ψ_0𝑖𝑋̆^𝑇_𝑖, (A.30) it is easy to show that Φ̆_𝑛−Υ^𝑇_𝑛Ξ_𝑛⁻¹Υ_𝑛+𝑜_𝑝(1)√𝑛𝛽−𝛽₀=1√𝑛^𝑛_𝑖=1𝑋̆_𝑖−Υ^𝑇_𝑛Ξ_𝑛⁻¹Ψ̆_0𝑖×𝜁_𝐾0𝑖+𝑅^∗𝑍_𝑖,𝑈_𝑖−Ψ_𝑖^∗𝑇Φ_𝑛⁻¹+𝑜_𝑝(1)Λ_𝑛=1√𝑛^𝑛_𝑖=1𝑋_𝑖𝜁_𝐾0𝑖−1√𝑛^𝑛_𝑖=1𝑋_𝑖Ψ_𝑖^∗𝑇Φ_𝑛⁻¹+𝑜_𝑝(1)Λ_𝑛+1√𝑛^𝑛_𝑖=1𝑋_𝑖𝑅^∗𝑍_𝑖,𝑈_𝑖≡𝐼₁+𝐼₂+𝐼₃, (A.31) where Φ̆_𝑛=𝑛⁻¹∑^𝑛_𝑖=1𝑋̆_𝑖𝑋̆^𝑇_𝑖, Υ_𝑛=𝑛⁻¹∑^𝑛_𝑖=1Ψ̆_0𝑖𝑋^𝑇_𝑖, 𝑋_𝑖=𝑋̆_𝑖−Υ^𝑇_𝑛Ξ_𝑛⁻¹Ψ̆_0𝑖.

Using the Central Limit Theorem, we can obtain 𝐼−−→₁^𝒟𝑁0,𝜎²_𝐾0Σ₀, (A.32) where “−−→^𝒟” means the convergence in distribution and Σ₀=𝐸𝑋𝑋^𝑇−𝐸𝑋Ψ^∗𝑇𝐸⁻¹Ψ^∗Ψ^∗𝑇𝐸Ψ^∗𝑋^𝑇−𝐸𝑋Ψ̆^𝑇₀𝐸⁻¹Ψ̆₀Ψ̆^𝑇₀𝐸Ψ̆₀𝑋^𝑇. (A.33) In addition, noting that ∑^𝑛_𝑖=1𝑋̆_𝑖Ψ_𝑖^∗𝑇=0 and ∑^𝑛_𝑖=1Ψ̆_0𝑖Ψ_𝑖^∗𝑇=0, we have 𝐼₂=0. Furthermore, we have 𝐼₃=1√𝑛^𝑛_𝑖=1𝑋_𝑖−𝐸Γ^𝑇_𝑛𝐸⁻¹Φ_𝑛Ψ^∗_𝑖𝑅^∗𝑍_𝑖,𝑈_𝑖+1√𝑛^𝑛_𝑖=1𝐸Γ^𝑇_𝑛𝐸⁻¹Φ_𝑛−Γ^𝑇_𝑛Φ_𝑛⁻¹Ψ^∗_𝑖𝑅^∗𝑍_𝑖,𝑈_𝑖−1√𝑛^𝑛_𝑖=1Υ^𝑇_𝑛Ξ_𝑛⁻¹Ψ̆_0𝑖𝑅^∗𝑍_𝑖,𝑈_𝑖≡𝐼₃₁+𝐼₃₂+𝐼₃₃. (A.34) Invoking 𝐸{[𝑋_𝑖−𝐸(Γ^𝑇_𝑛)𝐸⁻¹(Φ_𝑛)Ψ^∗_𝑖]Ψ_𝑖^∗𝑇}=0, then by Zhao and Xue [4], we have 1√𝑛^𝑛_𝑖=1𝑋_𝑖−𝐸Γ^𝑇_𝑛𝐸⁻¹Φ_𝑛Ψ^∗_𝑖Ψ_𝑖^∗𝑇=𝑂_𝑝(1). (A.35) This together with ‖Ψ_𝑗(𝑍^(𝑗))‖=𝑂(1) and ‖𝑅(𝑍,𝑈)‖=𝑜(1), we get 𝐼₃₁=𝑜_𝑝(1). Similarly, 𝐼₃₂=𝑜_𝑝(1). Noting that (1/√𝑛)∑^𝑛_𝑖=1Υ^𝑇_𝑛Ξ⁻¹Ψ^∗_0𝑖Ψ_𝑖^∗𝑇=0, so as above, we have 𝐼₃₃=𝑜_𝑝(1). Hence, we get that 𝐼₃=𝑜_𝑝(1).

By the law of large numbers, we have (1/𝑛)∑^𝑛_𝑖=1𝑋_𝑖𝑋−→^𝑇_𝑖^𝑃Σ₀, where “−→^𝑃” means the convergence in probability. Then using the Slutsky theorem, we get √𝑛(𝛽−𝛽₀)−−→^𝒟𝑁(0,𝜎²_𝐾0Σ₀⁻¹).

Acknowledgment

Lin and Zeng’s research are supported by NNSF projects (11171188, 10921101, and 11231005) of China, NSF and SRRF projects (ZR2010AZ001 and BS2011SF006) of Shandong Province of China and K C Wong-HKBU Fellowship Programme for Mainland China Scholars 2010-11. Wang’s research is supported by NSF project (ZR2011AQ007) of Shandong Province of China.