2.1. Additive Binary Response Models

In an additive binary response model independent and identically distributed (i.i.d.) tuples (Y_i,X_i,T_i) are observed (i = 1,…,n), where T_i = (T_i,1,…,T_i,d) is a random variable with components T_i,j in

, X_i is in

, and Y_i ∈ {0,1}. Conditionally given (X_i,T_i) the variable Y_i is distributed as a Bernoulli variable with parameter G{X_i^T β + α + m₁(T_i,1) + ··· + m_d(T_i,d)} where G is a known (link) function, β an unknown parameter in

, α an unknown scalar, and

unknown smooth functions. For identifiability we set E w_j(T_i,1)m_j(T_i,1) = 0 ∀j for some weight functions w_j. Given (X_i,T_i), the likelihood of Y_i is

where μ_i = G{X_i^Tβ + α + m₁(T_i,1) + ··· + m_d(T_i,d)}. The likelihood function is given by

where m⁺(t) is the additive function α + m₁(t₁) + ··· + m_d(t_d).

We now discuss how the additive nonparametric components can be estimated. Without loss of generality, we will do this for the first component m₁. We write r = q₁, s = q₂ + ··· + q_d and define the smoothed likelihood

where the vector T_i = (T_i,1,…,T_i,d) is a random variable with components T_i,j in

. For a vector u = (u₁,…,u_d) with components u_j in

we denote (u₂,…,u_d)^T by u₋₁; similarly, we write T_i,−1 = (T_i,2,…,T_i,d)^T. For a kernel function L defined on

put L_g(v) = (g₁·····g_s)⁻¹L(g₁⁻¹v₁,…,g_s⁻¹v_s) and for simplicity assume that L is a product kernel

. Similarly, define K_h(v) = h⁻¹K(h⁻¹v) for

and bandwidth vector

with product kernel

. The bandwidth vector g is related to smoothing in direction of the “nuisance” covariates. The relative speed of the elements of g to the elements of h and the choice of these bandwidths will be discussed subsequently.

Following Severini and Staniswalis (1994), we base our estimates on an iterative application of smoothed local and unsmoothed global likelihood functions. We define for β ∈ B

Equation (5) may be written as

. The result

is a multivariate kernel estimate of m⁺ that does not use the additive structure of m⁺. This

will be used in an additional step to obtain estimates

of the additive components α, m₁,…,m_d. The final additive estimate of m⁺(t) will then be given by

.

For the estimation of the nonparametric component m₁ the marginal integration method is applied. It is motivated by the fact that, up to a constant, m₁(t₁) is equal to

for a weight function w₋₁. Note that this method does not use iterations so that the explicit definition allows a detailed asymptotic analysis. A weight function w₋₁ is used for two reasons: it may be useful to avoid problems at the boundary, and it can be chosen to minimize the variance (compare Fan, Härdle, and Mammen, 1998). So we define

which estimates the function m₁ up to a constant. An estimate of the function m₁ is given by norming (with a weight function w₁)

The additive constant α can be estimated by

Again, the weight functions w₀ and w₁ may be useful to avoid problems at the boundary. After having estimated the remaining nonparametric components analogously, the final estimate of m is

2.2. Semiparametric Generalized Additive Models

We now come to the discussion of estimation in semiparametric generalized additive models. Suppose that we observe an independent sample (Y₁,X₁,T₁),…,(Y_n,X_n,T_n) with E [Y_i|X_i,T_i] = G{X_i^Tβ + m(T_i)}. Additional assumptions on the conditional distribution of Y_i will be given subsequently. For a positive function V the quasi-likelihood (QL) function is defined as

where μ is the (conditional) expectation of Y, i.e., μ = G{X^Tβ + m(T)}. The QL function has been introduced for the case that the conditional variance of Y is equal to σ²V(μ) where σ² is an unknown scale parameter. The function Q can be motivated by the following two considerations: clearly, Q(μ; y) is equal to −½(μ − y)²v⁻¹ where v⁻¹ is a weighted average of 1/V(s) for s between μ and y. Consequently, maximum QL estimates can be interpreted as a modification of weighted least squares. Another motivation comes from the fact that for exponential families the maximum QL estimate coincides with the maximum likelihood estimate. Note that the maximum likelihood estimate

, based on an i.i.d. sample Y₁,…,Y_n from an exponential family with mean μ(θ) and variance V{μ(θ)}, is given by

We consider three models.

Model A. (Y₁,X₁,T₁),…,(Y_n,X_n,T_n) is an i.i.d. sample with E [Y_i|X_i,T_i] = G{X_i^Tβ + m(T_i)}.

Model B. Model A holds, and the conditional variance of Y_i is equal to Var[Y_i|X_i,T_i] = σ²V(μ_i) where μ_i = G{X_i^Tβ + m(T_i)} and where σ² is an unknown scale parameter.

Model C. Model A holds, and the conditional distribution of Y_i belongs to an exponential family with mean μ_i and variance V(μ_i) with μ_i as in Model B.

The QL function is well motivated for Models B and C. The more general Model A is included for discussion of robustness issues, i.e., to discuss the case of a wrongly specified conditional variance in Models B and C. If not stated otherwise, all the following remarks and results treat the most general Model A. The QL function and the smoothed QL function are now defined as in (3) and (4) with (2) replaced by (12). The estimates

are defined as in (5)–(10). Asymptotics for

are presented in Section A.2 of the Appendix. In particular, Lemma A2.1 shows that

converges to a centered Gaussian variable where the bias δ_n¹(t₁) is of the form Ah₊² + Bg₊² + o_P(h₊² + g₊²), where h₊ = max_1≤j≤r h_j and g₊ = max_1≤j≤s g_j. For a definition of the terms Ah₊² and Bg₊² see Lemma A2.1. This lemma does not require that g₊ is of smaller order than h₊, an assumption that has been made in previous papers. Clearly, then the bias term Bg₊² would be asymptotically negligible, and therefore asymptotics suggests the choice g₊ = o(h₊). However, stochastic and numerical stability of the preestimator

demand that h₁ ×···× h_r·g₁ ×···× g_s is large. Otherwise too few observations would lie in the local support of the multidimensional kernel. Often, in practice even larger values for g_j than for h_l are needed for a satisfactory performance of

. The constant A in the bias depends on the value of m₁′ and m₁′′ at t₁, whereas the constant B depends on averages of powers of m_j′(t_j) and m_j′′(t_j) over t_j and over j ≠ 1. Typically the averaging leads to small values of B. For more discussion, especially on optimal rates and efficiency, we refer to Härdle, Huet, Mammen, and Sperlich (1998).

The remaining additive components m_j for j = 2,…,d are estimated in analogy to m₁. It can be checked that the estimates

are asymptotically independent. The variance of the estimate

can be consistently estimated (see Section A.2 of the Appendix). Consistency and asymptotic normality of β are shown in Lemma A2.2. It turns out that for asymptotic unbiasedness with rate

no undersmoothing is required in the nonparametric estimation. Further, an explicit expression for the asymptotic variance is given that, however, depends on unknown terms as, e.g., on the function m(·).

3.1. Bias Correction

Lemma A2.1 in the Appendix shows that if the elements of the bandwidth vectors h and g are of the same order, the bias of

depends on the shape of the other additive components m₂,…,m_d. This may lead to wrong interpretations of the estimate

. Bootstrap bias estimates will help to judge such effects.

In all three resampling schemes, one uses the data (X₁,T₁,Y₁*),…, (X_n,T_n,Y_n*) to calculate the estimate

. This is done with the same bandwidth h for the component t₁ and with the same g for the other d − 1 components. The bootstrap estimate of the mean of

is given by

, where E* denotes the conditional expectation given the sample (X₁,T₁,Y₁),…,(X_n,T_n,Y_n). The bias corrected estimate of m₁(t₁) is defined by

The theorem shows that the bias terms of order g₊² are removed by this construction.

THEOREM 3.1. Assume that Model A, Model B, or Model C holds and that the corresponding version of bootstrap is used. Suppose further that Assumptions (A1)–(A11) in the Appendix apply and that assumptions analogous to (A3) and (A4) hold for the estimation of the other additive components m_j for j = 2,…,d (h being always the bandwidth used for the estimated component m_j and g the bandwidth for the nuisance components). Furthermore, suppose that the elements of h and g tend to zero and that nh₁·····h_r g₁²·····g_s²(log n)⁻² tends to infinity. Then it holds that

where again h₊ = max_1≤j≤r h_j and g₊ = max_1≤j≤s g_j.

Bootstrap applications in nonparametric regression often use resampling from a modified estimate of the regression function. Suppose, e.g., that in the third step of the bootstrap algorithm

is replaced by

, where

is defined as

but with bandwidth vector h^O instead of h. Then if h_j^O/h₊ → ∞ (1 ≤ j ≤ r) one can show that the left-hand side of (13) is of order O_p{(h₊^O)⁴ + g₊⁴ + (nh₁^O···h_r^O)^−1/2}, where h₊^O is the maximal element of h^O. For appropriate choices of h^O, e.g., for h^O with (h₊^O)⁴ and (nh₁^O···h_r^O)^−1/2 of the same asymptotic order, this is of smaller order than the right-hand side of (13) with resampling from

.

3.2. Componentwise Hypothesis Testing

Interesting shape characteristics may be visible in plots of estimates of the additive components. The complicated nature of the model, though, makes it difficult to judge the statistical significance of such findings. Hypothesis tests in addition to uniform confidence bands are useful tools to analyze and interpret fitted components. We now discuss tests of the hypothesis that one component is linear:

Extensions to variable selection problems (H₀ : m₁ ≡ 0) or tests of polynomial forms are straightforward; see also the discussion that follows.

Our test is a modification of a general approach by Hastie and Tibshirani (1990). In semiparametric setups they propose to apply likelihood ratio tests and to use χ² approximations for the calculation of critical values. Approximate degrees of freedom are heuristically derived by calculating the expectation of asymptotic expansions of the test statistic under the null hypothesis. Here we propose more accurate distributional approximations. Furthermore, in the definition of the test statistic we correct for the bias of the nonparametric estimate. Our test statistic is asymptotically normal, but the convergence to the normal limit is very slow as mathematical arguments and simulations indicate. Therefore we propose the bootstrap for the calculation of critical values. Bias correction will be used in the test because otherwise it will have a nonnegligible effect on the power. For this reason,

is compared with a bootstrap estimate of its expectation under the hypothesis.

First, we calculate semiparametric estimates for the hypothetic model

Note that the α occurring in the preceding equation can be different from the α defined in Section 2.1 because X_i is replaced by (X_i,T_i,1). Estimation of the parametric components β, α, and γ₁ and of nonparametric components m₂,…,m_d can be done as in Section 2.1. This defines estimates

. Set

Second, for the bootstrap we proceed as follows: generate independent samples (Y₁*,…,Y_n*) (compare Section 3) but now with μ_i replaced by

. Then, using the data (X₁,T₁,Y₁*),…,(X_n,T_n,Y_n*) calculate the estimate

. The bootstrap estimate of the mean of

is given by

, where E* denotes the conditional expectation given the sample (X₁,T₁,Y₁),…,(X_n,T_n,Y_n). Third, we define the test statistic

with

. The weights [G′{…}]²/V(G{…}) in the summation of the test statistic are motivated by likelihood considerations (see Härdle et al., 1998) but could be replaced by some other weights. The test statistic R has an asymptotic normal distribution (see Lemma A3.1 in the Appendix). Mean and variance can be consistently estimated, and thus critical values for the test could be calculated using the normal approximation. But as mentioned before this approximation does not perform well. Again we recommend using bootstrap approximations. The bootstrap estimate of the distribution of R is given by the conditional distribution of the test statistic R*, defined by

The conditional distribution

(given the original data (X₁,T₁,Y₁),…,(X_n,T_n,Y_n)) is our bootstrap estimate of

(on the hypotheses (14)). Here,

denotes the distribution of R. The following theorem states consistency of the bootstrap.

THEOREM 3.2. Assume that Model A, Model B, or Model C holds and that the corresponding version of bootstrap is used. Furthermore suppose that assumptions (A1)–(A11) in the Appendix hold with X_i replaced by (X_i,T_i,1). Then, if additionally, n^1/2h₁·····h_r g₁²·····g_s²(log n)⁻¹ → ∞ and if all elements of h and g are of order o(n^−1/8), on the hypotheses (14), it holds that

where d_K denotes the Kolmogorov distance, which is defined for two probability measures μ and ν (on the real line) as

.

With similar arguments as in Härdle and Mammen (1993) one shows that the test R has nontrivial asymptotic power for deviations from the linear hypothesis of order n^−1/2(h₁·····h_r)^−1/4. This means that the test does not reject alternatives that have a distance of order n^−1/2. However, the test also detects local deviations (of order n^−1/2(h₁·····h_r)^−1/4) that are concentrated on shrinking intervals with length of order h. The test may be compared with overall tests that achieve nontrivial power for deviations of order n^−1/2. Typically, such tests have poorer power performance for deviations that are concentrated on shrinking intervals. For our test, the choice of the bandwidth determines how sensitively the test reacts on local deviations; i.e., for smaller h the test detects deviations that are more locally concentrated but at the cost of a poorer power performance for more global deviations. In particular, as an extreme case one can consider the case of a constant bandwidth h. This case is not covered by our theory. It can be shown that in that case R is an n^−1/2-consistent overall test.

Finally we want to emphasize that the same procedure works for any other linearly parameterized hypothesis

where θ₁,…,θ_q are unknown parameters but f₁,…,f_q are given. Moreover, the results of this section can be extended to tests of other parametric hypotheses on m₁:

where {m_θ : θ ∈ Θ} is a parametric family. This can be done similarly as in Härdle and Mammen (1993). However, this requires an asymptotic study of parametric estimates in the model (1) with parametric specification (17) for m₁.

Using an approach similar to the approach described earlier, one can construct F-type tests on the coefficients β. For testing H₀ : Hβ = c versus H₁ : Hβ ≠ c (with a k × p matrix H of rank k ≤ p and a constant

for a k ≥ 1) a natural test statistic is defined as

, where

is a consistent estimate of the matrix I, defined in Lemma A2.2. A natural estimate of I would be the bootstrap estimate. According to Lemma A2.2, on the hypothesis R_β has a central χ² distribution. This asymptotic result could be used for the approximate calculation of critical values. As before we recommend applying bootstrap. Then R_β will be compared with its bootstrap analog

. For simplicity, the same (bootstrap) covariance estimate has been used in the calculation of R_β and R_β*.

3.3. Testing Separability and Interactions

First note that our estimate of m₁ is robust against nonadditivity of the other components. In fact, in the construction of the estimate it is only used that m(x; t) is of the form

for an arbitrary function m_2,…,d. It is not assumed that the function m_2,…,d is additive, i.e., m_2,…,d(T₂,…,T_d) = m₂(T₂) + ··· + m_d(T_d). Also in the case that m(x; t) is not of the form (18), the estimate

makes sense because then it estimates the average (or marginal) effect of T₁. Nevertheless the hypothesis of additivity is of interest in its own right and an important step in a model choice procedure. Following the idea of Sperlich, Tjøstheim, and Yang (2002), we consider a split of the first covariate T₁ into two components T_1:1 and T_1:2 and consider the hypothesis

For other approaches to test additivity, see also Gozalo and Linton (2001). Estimates of m_1:1 and m_1:2 are constructed by marginal integration:

so that

is an estimate for the first-order interaction of T_1:1 and T_1:2.

For testing hypothesis (19) we proceed similarly as in Section 3.2. We define

where m_1:1,2* is an estimate based on a bootstrap sample. Bootstrap samples are generated as in Section 3.2 but now with

replaced by

The test statistic R_inter has an asymptotic normal distribution (see Lemma A3.2 in the Appendix). The bootstrap estimate of the distribution of R_inter is given by the conditional distribution of the test statistic R_inter*, with

where

is defined as

but now from a bootstrap sample instead of the original sample.

THEOREM 3.3. Under the assumptions of Theorem 3.2, on the hypotheses (19), it holds that

3.4. Testing the Link Function

Härdle, Mammen, and Proenca (2001) introduce a bootstrap test for the null hypothesis of a parametric generalized linear versus a single index model. We extend here their approach to test

with v(T,X,β) = β^TX + α + m₁(T₁) + ··· + m_d(T_d). We recommend a test statistic of the form

where h_L is an additional bandwidth and where

. For further details see also Section 4.

3.5. Uniform Bootstrap Confidence Bands

To construct uniform confidence bands we first define

where

is the estimate of the variance of

, defined in equation (A.2) in the Appendix. In the simulation study in Section 4 we also use a bootstrap estimate of σ₁(t₁). The distribution of S can be estimated by bootstrap as introduced in the beginning of Section 3. This defines the statistic

. In the definition of S* the norming

could be replaced by

. We write

. Here

is an estimate of the variance of

, that is defined similarly as

but that is calculated with a bootstrap resample instead of with the original sample. The first norming helps to save computation time; for the second choice bootstrap theory from other setups suggests higher order accuracy of bootstrap. Nevertheless, both bootstrap procedures can be used to construct valid uniform confidence bands:

THEOREM 3.4. Assume that Model A, Model B, or Model C holds and that the corresponding version of bootstrap is used. Furthermore suppose that assumptions (A1)–(A11) apply, that all elements of h and g are of order o(n^−1/8), and that nh₁·····h_r g₁²·····g_s²(log n)⁻² → ∞. Then it holds that

This gives uniform confidence intervals for m₁(t₁) − δ_n¹(t₁). For confidence bands of m₁ one needs a consistent estimate of δ_n¹(t₁). This could be done by plug-in or by bootstrap. Both approaches require oversmoothing, i.e., choice of a bandwidth vector h^O with h_j^O/h₊ → ∞; see also the discussion after Theorem 3.1. For related discussions in nonparametric estimation see Eubank and Speckman (1993) and Neumann and Polzehl (1998).