# Penalized Quadratic Inference Function-Based Variable Selection for Generalized Partially Linear Varying Coefficient Models with Longitudinal Data

Authors: Jinghua Zhang; Liugen Xue
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

Semiparametric generalized varying coefficient partially linear models with longitudinal data arise in contemporary biology, medicine, and life science. In this paper, we consider a variable selection procedure based on the combination of the basis function approximations and quadratic inference functions with SCAD penalty. The proposed procedure simultaneously selects significant variables in the parametric components and the nonparametric components. With appropriate selection of the tuning parameters, we establish the consistency, sparsity, and asymptotic normality of the resulting estimators. The finite sample performance of the proposed methods is evaluated through extensive simulation studies and a real data analysis.


---

## Body

## 1. Introduction

Identifying the significant variables is of great significance in all regression analysis. In practice, a number of variables are available for an initial analysis, but many of them may not be significant and should be excluded from the final model in order to increase the accuracy of prediction. Various procedures and criteria, such as stepwise selection and subset selection with Akaike information criterion (AIC), Mallows Cp, and Bayesian information criterion (BIC), have been developed. Nevertheless, these selection methods suffer from expensive computational costs. Many shrinkage methods have been developed for the purpose of computational efficiency, e.g., the nonnegative garrote [1], the LASSO [2], the bridge regression [3], the SCAD [4], and the one-step sparse estimator [5]. Among those, the SCAD possesses the virtues of continuity, unbiasedness, and sparsity. There are a number of works on the SCAD estimation methods in various regression models, e.g., [6–9]. Zhao and Xue [8] proposed a variable selection method to select significant variables in the parametric components and the nonparametric components simultaneously for the varying coefficient partially linear models (VCPLMs).On the other hand, longitudinal data occurs frequently in biology, medicine, and life science, in which it is often necessary to make repeated measurements of subjects over time. The responses from different subjects are independent, but the responses from the same subject are very likely to be correlated. This feature is called “within-cluster correlation”. Qu et al. [10] proposed a method of quadratic inference functions (QIFs) to treat the longitudinal data. The QIF can efficiently take the within-cluster correlation into account and is more efficient than the generalized estimating equation (GEE) [11] approach when the working correlation is misspecified. The QIF approach has been applied to many models, including varying coefficient models (VCM) [12, 13], partially linear models (PLM) [14], varying coefficient partially linear models (VCPLMs) [15], and generalized partially linear models (GPLM) [16]. Wang et al. [13] proposed a group SCAD procedure for variable selection of VCM with longitudinal data. More recently, Tian et al. [15] proposed a QIF-based SCAD penalty for the variable selection for VCPLM with longitudinal data.As introduced in Li and Liang [17], the generalized partially linear varying coefficient model (GPLVCM) possesses the great flexibility of a nonparametric regression model and provides the explanatory power of a generalized linear regression model, which arises naturally due to categorical covariates. Many models are the special case of GPLVCM, e.g., VCM, VCPLM, PLM, and GLM. Li and Liang [17] studied variable selection for GPLVCM, where the parametric components are identified via the SCAD but the nonparametric components are selected via a generalized likelihood ratio test instead of shrinkage. In this paper, we extend the QIF-based group SCAD variable selection procedure to GPLVCM with longitudinal data, and the B-spline methods are adopted to approximate the nonparametric component in the model. With suitable chosen tuning parameters, the proposed variable selection procedure is consistent, and the estimators of regression coefficients have oracle property, i.e., the estimators of the nonparametric components achieve the optimal convergence rate, and the estimators of the parametric components have the same asymptotic distribution as that based on the correct submodel.The rest of this paper is organized as follows. In Section2, we propose a variable selection procedure for the GPLVCM with longitudinal data. Asymptotic properties of the resulting estimators and an iteration algorithm are presented in Section 3. In Section 4, we carry out simulation studies to assess the finite sample performance of the method. A real data analysis is given in Section 5 to illustrate the proposed methodology. The details of proofs are provided in the appendix.


## 2. Methodology

### 2.1. GPLVCM with Longitudinal Data

In this article, we consider a longitudinal study withn subjects and mi observations over time for the ith subject (i=1,⋯,n) for a total of N=∑i=1nmi observations. Each observation consists of a response variable Yij and the predicator variables Xij,Zij,Uij, where Xij∈Rp,Zij∈Rq and Uij is a scalar. We assume that the observations from different subjects are independent, but those within the same subject are dependent. The generalized varying coefficient partially linear model (GPLVM) with longitudinal data takes the form
(1)μij=EYij∣Xij,Zij,Uij=hXijTβ+ZijTαUij,where μij is the expectation of Yij when Xij, Zij, and Uij are given, β=β1,⋯,βpT is an unknown p×1 regression coefficient vector, h· is a known smooth link function, and αu=α1u,α2u,⋯,αquT is a q×1 unknown monotonic smooth function vector. Without loss of generality, we assume U∼U0,1.We approximateα· by B-spline basis functions Bu=B1u,⋯,BLuT with the order of M, where L=K+M+1 and K is the number of interior knots, i.e.,
(2)αku≈αk∗u=BuTγk,k=1,⋯,q,where γk=γk1,⋯,γkLT is a L×1 vector of unknown regression coefficients. Accordingly, μij is approximated by
(3)μij=EYij∣Xij,Zij,Uij=hXijTβ+ZijT·Iq⊗BUijTγ,where γ=γ1T,⋯,γqTT and “⊗” is the Kronecker product. We use the B-spline basis functions because they are numerically stable and have bounded support [18]. The spline approach also treats a nonparametric function as a linear function with the basis functions as pseudodesign variables, and thus, any computational algorithm for the generalized linear models can be used for the GPLVCMs.To incorporate the within-cluster correlation, we apply the QIFs to estimateβ and γ, respectively. Denote θ=βT,γTT, we define the extended score gNθ as follows:
(4)gNθ=1n∑i=1ngiθ=1n∑i=1nμ̇iTAi−12M1Ai−12Yi−μi⋮μ̇iTAi−12MsAi−12Yi−μi,where μ̇i=∂μi/∂θ, Ai=diagVarYi1,⋯,VarYim is the marginal variance matrix of subject Yi, and M1,⋯,Ms are the base matrices to represent the inverse of the working correlation matrix R in GEE approach. Following Qu et al. [10], we define the quadratic inference functions to be
(5)Qnθ=ngNTθΩnθ−1gNθ,where Ωnθ=1/n∑i=1ngiθgiθT. Note that Ωn depends on θ. The QIF estimate θ~ is then given by
(6)θ~=argminθQnθ.


### 2.2. Penalized QIF

In real data analysis, the true regression model is always unknown. An overfitted model lowers the efficiency of estimation while an underfitted one leads to a biased estimator. A popular approach to identify the relevant predictors while estimating the nonzero parameters and functions in model (1) simultaneously is to exert some kind of “penalty” on the original objective function. Here, we choose the smoothly clipped absolute deviation (SCAD) penalty because it has several advantages such as unbiasedness, sparsity, and continuity. The SCAD-penalized quadratic inference function (PQIF) is defined as follows:
(7)Qnpθ=Qnθ+n∑l=1ppλ1βl+n∑k=1qpλ2γkH,where γkH=γkTHγk1/2,H=hijL×L,hij=∫01BiuBjTudu and pλ is the SCAD penalty function, where the derivative is defined as
(8)p′lω=λIω⩽l+aλ−ω+a−1λIω>λ,where a>2,ω>0,pλ0=0; here, we choose a=3.7 as in [4].Note that(9)γkH=∫01γkTBuBTUγkdu1/2=∫01α∗u2du1/2.This group-wised penalization ensures that the spline coefficient vector of the same nonparametric component is treated as an entire group in model selection.Denoteθ^ to be the penalized estimator obtained by minimizing the penalized objective function of (7). Then, β^=θ∧1,⋯,θ∧pT is the estimator of the parameter β and the estimator of the nonparametric function αu is calculated by α^u=BuTγ^, where γ^=γ∧1T,⋯,γ∧qTT=θ∧p+1,⋯,θ∧p+L,θ∧p+L+1,⋯,θ∧p+qLT.


## 2.1. GPLVCM with Longitudinal Data

In this article, we consider a longitudinal study withn subjects and mi observations over time for the ith subject (i=1,⋯,n) for a total of N=∑i=1nmi observations. Each observation consists of a response variable Yij and the predicator variables Xij,Zij,Uij, where Xij∈Rp,Zij∈Rq and Uij is a scalar. We assume that the observations from different subjects are independent, but those within the same subject are dependent. The generalized varying coefficient partially linear model (GPLVM) with longitudinal data takes the form
(1)μij=EYij∣Xij,Zij,Uij=hXijTβ+ZijTαUij,where μij is the expectation of Yij when Xij, Zij, and Uij are given, β=β1,⋯,βpT is an unknown p×1 regression coefficient vector, h· is a known smooth link function, and αu=α1u,α2u,⋯,αquT is a q×1 unknown monotonic smooth function vector. Without loss of generality, we assume U∼U0,1.We approximateα· by B-spline basis functions Bu=B1u,⋯,BLuT with the order of M, where L=K+M+1 and K is the number of interior knots, i.e.,
(2)αku≈αk∗u=BuTγk,k=1,⋯,q,where γk=γk1,⋯,γkLT is a L×1 vector of unknown regression coefficients. Accordingly, μij is approximated by
(3)μij=EYij∣Xij,Zij,Uij=hXijTβ+ZijT·Iq⊗BUijTγ,where γ=γ1T,⋯,γqTT and “⊗” is the Kronecker product. We use the B-spline basis functions because they are numerically stable and have bounded support [18]. The spline approach also treats a nonparametric function as a linear function with the basis functions as pseudodesign variables, and thus, any computational algorithm for the generalized linear models can be used for the GPLVCMs.To incorporate the within-cluster correlation, we apply the QIFs to estimateβ and γ, respectively. Denote θ=βT,γTT, we define the extended score gNθ as follows:
(4)gNθ=1n∑i=1ngiθ=1n∑i=1nμ̇iTAi−12M1Ai−12Yi−μi⋮μ̇iTAi−12MsAi−12Yi−μi,where μ̇i=∂μi/∂θ, Ai=diagVarYi1,⋯,VarYim is the marginal variance matrix of subject Yi, and M1,⋯,Ms are the base matrices to represent the inverse of the working correlation matrix R in GEE approach. Following Qu et al. [10], we define the quadratic inference functions to be
(5)Qnθ=ngNTθΩnθ−1gNθ,where Ωnθ=1/n∑i=1ngiθgiθT. Note that Ωn depends on θ. The QIF estimate θ~ is then given by
(6)θ~=argminθQnθ.


## 2.2. Penalized QIF

In real data analysis, the true regression model is always unknown. An overfitted model lowers the efficiency of estimation while an underfitted one leads to a biased estimator. A popular approach to identify the relevant predictors while estimating the nonzero parameters and functions in model (1) simultaneously is to exert some kind of “penalty” on the original objective function. Here, we choose the smoothly clipped absolute deviation (SCAD) penalty because it has several advantages such as unbiasedness, sparsity, and continuity. The SCAD-penalized quadratic inference function (PQIF) is defined as follows:
(7)Qnpθ=Qnθ+n∑l=1ppλ1βl+n∑k=1qpλ2γkH,where γkH=γkTHγk1/2,H=hijL×L,hij=∫01BiuBjTudu and pλ is the SCAD penalty function, where the derivative is defined as
(8)p′lω=λIω⩽l+aλ−ω+a−1λIω>λ,where a>2,ω>0,pλ0=0; here, we choose a=3.7 as in [4].Note that(9)γkH=∫01γkTBuBTUγkdu1/2=∫01α∗u2du1/2.This group-wised penalization ensures that the spline coefficient vector of the same nonparametric component is treated as an entire group in model selection.Denoteθ^ to be the penalized estimator obtained by minimizing the penalized objective function of (7). Then, β^=θ∧1,⋯,θ∧pT is the estimator of the parameter β and the estimator of the nonparametric function αu is calculated by α^u=BuTγ^, where γ^=γ∧1T,⋯,γ∧qTT=θ∧p+1,⋯,θ∧p+L,θ∧p+L+1,⋯,θ∧p+qLT.


## 3. Asymptotic Properties

### 3.1. Oracle Property

We next establish the asymptotic properties of the resulting penalized QIF estimators. We first introduce some notations. Letβ0 and α0· denote the true values of β· and α·. In addition, γ0 is the spline coefficient vector from the spline approximation to α0·. Without loss of generality, we assume that β0l≠0,l=1,⋯,p1 and β0l=0,l=p1+1,⋯,p, i.e., only the first p1 component of β0 is nonzero. Similarly, we assume that α0k·≠0,k=1,⋯,q1 and α0k·=0,k=q1+1,⋯,q, i.e., only the first q1 component of α0· is nonzero. For convenience and simplicity, let C denote a positive constant that may have different values at each appearance throughout this paper and ∥A∥ denote the modulus of the largest singular value of matrix or vector A. Before the proof of our main theorems, we list some regularity conditions used in this paper.Assumption (A1).
The spline regression parameterγ is identifiable, that is, γ0 is the spline coefficient vector from the spline approximation to α0·. In addition, there is a unique θ0=β0,γ0∈S satisfying EgNθ0=0, where S is the parameter space.Assumption (A2).
The weight matrixΩn=1/n∑i=1ngiθgiTθ converges almost surely to a constant matrix Ω0, where Ω0 is invertible.Assumption (A3).
The covariate matricesXi and Zi, i=1,⋯,n, satisfy supiEXi4<∞ and supiEZi4<∞.Assumption (A4).
The errorεi=Yi−μi satisfies EεiεiT=Vi,supiVi<∞, and there exists a positive constant δ such that supiEεi2+δ<∞.Assumption (A5).
All marginal variancesAi≥0 and supiAi<∞.Assumption (A6).
mi is a bounded sequence of positive integers.Assumption (A7).
αiu,i=1,2,⋯,q is rth continuous differentiable on 0,1, where r≥2.Assumption (A8).
The inner knotsci,i=1,⋯,K satisfy
(10)max1≤i≤Khi+1−hi=oK−1,maxhiminhi≤C0,where hi=ci−ci−1.Assumption (A9).
The link functionh· is 2th continuous differentiable and Eh2+δ<∞ for some δ>2.Assumption (A10).
an=On−1/2;bn⟶0 as n⟶∞, where
(11)an=maxk,lp′λ1β0l,p′λ2γ0kH:βol≠0,γ0k≠0,bn=maxk,lp′′λ1β0l,p′′λ2γ0kH:βol≠0,γ0k≠0.Theorem1 indicates that the estimator of nonparametric components achieve the optimal convergence rate.Theorem 1.
Assume that Assumptions (A.1)–(A.10) hold and the number of knotsK=ON1/2r+1, then
(12)α^k·−α0k·=Opn−r/2r+1,k=1,⋯,q.
Furthermore, under suitable condition, Theorem1 shows that the penalized QIF estimator has the sparsity property.Theorem 2.
Assume that the conditions in Theorem1 hold and λmax⟶0,nλmin⟶∞ as n⟶∞, with probability approaching 1,
(13)β^l=0,l=p1+1,⋯,p,α^k·≡0,k=q1+1,⋯,q,where λmax=maxλ1,λ2,λmin=minλ1,λ2.
Theorems1 and 2 indicate that with the tune parameter λ being suitably chosen, the proposed selection method possesses model selection consistency. Next, we establish the asymptotic property for the estimator of the nonzero parametric components. Let β∗=β1,⋯,βp1T,α∗·=α1∗·,⋯,αq1∗·T and let β0∗ and α0∗· denote their true value, respectively. In addition, let γ∗=γ1T,⋯,γq1TT and γ0∗=γ01T,⋯,γ0q1TT denote the spline coefficient vector of α∗· and α0∗·, respectively, and let Xi∗ and Zi∗,i=1,⋯,n denote their correspondent covariate. Let X~i=H′ηiXi∗,X~=X~1T,⋯,X~nT,W~i=H′ηiWi∗,W~=W~1T,⋯,W~nT, and
(14)Γ=EX~TτX~−EX~TτW~uEW~TτW~u−1W~TτX~,Δ=Eτ−EX~TτW~uEW~TτW~u−1EW~Tτuε⊗2,where Δ⊗2=ΔΔT, τ=τijn×n is a n×n block matrix with its i,j block taking the form
(15)τij=∑k=1s∑l=1sAi−1/2MkAi−1/2H′ηiPi∗Ωlk−1Pj∗TH′ηjAj−1/2MlAj−1/2.
Theorem3 states that β∧∗ is asymptotically normally distributed.Theorem 3.
Suppose that Assumptions (A.1)–(A.9) hold and the number of knotsK=ON1/2r+1, then
(16)nβ∧∗−β0∗⟶LN0,Σ,where Σ=ΓΔ−1Γ−1 and ⟶L represents the convergence in distribution.


### 3.2. Selection of Tuning Parameters

Theorems1–3 imply that the proposed variable selection procedure possessed the oracle property. However, this attractive feature relies on the choice of tuning parameters λi. The popular criteria to choose λi include cross-validation, generalized cross-validation, AIC, and BIC. Wang et al. [19] suggested using BIC for the SCAD estimator in linear models and partially linear models and proved its model selection consistency property, i.e., the optimal parameter chosen by BIC can identify the true model with probability tending to one. Tian proved that for partially linear models. Hence, we adopt BIC to choose the optimal λ1,λ2. Following [19–21], we simplify the tuning parameters as
(17)λ1=λ0γ~k0H,λ2=λ0β~k0,where β~k0 and γ~k0 are the unpenalized QIF estimates. Consequently, the original two-dimensional problem becomes a univariate problem about λ0, which can be selected according to the following BIC-type criterion:
(18)BICλ=Qnpθ^λ+dfλ×logn,where θ^λ=β^λ,γ^1λT,⋯,γ^qλT is the regression coefficient estimated by minimizing the penalized QIF in (2.8) for a given λ and dfλ is the number of nonzero coefficients of β^λ and γ^1λH,⋯,γ^qλH. Thus, the tuning parameter λ is obtained by
(19)λ^=argminλBICλ.From Theorem 4 of Tian et al. [15], the BIC tuning parameter selector enables us to select the true model consistently.


### 3.3. An Algorithm Using Local Quadratic Approximation

Based on Fan and Li’s local quadratic approximating approach [4], we propose an iterative algorithm to minimize the PQIF (7). Similar with Tian et al. [15], we choose the unpenalized QIF estimator θ~ as the initial estimator. Let θk=β1k,⋯,βpk,γ1kT,⋯,γqkTT be the value of θ at the kth iteration. If βlk (or γlk) is close to 0 (or 0), i.e., βlk⩽ϵ (or γlkH⩽ϵ) with some small threshold value ϵ, then we set βlk=0 (or γlk=0). We consider ϵ=10−6 in our simulations.Supposeβlk+1=0, for l=pk+1,⋯,p, and γlk+1=0, for l=qk+1,⋯,q, and βk+1=β1k+1,⋯,βpkk+1,βpk+1k+1,⋯,βpk+1T=βNk+1T,βZk+1TT, where βNk+1=β1k+1,⋯,βpkk+1T are the nonzero parametric components and βZk+1=βpk+1k+1,⋯,βpk+1T=0. Similarly, let γk+1=γ1k+1T,⋯,γqkk+1T,γqk+1k+1T,⋯,γqk+1TT=γNk+1T,γZk+1TT, where γNk+1=γ1k+1T,⋯,γqkk+1TT and γZk+1=γqk+1k+1T,⋯,γqk+1TT correspond to qk zero functions and q−qK zero functions, respectively. Let θ=βNT,βZT,γNT,γZTT denote a vector which has the same length and same partition with θk+1.For the parametric term, ifβlk>ϵ, the penalty function at βl≈βlk is approximated by
(20)pλβl≈pλβlk+12pλ′βlkβlkβl2−βlk2.Similarly, to the nonparametric component, ifγlkH>ϵ, the penalty function at γl≈γlk is approximated by
(21)pλγlH≈pλγlkH+12pλ′γlkHγlkHγlH2−γlkH2=pλγlkH+12p′λγlkHγlkHβlTHβl−βlkTHβlk,where p′λ is the first-order derivative of the penalty function pλ. This leads to the local approximation of the PQIF Qnpθ by a quadratic function:
(22)Qnθk+Q̇nθkTω11−ω11k+12ω11−ω11kTQ¨nθkω11−ω11k+n2ω11tΛθkω11,where Q̇nθk=∂Qnθk/∂ω11, Q¨nθk=∂2Qnθk/∂ω11∂ω11T, with ω11=βNT,γZTT, and
(23)Λθk=diagp′λ2β1kβ1k,⋯,p′λ2βpkkβpkk,p′λ1γ1kHγ1kHH,⋯,p′λ1γqkkHγqkkHH.Minimizing the quadratic function (22), we obtain ω11k+1. The Newton-Raphson method then iterates the following process to convergence:
(24)ω11k+1=ω11k−Q¨nω11k+nΛω11k−1Q̇nω11k+nΛω11kω11k.


## 3.1. Oracle Property

We next establish the asymptotic properties of the resulting penalized QIF estimators. We first introduce some notations. Letβ0 and α0· denote the true values of β· and α·. In addition, γ0 is the spline coefficient vector from the spline approximation to α0·. Without loss of generality, we assume that β0l≠0,l=1,⋯,p1 and β0l=0,l=p1+1,⋯,p, i.e., only the first p1 component of β0 is nonzero. Similarly, we assume that α0k·≠0,k=1,⋯,q1 and α0k·=0,k=q1+1,⋯,q, i.e., only the first q1 component of α0· is nonzero. For convenience and simplicity, let C denote a positive constant that may have different values at each appearance throughout this paper and ∥A∥ denote the modulus of the largest singular value of matrix or vector A. Before the proof of our main theorems, we list some regularity conditions used in this paper.Assumption (A1).
The spline regression parameterγ is identifiable, that is, γ0 is the spline coefficient vector from the spline approximation to α0·. In addition, there is a unique θ0=β0,γ0∈S satisfying EgNθ0=0, where S is the parameter space.Assumption (A2).
The weight matrixΩn=1/n∑i=1ngiθgiTθ converges almost surely to a constant matrix Ω0, where Ω0 is invertible.Assumption (A3).
The covariate matricesXi and Zi, i=1,⋯,n, satisfy supiEXi4<∞ and supiEZi4<∞.Assumption (A4).
The errorεi=Yi−μi satisfies EεiεiT=Vi,supiVi<∞, and there exists a positive constant δ such that supiEεi2+δ<∞.Assumption (A5).
All marginal variancesAi≥0 and supiAi<∞.Assumption (A6).
mi is a bounded sequence of positive integers.Assumption (A7).
αiu,i=1,2,⋯,q is rth continuous differentiable on 0,1, where r≥2.Assumption (A8).
The inner knotsci,i=1,⋯,K satisfy
(10)max1≤i≤Khi+1−hi=oK−1,maxhiminhi≤C0,where hi=ci−ci−1.Assumption (A9).
The link functionh· is 2th continuous differentiable and Eh2+δ<∞ for some δ>2.Assumption (A10).
an=On−1/2;bn⟶0 as n⟶∞, where
(11)an=maxk,lp′λ1β0l,p′λ2γ0kH:βol≠0,γ0k≠0,bn=maxk,lp′′λ1β0l,p′′λ2γ0kH:βol≠0,γ0k≠0.Theorem1 indicates that the estimator of nonparametric components achieve the optimal convergence rate.Theorem 1.
Assume that Assumptions (A.1)–(A.10) hold and the number of knotsK=ON1/2r+1, then
(12)α^k·−α0k·=Opn−r/2r+1,k=1,⋯,q.
Furthermore, under suitable condition, Theorem1 shows that the penalized QIF estimator has the sparsity property.Theorem 2.
Assume that the conditions in Theorem1 hold and λmax⟶0,nλmin⟶∞ as n⟶∞, with probability approaching 1,
(13)β^l=0,l=p1+1,⋯,p,α^k·≡0,k=q1+1,⋯,q,where λmax=maxλ1,λ2,λmin=minλ1,λ2.
Theorems1 and 2 indicate that with the tune parameter λ being suitably chosen, the proposed selection method possesses model selection consistency. Next, we establish the asymptotic property for the estimator of the nonzero parametric components. Let β∗=β1,⋯,βp1T,α∗·=α1∗·,⋯,αq1∗·T and let β0∗ and α0∗· denote their true value, respectively. In addition, let γ∗=γ1T,⋯,γq1TT and γ0∗=γ01T,⋯,γ0q1TT denote the spline coefficient vector of α∗· and α0∗·, respectively, and let Xi∗ and Zi∗,i=1,⋯,n denote their correspondent covariate. Let X~i=H′ηiXi∗,X~=X~1T,⋯,X~nT,W~i=H′ηiWi∗,W~=W~1T,⋯,W~nT, and
(14)Γ=EX~TτX~−EX~TτW~uEW~TτW~u−1W~TτX~,Δ=Eτ−EX~TτW~uEW~TτW~u−1EW~Tτuε⊗2,where Δ⊗2=ΔΔT, τ=τijn×n is a n×n block matrix with its i,j block taking the form
(15)τij=∑k=1s∑l=1sAi−1/2MkAi−1/2H′ηiPi∗Ωlk−1Pj∗TH′ηjAj−1/2MlAj−1/2.
Theorem3 states that β∧∗ is asymptotically normally distributed.Theorem 3.
Suppose that Assumptions (A.1)–(A.9) hold and the number of knotsK=ON1/2r+1, then
(16)nβ∧∗−β0∗⟶LN0,Σ,where Σ=ΓΔ−1Γ−1 and ⟶L represents the convergence in distribution.


## 3.2. Selection of Tuning Parameters

Theorems1–3 imply that the proposed variable selection procedure possessed the oracle property. However, this attractive feature relies on the choice of tuning parameters λi. The popular criteria to choose λi include cross-validation, generalized cross-validation, AIC, and BIC. Wang et al. [19] suggested using BIC for the SCAD estimator in linear models and partially linear models and proved its model selection consistency property, i.e., the optimal parameter chosen by BIC can identify the true model with probability tending to one. Tian proved that for partially linear models. Hence, we adopt BIC to choose the optimal λ1,λ2. Following [19–21], we simplify the tuning parameters as
(17)λ1=λ0γ~k0H,λ2=λ0β~k0,where β~k0 and γ~k0 are the unpenalized QIF estimates. Consequently, the original two-dimensional problem becomes a univariate problem about λ0, which can be selected according to the following BIC-type criterion:
(18)BICλ=Qnpθ^λ+dfλ×logn,where θ^λ=β^λ,γ^1λT,⋯,γ^qλT is the regression coefficient estimated by minimizing the penalized QIF in (2.8) for a given λ and dfλ is the number of nonzero coefficients of β^λ and γ^1λH,⋯,γ^qλH. Thus, the tuning parameter λ is obtained by
(19)λ^=argminλBICλ.From Theorem 4 of Tian et al. [15], the BIC tuning parameter selector enables us to select the true model consistently.


## 3.3. An Algorithm Using Local Quadratic Approximation

Based on Fan and Li’s local quadratic approximating approach [4], we propose an iterative algorithm to minimize the PQIF (7). Similar with Tian et al. [15], we choose the unpenalized QIF estimator θ~ as the initial estimator. Let θk=β1k,⋯,βpk,γ1kT,⋯,γqkTT be the value of θ at the kth iteration. If βlk (or γlk) is close to 0 (or 0), i.e., βlk⩽ϵ (or γlkH⩽ϵ) with some small threshold value ϵ, then we set βlk=0 (or γlk=0). We consider ϵ=10−6 in our simulations.Supposeβlk+1=0, for l=pk+1,⋯,p, and γlk+1=0, for l=qk+1,⋯,q, and βk+1=β1k+1,⋯,βpkk+1,βpk+1k+1,⋯,βpk+1T=βNk+1T,βZk+1TT, where βNk+1=β1k+1,⋯,βpkk+1T are the nonzero parametric components and βZk+1=βpk+1k+1,⋯,βpk+1T=0. Similarly, let γk+1=γ1k+1T,⋯,γqkk+1T,γqk+1k+1T,⋯,γqk+1TT=γNk+1T,γZk+1TT, where γNk+1=γ1k+1T,⋯,γqkk+1TT and γZk+1=γqk+1k+1T,⋯,γqk+1TT correspond to qk zero functions and q−qK zero functions, respectively. Let θ=βNT,βZT,γNT,γZTT denote a vector which has the same length and same partition with θk+1.For the parametric term, ifβlk>ϵ, the penalty function at βl≈βlk is approximated by
(20)pλβl≈pλβlk+12pλ′βlkβlkβl2−βlk2.Similarly, to the nonparametric component, ifγlkH>ϵ, the penalty function at γl≈γlk is approximated by
(21)pλγlH≈pλγlkH+12pλ′γlkHγlkHγlH2−γlkH2=pλγlkH+12p′λγlkHγlkHβlTHβl−βlkTHβlk,where p′λ is the first-order derivative of the penalty function pλ. This leads to the local approximation of the PQIF Qnpθ by a quadratic function:
(22)Qnθk+Q̇nθkTω11−ω11k+12ω11−ω11kTQ¨nθkω11−ω11k+n2ω11tΛθkω11,where Q̇nθk=∂Qnθk/∂ω11, Q¨nθk=∂2Qnθk/∂ω11∂ω11T, with ω11=βNT,γZTT, and
(23)Λθk=diagp′λ2β1kβ1k,⋯,p′λ2βpkkβpkk,p′λ1γ1kHγ1kHH,⋯,p′λ1γqkkHγqkkHH.Minimizing the quadratic function (22), we obtain ω11k+1. The Newton-Raphson method then iterates the following process to convergence:
(24)ω11k+1=ω11k−Q¨nω11k+nΛω11k−1Q̇nω11k+nΛω11kω11k.


## 4. Simulation Studies

### 4.1. Assessing Rule

In this section, we conduct a simulation study to assess the finite sample performance of the proposed procedures. Following [17], the performance of estimator β^ will be assessed by the generalized mean square error (GMSE), which is defined as
(25)GMSE=1n∑i=1nβ^−βXi∗Xi∗Tβ^−β.The performance of estimatorα^· will be assessed by the square root of average square errors (RASE)
(26)RASE=1M∑v=1M∑k=1qα∧kuv−αkuv21/2,where uv,v=1,⋯,M are the grid points where the function α^u is evaluated. In our simulation, M=300 is used.To assess the performance of the variable selection, we use “C” to denote the average number of zero regression coefficients that are correctly estimated as zero and use “IC” to denote the average number of nonzero regression coefficients that are erroneously set to zero. The more closer the value of “C” to the number of true zero coefficient in the model and the more closer the value of “IC” to zero, the better the performance of the variable selection procedure is.In our simulations, we use the sample quantiles ofUij as knots and take the number of internal knots to be 3, that is, ON1/5. This particular choice is consistent with the asymptotic theory in Section 3 and performs well in the simulations. For each simulated dataset, the proposed estimation procedures for finding out penalized QIF estimators with SCAD and LASSO penalty functions are considered. The tuning parameters λ1,λ2 for the penalty functions are chosen by BIC from 50 equispaced grid points in −15,5. For each of these methods, the average of zero coefficients over the 500 simulated datasets is reported.


### 4.2. Study 1 (Partial Penalty)

Consider a Bernoulli response(27)logitYij=XijTβ+αUij,where β=2,1.5,0.7,017TT,m=6,Xij∼N0,I20,αUij=0.4cosπ/2Uij, and Uij are drawn independently from U0,1. Response variable Yij with compound symmetry correlation structure (CS) is generated according to Oman [22]. In our simulation study, we consider ρ=0.25 and 0.75, representing weak and strong correlations, respectively. In some situations, we prefer not to shrink some certain components in the variable selection procedure when some kind of prior information is available. Partial penalty arises naturally for such case. In this example, we only exert penalty on the parametric component, i.e., coefficient β. In this situation, the PQIF (7) becomes
(28)Qnpθ=Qnθ+n∑l=1ppλ1βl.The variable selection result is reported in Tables1 and 2.Table 1
Variable selection for the parametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.001115.8300.000616.24600.000516.7460LASSO0.000614.8100.000515.34600.000415.5740ρ=0.25SCAD0.001115.7500.000616.7000.000416.8460LASSO0.000714.8200.000614.9600.000515.350Table 2
RASE ofα^u under different methods.
Methodn=150n=200n=300ρ=0.75SCAD0.19200.20510.1054LASSO0.09990.08400.1064ρ=0.25SCAD0.24490.24600.0694LASSO0.13990.12050.1033Tables1 and 2 show that the performance of the proposed variable selection approach improves as n increases, e.g., the number of correctly recognized zero coefficient increases to the number of true zero coefficient in the model and the GMSE of β^ decreases as n increases. In addition, the RASE of α^u also decreases as n increases, which means the estimated curve of α^u fits better to the true line of αu when the sample size increases. Moreover, the SCAD penalty method outperforms the LASSO penalty ones in the sense of correct variable selection rate, which significantly reduces the model uncertainty and complexity.


### 4.3. Study 2 (Fixed-Dimensional Setup)

In this example, we generate data from the following model:(29)logitYij=1∣Xij,Uij=XijTβ+ZijTαUij,where β=2,1.5,0.7,07T and αu=α1u,α2u,05TT with α1u=0.8cosπ/2u,α2u=1.5+u2, Xij and Zijj=1,⋯,6 come from a multivariate normal distribution with mean zero, marginal variance 1 and correlation coefficient 0.5, and u∼U0,1. Response variable Yij with compound symmetry correlation structure (CS) is generated by the same method as study 1 and we also consider ρ=0.25 and 0.75, representing weak and strong correlations, respectively. We generated 500 datasets for each pair of N,ρ. The results are also reported in Tables 3 and 4.Table 3
Variable selection for the parametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.00486.7600.00366.84600.00306.8640LASSO0.00394.69400.00334.76600.00285.0740ρ=0.25SCAD0.00476.7600.00356.71800.00286.8460LASSO0.00384.81400.00354.9800.00295.0480Table 4
Variable selection for the nonparametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.16964.3500.12214.6600.08124.830LASSO0.19324.3800.15404.3600.12354.570ρ=0.25SCAD0.16364.4200.10764.7200.03444.850LASSO0.19824.4000.11604.6800.03984.760Table3 reports the variable selection for the parametric components; it shows that the performances become better and better as n increases, e.g., the number of correctly recognized zero coefficients, which is denoted as values in the column labeled “C,” becomes more and more closer to the true number of zero regression coefficients in the model. At the same time, the GMSE decreases steadily as n increases. Table 4 shows that, for the nonparametric components, the performances of the proposed variable selection method are similar to those of the method for the parametric components. As n increases, the RASE of the estimated nonparametric function also becomes smaller and smaller. This reflects that the estimate curves fit better to the corresponding true line as the sample size increases. Moreover, the SCAD penalty method outperforms the LASSO penalty ones in the sense of correct variable selection rate, which significantly reduces the model uncertainty and complexity.To study the influence of misspecified correlation structure to the proposed approach, we perform variable selection when the working correlation structure is specified to be CS and first-order autoregressive (AR-1), respectively. The result is listed in Table5. It is known that the QIF estimator is insensitive to misspecification in correlation structure. Table 5 shows that the proposed variable selection procedure gives similar results even when the correlation structure is misspecified. This indicates that our method is robust.Table 5
Variable selection when the trueR is CS when n=300.
WorkingRMethodβα·GMSECICRASECICρ=0.75CSSCAD0.00306.86400.08124.830LASSO0.00285.07400.12354.570AR-1SCAD0.00336.85600.09354.820LASSO0.00344.92400.12304.570ρ=0.25CSSCAD0.00286.84600.03444.850LASSO0.00295.04800.03984.760AR-1SCAD0.00306.84600.03544.860LASSO0.00315.04800.04114.750


### 4.4. Study 3 (High-Dimensional Setup)

In this example, we discuss how the proposed variable selection procedure can be applied to the “largen, diverging p/q” setup for longitudinal models. We consider the high-dimensional setup of study 2. In this simulation, we take n=300,m=6,p=20=ON1/4,q=10=ON1/4. The true coefficient vector is β=2,1.5,0.7,017TT,αu=α1u,α2u,010TT, where α1u and α2u are defined in study 2. The other settings are the same with study 2. The results are reported in Table 6. It is easy to see that the proposed variable selection procedure is able to correctly identify the true model and works well in the “large n, diverging p/q” setup.Table 6
Variable selection under high-dimensional setup.
Methodβα·GMSECICRASECICρ=0.75SCAD0.003616.66400.11489.6560LASSO0.003315.57400.12399.5460ρ=0.25SCAD0.003416.84600.10479.8750LASSO0.003915.3500.11389.8020


## 4.1. Assessing Rule

In this section, we conduct a simulation study to assess the finite sample performance of the proposed procedures. Following [17], the performance of estimator β^ will be assessed by the generalized mean square error (GMSE), which is defined as
(25)GMSE=1n∑i=1nβ^−βXi∗Xi∗Tβ^−β.The performance of estimatorα^· will be assessed by the square root of average square errors (RASE)
(26)RASE=1M∑v=1M∑k=1qα∧kuv−αkuv21/2,where uv,v=1,⋯,M are the grid points where the function α^u is evaluated. In our simulation, M=300 is used.To assess the performance of the variable selection, we use “C” to denote the average number of zero regression coefficients that are correctly estimated as zero and use “IC” to denote the average number of nonzero regression coefficients that are erroneously set to zero. The more closer the value of “C” to the number of true zero coefficient in the model and the more closer the value of “IC” to zero, the better the performance of the variable selection procedure is.In our simulations, we use the sample quantiles ofUij as knots and take the number of internal knots to be 3, that is, ON1/5. This particular choice is consistent with the asymptotic theory in Section 3 and performs well in the simulations. For each simulated dataset, the proposed estimation procedures for finding out penalized QIF estimators with SCAD and LASSO penalty functions are considered. The tuning parameters λ1,λ2 for the penalty functions are chosen by BIC from 50 equispaced grid points in −15,5. For each of these methods, the average of zero coefficients over the 500 simulated datasets is reported.


## 4.2. Study 1 (Partial Penalty)

Consider a Bernoulli response(27)logitYij=XijTβ+αUij,where β=2,1.5,0.7,017TT,m=6,Xij∼N0,I20,αUij=0.4cosπ/2Uij, and Uij are drawn independently from U0,1. Response variable Yij with compound symmetry correlation structure (CS) is generated according to Oman [22]. In our simulation study, we consider ρ=0.25 and 0.75, representing weak and strong correlations, respectively. In some situations, we prefer not to shrink some certain components in the variable selection procedure when some kind of prior information is available. Partial penalty arises naturally for such case. In this example, we only exert penalty on the parametric component, i.e., coefficient β. In this situation, the PQIF (7) becomes
(28)Qnpθ=Qnθ+n∑l=1ppλ1βl.The variable selection result is reported in Tables1 and 2.Table 1
Variable selection for the parametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.001115.8300.000616.24600.000516.7460LASSO0.000614.8100.000515.34600.000415.5740ρ=0.25SCAD0.001115.7500.000616.7000.000416.8460LASSO0.000714.8200.000614.9600.000515.350Table 2
RASE ofα^u under different methods.
Methodn=150n=200n=300ρ=0.75SCAD0.19200.20510.1054LASSO0.09990.08400.1064ρ=0.25SCAD0.24490.24600.0694LASSO0.13990.12050.1033Tables1 and 2 show that the performance of the proposed variable selection approach improves as n increases, e.g., the number of correctly recognized zero coefficient increases to the number of true zero coefficient in the model and the GMSE of β^ decreases as n increases. In addition, the RASE of α^u also decreases as n increases, which means the estimated curve of α^u fits better to the true line of αu when the sample size increases. Moreover, the SCAD penalty method outperforms the LASSO penalty ones in the sense of correct variable selection rate, which significantly reduces the model uncertainty and complexity.


## 4.3. Study 2 (Fixed-Dimensional Setup)

In this example, we generate data from the following model:(29)logitYij=1∣Xij,Uij=XijTβ+ZijTαUij,where β=2,1.5,0.7,07T and αu=α1u,α2u,05TT with α1u=0.8cosπ/2u,α2u=1.5+u2, Xij and Zijj=1,⋯,6 come from a multivariate normal distribution with mean zero, marginal variance 1 and correlation coefficient 0.5, and u∼U0,1. Response variable Yij with compound symmetry correlation structure (CS) is generated by the same method as study 1 and we also consider ρ=0.25 and 0.75, representing weak and strong correlations, respectively. We generated 500 datasets for each pair of N,ρ. The results are also reported in Tables 3 and 4.Table 3
Variable selection for the parametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.00486.7600.00366.84600.00306.8640LASSO0.00394.69400.00334.76600.00285.0740ρ=0.25SCAD0.00476.7600.00356.71800.00286.8460LASSO0.00384.81400.00354.9800.00295.0480Table 4
Variable selection for the nonparametric components under different methods.
Methodn=150n=200n=300GMSECICGMSECICGMSECICρ=0.75SCAD0.16964.3500.12214.6600.08124.830LASSO0.19324.3800.15404.3600.12354.570ρ=0.25SCAD0.16364.4200.10764.7200.03444.850LASSO0.19824.4000.11604.6800.03984.760Table3 reports the variable selection for the parametric components; it shows that the performances become better and better as n increases, e.g., the number of correctly recognized zero coefficients, which is denoted as values in the column labeled “C,” becomes more and more closer to the true number of zero regression coefficients in the model. At the same time, the GMSE decreases steadily as n increases. Table 4 shows that, for the nonparametric components, the performances of the proposed variable selection method are similar to those of the method for the parametric components. As n increases, the RASE of the estimated nonparametric function also becomes smaller and smaller. This reflects that the estimate curves fit better to the corresponding true line as the sample size increases. Moreover, the SCAD penalty method outperforms the LASSO penalty ones in the sense of correct variable selection rate, which significantly reduces the model uncertainty and complexity.To study the influence of misspecified correlation structure to the proposed approach, we perform variable selection when the working correlation structure is specified to be CS and first-order autoregressive (AR-1), respectively. The result is listed in Table5. It is known that the QIF estimator is insensitive to misspecification in correlation structure. Table 5 shows that the proposed variable selection procedure gives similar results even when the correlation structure is misspecified. This indicates that our method is robust.Table 5
Variable selection when the trueR is CS when n=300.
WorkingRMethodβα·GMSECICRASECICρ=0.75CSSCAD0.00306.86400.08124.830LASSO0.00285.07400.12354.570AR-1SCAD0.00336.85600.09354.820LASSO0.00344.92400.12304.570ρ=0.25CSSCAD0.00286.84600.03444.850LASSO0.00295.04800.03984.760AR-1SCAD0.00306.84600.03544.860LASSO0.00315.04800.04114.750


## 4.4. Study 3 (High-Dimensional Setup)

In this example, we discuss how the proposed variable selection procedure can be applied to the “largen, diverging p/q” setup for longitudinal models. We consider the high-dimensional setup of study 2. In this simulation, we take n=300,m=6,p=20=ON1/4,q=10=ON1/4. The true coefficient vector is β=2,1.5,0.7,017TT,αu=α1u,α2u,010TT, where α1u and α2u are defined in study 2. The other settings are the same with study 2. The results are reported in Table 6. It is easy to see that the proposed variable selection procedure is able to correctly identify the true model and works well in the “large n, diverging p/q” setup.Table 6
Variable selection under high-dimensional setup.
Methodβα·GMSECICRASECICρ=0.75SCAD0.003616.66400.11489.6560LASSO0.003315.57400.12399.5460ρ=0.25SCAD0.003416.84600.10479.8750LASSO0.003915.3500.11389.8020


## 5. Application to Infectious Disease Data

We apply the proposed method to analyze an infectious disease data (indon.data), which has been well analyzed by many authors, such as [16, 23–27]. In this study, a total of 275 preschool children were examined every three months for 18 months. The response is the presence of respiratory infection (1=yes, 0=no). The primary interest is in studying the relationship between the risk of respiratory infection and vitamin A deficiency (1=yes, 0=no).In our study, we consider the following GPLVCM model(30)logitμij∣Xij,tij=∑k=16βixij+α0tij+zijα1tij,where t is age, X1 is vitamin A deficiency, X2,X3 are the seasonal cosine and seasonal sine variables, respectively, which indicate the season when those examinations took place, X4 is gender (1=female, 0=male), X5 is height, X6 is stunting status (1=yes, 0=no), and Z1=X52 is the square of height. The with-cluster correlation structure is assumed to be exchangeable, i.e., compound symmetric. This structure is also used in [16, 26, 27].We apply the proposed QIF-based group SCAD variable selection procedure to the above model and recognize five nonzero coefficients and one nonzero functionα0t, where β1=0.842, β2=−0.685, β3=−0.309, β4=−0.554, and β6=0.966. The results are generally consistent with those previous studies, but our results show that the height has no significant impact on the infectious rate and can be removed from the model. Figure 1 reports the curve of baseline age function α0t estimated by QIF-based group SCAD that is estimated by QIF and that is estimated by QIF-based SCAD partial penalty to β in [16], where the GPLM without the varying coefficient term is used. Figure 1 implies that the probability of having respiratory infection increases at the very early stage, then decreases steadily, and declines dramatically when the age is over 5.5 years old. This also coincides with previous results [16, 26, 27].Figure 1
The estimated function on age for the infectious disease data.


## 6. Conclusion and Discussion

We proposed a QIF-based group SCAD variable selection procedure for the generalized partially linear varying coefficient models with longitudinal data. This procedure can select significant variables in the parametric components and nonparametric components simultaneously. Under mild conditions, the estimators of regression coefficients have oracle property. Simulation studies indicate that the proposed procedure is very effective in selecting significant variables and estimating the regression coefficients.In this paper, we assume that the dimensions of the covariatesX and Z are fixed. Study 3 in simulations shows that the proposed approach still have desired results when the dimensions p and q go to infinity as n⟶∞. However, when in ultrahigh-dimensional case, the proposed variable selection procedure may not work well anymore. As a future research topic, it is interesting to consider the variable selection for the generalized partially linear varying coefficient models with ultrahigh-dimensional covariates.



---

*Source: 3505306-2020-10-05.xml*