Model Detection for Additive Models with Longitudinal Data

In this paper, we consider the problem of variable selection and model detection in additive models with longitudinal data. Our approach is based on spline approximation for the components aided by two Smoothly Clipped Absolute Deviation (SCAD) penalty terms. It can perform model selection (finding both zero and linear components) and estimation simultaneously. With appropriate selection of the tuning parameters, we show that the proposed procedure is consistent in both variable selection and linear components selection. Besides, being theoretically justified, the proposed method is easy to understand and straightforward to implement. Extensive simulation studies as well as a real dataset are used to illustrate the performances.


Introduction
Longitudinal data arise frequently in biological and economic applications.The challenge in analyzing longitudinal data is that the likelihood function is difficult to specify or formulate for non-normal responses with large cluster size.To allow richer and more flexible model structures, an effective semi-parametric regression tool is the additive model introduced by [1], which stipulates that where Y is a varaible of interest and which is required for identifiability of model (1.1), 1, , l d =  .We propose a penalized method for variable selection and model detection in model (1.1) and show that the proposed method can correctly select the nonzero components with probability approaching one as the sample size goes to infinity.
Statistical inference of additive models with longitudinal data has also been considered by some authors.By extending the generalized estimating equations approach, [2] studied the estimation of additive model with longitudinal data.[3] focuses on a nonparametric additive time-varying regression model for longitudinal data.[4] considered the generalized additive model when responses from the same cluster are correlated.However, in semiparametric regression modeling, it is generally difficult to determine which covariates should enter as nonparametric components and which should enter as linear components.The commonly adopted strategy in practice is just to consider continuous entering as nonparametric components and discrete covariates entering as parametric.Traditional method uses hypothesis testing to identify the linear and zero component.But this might be cumbersome to perform in practice whether there are more than just a few predictor to test.[5] proposed a penalized procedure via the LASSO penalty; [6] presented a unified variable selection method via the adaptive LASSO.But these methods are for the varying coefficient models.[7] established a model selection and semiparametric estimation method for additive quantile regression models by two-fold penalty.To our knowledge, the model selection and variable selection simultaneously with longitudinal data have not been investigated.We make several novel contributions: 1) We develop a new strategies for model selection and variable selection in additive model with longitudinal data; 2) We develop theoretical properties for our procedure.
In the next section, we will propose the two-fold SCAD penalization procedure based on QIF and computational algorithm; furthermore we present its theoretical properties.In particular, we show that the procedure can select the true model with probability approaching one, and show that newly proposed method estimates the non-zero function components in the model with the same optimal mean square convergence rate as the standard spline estimators.Simulation studies and an application of proposed methods in a real data example are included in Sections 3 and 4, respectively.Technical lemmas and proofs are given in Appendix.

Additive Models with Two Fold Penalized Splines
Consider a longitudinal study with n subjects and i n observations over time for the ith subject ( ) observation.Each observation consists of a response variable ij Y and a covariate vector taken from the ith subject at time ij t .We assume that the full data set ,  is observed and can be modelled as where ij ε is random error with ( ) At the start of the analysis, we do not know which component functions in model (1.1) are linear or actually zero.We adopt the centered B-spline basis, where . Equally-spaced knots are used in this article for simplicity of proof.Other regular knot sequences can also be used, with similar asymptotic results.Suppose that ( ) l m ⋅ can be approximated well by a spline function, so that .
To simplify notation, we first assume equal cluster size i n m = < ∞ , and let ( ) , , , d Jd β β β β + =  be the collection of the coefficients in (2.3), and , , 1, , , , then we have an approximation We can also write the approximation of (2.1) in matrix notation as where , , , , , , , , . [8] introduced the QIF that approximates the inverse of R by a linear combination of some basis matrixes, i.e.
where I is the identity and i M are known symmetric basis matrices and 0 1 , , , K a a a  are unknown constants.The advantage of the QIF approach is that it does not require the estimation of the linear coefficients i a 's associated with the working correlation matrix, which are treated as nuisance parameters here.
The vector

( )
n β G contains more estimating equations than parameters, but these estimating equations can be combined optimally using the generalised method of the moment.So according to [8], the QIF approach estimates β by setting n G as close to zero as possible, in the sense of minimizing the quadratic inference function ( ) where ( ) ( ) ( )

C g g
Our main goal is to find both zero components (i.e., 0 j m ≡ ) and linear compoents (i.e., j m is a linear function).The former can be achieved by shrinking sp j m to zero.For the latter, we want to shrink the second derivative ( ) sp '' j m to zero instead.This suggests the following minimization problem where ( ) p λ ⋅ are two penalties used to find zero and linear coefficients respectively, with two regularization parameters 1 λ and 2 λ , and

Asymptotic Properties
To study the rate of convergence for μ and β , we first introduce some notations and present regularity conditions.We assume equal cluster sizes ( ) For convenience, we assume that ( ) . The asymptotic result still hold for data with unequal cluster sizes i m using a cluster-specific transformation as discuss in [4].For any matrix A , A denotes the modulus of the largest singular value of A .To prove the theoretical arguments, we need the following assumptions:  are compactly supported, and without loss of generality, we assume that each x , is absolutely conti- nuous and there exist constants 1 C and 2 C such that ( ) . Assume the modular of the singular value of M is bounded away from 0 and infinity.(A6) The matrix A defined in Theorem 3 is positive definite.
Theorem 1. Suppose that the regularity conditions A1-A5 hold and the number of knots ( ) ( ) , 0 λ λ → .Then there exists a local minimizer of (2.7) such that ( ) ( ) , it reduces to a special case where the responses are i.i.d.The rate of convergence given here is the same optimal rate as that obtain for polynomial spline regression for independent data [9] [10].The main advantage of the QIF approach is that it incorporates within-cluster correlation by optimally combing estimating equations without estimating the correlation parameters.the estimator of two fold penalized QIF achieve the same rate of convergence as un-penalized estimator.Furthermore, we prove that the penalized estimators in Theorem 1 possess the sparsity property, ˆ0 l m = almost surely for 1, , l s d = +  .The sparsity property ensures that the proposed model selection is consistent, that is, it selects the correct variables with probability tending to 1 as the sample size goes to infinity.Theorem 2. Under the same assumptions of Theorem 1, and if the tuning parameter Then with probability approaching 1. a) ˆ0, 1 m is a linear function for 1 1 d j s + ≤ ≤ Theorem 2 also implies that above additive model selection possesses the consistency property.The results in Theorems 2 are similar to semiparametric estimation of additive quantile regression model in [7].However, the theoretical proof is very different from the penalized quantile loss function due to the two fold penalty and longitudinal data.
Finally, in the same spirit of that [11], we come to the question of whether the SIC can identify the true model in our setting.Theorem 3. Suppose that the regularity conditions A1-A5 hold and the number of knots

Simulation Study
In this section, we conducted Monte Carlo studies for the following longitudinal data and additive model.the   ( )  are generated independently from uniform.The error ( ) , , tion with mean 0, a common marginal variance 2 1 σ = , it has first-order autoregressive (AR-1) and an com- The predictors  , where Φ is the standard normal c.d.f. and To illustrate the effect on estimation efficiency, we compare the penalized QIF approach in [4] (PQIF) and an Oracle model (ORACLE).here the full model consists of all ten variable, and oracle model only contains the first five relevant variables and we know it's a partial additive model.The oracle model is only available in simulation studies where the true information is known.In all simulation, the number of replications is 100 and the result are summarized in Table 1 and Table 2 with the one penalty QIF when the error are Gaussian, and we also list the oracle model as a benchmark, the oracle model is only available in simulation studies where the true information is known in Table 1, in which the column labeled "NNC" presents the average number of nonparametric components selected, the column "NNT" depicts the average number of nonparametric components selected that are truly nonparametric (truly nonzero for one penalty QIF), "NLC" presents the average number of linear components, "NLT" depicts the average number of linear components selected that are truly linear.
In Table 2, we conduct some simulations to evaluate finite sample performance of the proposed method.Let ˆj m ⋅ will be assessed by using the square root of average square errors(RASE), we compare the performance of above estimators.On the nonparametric coponents, the errors for estimators with a single penalty and our procedure are similar, and both are qualitatively close to those of the oracle estimator.However, for the parametric components, our estimator is obviously more efficient,leading to about 40% -50% reduction in RASE.

Real Data Analysis
In this subsection, we analyze data from the Multi-Center AIDS Cohort Study.The dataset contains the human immunodeficiency virus, HIV, status of 283 homosexual men who were infected with HIV during the follow-up period between 1984 and 1991.All individuals were scheduled to have their measurements made during semiannual visits.Here , 1, , , 1, , denotes the time length in years between seroconversion and the j-th measurement of the i-th individual after the infection.[12] analyzed the dataset using partial linear models.The primary interest was to describe the trend of the mean CD4 percentage depletion over time and to evaluate the effects of cigarette smoking, pre-HIV infection CD4 percentage, and age at infection on the mean CD4 cell percentage after the infection.
In our analysis, the response variable is the CD4 cell percentage of a subject at distinct time points after HIV infection.We take four covariates for this study: 1 X , the CD4 cell percentage level before HIV infection; and 2 X , age at HIV infection; 3 X the individual's smoking status, which takes binary values 1 or 0, according to whether a individual is a smoker or nonsmoker; ij T denote 1, , , 1, , i n j m = =   , denotes the time length in years between seroconversion and the j -th measurement of the i -th individual after the infection.We construct the following additive model; the partially linear additive models instead of additive model because of the binaray variable X , but we not select the linear component.using our procedure, we wang to ensure which is linear component and which is zero in the non-parametirc function.For implement our procedure, linear transformation be used to the variable , , X X T .The result of our procedure select the 1 m is zero function and select the 2 m is a linear function, depletes rather quickly at the beginning of HIV infection, but the rate of depletion appears to be slowing down at four years after the infection.This result is the same as before [13].

Concluding Remark
In summary, we present a two-fold penalty variable selection procedure in this paper, which can select linear component and significant covariate and estimate unknown coefficient function simultaneously.The simulation study shows that the proposed model selection method is consistent with both variable selection and linear components selection.Besides, being theoretically justified, the proposed method is easy to understand and straightforward to implement.Further study of the problem is how to use the multi-fold penalty to solve the model selection and variable selection in generalized additive partial linear models with longitudinal data.

Appendix: Proofs of Theorems
For convenience and simplicity, let C denote a positive constant that may be different at each appearance throughout this paper.Before we prove our main theorems, we list some regularity conditions that are used in this paper.Lemma 1.Under the conditions (A1)-(A6), minimizing the no penalty QIF ( ) satisfying the condition (4).
There exists a constant 0 C > and a spline function . Using the triangular in equality According to [8] entail , we want to show that for large n and any 0 ε > , there exist a constant C large enough such that ( ) ( ) ( ) ( ) As a result, this implies that , 1 o λ λ = .We have By the definition of SCAD penalty function, removing the regularizing terms in (A2) with n Q  and n Q  being the gradient vector and hessian matrix n Q , respectively.Following [8], and Lemma A1 in supplement, for any β with ( ) ( ) where n G  is the first order derivative of n G .Therefore, by choosing C large enough, the second term on (A3) dominates its first term.therefore (A1) holds when C and n are large enough.This completes the proof of Theorem 1.  Proof of Theorem 2. We only show part (b), as an illustration and part (a) is similar.Suppose for some ˆˆˆˆ.
As in the proof of Theorem 1, we have and thus with probability approaching   , we denote the estimator of two fold penalty ˆλ β , and denote by β the minimizer when the optimal sequence of regularization parameters is chosen.There are four separate cases to consider CASE 1: j j B λ β represents a linear component for some   β is nonzero for j s ≥ .The case is similar to case 3. Thus the proof is omitted.
of predictor variables, µ is a unknown constant, and functions.As in most work on nonparametric smoothing, estimation of the non-parametric functions SIC can select the true model with pro- bability tending to 1.
− .Thus the last 5 variables in this model are null variables and do not contribute the model.The covariates ( ) ( ) pound symmetry (CS) correlation (i.e.exchangeable correlation) structure with different within correlation coefficient, and weak within correlation structure.

Figure 1 .
Figure 1.The estimator of 1. by the definition of SCAD penalty.

Proof of Theorem 3 .
to the proof of Theorem 1, by choosing C large enough, the second term on the right had side of (A4) dominates its first term. For any regularization parameters

.
In Table1, the model selection result for both our procedure

Table 1 .
The estimation results for our estimator (TFPQIF) and sparse additive estimator (PQIF) and ORACLE esitmator.
the SIC cannot select such λ .CASE 2: ˆj λ β is zero for some 1 j s ≤ ≤ .The proof is very similar with CASE 1 and therefore omitted..Here when considering CASE 3, we implicity exclude all previous cases that no underfitting cases.β~ is the estimator of minimizing the no penalty QIF(2.6)