Copyright © Queen’s Printer and Controller of HMSO 2017. This work was produced by Hemingway et al.![Tutorial Tutorial](/uploads/1/2/5/7/125719167/391712571.jpg)
OpenSCAD is a powerful, free and open source modelling tool that lets you make precise 3D models with just a few lines of code. With the right setup, you can even make adaptable, 'parametric' models. In this OpenSCAD tutorial, we'll get acquainted with the software, review a simple example and then discuss how to make parametric designs. Creating a design table to drive Parametric Blocks parameters Create a design table. First, you must create the design table as a CSV (comma separated) file. You can use a spreadsheet app or a simple text editor. Below are the steps to use Wordpad, but you can use any app that allows you to save it as plain text with a.CSV extension.
under the terms of a commissioning contract issued by the Secretary of State for Health. This issue may be freely reproduced for the purposes of private research and study and extracts (or indeed, the full report) may be included in professional journals provided that suitable acknowledgement is made and the reproduction is not associated with any form of advertising. Applications for commercial reproduction should be addressed to: NIHR Journals Library, National Institute for Health Research, Evaluation, Trials and Studies Coordinating Centre, Alpha House, University of Southampton Science Park, Southampton SO16 7NS, UK.Included under terms of UK Non-commercial Government License.
Ann Stat. Author manuscript; available in PMC 2012 Dec 12.
Published in final edited form as:
doi: 10.1214/11-AOS885SUPP
NIHMSID: NIHMS355290
See other articles in PMC that cite the published article.
Abstract
We study generalized additive partial linear models, proposing the use of polynomial spline smoothing for estimation of nonparametric functions, and deriving quasi-likelihood based estimators for the linear parameters. We establish asymptotic normality for the estimators of the parametric components. The procedure avoids solving large systems of equations as in kernel-based procedures and thus results in gains in computational simplicity. We further develop a class of variable selection procedures for the linear parameters by employing a nonconcave penalized quasi-likelihood, which is shown to have an asymptotic oracle property. Monte Carlo simulations and an empirical example are presented for illustration.
Key words and phrases: Backfitting, generalized additive models, generalized partially linear models, LASSO, nonconcave penalized likelihood, penalty-based variable selection, polynomial spline, quasi-likelihood, SCAD, shrinkage methods
1. Introduction
Generalized linear models (GLM), introduced by Nelder and Wedderburn (1972) and systematically summarized by McCullagh and Nelder (1989), are a powerful tool to analyze the relationship between a discrete response variable and covariates. Given a link function, the GLM expresses the relationship between the dependent and independent variables through a linear functional form. However, the GLM and associated methods may not be flexible enough when analyzing complicated data generated from biological and biomedical research. The generalized additive model (GAM), a generalization of the GLM that replaces linear components by a sum of smooth unknown functions of predictor variables, has been proposed as an alternative and has been used widely [Hastie and Tibshirani (1990), Wood (2006)]. The generalized additive partially linear model (GAPLM) is a realistic, parsimonious candidate when one believes that the relationship between the dependent variable and some of the covariates has a parametric form, while the relationship between the dependent variable and the remaining covariates may not be linear. GAPLM enjoys the simplicity of the GLM and the flexibility of the GAM because it combines both parametric and nonparametric components.
There are two possible approaches for estimating the parametric component and the nonparametric components in a GAPLM. The first is a combination of kernel-based backfitting and local scoring, proposed by Buja, Hastie and Tibshirani (1989) and detailed by Hastie and Tibshirani (1990). This method may need to solve a large system of equations [Yu, Park and Mammen (2008)], and also makes it difficult to introduce a penalized function for variable selection as given in Section 4. The second is an application of the marginal integration approach [Linton and Nielsen (1995)] to the nonparametric component of the generalized partial linear models. They treated the summand of additive terms as a nonparametric component, which is then estimated as a multivariate nonparametric function. This strategy may still suffer from the “curse of dimensionality” when the number of additive terms is not small [Härdle et al. (2004)].
The kernel-based backfitting and marginal integration approaches are computationally expensive. Marx and Eilers (1998), Ruppert, Wand and Carroll (2003) and Wood (2004) studied penalized regression splines, which share most of the practical benefits of smoothing spline methods, combined with ease of use and reduction of the computational cost of backfitting GAMs. Widely used R/Splus packages gam and mgcv provide a convenient implementation in practice. However, no theoretical justifications are available for these procedures in the additive case. See Li and Ruppert (2008) for recent work in the one-dimensional case.
In this paper, we will use polynomial splines to estimate the nonparametric components. Besides asymptotic theory, we develop a flexible and convenient estimation procedure for GAPLM. The use of polynomial spline smoothing in generalized nonparametric models goes back to Stone (1986), who first obtained the rate of convergence of the polynomial spline estimates for the generalized additive model. Stone (1994) and Huang (1998) investigated polynomial spline estimation for the generalized functional ANOVA model. More recently, Xue and Yang (2006) studied estimation of the additive coefficient model for a continuous response variable using polynomial spline methods. Our models emphasize possibly non-Gaussian responses, and combine both parametric and nonparametric components through a link function. Estimation is achieved through maximizing the quasi-likelihood with polynomial spline smoothing for the nonparametric functions. The convergence results of the maximum likelihood estimates for the nonparametric parts in this article are similar to those for regression established by Xue and Yang (2006). However, it is very challenging to establish asymptotic normality in our general context, since it cannot be viewed simply as an orthogonal projection, due to its nonlinear structure. To the best of our knowledge, this is the first attempt to establish asymptotic normality of the estimators for the parametric components in GAPLM. Moreover, polynomial spline smoothing is a global smoothing method, which approximates the unknown functions via polynomial splines characterized by a linear combination of spline basis. After the spline basis is chosen, the coefficients can be estimated by an efficient one-step procedure of maximizing the quasi-likelihood function. In contrast, kernel-based methods, such as those reviewed above, in which the maximization must be conducted repeatedly at every data point or a grid of values, are more time-consuming. Thus, the application of polynomial spline smoothing in the current context is particularly computationally efficient compared to some of its counterparts.
In practice, a large number of variables may be collected and some of the insignificant ones should be excluded before forming a final model. It is an important issue to select significant variables for both parametric and nonparametric regression models; see Fan and Li (2006) for a comprehensive overview of variable selection. Traditional variable selection procedures such as stepwise deletion and subset selection may be extended to the GAPLM. However, these are also computationally expensive because, for each submodel, we encounter the challenges mentioned above.
To select significant variables in semiparametric models, adopted Fan and Li’s (2001) variable selection procedures for parametric models via nonconcave penalized quasi-likelihood, but their models do not cover the GAPLM. Of course, before developing justifiable variable selection for the GAPLM, it is important to establish asymptotic properties for the parametric components. In this article, we propose a class of variable selection procedures for the parametric component of the GAPLM and study the asymptotic properties of the resulting estimator. We demonstrate how the rate of convergence of the resulting estimate depends on the regularization parameters, and further show that the penalized quasi-likelihood estimators perform asymptotically as an oracle procedure for selecting the model.
The rest of the article is organized as follows. In Section 2, we introduce the GAPLM model. In Section 3, we propose polynomial spline estimators via a quasi-likelihood approach, and study the asymptotic properties of the proposed estimators. In Section 4, we describe the variable selection procedures for the parametric component, and then prove their statistical properties. Simulation studies and an empirical example are presented in Section 5. Regularity conditions and the proofs of the main results are presented in the Appendix.
2. The models
Let Y be the response variable, X = (X1, …, Xd1)T ∈ Rd1 and Z = (Z1, …, Zd2)T ∈ Rd2 be the covariates. We assume the conditional density of Y given (X, Z) = (x, z) belongs to the exponential family
(1)
for known functions ℬ and ?, where ξ is the so-called natural parameter in parametric generalized linear models (GLM), is related to the unknown mean response by
In parametric GLM, the mean function μ is defined via a known link function g by g{μ(x, z)} = xTα + zTβ, where α and β are parametric vectors to be estimated. In this article, g(μ) is modeled as an additive partial linear function
(2)
where β is a d2-dimensional regression parameter, are unknown and smooth functions and E{ηk(Xk)} = 0 for 1 ≤ k ≤ d1 for identifiability.
If the conditional variance function var(Y|X = x, Z = z) = σ2V {μ(x, z)} for some known positive function V, then estimation of the mean can be achieved by replacing the conditional loglikelihood function log{fY|X,Z(y|x, z)} in (1) by a quasi-likelihood function Q(m, y), which satisfies
The first goal of this article is to provide a simple method of estimating β and in model (2) based on a quasi-likelihood procedure [Severini and Staniswalis (1994)] with polynomial splines. The second goal is to discuss how to select significant parametric variables in this semiparametric framework.
3. Estimation method
3.1. Maximum quasi-likelihood
Let (Yi, Xi, Zi), i = 1, …, n, be independent copies of (Y, X, Z). To avoid confusion, let and β0 be the true additive function and the true parameter values, respectively. For simplicity, we assume that the covariate Xk is distributed on a compact interval [ak, bk], k = 1, …, d1, and without loss of generality, we take all intervals [ak, bk] = [0, 1], k = 1, …, d1. Under smoothness assumptions, the η0k’s can be well approximated by spline functions. Let ?n be the space of polynomial splines on [0, 1] of order r ≥ 1. We introduce a knot sequence with J interior knots
ξ−r+1 = ⋯ = ξ−1 = ξ0 = 0 < ξ1 < ⋯ < ξJ < 1 = ξJ+1 = ⋯ = ξJ+r,
where J ≡ Jn increases when sample size n increases, where the precise order is given in condition (C5) in Section 3.2. According to Stone (1985), ?n consists of functions ћ satisfying:
- ћ is a polynomial of degree r − 1 on each of the subintervals Ij = [ξj, ξj + 1), j = 0, …, Jn − 1, IJn = [ξJn, 1];
- for r ≥ 2, ћ is r − 2 times continuously differentiable on [0, 1].
Equally-spaced knots are used in this article for simplicity of proof. However other regular knot sequences can also be used, with similar asymptotic results.
We will consider additive spline estimates of η0. Let ?n be the collection of functions η with the additive form , where each component function . We seek a function η ∈ ?n and a value of β that maximize the quasi-likelihood function
(3)
For the kth covariate xk, let bj,k(xk) be the B-spline basis functions of order r. For any η ∈ ?n, write η(x) = γTb(x), where b(x) = {bj,k(xk), j = 1, …, Jn + r, k = 1, …, d1}T is the collection of the spline basis functions, and γ = {γj,k, j = 1, …, Jn + r, k = 1, …, d1}T is the spline coefficient vector. Thus, the maximization problem in (3) is equivalent to finding β and γ to maximize
(4)
We denote the maximizer as and = {j,k, j = 1, …, Jn + r, k = 1, …, d1}T. Then the spline estimator of η0 is (x) = Tb(x), and the centered spline component function estimators are
The above estimation approach can be easily implemented because this approximation results in a generalized linear model. However, theoretical justification for this estimation approach is very challenging [Huang (1998)].
Let Nn = Jn + r − 1. We adopt the normalized B-spline space introduced in Xue and Yang (2006) with the following normalized basis
which is convenient for asymptotic analysis. Let B(x) = {Bj,k(xk), j = 1, …, Nn, k = 1, …, d1}T and Bi = B(Xi). Finding (γ, β) that maximizes (4) is mathematically equivalent to finding (γ, β) which maximizes
Then the spline estimator of η0 is (x) = TB(x), and the centered spline estimators of the component functions are
We show next that estimators of both the parametric and nonparametric components have nice asymptotic properties.
3.2. Assumptions and asymptotic results
Let υ be a positive integer and α ∈ (0, 1] such that p = υ + α > 2. Let ℋ(p) be the collection of functions g on [0, 1] whose υth derivative, g(υ), exists and satisfies a Lipschitz condition of order α, |g(υ)(m*) − g(υ)(m)| ≤ C|m* − m|α, for 0 ≤ m*, m ≤ 1, where C is a positive constant. Following the notation of Carroll et al. (1997), let ρℓ(m) = {dg−1(m)/dm}ℓ/V {g−1(m)} and qℓ(m, y) = ∂ℓ/∂mℓQ{g−1(m), y}, so that
For simplicity of notation, write T = (X, Z) and A⊗2 = AAT for any matrix or vector A. We make the following assumptions:
- The function is continuous and each component function η0k(·) ∈ ℋ(p), k = 1, …, d1.
- (C2)The function q2(m, y) < 0 and for m ∈ R and y in the range of the response variable.
- The distribution of X is absolutely continuous and its density f is bounded away from zero and infinity on [0, 1]d1.
- (C4)The random vector Z satisfies that for any unit vector ω ∈ Rd2
- (C5)The number of knots n1/(2p) ≪ Nn ≪ n1/4.
Remark 1. The smoothness condition in (C1) describes a requirement on the best rate of convergence that the functions η0k(·)’s can be approximated by functions in the spline spaces. Condition (C2) is imposed to ensure the uniqueness of the solution; see, for example, Condition 1a of Carroll et al. (1997) and Condition (i) of . Condition (C3) requires a boundedness condition on the covariates, which is often assumed in asymptotic analysis of nonparametric regression problems; see Condition 1 of Stone (1985), Assumption (B3)(ii) of Huang (1999) and Assumption (C1) of Xue and Yang (2006). The boundedness assumption on the support can be replaced by a finite third moment assumption, but this will add much extra complexity to the proofs. Condition (C4) implies that the eigenvalues of E(Z⊗2|X = x) are bounded away from 0 and ∞. Condition (C5) gives the rate of growth of the dimension of the spline spaces relative to the sample size.
For measurable functions ϕ1, ϕ2 on [0, 1]d1, define the empirical inner product and the corresponding norm as
If ϕ1 and ϕ2 are L2-integrable, define the theoretical inner product and corresponding norm as
Let be the empirical and theoretical norm of ϕ on [0, 1], defined by
where fk(·) is the density function of Xk.
Theorem 1 describes the rates of convergence of the nonparametric parts.
Theorem 1. Under conditions (C1)–(C5), for k = 1, …, d1, .
Let m0(T) = η0(X) + ZTβ0 and define
where
(7)
is the projection of Γ onto the Hilbert space of theoretically centered additive functions with a norm . To obtain asymptotic normality of the estimators in the linear part, we further impose the conditions:
- The additive components in (7) satisfy that Γk(·) ∈ ℋ(p), k = 1, …, d1.
- (C7)For ρℓ, we have
- There exists a positive constant C0, such that E[{Y − g−1(m0(T))}2|T] ≤ C0, almost surely.
The next theorem shows that the maximum quasi-likelihood estimator of β0 is root-n consistent and asymptotically normal, although the convergence rate of the nonparametric component η0 is of course slower than root-n.
Theorem 2. Under conditions (C1)–(C8), , whereΩ = E[ρ2{m0(T)}⊗2].
The proofs of these theorems are given in the Appendix.
It is worthwhile pointing out that taking the additive structure of the nuisance parameter into account leads to a smaller asymptotic variance than that of the estimators which ignore the additivity [Yu and Lee (2010)]. Carroll et al. (2009) had the same observation for a special case with repeated measurement data when g is the identity function.
4. Selection of significant parametric variables
In this section, we develop variable selection procedures for the parametric component of the GAPLM. We study the asymptotic properties of the resulting estimator, illustrate how the rate of convergence of the resulting estimate depends on the regularization parameters, and further establish the oracle properties of the resulting estimate.
4.1. Penalized likelihood
Building upon the quasi-likelihood given in (3), we define the penalized quasi-likelihood as
where pλj(·) is a prespecified penalty function with a regularization parameter λj. The penalty functions and regularization parameters in (8) are not necessarily the same for all j. For example, we may wish to keep scientifically important variables in the final model, and therefore do not want to penalize their coefficients. In practice, λj can be chosen by a overflow='scroll'>p λ j ( | β | ) = 0.5 λ j 2 I ( | β | ≠ 0 ) . The traditional best-subset variable selection can be viewed as a penalized least squares with the L0 penalty because is essentially the number of nonzero regression coefficients in the model. Of course, this procedure has two well known and severe problems. First, when the number of covariates is large, it is computationally infeasible to do subset selection. Second, best subset variable selection suffers from high variability and instability [Breiman (1996), Fan and Li (2001)].
The Lasso is a regularization technique for simultaneous estimation and variable selection [Tibshirani (1996), Zou (2006)] that avoids the drawbacks of the best subset selection. It can be viewed as a penalized least squares estimator with the L1 penalty, defined by pλj(|β|) = λj|β|. Frank and Friedman (1993) considered bridge regression with an Lq penalty, in which pλj(|β|) = λj|β|q (0 < q < 1). The issue of selection of the penalty function has been studied in depth by a variety of authors. For example, Fan and Li (2001) suggested using the SCAD penalty, defined by
where pλj(0) = 0, and λj and a are two tuning parameters. Fan and Li (2001) suggested using a = 3.7, which will be used in Section 5.
Substituting η by its estimate in (8), we obtain a penalized likelihood
Maximizing ℒP(β) in (9) yields a maximum penalized likelihood estimator MPL. The theorems established below demonstrate that MPL performs asymptotically as well as an oracle estimator.
4.2. Sampling properties
We next show that with a proper choice of λj, the maximum penalized likelihood estimator MPL has an asymptotic oracle property. Let , where β10 is assumed to consist of all nonzero components of β0 and β20 = 0 without loss of generality. Similarly we write Denote and
(10)
Theorem 3. Under the regularity conditions given in Section 3.2, and if an → 0 and wn → 0 as n → ∞, then there exists a local maximizerMPLof ℒP(β) defined in(9)such that its rate of convergence is OP(n−1/2 + an), where an is given in(10).
Next, define and a diagonal matrix , where s is the number of nonzero components of β0. Define T1 = (X, Z1) and , and further let
where is the projection of Γ1 onto the Hilbert space of theoretically centered additive functions with the norm .
Theorem 4. Suppose that the regularity conditions given in Section 3.2 hold, and that. If, then the root-n consistent estimatorMPLin Theorem 3 satisfies.
4.3. Implementation
As pointed out by , many penalty functions, including the L1 penalty and the SCAD penalty, are irregular at the origin and may not have a second derivative at some points. Thus, it is often difficult to implement the Newton–Raphson algorithm directly. As in Fan and Li (2001), , we approximate the penalty function locally by a quadratic function at every step in the iteration such that the Newton–Raphson algorithm can be modified for finding the solution of the penalized likelihood. Specifically, given an initial value β(0) that is close to the maximizer of the penalized likelihood function, the penalty pλj(|βj|) can be locally approximated by the quadratic function as when is not very close to 0; otherwise, set j = 0. In other words, for . For instance, this local quadratic approximation for the L1 penalty yields
Standard error formula forMPL. We follow the approach in to derive a sandwich formula for the estimator MPL. Let
A sandwich formula is given by
Following conventional techniques that arise in the likelihood setting, the above sandwich formula can be shown to be a consistent estimator and will be shown in our simulation study to have good accuracy for moderate sample sizes.
Choice of λj’s. The unknown parameters (λj) can be selected using>β ̂ MPL, we maximize ℓ(γ, MPL) with respect to γ. The solution is denoted by MPL, and the corresponding estimator of η0 is defined as
(11)
Here the GCV statistic is defined by
where e(λ1, …, λd2) = tr[{ℓ″(MPL) − nΣλ(MPL)}−1 ℓ″(MPL)] is the effective number of parameters and D(Y, μ) is the deviance of Y corresponding to fitting with λ. The minimization problem over a d2-dimensional space is difficult. However, conjectured that the magnitude of λj should be proportional to the standard error of the unpenalized maximum pseudo-partial likelihood estimator of βj. Thus, we suggest taking λj = λSE(j) in practice, where SE(j) is the estimated standard error of j, the unpenalized likelihood estimate defined in Section 3. Then the minimization problem can be reduced to a one-dimensional problem, and the tuning parameter can be estimated by a grid search.
![Tutorial Tutorial](/uploads/1/2/5/7/125719167/391712571.jpg)
5. Numerical studies
5.1. A simulation study
We simulated 100 data sets consisting of n = 100, 200 and 400 observations, respectively, from the GAPLM:
(12)
where
and the true parameters β = (3, 1.5, 0, 0, 0, 0, 2, 0)T. X1 and X2 are independently uniformly distributed on [0, 1]. Z1 and Z2 are normally distributed with mean 0.5 and variance 0.09. The random vector (Z1, …, Z6, X1, X2) has an autoregressive structure with correlation coefficient ρ = 0.5.
In order to determine the number of knots in the approximation, we performed a simulation with 1,000 runs for each sample size. In each run, we fit, without any variable selection procedure, all possible spline approximations with 0–7 internal knots for each nonparametric component. The internal knots were equally spaced quantiles from the simulated data. We recorded the combination of the numbers of knots used by the best approximation, which had the smallest prediction error (PE), defined as
(2, 2) and (5, 3) are most frequently chosen for sample sizes 100 and 400, respectively. These combinations were used in the simulations for the variable selection procedures.
The proposed selection procedures were applied to this model and B-splines were used to approximate the two nonparametric functions. In the simulation and also the empirical example in Section 5.2, the estimates from ordinary logistic regression were used as the starting values in the fitting procedure.
To study model fit, we also defined model error (ME) for the parametric part by
(14)
The relative model error is defined as the ratio of the model error between the fitted model using variable selection methods and using ordinary logistic regression.
The simulation results are reported in Table 1, in which the columns labeled with “C” give the average number of the five zero coefficients correctly set to 0, the columns labeled with “I” give the average number of the three nonzero coefficients incorrectly set to 0, and the columns labeled with “MRME” give the median of the relative model errors.
TABLE 1
Results from the simulation study in Section 5.1. C, I and MRME stand for the average number of the five zero coefficients correctly set to 0, the average number of the three nonzero coefficients incorrectly set to 0, and the median of the relative model errors. The model errors are defined in (14)
n | Method | C | I | MRME |
---|---|---|---|---|
100 | ORACLE | 5 | 0 | 0.27 |
SCAD | 4.29 | 0.93 | 0.60 | |
Lasso | 3.83 | 0.67 | 0.51 | |
BIC | 4.53 | 0.95 | 0.54 | |
400 | ORACLE | 5 | 0 | 0.33 |
SCAD | 4.81 | 0.27 | 0.49 | |
Lasso | 3.89 | 0.10 | 0.67 | |
BIC | 4.90 | 0.35 | 0.46 |
Summarizing Table 1, we conclude that BIC performs the best in terms of correctly identifying zero coefficients, followed by SCAD and LASSO. On the other hand, BIC is also more likely to set nonzero coefficients to zero, followed by SCAD and LASSO. This indicates that BIC most aggressively reduce the model complexity, while LASSO tends to include more variables in the models. SCAD is a useful compromise between these two procedures. With an increase of sample sizes, both SCAD and BIC nearly perform as if they had Oracle property. The MRME values of the three procedures are comparable. Results of the cases not depicted here have characteristics similar to those shown in Table 1. Readers may refer to the online supplemental materials.
We also performed a simulation with correlated covariates. We generated the response Y from model (12) again but with β = (3.00, 1.50, 2.00). The covariates Z1, Z2, X1 and X2 were marginally normal with mean zero and variance 0.09. In order, (Z1, Z2, X1, X2) had autoregressive correlation coefficient ρ, while Z3 is Bernoulli with success probability 0.5. We considered two scenarios: (i) moderately correlated covariates (ρ = 0.5) and (ii) highly correlated (ρ = 0.7) covariates. We did 1,000 simulation runs for each case with sample sizes n = 100, 200 and 400. From our simulation, we observe that the estimator becomes more unstable when the correlation among covariates is higher. In scenario (i), all simulation runs converged. However, there were 6, 3 and 7 cases of nonconvergence over the 1,000 simulation runs for sample sizes 100, 200 and 400, respectively, in scenario (ii). In addition, the variance and bias of the fitted functions in scenario (ii) were much larger than those in scenario (i), especially on the boundaries of the covariates’ support. This can be observed in Figures 1 and and2,2, which present the mean, absolute value of bias and variance of the fitted nonparametric functions for ρ = 0.5 and ρ = 0.7 with sample size n = 100. Similar results are obtained for sample sizes n = 200 and 400, but are not given here.
The mean, absolute value of the bias and variance of the fitted nonparametric functions when n = 100 and ρ = 0.5 [the left panel for η1(x1) and the right for η2(x2)]. 95% CB stands for the 95% confidence band.
The mean, absolute value of the bias and variance of the fitted nonparametric functions when n = 100 and ρ = 0.7. The left panel is for η1(x1) and the right panel is for η2(x2). Here 95% CB stands for the 95% confidence band.
5.2. An empirical example
We now apply the GAPLM and our variable selection procedure to a data set from the Pima Indian diabetes study [Smith et al. (1988)]. This data set is obtained from the UCI Repository of Machine Learning Databases, and is selected from a larger data set held by the National Institutes of Diabetes and Digestive and Kidney Diseases. All patients in this database are Pima Indian women at least 21 years old and living near Phoenix, Arizona. The response Y is the indicator of a positive test for diabetes. Independent variables from this data set include: NumPreg, the number of pregnancies; DBP, diastolic blood pressure (mmHg); DPF, diabetes pedigree function; PGC, the plasma glucose concentration after two hours in an oral glucose tolerance test; BMI, body mass index [weight in kg/(height in m)2]; and AGE (years). There are in total 724 complete observations in this data set.
In this example, we explore the impact of these covariates on the probability of a positive test. We first fit the data set using a linear logistic regression model: the estimated results are listed in the left panel of Table 2. These results indicate that NumPreg, DPF, PGC and BMI are statistically significant, while DBP and AGE are not statistically significant.
TABLE 2
Results for the Pima study. Left panel: estimated values, associated standard errors and P -values by using GLM. Right panel: Estimates, associated standard errors using the GAPLM with the proposed variable selection procedures
GLM | GAPLM | ||||||
---|---|---|---|---|---|---|---|
Est. | s.e. | z value | Pr(>|z|) | SCAD (s.e.) | LASSO (s.e.) | BIC (s.e.) | |
NumPreg | 0.118 | 0.033 | 3.527 | 0 | 0 (0) | 0.021 (0.019) | 0 (0) |
DBP | −0.009 | 0.009 | −1.035 | 0.301 | 0 (0) | −0.006 (0.005) | 0 (0) |
DPF | 0.961 | 0.306 | 3.135 | 0.002 | 0.958 (0.312) | 0.813 (0.262) | 0.958 (0.312) |
PGC | 0.035 | 0.004 | 9.763 | 0 | 0.036 (0.004) | 0.034 (0.003) | 0.036 (0.004) |
BMI | 0.091 | 0.016 | 5.777 | 0 | |||
AGE | 0.017 | 0.01 | 1.723 | 0.085 |
However, a closer investigation shows that the effect of AGE and BMI on the logit transformation of the probability of a positive test may be nonlinear, see Figure 3. Thus, we employ the following GAPLM for this data analysis,
logit{P(Y = 1)} = η0 + β1NumPreg + β2DBP + β3DPF + β4PGC + η1(BMI) + η2(AGE).
The patterns of the nonparametric functions of BMI and Age (solid lines) with ± s.e. (shaded areas) using the R function, gam, for the Pima study.
Using B-splines to approximate η1(BMI) and η2(AGE), we adopt 5-fold cross-validation to select knots and find that the approximation with no internal knots performs well for the both nonparametric components.
We applied the proposed variable selection procedures to the model (15), and the estimated coefficients and their standard errors are listed in the right panel of Table 2. Both SCAD and BIC suggest that DPF and PGC enter the model, whereas NumPreg and DBP are suggested not to enter. However, the LASSO suggests an inclusion of NumPreg and DBP. This may be because LASSO admits many variables in general, as we observed in the simulation studies. The nonparametric estimators of η1(BMI) and η2(AGE), which are obtained by using the SCAD-based procedure, are similar to the solid lines in Figure 3. It is worth pointing that the effect of AGE on the probability of a positive test shows a concave pattern, and women whose age is around 50 have the highest probability of developing diabetes. Importantly, the linear logistic regression model does not reveal this significant effect.
It is interesting that the variable NumPreg is statistically insignificant when we fit the data using GAPLM with the proposed variable selection procedure, but shows a statistically significant impact when we use GLM. One might reasonably conjecture that this phenomenon might be due to model misspecification. To test this, we conducted a simulation as follows. We generated the response variables using the estimates and functions obtained by GAPLM with the SCAD. Then we fit a GLM for the generated data set. We repeated the generation and fitting procedures 5,000 times and found that NumPreg is identified positively significant 67.42% percent of the time at level 0.05 in the GLMs. For DBP, DPF, PGC, BMI and AGE, the percentages that they are identified as statistically significant at the level 0.05 are 4.52%, 90.36%, 100% and 99.98% and 56.58%, respectively. This means that NumPreg can incorrectly enter the model, with more than 65% probability, when a wrong model is used, while DBP, DPF, PGC, BMI and AGE seem correctly to be classified as insignificant and significant covariates even with this wrong GLM model.
6. Concluding remarks
We have proposed an effective polynomial spline technique for the GAPLM, then developed variable selection procedures to identify which linear predictors should be included in the final model fitting. The contributions we made to the existing literature can be summarized in three ways: (i) the procedures are computationally efficient, theoretically reliable, and intuitively appealing; (ii) the estimators of the linear components, which are often of primary interest, are asymptotically normal; and (iii) the variable selection procedure for the linear components has an asymptotic oracle property. We believe that our approach can be extended to the case of longitudinal data [Lin and Carroll (2006)], although the technical details are by no means straightforward.
An important question in using GAPLM in practice is which covariates should be included in the linear component. We suggest proceeding as follows. The continuous covariates are put in the nonparametric part and the discrete covariates in the parametric part. If the estimation results show that some of the continuous covariate effects can be described by certain parametric forms such as a linear form, either by formal testing or by visualization, then a new model can be fit with those continuous covariate effects moved to the parametric part. The procedure can be iterated several times if needed. In this way, one can take full advantage of the flexible exploratory analysis provided by the proposed method. However, developing a more efficient and automatic criterion warrants future study. It is worth pointing out the proposed procedure may be instable for high-dimensional data, and may encounter collinear problems. Addressing these challenging questions is part of ongoing work.
Acknowledgments
The authors would like to thank the Co-Editor, Professor Tony Cai, an Associate Editor and two referees for their constructive comments that greatly improved an earlier version of this paper.
APPENDIX
Throughout the article, let ‖ · ‖ be the Euclidean norm and ‖φ‖∞ = supm |φ(m)| be the supremum norm of a function φ on [0, 1]. For any matrix A, denote its L2 norm as ‖A‖2 = sup‖x‖≠0 ‖Ax‖/‖x‖, the largest eigenvalue.
A.1. Technical lemmas
In the following, let ℱ be a class of measurable functions. For probability measure Q, the L2(Q)-norm of a function f ∈ ℱ is defined by (∫|f|2dQ)1/2. According to van der Vaart and Wellner (1996), the δ-covering number ?(δ, ℱ, L2(Q)) is the smallest value of ? for which there exist functions f1, …, f?, such that for each f ∈ ℱ, ‖f − fj‖ ≤ δ for some j ∈ {1, …, ?}. The δ-covering number with bracketing ?[·](δ, ℱ, L2(Q)) is the smallest value of ? for which there exist pairs of functions , such that for each f ∈ ℱ, there is a j ∈ {1, …, ?} such that . The δ-entropy with bracketing is defined as log ?[·](δ, ℱ, L2(Q)). Denote . Let Qn be the empirical measure of Q. Denote for any measurable class of functions ℱ.
We state several preliminary lemmas first, whose proofs are included in the supplemental materials. Lemmas 1–3 will be used to prove the remaining lemmas and the main results. Lemmas 4 and 5 are used to prove Theorems 1–3.
Lemma 1 [Lemma 3.4.2 of van der Vaart and Wellner (1996)]. Let M0be a finite positive constant. Let ℱ be a uniformly bounded class of measurable functions such that Qf2 < δ2and ‖f‖∞ < M0. Then
where C0is a finite constant independent of n.
Lemma 2 [Lemma A.2 of Huang (1999)]. For any δ > 0, let
Θn = {η(x) + zTβ; ‖β − β0‖ ≤ δ, η ∈ ?n, ‖η−η0‖2 ≤ δ}.
Then, for any ε ≤ δ, log ?[·](δ, Θn, L2(P)) ≤ cNn log (δ/ε).
For simplicity, let
(16)
Lemma 3. Under conditions (C1)–(C5), for the above random matrixWn, there exists a positive constant C such that, a.s.
According to a result of de Boor [(2001), page 149], for any function g ∈ ℋ(p) with p < r − 1, there exists a function , such that , where C is some fixed positive constant. For η0 satisfying (C1), we can find = {j,k, j = 1, …, Nn, k = 1, …, d1}T and an additive spline function = TB(x) ∈ ?n, such that
(17)
Let
(18)
In the following, let . Further let
Lemma 4. Under conditions (C1)–(C5), , whereis in(18), .
In the following, denote = (T, T)T, = (T, T)T and
(19)
Lemma 5. Under conditions (C1)–(C5),
A.2. Proof of Theorem 1
According to Lemma 5,
thus and
By Lemma 1 of Stone (1985), , for each 1 ≤ k ≤ d1. Equation (17) implies that . Then
Similarly,
and , for any k = 1, …, d1.
A.3. Proof of Theorem 2
We first verify that
(20)
(21)
where is defined in (6).
Define
(22)
Noting that ρ2 is a fixed bounded function under (C7), we have , for l = 1, …, d2. By Lemma 2, the logarithm of the ε-bracketing number of the class of functions
?1(δ) = {ρ2{m(x, z)}{z − Γ(x)}:m ∈ ℳn, ‖m − m0‖ ≤ δ}
is c{Nn log (δ/ε) + log(δ−1)}, so the corresponding entropy integral
According to Lemmas 4 and 5 and Theorem 1, . By Lemma 7 of Stone (1986), we have , thus
(23)
Thus by Lemma 1 and Theorem 1, for ,
where according to condition (C5). By the definition of , for any measurable function ϕ, E[ϕ(X)ρ2{m0(T)}] = 0. Hence (21) holds. Similarly, (20) follows from Lemmas 1 and 5.
According to condition (C6), the projection function , where the theoretically centered function Γk ∈ ℋ(p). By the result of de Boor [(2001), page 149], there exists an empirically centered function , such that , k = 1, …, d1. Denote and clearly add ∈ ?n. For any ν ∈ Rd2, define ν = (x, z) + νT{z − add(x)} = {(x) − νTadd(x)} + ( + ν)Tz ∈ ℳn, where ℳn is given in (22). Note that ν maximizes the function for all m ∈ ℳn when ν = 0, thus . For simplicity, denote i ≡ (Ti), and we have
(24)
For the first term in (24), we get
We decompose II into two terms II1 and II2 as follows:
We next show that
(26)
where . Using an argument similar to the proof of Lemma 5, we have
where K = (INnd1, 0(Nnd1) × d2) and INnd1 is a diagonal matrix. Note that the expectation of the square of the sth column of is
Thus, (26) holds by Markov’s inequality. Based on (21), we have . Using similar arguments and (20) and (21), we can show that
According to (23) and condition (C5), we have
Combining (24) and (25), we have
Note that
Thus the desired distribution of follows.
A.4. Proof of Theorem 3
Let τn = n−1/2 + an. It suffices to show that for any given ζ > 0, there exists a large constant C such that
(27)
Denote
and , where s is the number of components of β10. Note that pλn (0) = 0 and pλn (|β|) ≥ 0 for all β. Thus, ℒP(β0 + τnu) − ℒP(β0) ≤ Un,1 + Un,2. Let . For Un,1, note that
Mimicking the proof for Theorem 2 indicates that
(28)
where the orders of the first term and the second term are , respectively. For Un,2, by a Taylor expansion and the Cauchy–Schwarz inequality, n−1Un,2 is bounded by . As wn → 0, both the first and second terms on the right-hand side of (28) dominate Un,2, by taking C sufficiently large. Hence, (27) holds for sufficiently large C.
A.5. Proof of Theorem 4
The proof of is similar to that of Lemma 3 in . We therefore omit the details and refer to the proof of that lemma.
Let , for MPL in (11), and . Define . For any ν1 ∈ Rs, where s is the dimension of β10, define
Note that for all m ∈ ℳ1n when ν1 = 0. Mimicking the proof for Theorem 2 indicates that
Thus, asymptotic normality follows because
Footnotes
1Supported by NSF Grant DMS-09-05730.
2Supported by a Merck Quantitative Sciences Fellowship Program.
3Supported by NSF Grants DMS-08-06097 and DMS-10-07167.
4Supported by a grant from the National Cancer Institute (CA57030) and by Award Number KUS-CI-016-04, made by King Abdullah University of Science and Technology (KAUST).
SUPPLEMENTARY MATERIAL
Detailed proofs and additional simulation results of: Estimation and variable selection for generalized additive partial linear models (DOI: 10.1214/11-AOS885SUPP;.pdf). The supplemental materials contain detailed proofs and additional simulation results.
Contributor Information
Li Wang, Department of Statistics, University of Georgia, Athens, Georgia 30602, USA, ude.agu@gnawylil.
Xiang Liu, Department of Biostatistics, and Computational Biology, University of Rochester, Rochester, New York 14642, USA, ude.retsehcor.tsb@uilx.
Hua Liang, Department of Biostatistics, and Computational Biology, University of Rochester, Rochester, New York 14642, USA, ude.retsehcor.tsb@gnailh.
Raymond J. Carroll, Department of Statistics, Texas A&M University, College Station, Texas 77843-3143, USA, ude.umat.tats@llorrac..
REFERENCES
- Breiman L. Heuristics of instability and stabilization in model selection. Ann. Statist. 1996;24:2350–2383. MR1425957. [Google Scholar]
- Buja A, Hastie T, Tibshirani R. Linear smoothers and additive models (with discussion) Ann. Statist. 1989;17:453–555. MR0994249. [Google Scholar]
- Carroll RJ, Fan JQ, Gijbels I, Wand MP. Generalized partially linear single-index models. J. Amer. Statist. Assoc. 1997;92:477–489. MR1467842. [Google Scholar]
- Carroll RJ, Maity A, Mammen E, YU K. Nonparametric additive regression for repeatedly measured data. Biometrika. 2009;96:383–398. MR2507150. [Google Scholar]
- Craven P, Wahba G. Smoothing noisy data with spline functions. Numer. Math. 1979;31:377–403. MR0516581. [Google Scholar]
- De Boor C. Applied Mathematical Sciences. Vol. 27. New York: Springer; 2001. A Practical Guide to Splines. revised ed. MR1900298. [Google Scholar]
- Fan J, LI R. Variable selection via nonconcave penalized likelihood and its oracle properties. J. Amer. Statist. Assoc. 2001;96:1348–1360. MR1946581. [Google Scholar]
- Fan J, LI R. Statistical challenges with high dimensionality: Feature selection in knowledge discovery. International Congress of Mathematicians. 2006;Vol. III:595–622. Eur. Math. Soc. Zürich. MR2275698. [Google Scholar]
- Frank IE, Friedman JH. A statistical view of some chemometrics regression tools (with discussion) Technometrics. 1993;35:109–148.[Google Scholar]
- HäRDLE W, Huet S, Mammen E, Sperlich S. Bootstrap inference in semiparametric generalized additive models. Econometric Theory. 2004;20:265–300. MR2044272. [Google Scholar]
- Hastie TJ, Tibshirani RJ. Generalized Additive Models. Monographs on Statistics and Applied Probability. Vol. 43. London: Chapman and Hall; 1990. MR1082147. [Google Scholar]
- Huang J. Functional ANOVA models for generalized regression. J. Multivariate Anal. 1998;67:49–71. MR1659096. [Google Scholar]
- Huang J. Efficient estimation of the partially linear additive Cox model. Ann. Statist. 1999;27:1536–1563. MR1742499. [Google Scholar]
- Hunter D, LI R. Variable selection using MM algorithms. Ann. Statist. 2005;33:1617–1642. MR2166557. [PMC free article] [PubMed] [Google Scholar]
- Li R, Liang H. Variable selection in semiparametric regression modeling. Ann. Statist. 2008;36:261–286. MR2387971. [PMC free article] [PubMed] [Google Scholar]
- Li Y, Ruppert D. On the asymptotics of penalized splines. Biometrika. 2008;95:415–436. MR2521591. [Google Scholar]
- Lin X, Carroll RJ. Semiparametric estimation in general repeated measures problems. J. R. Stat. Soc. Ser. B Stat. Methodol. 2006;68:69–88. MR2212575. [Google Scholar]
- Linton O, Nielsen JP. A kernel method of estimating structured nonparametric regression based on marginal integration. Biometrika. 1995;82:93–101. MR1332841. [Google Scholar]
- Marx BD, Eilers PHC. Direct generalized additive modeling with penalized likelihood. Comput. Statist. Data Anal. 1998;28:193–209.[Google Scholar]
- Mccullagh P, Nelder JA. Monographs on Statistics and Applied Probability. 2nd ed. Vol. 37. London: Chapman and Hall; 1989. Generalized Linear Models. MR0727836. [Google Scholar]
- Nelder JA, Wedderburn RWM. Generalized linear models. J. Roy. Statist. Soc. Ser. A. 1972;135:370–384.[Google Scholar]
- Ruppert D, Wand M, Carroll R. Semiparametric Regression. Cambridge: Cambridge Univ. Press; 2003. MR1998720. [Google Scholar]
- Severini TA, Staniswalis JG. Quasi-likelihood estimation in semiparametric models. J. Amer. Statist. Assoc. 1994;89:501–511. MR1294076. [Google Scholar]
- Smith JW, Everhart JE, Dickson WC, Knowler WC, Johannes RS. Proc. Annu. Symp. Comput. Appl. Med. Care. Washington, DC: IEEE Computer Society Press; 1988. Using the ADAP learning algorithm to forecast the onset of diabetes mellitus; pp. 261–265. [Google Scholar]
- Stone CJ. Additive regression and other nonparametric models. Ann. Statist. 1985;13:689–705. MR0790566. [Google Scholar]
- Stone CJ. The dimensionality reduction principle for generalized additive models. Ann. Statist. 1986;14:590–606. MR0840516. [Google Scholar]
- Stone CJ. The use of polynomial splines and their tensor products in multivariate function estimation. Ann. Statist. 1994;22:118–184. MR1272079. [Google Scholar]
- Tibshirani R. Regression shrinkage and selection via the Lasso. J. R. Stat. Soc. Ser. B Stat. Methodol. 1996;58:267–288. MR1379242. [Google Scholar]
- Van Der Vaart AW, Wellner JA. Weak Convergence and Empirical Processes with Applications to Statistics. New York: Springer; 1996. MR1385671. [Google Scholar]
- Wood SN. Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Assoc. 2004;99:673–686. MR2090902. [Google Scholar]
- Wood SN. Generalized Additive Models. Boca Raton, FL: Chapman & Hall/CRC Press; 2006. MR2206355. [Google Scholar]
- Xue L, Yang L. Additive coefficient modeling via polynomial spline. Statist. Sinica. 2006;16:1423–1446. MR2327498. [Google Scholar]
- Yu K, Lee YK. Efficient semiparametric estimation in generalized partially linear additive models. J. Korean Statist. Soc. 2010;39:299–304. MR2730081. [Google Scholar]
- Yu K, Park BU, MAMMEN E. Smooth backfitting in generalized additive models. Ann. Statist. 2008;36:228–260. MR2387970. [Google Scholar]
- Zou H. The adaptive Lasso and its oracle properties. J. Amer. Statist. Assoc. 2006;101:1418–1429. MR2279469. [Google Scholar]