{'fraction_non_alphanumeric': 0.03676556957267367, 'fraction_numerical': 0.015340522000610802, 'mean_word_length': 4.210919329171257, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularization in Generalized Semiparametric Mixed-Effects Model for Longitudinal Data

This paper considers the problem of simultaneous variable selection and estimation in the generalized semiparametric mixed-effects model for longitudinal data when the number of parameters diverges with the sample size. A penalization type of generalized estimating equation method is proposed while using regression spline to approximate the nonparametric component. Our approach applies SCAD to the estimating equation objective function in order to simultaneously estimate parameters and select the important variables. The proposed procedure involves the specification of the posterior distribution of the random effects, which cannot be evaluated in a closed form. However, it is possible to approximate this posterior distribution by producing random draws from the distribution using a Metropolis algorithm, which does not require the specification of the posterior distribution. For practical implementation, we develop an appropriate iterative algorithm to select the significant variables and estimate the nonzero coefficient functions. Under some regularity conditions, the resulting estimators enjoy the oracle properties, under a high-dimensional regime. Simulation studies are carried out to assess the performance of our proposed method, and a real data set is analyzed to illustrate the procedure.


Introduction
Longitudinal studies are often conducted in epidemiology, social science and other biomedical research areas. A challenge in the analysis of longitudinal data is that the repeated measurements from the same subjects are correlated over time. A popular way for incorporating this correlation within the likelihood framework is to use the linear mixed effects model (LMM; Laird and Ware (1982)) to analyze continuous longitudinal data and the generalized linear mixed effects model (GLMM; (Zeger and Karim , 1991)) to analyse discrete longitudinal data, where the random component takes care the correlation among observations from the same subjects. However, the traditional GLMM assumes parametric fixed effects that may be too restrictive to account complex covariate effects, especially when the variety of response over time is in a complicated manner. To eliminate the limitation of the GLMMs for modeling non linear time trend, a generalized semiparametric mixed effects model (GSMM), a natural extension of the GLMMs and semiparametric mixed models (SMM; (Zeger and Diggle , 1994)), is widely used to analyze longitudinal data by incorporating the within subject correlation using random effects and an arbitrary smooth function to model the time effect. Further developments along this line in the framework of GSMM can be found in Fan et al. (2007), Qin and Zhu (2007), Zhu (2009), Liang (2009), and Kurum et al. (2016) to mention afew.
High dimensional longitudinal data, which consist of repeated measurements on a large number of covariates, have become increasingly common. Inclusion of redundant variables may hinder accuracy and efficiency for both estimation and inference.
Thus, it is important to develop new statistical methodology and theory of variable selection and estimation for high-dimensional longitudinal data. There is a large body of variable selection methods for cross-sectional data. Among all, we refer to bridge regression (Frank and Friedman , 1993), Lasso (Tibshirani , 1996), adaptive Lasso (Zou , 2006), Elastic-net (Zou and Hastie , 2005), and SCAD (Fan and Li , 2001). The literature on variable selection for longitudinal data is rather limited due to the challenges imposed by incorporating the intracluster correlation. Fan and Li (2004) extended the SCAD procedure to the semiparametric model for longitudinal data. Bondell et al. (2010) proposed simultaneous selection of the fixed and random factors using a penalized joint log likelihood for the LMM. Ni et al. (2010) proposed a double-penalized likelihood approach for simultaneous model selection and estimation for the SMM. Ma et al. (2013) applied proper penalty functions in the additive semiparametric model. Chu et al. (2016) developed a screening procedure for ultrahigh dimensional longitudinal data. In contrast to extensive attention on model selection for Gaussian longitudinal data, research on model selection for non-Gaussian longitudinal data in the framework of the GLM remains largely unexplored.
To do variable selection, Pan (2001) developed a quasi-likelihood information criterion (QIC) which is analogous to AIC; Cantoni et al. (2005) generalized Mallow's C p criterion, and Wang and Qu (2009) proposed a BIC criterion based on the quadratic inference function. These are best subset type model selection procedures which become computationally intensive when the number of parameters is moderately large.
Regarding regularization methods for longitudinal data, Fu (2003) proposed a generalization of the bridge and Lasso penalties to the generalized estimating equations (GEE) model. Xu and Zhu (2010) extended the independence screening method to deal with the high dimensional longitudinal GLMs. Dziak (2006) generalized the Lasso and SCAD methods to the longitudinal GLMs. The SCAD-penalized selection procedures were illustrated in Xue et al (2010) for the generalized additive model with correlated data. In all aforementioned studies on the penalized GLM for longitudinal data, the dimension of predictors is fixed. Xu et al. (2012) proposed a weighted least squares type function to study the longitudinal GLMs with a diverg-ing number of parameters. For correlated discrete outcome data, the joint likelihood function does not have a closed form if the correlation information is taken into account. When the dimension of parameters diverges, numerical approximation to the joint likelihood function tends to be computationally infeasible as it often involves high-dimensional integration. This motivated Liang and Zeger (1986) to develop an approach of the GEE which is a multivariate analogue of the quasi-likelihood. Johnson et al. (2008) recently derived the asymptotic theory for the penalized estimating equations for independent data. Wang et al. (2012) employed rather different techniques than those in Johnson et al. (2008) and proposed the SCAD-penalized GEE for analyzing longitudinal data with high dimensional covariates. To the best of our knowledge, regularization in the GSMM is neglected.
In this paper, we focus on the GSMM with longitudinal data by allowing for non-Gaussian data and nonlinear link functions. We consider the case where the number of variables p is allowed to increase with the number of sample size n under a high dimensional regime. Similar to the work of Wang et al. (2012), we apply the penalty functions to the estimating equation objective function such that the proposed procedure can simultaneously estimate parameters and select the important variables.
Our method is rather different from their work because of including random effects and a nonparametric component in the model. We adopt spline regression to estimate the nonparametric components. The proposed penalized estimation involves the specification of the posterior distribution of the random effects, which cannot be evaluated in a closed form. However, it is possible to approximate this posterior distribution by producing random draws from a distribution using the Metropolis algorithm (Tanner , 1993), which does not require the specification of the posterior distribution. We establish the asymptotic theory for the proposed method in a highdimensional framework where the number of covariates increases with the sample size. To estimate the parameters, a computationally flexible iterative algorithm is developed. Furthermore, we propose a sandwich formula to estimate the asymptotic covariance matrix.
The rest of this paper is organized as follows. Section 2, formulates the model, discusses the approximation of the nonparametric function using splines and considers the estimation under the GEE framework. Section 3 include selection of the regularization parameters and the model selection procedure for assessing the fits of candidate models is also addressed. Moreover, asymptotic properties of the estimators and a Monte Carlo Newton-Raphson algorithm to implement the procedures is given. In Section 4, we apply a number of simulations to assess the finite sample performance of the proposed estimation method in the GSMM. A real data analysis is also presented in this section to augment the theoretical results. Some concluding remarks are given in Section 5.
Further, the proofs of the main results as well as some instrumental lemmas are provided in a separate supplementary file.

Model specification
Consider a longitudinal study with n subjects and n i observations over time for the ith subject (i " 1, . . . , n). Let u i be a qˆ1 vector of random effects corresponding to the ith subject, and y ij be an observation of the ith subject measured at time t ij for i " 1, . . . , n and j " 1, . . . , n i . Suppose that y i1 , . . . , y in i given u i are conditionally independent and each y ij |u i is distributed as an exponential family distribution whose probability density function is given by where φ is a scale parameter, cp., .q is a function only depending on y ij and φ, and θ ij is the (scalar) canonical parameter. The conditional expectations and variances of y ij given u i are given by µ ij " Epy ij |u i q " b . pθ ij q and ν ij " varpy ij |u i q " φb .. pθ ij q, respectively, where b . pθq " Bbpθq Bθ and b .. pθq " B 2 bpθq Bθ 2 . In this paper, we assume that the conditional mean µ ij satisfies where gp.q is a known monotonic link function, X J ij is a p nˆ1 vector of explanatory variables, β n is a p nˆ1 vector of unknown parameters of the fixed effects, Z J ij is a qˆ1 vector of explanatory variables relating to the random effects, f p.q is an unknown smooth function which is continuous and twice differentiable function on some finite interval. The dimension of the covariates p n is allowed to depend on the number of subjects n. To complete the specification, assume that the random effects u " tu 1 , . . . , u q u independently follow the same distribution, depending on parameters Σ as The model defined in Eqs.
(1)-(3) is referred to as generalized semiparametric mixed model (GSMM). Specific assumptions will be considered for the number of variables p n in section 3.3.

Smoothing spline approximation
Following the most literature such as He et al. (2005) and Qin and Zhu (2009), the unspecified smooth function can be approximated sufficiently well by the following polynomial spline where d is the degree of the polynomial component, L n is the number of interior knots (rate of L n will be specified in Section 3.3), t plq i is referred as knots of the ith subject, Bpt ij q "´1, t ij , . . . , t d ij ,`t ij´t p1q i˘d , . . . ,`t ij´t pLnq i˘d¯i s a h nˆ1 vector of basis functions, h n is the number of basis functions used to approximate f pt ij q, h n " d`1`L n , paq`" maxp0, aq, and α n " pα 0 , . . . , α d , α d`1 , . . . , α d`1`Ln q J is the spline coefficients vector of dimension h. Regression splines have some desirable properties in approximating a smooth function. It often provides good approximations with a small number of knots. The spline approach also treats a nonparametric function as a linear function with the basis functions as pseudo design variables, and thus this linearizes our regression model (2) so that our regression problem becomes gpµ ij q fi η ij " X J ij β n`Z J ij u i`B pt ij qα n , i " 1, . . . , n; j " 1, . . . , n i .
For convenience, model (4) can take the form where D ij "`X J ij , B j pt i q J˘J being a pp n`hn qˆ1 design matrix combining the fixed-effects and spline-effects design matrices for the jth outcome of the ith subject, and θ n " pβ J n , α J n q J is a pp n`hn qˆ1 combined regression parameters vector that must be estimated.
According to the linearization of the GSMM (2) using spline approach, any computational algorithm developed for the GMM can be used for the GSMM. We now formulate the linearization of the GSMM in the seamless form

Estimation procedure
For linearization of the GSMM defined in (5), the classical likelihood function can be defined as where y i " py i1 , . . . , y in i q J , u " pu 1 , . . . , u n q, and p y i |u i py i |u i , θ n , φq " n i ź j"1 ppy ij |u i , θ n , φq.
For the maximum likelihood (ML) estimates of the parameters θ n , φ and Σ, one can maximize this likelihood function by using suitable numerical techniques. The EM algorithm is an attractive method to obtain the ML estimates, in presence of incomplete data, which avoids explicit calculation of the observed data log-likelihood.
To set up the EM algorithm we consider the random effects, u i , to be the missing data. The complete data, is then (y i , u i ) and the complete data log-likelihood is given by Considering u i to be the missing has an advantage that in the M-step, maximization can be accomplished with respect to the parameters θ n and φ only in the first term of (7). Thus, the M-step with respect to θ n and φ uses only the GLM part of the likelihood function. Hence, the procedure is similar to a standard GLM computation assuming u i is known. Therefore, maximizing with respect to Σ, in the second term, can be handled by the maximum likelihood using the distribution of p u i pu i |Σq after replacing sufficient statistics with the conditional expected values. It thus makes sense to develop a simulation analogous to the Newton-Raphson approach for fitting the GPLMM. Using this separation as in (7), the ML equations for θ n and Σ take the following forms E " Blnp y ij |u i py ij |u i , θ n q Bθ n |y ij ı " 0, The ML estimates of θ n and Σ can be obtained by solving the preceding equations numerically. McCulloch (1997) developed a Monte Carlo Newton-Raphson (MCNR) algorithm for solving these estimating equations, and obtained approximate ML estimates of the parameters for the GMM. Motivated by his work and combining it with the GEE of Liang and Zeger (1986), the optimal estimating equation for θ n is given by where µ i pθ n , u i q " pµ i1 , . . . , µ in i q J and V i pθ n , u i q is the covariance matrix of y i |u i .
In real applications the true intracluster covariance structure is often unknown. The GEE procedure adopts a working covariance matrix, which is specified through a working correlation matrix Rpρq : V i pθ n , u i q " A 1 2 i pθ n , u i qRpρqA 1 2 i pθ n , u i q, where ρ is a finite dimensional parameter and A i pθ n , u i q " diagpν i1 , . . . , ν in i q. Some commonly used working correlation structures include independence, autocorrelation (AR)-1, equally correlated (also called compound symmetry), or unstructured correlation, among others. For a given working correlation structure, ρ can be estimated using the residual-based method of moments. With the estimated working correlation matrix p R " Rpp ρq, the estimating equations in (8) reduces to where D i " pD J i1 , . . . , D J in i q J . We formally define the estimator as the solution p θ n of the above estimating equations. For ease of exposition, we assume φ " 1 and n i " m ă 8 in the rest of the article. Extension of the methodology to the cases of unequal n i is straightforward. We vary the dimension of A i and replace p R by p R i , which is the n iˆni matrix using the specified working correlation structure and the corresponding initial parameter ρ estimator.

Regularization in the GPLMM
In order to select important covariates and estimate them simultaneously, the log likelihood (8) is expanded to include the penalty term ř pn k"1 p λn p|β nk |q which yields the following penalized log likelihood p pβ n , α n , D, φq " where p λ p|β nk |q is any penalty function and λ n is a tuning parameter that determines the amount of shrinkage. Since the coefficients θ n depends to the first and third terms of (10), we propose the penalized estimating equation U n pθ n q " S n pθ n q´q λn p|β n |qsignpβ n q, with q λn p|β n |q "`q λn p|β n1 |q, . . . , q λn p|β npn |q˘J is a p nˆ1 vector of penalty functions, signpβ n q "`signpβ n1 q, . . . , signpβ npn q˘J with signpaq " Ipa ą 0q´Ipa ă 0q and q λn p|β nk |q " p 1 λn p|β nk |q. Note that we assume the semiparametric part contains significant contribution in the model and the proposed penalized estimating equation has been defined to shrink small components of the coefficient β n to zero not α n . Thus, the method performing variable selection for fixed effects, produces estimators of the nonzero components and the nonparametric component.
Among all penalty functions, the smoothing clipped absolute deviation (SCAD) penalty proposed by Fan and Li (2001) can be used to retain the good features of both subset selection and ridge regression, for producing sparse solutions, and to ensure continuity of the selected models. Therefore, we will use the SCAD penalty in our simulation and application studies. The SCAD penalty function is defined by ; a ą 2, where the notation p.q`stands for the positive part of p.q.
Our proposed estimator for θ n is the solution of U n pθ n q " 0. Because U n pθ n q has discontinuous points, an exact solution to U n pθ n q " 0 may not exist. We formally define the estimator p θ n to be an approximate solution, i.e., U n p p θ n q " opa n q for a sequence a n Ñ 0. Alternatively, since the penalty function is singular at the origin, it is challenging to obtain the estimator of θ n by solving U n pθ n q " 0. Following Fan and Li (2001) we locally approximate the penalty function by a quadratic function.
In the neighborhoods of the true parameter values β n0k , |β n0k | ą 0, the derivative of the penalty function is well approximated by for a small numbers e.g. " 10´6. Here, and 0 hn denotes a zero vector of dimension h n .
Note that, in general, the expectations in (11) cannot be computed in a closed form as the conditional distribution of u i |y i involves the marginal distribution of y i , which is not easy to be computed explicitly. Similar to McCulloch (1997), here, we use an alternative method that produces random observations from the conditional distribution of u i |y i by using a Metropolis algorithm where the specification of the density of y i is not required.
In the Metropolis algorithm, p u is chosen as the candidate distribution from which potential new draws are made; then we specify the acceptance function that provides the probability of accepting the new value (as opposed to retaining the previous value). In the forthcoming section we outline the computational procedure used for sample generation.

MCNR algorithm
Let U denote the previous draw from the conditional distribution of U |y, and generate a new value uk for the jth component of U˚" pu 1 , . . . , u k´1 , uk, u k`1 , . . . , u nq q by using the candidate distribution p u , accept U˚as the new value with probability α k pU , U˚q " min ! 1, p u|y pU˚|y, θ n , Dqp u pU |Dq p u|y pU |y, θ n , Dqp u pU˚|Dq ) .
otherwise, reject it and retain the previous value U . The second term in brace in (11) can be simplified to p u|y pU˚|y, θ n , Dqp u pU |Dq p u|y pU |y, θ n , Dqp u pU˚|Dq " Note that, the calculation of the acceptance function α k pU , U˚q here involves only the specification of the conditional distribution of y|u which can be computed in a closed form.
Algorithm 1 describes the Metropolis step into the NewtonRaphson iterative equation (11) for the Monte Carlo estimates of expected values.
Generate N observations U p1q , . . . , U pN q from the distribution p u|y pu|y, θ pm k q n , Σ pm k q q using the Metropolis algorithm. Use these observations to find the Monte Carlo estimates of the expectations. Specially, a) Compute θ pm k`1 q n from the expression step 3. Go to step 2 until convergence is achieved. Choose θ pm k`1 q n and Σ pm k`1 q to be the MCNR estimates of θ n and Σ.

Choice of regularization parameters
To implement the proposed method, several parameters need to be chosen appropriately. One needs to choose the knot sequence in the polynomial spline approximation, λ n , and a in the SCAD penalty function. It is important that the number of distinct knots h, must increase with the sample size n. On the other hand, too many knots would increase the variance of our estimators. Therefore, the number of knots must be properly chosen to balance the bias and variance. For computational convenience, we use equally spaced knots with the number of interior knots L n « n 1{p2r`1q . A similar strategy for knot selection can also be found in He et al. (2002); Qin and Zhu (2007); Sinha and Sattar (2015). To reduce the computational burden, we we follow Fan and Li (2001) and set a " 3.7. Finally we need to choose λ n .
Due to the lack of joint likelihood in the generalized model, to select the tuning parameter λ n we use generalized cross validation (GCV) suggested by Fan and Li (2001) as defined by is the residual sum of squares, and is the effective number of parameters. Then, λ opt is the minimizer of the GCV λn .
Note that W i in (13) is an n iˆni covariance matrix of y i , that can be computed

Asymptotic properties
In this section, we study the asymptotic properties of the estimators p β under the proposed penalized estimating equations. We assume that n i is bounded but the number of subjects n goes to infinity in our asymptotic study. Also consider the framework where the true value of β 0 is partitioned β 0 " pβ J 01 , β J 02 q J and the corresponding design matrix into X i "`X ip1q , X ip2q˘. In our study, the true regression ps " s˚`h n q-dimensional vector of true values that the elements are all nonzero and θ 02 " β 02 " 0. Here, s˚is the dimension of θ 01 and assume that only a small number of covariates contribute to the response i.e. ∫˚" t1 ď j ď p; β j ‰ 0u has cardinal |s| " s˚ă p. Consequently, estimated values and the design matrix is re- Meanwhile, If Eq. (9) has multiple solutions, only a sequence of consistent estimator p θ n is considered. A sequence p θ n is said to be a consistent sequence, if p β n´βn0 Ñ 0 and sup t |B J ptq p α n´f0 ptq| Ñ 0 in probability as n Ñ 8.
The following regularity conditions are require for the main results.
(A.1) n i is a bounded sequence of positive integers, and the distinct values of t ij form a quasi-uniform sequence that grows dense on r0, 1s, and the kth derivative The unknown parameter β n belongs to a compact subset B Ď R p , the true parameter value β n0 lies in the interior of B; (A.4) There exist two positive constants, b 1 and b 2 , such that where λ min (resp. λ max ) denotes the minimum (resp. maximum) eigenvalue of a matrix; (A.5) The common true correlation matrix R 0 has eigenvalues bounded away from zero and`8; the estimated working correlation matrix R satisfies } p R´1Ŕ´1 where R is a constant positive definite matrix with eigenvalues bounded away from zero and`8; we do not require R to be the true There exists a finite constant M 1 ą 0 such that Ep} i pθ n0 , u i q} 2`δ q ď M 1 , for all i and some δ ą 0; and there exist positive constants M 2 and M 3 such that Let B n " tθ n : }θ n´θn0 } ď ∆ a p n {n, then µ . pD J ij θ n q, 1 ď i ď n, 1 ď j ď m, are uniformly bounded away from 0 and 8 on B n ; µ .. pD J ij θ n q and µ p3q pD J ij θ n q, 1 ď i ď n, 1 ď j ď m, are uniformly bounded by a finite positive constant M 2 on B n ; (A.8) Assuming min 1ďkďsn |θ n0k |{λ n Ñ 8 as n Ñ 8 and s 3 n n´1 " op1q, λ n Ñ 0, s n plog nq 2 " opnλ 2 n q, log p n " opnλ 2 n { log nq, p n s 4 n plog nq 6 " opn 2 λ 2 n q, and p n s 3 n plog nq 8 " opn 2 λ 4 n q.
Under (A.1), the total sample size n is of the same order as the number of subjects for independent data. Condition (A.5) is a similar assumption in the Liang and Zeger (1986), assumes that the estimator of the working correlation matrix parameter p τ satisfies ? npp τ´τ 0 q " O p p1q for some τ 0 . Note that when a nonparametric moment estimator is used for the working correlation matrix, we have R " R 0 . The first part of condition (A.6) is similar to the condition in Lemma 2 of Xie and Yang (2003) and condition pÑ δ q in Balan and Schiopu-Kratina (2005); the second part is satisfied for Gaussian distribution, sub-Gaussian distribution, and Poisson distribution, etc. Condition (A.7) requires µ pkq ij pD J ij θ n q, which denotes the kth derivative of µ pkq ij ptq evaluated at D J ij θ n , to be uniformly bounded when θ n is in a local neighborhood around θ n0 , k " 1, 2, 3. This condition is generally satisfied for the GEE. For example, when the marginal model follows a Poisson distribution, µptq " expptq; thus µ pkq ij pD J ij θ n q " exppD J ij θ n q, k " 1, 2, 3, are uniformly bounded around θ n0 on B n . Now, consider the following estimating equation Let M n pθ n q to be the covariance matrix of S n pθq, then However, for the estimating equation approach, alternative techniques are necessary to establish the asymptotic theory, in the supplementary file we present some useful lemmas. By Lemma 1, we approximate f 0 ptq by Bptqα 0 , then have Theorems 1-3 below characterize the existency, consistency and normality of the proposed penalized estimator when p n Ñ 8. P n´| U nk p p θ n q| " 0, k " 1, . . . , sn, psn`1q, . . . , ps n " sn`h n q¯Ñ 1, (ii) P n´| U nk p p θ n q| ď λ n log n , k " psn`h n`1 q, . . . , p n¯Ñ 1, . . , s n , S nk p p θ n q k " ps n`1 q, . . . , p n , and S nk p p θ n q denotes the kth element of S n p p θ n q.
Theorem 2 provides a more precise characterization for the approximate solution of the penalized estimating equation procedure.
Theorem 2. (Consistency). Assume conditions (A1)-(A8) and that n´1p 2 n " op1q. Then, U n pθ n q " op1q has a root p θ n such that The following theorem shows that under regularity conditions, all the covariates with zero coefficients can be detected simultaneously with probability tending to 1, and the estimators of all the non-zero coefficients are asymptotically normally distributed.
Theorem 3. (Asymptotic normality). Assume (A.1)-(A.8). If the number of knots L n " O p pn 1{p2r`1q q, and n´1p 3 n " op1q, then @ξ n P R pn such that }ξ n } " 1, we have Theorem 3 is often referred to as the oracle properties of variable selection, that is, the procedure estimates the true zero coefficient as zero with probability approaching one and estimates the nonzero coefficients as efficiently as if the true model is known in advance.
From the Algorithm 1, we obtain the following sandwich formula to estimate the asymptotic covariance matrix of p θ n : Covp p θ n q « rH n p p θ n , u i q`nE n p p θ n qs´1M n p p θ n , u i qrH n p p θ n , u i q`nE n p p θ n , u i qs´1, where H n and E n are defined in Section 3, and

Numerical Studies
In this section, we first conduct simulation study to illustrate the consistency and the sure screening property of the proposed procedure empirically, and compare its finite sample performance with some other different model settings. We further apply our proposed method for analyzing two real data sets.

Simulation studies
We simulated 100 data sets from the GSMM for longitudinal data. The underlying model is the random intercept Poisson model so that the nonlinear function were set to a sinusoidal function y ij |b i ∼ Poispµ ij q, i " 1, . . . , n, j " 1, . . . , n i , where i " 1, . . . , n (n " 50, 100 and 150), and j " 1, . . . , n i which the number of observations per subjects assumed to be fixed at n i " 5. The true regression coefficients are β " p´1,´1, 2, 0, . . . , 0q with the mutually independent covariates X J ij " px p1q ij , . . . , x ppq ij q are drawn independently from uniform distribution on p´1, 1q. The measurement time points t ij are drawn from uniform distribution on p0, 1q. The random effect process b i is taken to be a Gaussian process with mean 0, variance σ 2 " 0.25. The predictor dimension p n is diverging but the dimension of the true model is fixed to be 3.
Regarding the choice of the dimensionality of the parametric component, p n , authors recommended many suggestions as a sensible choice. For example p n " r n 2 s, p n " r4.5n 1{4 s, and p n " r n b logpnq s, where b ą 1 and ras stands for the largest integer no larger than a. These only discuss the situation p Ñ 8 as n Ñ 8 with p n ă n.
For case p n ąą n, we can mention to logpp n q " o p pn b q, where 0 ă b ă 1. Of course challenges arise when p is much larger than n, choosing a larger value of p n increases the probability that variable selection methods will include all of the correct variables, but including more inactive variables will tend to have a slight detrimental effect on the performance of the final variable selection and parameter estimation method. We have found that this latter effect is most noticeable in models where the response provides less information. In particular, the binary response of a logistic regression model and, to a lesser extent, the integer valued response in a Poisson regression model are less informative than the real valued response in a linear model. We therefore used the pairs of pn, p n q as p50, 11q, p100, 14q, p150, 16q and p30, 100q, p100, 500q, p200, 2000q respectively for cases p n ă n and p n ąą n. The choices seem to be satisfactory and consist of all suggestions, as the performance of the procedures is quite robust to different choices of p n .
The performance of estimators was evaluated separately for the structural com-ponents and the variance. By averaging Performance of the proposed penalized procedure compared with the unpenalized one and the penalized GLMM where each simulated data set was fitted under these three methods. For evaluating estimation accuracy, we report the empirical mean square error (MSE), defined as ř 100 k"1 } p β k n´βn0 }{100 where p β k n is the estimator of β n0 obtained using the kth generated data set. The performance in variable selection is The results of Table 1 summarize the estimation accuracy and model selection properties of the penalized GSMM (P-GSMM), the unpenalized GSMM, and the penalized GLMM (P-GLMM) for the different values of pn, p n q. In terms of estimation accuracy the penalized GSMM procedure performs closely to the penalized GLMM, whereas our proposed approach gives the smallest MSE, and consistently outperforms its penalized GLMM counterpart. In terms of model selection we observe that the unpenalized GSMM generally does not lead to a sparse model. Furthermore, the penalized GSMM and the penalized GLMM successfully selects all covariates with nonzero coefficients (i.e., I rates are zero), but it is obvious that the proposed approach has slightly stronger sparsity (i.e., a fairly higher number of Cs) than the penalized GLMM. For penalized GSMM, The probability of identifying the exact underlying model is about 80% and this rate grows by increasing the sample size, confirming the good asymptotic properties of the penalized estimators. The results are the same in both cases of p n ă n and p n ąą n, but when p n ąą n zero coefficients tends to increasingly included in the model.
To further investigate the performance of the proposed method, Table 2 reports its bias, the estimated standard deviation (calculated from the sandwich variance formula), the empirical standard deviation, and the empirical coverage probability of 95% confidence interval for estimating β 1 , β 2 , and β 3 . The estimated standard deviation is close to the empirical standard deviation, and the empirical coverage probability is close to 95%. This indicated good performance of the sandwich variance formula.
These observations suggest that considering partial part is important to modify the estimation accuracy and model selection when the growth curves of the data exhibit a nonlinear fashion over time, especially in a complicated manner. On the other hand, penalized GSMM allows us to make systematic inference on all model parameters by representing a partially model as a modified penalized GLMM.  For proposed method, the estimated baseline function f ptq is also evaluated through visualization. We plot and compare the estimated f ptq and pointwise biases, for the cases of p n ă n and p n ąą n by two sample size n " 50 and 100. We also plot the pointwise standard deviations (calculated from the sandwich variance formula), and coverage probability of 95% confidence intervals. Figures 1 shows that the p n ă n case, our approach yields smaller overall biases and standard deviations than the p n ąą n case. Also, it can be seen that the empirical coverage probability for f ptq is close to 95% for two cases. Figure 2 depicts the results for n " 100. As shown, larger sample size modified the biases and the differences between two cases.
Nevertheless, the case of p n ă n has smaller standard deviation.

Real data analyses 4.2.1 AIDS data
In this section, to illustrate our method, we considered the longitudinal CD4 cell count data among HIV seroconverters. This dataset contains 2376 observations of CD4 cell counts on 369 men infected with the HIV virus; see (Zeger and Diggle , 1994) for a detailed description of this dataset. Figure 3  CESD scale (larger values indicate increased depressive symptoms) and YEAR (the effect of time since seroconversion). Since there seems to exist a positive correlation among responses from the same patient, we need to incorporate a correlation structure into the estimation scheme. Zeger and Diggle (1994) found that the compound symmetry covariance matrix fitted the data reasonably well. This data analysed by many authors such as Wang et al. (2005), Huang et al. (2007) and Ma et al. (2013).
Their analysis was conducted on square root transformed CD4 numbers whose distribution is more nearly Gaussian. In our analysis, we fit the data using an GSMM, without transforming the CD4 by adopting the Poisson regression. To take advantage of flexibility of partially linear models, we let YEAR be modeled nonpara-metrically, the remaining parametrically. It is of interest to examine whether there are any interaction effects between the parametric covariates, so we included all these interactions in the parametric part. We further applied the proposed approach to select significant variables. We used the SCAD penalty, and the tuning parameter λ " 0.45. To compare the performance of our proposed method (P-SMM) with other two existing scenarios, including the unpenalized GSMM, and the penalized GLMM (P-GLMM), we use the standard errors (SE) were all calculated using the sandwich method. To best identify a model supported by the data, we adopt the Akaike information criterion (AIC; Akaike (1973)) and the Bayesian information criterion (BIC; Schwarz (1978)). They are defined as where max is the maximized log-likelihood value, m is the number of free parameters in the model. Table 3 presents the summary of the fitting results including the values of standard errors, together with max , AIC, and BIC under the three models.
Judging from Table 3, the P-GSMM tends to exhibit slightly standard errors compared to GSMM and P-GLMM, nevertheless this difference is not more dramatic.
Meanwhile, the values of AIC, BIC of our proposed model are smaller than those for the other two competing models, revealing that the P-GSMM can provide better fitting performance. Under P-GSMM, SMOKE, DRUGS, SEXP, SOM KE˚DRU G and DRU G˚SEXP are identifies as significant covariates. One notes some slight selection difference when P-GLMM is used, which suggests that AGE˚SM OKE, AGE˚DRU G, and SM OKE˚SEXP may also be significant. We also find some significant interactions among some covariates which may be ignored by Wang et al. (2005) and Huang et al. (2007). The Results for nonparametric curve estimates using the P-GSMM estimators are plotted in Figure 3 for Y EAR. It shows the Table 3: Summary of parameter estimates along with standard errors (in parentheses) under the three fitted models for the AIDS data. estimated nonparametric function f ptq, its 95% pointwise confidence bands, standard deviation, and 95% coverage probability given by the empirical and sandwich formula variance. We can see that the baseline function f ptq has decreasing effect as time passing. Therefore, one can see that it is more reasonable to put it as a nonparametric component. We notice the disparity between the empirical and the sandwich formula standard deviation in the boundary positions and the sandwich formula standard deviations are smaller in which case the coverage probability recede from 95%.

Yeast Cell-Cycle Gene Expression data
A yeast cell-cycle gene expression data collected in the CDC15 experiment of Spellman et al. (1998) where genome-wide mRNA levels of 6178 yeast ORFs (abbreviation for open reading frames, which are DNA sequences that can determine which amino acids will be encoded by a gene) at 7 minute intervals for 119 minutes, which covers two cell-cycle periods for a total of 18 time points, measured. The cell cycle is a tightly regulated life process where cells grow, replicate their DNA, segregate their chromosomes, and divide into as many daughter cells as the environment allows.
The cell-cycle process is commonly divided into M/G1-G1-S-G2-M stages. Refer to Wang et al. (2012), for more detailed description of this data set.
Transcription factors (TFs) have been observed to play critical roles in gene expression regulation. A TF (sometimes called a sequence-specific DNA-binding factor) is a protein that binds to specific DNA sequences, thereby controlling the flow (or transcription) of genetic information from DNA to mRNA. To better understand the phenomenon underlying cell-cycle process, it is important to identify TFs that regulate the gene expression levels of cell cycle-regulated genes. It is not clear where these TFs regulate all cell cycle genes, however.
We applied our methods to identify the key TFs. The dataset that we use present a subset of 283 cell-cycled-regularized genes observed over 4 time points at G1 stage.
The response variable Y ij is the log-transformed gene expression level of gene i measured at time point j, for i " 1, . . . , 283 and t " 1, . . . , 4. We use the following semiparametric mixed model where the covariates x pkq ij , k " 1, . . . , 96, is the matching score of the binding probability of the kth TF on the promoter region of the ith gene. The binding probability is computed using a mixture modeling approach based on data from a ChIP binding experiment; see Wang et al. (2007) for details. Covariates x pkq ij is standardized to have mean zero and variance 1. t ij denotes time, f pt ij q models the nonparametric time effect, and b i is the random intercept. Our goal is to identify the TFs that might be related to the expression patterns of these 283 cellcycleregulated genes. Therefore we apply a penalization procedure by the proposal P-GSMM and also by ignoring the nonparametric component f pt ij q using P-GLMM. Table 4: Summary of parameter estimates along with standard errors (in parentheses) under the P-GSMM and P-GLMM for the Yeast Cell-Cycle Gene Expression data.   Table 4 summarizes the TFs identified when p-GSMM and p-GLMM are adopted.
Our analysis reveals that a total of 13 and 16 TFs related to yeast cell-cycle processes are identified respectively by the P-GSMM and P-GLMM. The sets of TFs selected at different methods have only small overlaps. These common Tfs are GAT3, MBP1, MSN4, PHD1, SMP1, SWI4, and SWI6. For stage G1, MBP1, SWI4, and SWI6 are three TFs that have been proved important in the aforementioned biological experiments and they have been selected by the two methods. However, model selection criteria, including the values of standard errors, together with max , AIC, and BIC confirm the superiority of our proposed model.

Concluding Remarks
In general, when the number of covariates is large, identifying the exact underlying model is a challenging task, in particular when some of the nonzero signals are relatively weak. Here, We have developed a general methodology for simultaneously selecting variables and estimating the unknown components in semiparametric mixed effects model for non Gaussian longitudinal data when the number of parameters diverges with the sample size. We approximate the nonparametric components based on polynomial B-spline smoothing, propose a penalized estimating equation approach to do variable selection, and the estimation of the both parametric and nonparametric components obtained in the framework of maximum likelihood estimation. We apply the penalty functions to the estimating equation objective function such that the proposed procedure can simultaneously estimate parameters and select the important variables. The procedure involves the specification of the posterior distribution of the random effects, which cannot be evaluated in a closed form, and we use a Metropolis algorithm, which does not require the specification of the posterior distribution. To implement the procedure in practice, a computationally flexible iterative algorithm is developed. We obtain asymptotic normality of the resulting estimators and demonstrate that the asymptotic normality of the estimated coefficients for the linear part is retained. To investigate the performance of our approach, we compared it with the unpenalized generalized semiparametric mixed effects model and penalized generalized linear mixed effects model throw a simulation study and the analysis of two data set including AIDS data and Yeast Cell-Cycle Gene Expression data. Results show that the proposal works well, and can estimate the nonzero coefficients and select model simultaneously, the model outperforms the penalized generalized linear mixed effects counterparts on the provision of likelihood-based model selection and estimating the regression coefficients. In addition, we can find that estimation is more efficient when the partially part is taken into consideration. The results are consistent in both cases of p n ă n and p n ąą n. Our primary interest is the linear components, and we treat the nonparametric functions as nuisance components; thus we limit our discussions to estimation and variable selection for the linear part.
Nonetheless, this may be extended to the nonparametric components using techniques similar to those in Xue (2009). However, the technical details deserve careful consideration, and this is an interesting topic of future research.