Mem. Fac. Sci. Eng. Shimane Univ.
Series B: Mathematical Science
44 (2011), pp. 33–53

KERNEL REGRESSION FOR BINARY RESPONSE DATAHIDENORI OKUMURA
Communicated by Kanta Naito(Received: February 7, 2011)
Abstract. This paper is based on the author’s thesis, “Kernel Regression for
Binary Response Data”. We consider kernel-based estimators with additional
weights of the regression functions in nonparametric binomial and binary regressions. Firstly, in the binomial regression with a single covariate and a fixedcovariate design, we introduce a weighted Nadaraya-Watson estimator and its
bias adjusted estimator discussed in Okumura and Naito [13]. Secondly, in thebinary regression with multiple covariates and a random covariate design, we
propose weighted local linear estimators.

1. Introduction
In various ﬁelds such as pharmacology, toxicology, econometrics etc, researchers
are often interested in understanding the relationship between a binary response
and a single covariate/covariate vector. The formal description of the relationshipis as follows. Suppose that the response Y of a subject, which is encoded as Y = 1 if
it exhibits interest and Y = 0 if it does not, is observed at each covariate/covariate
vector x. The relationship can be expressed as
Pr(Y = 1|x) = p(x),where p(x) is a function of x. The function p(x) is referred to as the dose response
curve in biostatistics. The aim of this paper is to estimate p(x) with respect to x,
which has various applications to practical situations.
This study focuses on the kernel smoothing methods used in nonparametricbinary/binomial regression. The Nadaraya-Watson estimator (NWE) is one of themost popular and simplest estimators, which is referred to as the local constantestimator. Copas [2] discussed the application of NWE in the binary settings witha single covariate (Staniswalis and Cooper [17] and Müller and Schmitt [10]). The
ordinary local linear estimator (OLLE) is superior to NWE in terms of the minimax eﬃciency and the boundary eﬀect when the response variable is continuous ingeneral (Fan [3] and Fan and Gijbels [4]). NWE is bounded on [0, 1] while OLLE is2000 Mathematics Subject Classiﬁcation. 62G08, 62G20.
Key words and phrases. kernel regression, binary response.
33

34

H. OKUMURA

not necessarily bounded. Aragaki and Altman [1] used the truncated LLE boundedin binary settings with a single covariate. To avoid the boundary problem, the local
linear logistic estimator (LLLE) is considered, which is a special case of the local
polynomial quasi-likelihood estimators discussed by Fan et al. [5]. Since LLLE
is deﬁned through the optimal solution of an optimization problem, an iterative
calculation is required to derive the solution in practice. Moreover, in the case of
single covariate settings, Signorini and Jones [16] proposed an NW-type estimator
involving two bandwidths and Hazelton [7] proposed bias reduction for NWE.
We propose the kernel estimators of p(x) with additional weights for binary andbinomial regression data. In most studies on binomial settings, the dimension of
the covariate is assumed to be one. On the other hand, in binary settings, it is
not necessarily so. Hence, we consider the settings with multiple covariates, i.e.,a covariate vector. This paper deals with the binomial regression problem with asingle covariate and the binary regression problem with multiple covariates.
The remainder of this paper is organized as follows. In Section 2, the standard
kernel smoothing methods are summarized from a general viewpoint, and NWE,
OLLE and LLLE are given.
Sections 3 and 4 investigate the use of kernel estimators in nonparametric binomial regression with a single covariate. In Section 3, for the binomial data, aweighted NWE is derived by considering the variance of the proportion of responses
at each covariate. Moreover, the eﬃcient bias-adjusted estimator proposed in Okumura and Naito [13] is reviewed. In addition, an application of the bias-adjusted
estimator to the quantal bioassay, discussed in Okumura and Naito [12], is also
given. In Section 4, two methods for data-driven bandwidth selection for the biasadjusted estimator discussed in Okumura and Naito [14] are described: the rule ofthumb (ROT) method and the plug-in (PI) method. A scale adjustment and aneﬃcient ROT method are discussed to improve the performance of the proposed
method from practical viewpoints. The asymptotic properties of the proposed PI
bandwidth selector are given.
Sections 5 and 6 are concerned with the kernel estimators in the nonparametricbinary regression with multiple covariates. In particular, we will discuss the eﬃciency of weighting for the local linear type estimators (LLTEs). In Section 5, the
LLTEs with weights are proposed, and their asymptotic properties are given. In
Section 6, the ROT method and the PI method for data-driven bandwidth selection are discussed as in Section 3. Finally, we conclude the paper. Proofs of all
theorems are omitted.
2. Regression with binary data
2.1. Binomial regression models. In bioassay, binomial responses are often observed at K levels designed through a single covariate, for example, a drug doseor a toxic dose. The sample can be expressed as independent random variables
Yi (i = 1, . . . , K) having the binomial distribution Bi(Ni , p(xi )) corresponding tocovariates xi , where p is unknown. Then, the binomial data is generally expressed
as B1 = {(xi , Yi , Ni ) : i = 1, . . . , K}. We consider the equispaced ﬁxed design.

KERNEL BINARY REGRESSION

35

However, this assumption is not essential and can be relaxed by the manipulation
described in Müller and Schmitt [10]. Without lack of generality, the covariates xi
are assumed to be
i−1
xi =
,
i = 1, ..., K.
K −1
Then, a regression model to estimate p(x) is given as Ȳi = p(xi ) + ϵ̄i , i = 1, . . . , K,
where Ȳi = Yi /Ni and ϵ̄i = Ȳi − p(xi ). The data set in this setting can be expressed∑ i
as B2 = {(xi , Ȳi ) : i = 1, . . . , K}. If Yi can be expressed as Yi = N
j=1 Yij , where Yij
has the Bernoulli distribution Bi(1, p(xi )) for j = 1, . . . , Ni , then another regression
model can be given as Yij = p(xi ) + ϵi , i = 1, . . . , K, j = 1, . . . , Ni , where ϵi =Yij − p(xi ). Then, B3 = {(xi , Yij ) : i = 1, . . . , K, j = 1, . . . , Ni } is a general
expression of the binary data.
2.2. Binary regression models. We consider the random design of multiple
covariates. In this case, the outcome at each covariate vector is mostly observed as a binary scale. The binary data can be expressed as follows. Let(X1 , Y1 ), . . . , (Xn , Yn ) be independent random vectors, where Xi = (Xi1 , . . . , Xid )T
is a d-dimensional random vector with a probability density function φ(x) and Yi
is a binary random variable. In addition, assume that Pr(Yi = 1|Xi = x) = p(x)for all i, where p(x) is an unknown function of x. The function p(x) is to be estimated. A binary regression model is given as Yi = p(Xi ) + εi , i = 1, ..., n, where
εi = Yi − p(Xi ). Then, B4 = {(Xi , Yi ) : i = 1, . . . , n} is a general expression of thebinary data with multiple covariates.
2.3. Kernel regression approaches. We will introduce the standard kernel regression estimators. These can be derived from a general viewpoint. Assume that
(X1 , Y1 ), . . . , (Xn , Yn ) are observed independently, where Xi is a d-dimensional random/nonrandom vector and Yi has the distribution Bi(1, p(Xi )). The kernel regression approaches are the methods to estimate p(x) for a ﬁxed x by ﬁtting aparametric approximation to the observations (Xi , Yi ) using the weighted sum of
divergence measures. In general, the standard kernel regression approaches can be
described as follows. Let D(s, t) be a measure that speciﬁes the divergence of sand t and Zi = (1, (Xi − x)T )T for i = 1, . . . , n. We consider the following criterion
for estimation:n∑ℓ(β r : g, D) =KH (Xi − x)D(Yi , g −1 (β Tr Zi )),
i=1

where H is a symmetric positive deﬁnite d × d-matrix called the bandwidth matrix,
KH (x) = |H|−1/2 ψ(H−1/2 x),∫
ψ(x) is the nonnegative d-variate kernel function with ψ(x)dx = 1, β r =(β0 , rβ1 , . . . , rβd )T for r = 0, 1 and g −1 (β Tr Zi ) is a parametric approximation top(Xi ) in a neighborhood of x. The kernel estimators are obtained by minimizing (2.1). Let ℓ0 (β0 ; g, D) = ℓ(β 0 ; g, D) and ℓ1 (β; g, D) = ℓ(β 1 ; g, D), where

36

H. OKUMURA

β = (β0 , β1 , . . . , βd )T . The simplest kernel estimator is NWE that is deﬁned as∑nwi Yip̂N W (x; H) = ∑i=1,n
i=1 wiwhere wi = KH (Xi − x). By setting g1 (t) = t and D1 (s, t) = (s − t)2 , we see thatℓ0 (β0 ; g1 , D1 ) =n∑wi (Yi − β0 )2 .

i=1

Then, p̂N W (x; H) is given as β̂0 that minimizes ℓ0 (β0 ; g1 , D1 ). OLLE can be expressed asp̄OLL (x; H) = eT1(d+1) (XTx Kx Xx )−1 XTx Kx Y,
where ek(d+1) is the (d + 1) column vector having 1 in the kth entry and 0 in allother entries,
1 (X1 − x)T..
,Xx =  ...
.
1 (Xn − x)T
Y = (Y1 , . . . , Yn )T and Kx = diag(w1 , . . . , wn ). Setℓ1 (β; g1 , D1 ) =n∑wi (Yi − β T Zi )2 .

i=1

Then, p̄OLL (x; H) is given as β̂0 of β̂ = (β̂0 , . . . , β̂d )T that minimizes ℓ1 (β; g1 , D1 ).
Fan et al. [5] developed local quasi-likelihood estimators (LQLEs) for generalized
linear models with one-parameter exponential families, which can be constituted
naturally for binomial regression models. Their method estimates η(x) locallyunder the assumption that p(x) = g −1 (η(x)). We assume that g −1 (t) is strictlyincreasing for any t satisfying 0 < g −1 (t) < 1 and − log g −1 (t) and − log(1 − g −1 (t))
are convex. Because of its mathematical tractability and eﬃciency, we only consider
the local linear estimator (LLE) derived in their theory, which corresponds to thecase p = 1 and r = 0 in their notations. Set D2 (s, t) = −{s log t+(1−s) log(1−t)},
then,
(2.1)

ℓ1 (β; g, D2 ) = −n∑

wi {Yi log g −1 (β T Zi ) + (1 − Yi ) log(1 − g −1 (β T Zi ))}.

i=1

The LLE η̂LL (x; H) of η(x) is given by β̂0 of β̂ that minimizes this ℓ1 (β; g2 , D2 ).
The local maximum likelihood estimator (LMLE) of p(x) is given asp̂LM L (x; H) = g −1 (η̂LL (x; H)).
Note that β̂ = η̂LL (x; H) cannot be expressed explicitly. This derivation requires
iterative calculations. A method to identify the existence of β̂ is proposed in
Okumura [11]. In particular, if g −1 (t) is the logistic function G(t) = et /(1 + et ),LMLE is referred to as LLLE, and is denoted as p̂LLL (x; H).

KERNEL BINARY REGRESSION

37

3. Nonparametric binomial regression
We will discuss estimators for the binomial regression models with a single covariate described in Section 2.1.
3.1. Nadaraya-Watson type estimators. The NWE of p(x) based on data B3
can be written as∑K
wi Yip̂N W (x; h) = ∑Ki=1
i=1 wi Ni(cf. Lloyd [8]), where wi = ϕh (x) = h−1 ϕ(h−1 (xi − x)), ϕ(x) is a kernel function
and h is the bandwidth. Note that p̂N W (x; h) can be constructed based on thebinomial data B1 . Throughout this paper, we adopt as ϕ(x) a symmetricunivariate∫
probability
densityfunction
whose
support
is
[−1,
1]
with
xϕ(x)dx= 0 and∫ 2
x ϕ(x)dx > 0.
In addition, the NWE of p(x) based on data B2 is given as∑K
wi Ȳip̂M (x; h) = ∑i=1.
K
i=1 wiThis estimator was proposed in Müller [9], and utilized in Staniswalis and Cooper[17] in the multivariate covariate setting. Note that p̂M (x; h) is equal to p̂N W (x; h)if Ni s are equal. From the viewpoint of local polynomial smoothing (Ruppert andWand [15]), p̂N W (x; h) based on B3 is characterized as the minimizer of(3.1)ℓ0 (β0 ; g1 , D1 ) =

K
∑

wi

i=1

Ni
∑

(Yij − β0 )2

j=1

with respect to β0 . Similarly, p̂M (x; h) based on B2 is characterized as the minimizer
of
K
∑
(3.2)
ℓ0 (β0 ; g1 , D1 ) =wi (Ȳi − β0 )2i=1

with respect to β0 . In both (3.1) and (3.2), the diﬀerences in the variations of Yij s
are not considered in the weights wi s. However, it seems natural that the weights
should be related to the variations of the responses. Since Ȳi has the variance
V [Ȳi ] = vi /Ni at each xi , where vi = pi (1 − pi ), our approach for constructing theestimator of p(x) is to use the weighted kernels deﬁned by Ni wi /vi , i = 1, . . . , K,
which provides the criterion
K
∑
Ni wi
i=1

vi

(Ȳi − β0 )2

with respect to β0 . Through the minimization, we have the ideal estimator p∗ (x; h)as∑K
Yi wi /vi∗.
p (x; h) = ∑Ki=1
i=1 Ni wi /vi

38

H. OKUMURA

The phrase ‘ideal ’ means that p∗ (x; h) contains the vi s that are unknown. Therefore, vi must be estimated and a reliable estimator should be built up by substituting vi s with the appropriate estimators v̂i s. To obtain the estimator of vi , weadopt a pilot estimator of pi deﬁned by√
Yi + Ni /2
√
p̂i =
,
Ni + Ni
which is known as the Bayes
estimator with respect to the Beta distribution√
Beta(s, t) with s = t = Ni /2 that has the minimax property. Hence, the estimator of p(x) is given by∑K
Yi wi /v̂ip̂(x; h) = ∑Ki=1,
i=1 Ni wi /v̂iwhere v̂i = p̂i (1 − p̂i ), i = 1, . . . , K. We call p̂(x; h) the weighted NWE (WNWE).
3.2. Theoretical performance. In this section, we develop a theory for the behavior of the proposed estimator and other estimators. In subsection 3.3.1, we
consider their performance in the near exact situation. Subsection 3.3.2 considers
the performance under a situation familiar in nonparametric smoothing.
3.2.1. Near exact performance. It is known that the consistent property of kernel
estimators is guaranteed under the situation where h tends to zero as K and Ni s
go to inﬁnity. However, their consistency is not held under such a near exact
situation where h and K are ﬁxed and Ni s go to inﬁnity. We assess the behaviorof the estimators through their mean squared errors (MSEs) under the near exact
situation. First, under the situation where h, K and Ni s are ﬁxed, the MSEs ofbasic estimators p̂N W (x; h) and p̂M (x; h) can be easily obtained as follows:
[∑
(3.3)

MSE[p̂N W (x; h)] =

]2 ∑
K
K
2
N
p
w
i i i
i=1 Ni vi wi
,
−
p(x)
+
∑
∑i=1
K
2
( K
i=1 Ni wi )i=1 Ni wi

[∑
(3.4)

MSE[p̂M (x; h)] =

]2 ∑
K
K
2
p
w
i
i
i=1
i=1 vi wi /Ni− p(x) + ∑
.
∑K
2
( K
i=1 wii=1 wi )Note that (3.4) and (3.5) can be derived without any asymptotic manipulations.
In this sense, these are exact evaluations. On the other side, one assumption isneeded to obtain the MSE of p̂(x; h).
∑Assumption 1 K and h are ﬁxed, and N = K
i=1 Ni → ∞, Ni /N → 1/K for
i = 1, 2, ..., K.
We explain the reason for making this assumption. In fact, there is often a casewhere the size of K is limited. In such a situation, one is interested in the behaviorof the estimators when the number of subjects increases. Assumption 1 expresses
the minimum situation wherein the MSE of p̂(x; h) is calculable. Thus, we canKERNEL BINARY REGRESSION

39

evaluate the performance of the estimators under the near exact situation. Direct
calculations including Taylor expansion yield the following.
Theorem 3.1. Under Assumption 1, we have
[∑
]2
K
N
w
/(1
−
p
)
i
i=1 i iMSE[p̂(x; h)] =− p(x)
∑K
N
w
/v
i
i
i
i=1
]
{
} [ ∑K
N
w
/(1
−
p
)
G21 G22i
i=1 i i+2 √− p(x)
∑K
N
N
N
w
/v
i
i
i
i=1
∑K
2
2
Ni wi /vi
GH
1
+ 21 + ∑Ki=1++ o( ),
NN( i=1 Ni wi /vi )2 N
where}
1 ∑ Nj wi wj { 31 ∑ Ni Nj wi2 wj
p
−
p
(1
−
3v
)
−
(pi − pj )(1 − 2pi )2 ,ji
i
B 2 i̸=j N vi2 vj
B 3 i,j N 2 vi3 vj√
1 ∑ Nj Ni wi wj√(pi − pj )(1 − 2pi )2 ,= − 2
2
2B i,j N N vi vj
1 ∑ Ni wi wj= G1 +(pi − pj )(1 − 2pi )2 (1 − vi )4B 2 i,j N vi3 vj√
1 ∑ Nj Ni Nk wi wj wk− 3(pi − pj )(1 − 2pi )2 (1 − 2pk )2 .
4B i,j,k
N 2 vi2 vj vk2

G1 =
G21
G22

2 ∑ Ni Nj wi2 wj
H = − 3(pi − pj )(1 − 2pi )B i,j N 2 vi2 vj+

B =1 ∑ Ni Nj Nk wi2 wj wk(1 − 2pi )2 (pi − pj )(pi − pk ),B 4 i,j,k
N 3 vi3 vj vk

K∑
N i wi
i=1

N vi

.

Theorem 3.1 evaluates the behavior of p̂(x; h) among the covariates, and plays a keyrole in deriving the asymptotic properties of the estimators discussed in subsection
3.3.2. The proof of Theorem 3.1 is given in Okumura and Naito [13].
3.2.2. Asymptotic performance. We focus on the estimators discussed above fromthe viewpoint of the estimation of p(x). Then, it is natural to consider a nonparametric smoothing situation such as the smoothing parameter tends to zero as the
sample size grows. Indeed, no estimator including the bandwidth h is endowed withconsistency without such assumptions. The necessary assumptions are as follows.
Assumption 2 h → 0 as K → ∞ and Ni = N1 → ∞ for i = 2, . . . , K in such amanner that Kh3+ε = O(1) and N1 h2−ε = O(1) for some 0 < ε < 1.

40

H. OKUMURA

Assumption 3 The support of the kernel ϕ(x) is [−1, 1] and ϕ(x) has continuousand bounded derivatives of order m for any x in [−1, 1] with∫ 1
ℓ < k or ℓ = k + 1,
 0,
ℓ (k)
k
ℓ!,
ℓ = k,x ϕ (x)dx =
(−1)
 c ,
−1
otherwise,ℓ,kwhere 0 ≤ k ≤ m and the cℓ,k s are non-positive constants.
Note that the symmetric beta densities on the interval [−1, 1] satisfy Assumption
3 (cf. Wand and Jones [18]):
ϕ(x) = ϕB (x; r) = {22r+1 B(r + 1, r + 1)}−1 (1 − x2 )r 1{|x|<1} (x),r = 0, 1, . . . ,where B(·, ·) is the beta function. Particular cases are the uniform, Epanechnikov,
biweight and triweight kernels for r = 0, 1, 2 and 3, respectively.
Assumption 4 The curve p(x) has continuous and bounded derivatives of orderm + 2 for any x in [0,1], and satisﬁes 0 < p(x) < 1 for any x in [0, 1].
Note again that p̂N W (x; h) is equal to p̂M (x; h) under these assumptions, sinceNi = N1 for i = 2, . . . , K. Under Assumptions 2–4, the MSE of p̂N W (x; h) is given
by
(
)
h
6
MSE[p̂N W (x; h)] = AMSE[p̂N W (x; h)] + O
+h ,
K
where
h4 µ2 (ϕ)2 (2) 2 v(x)R(ϕ)
AMSE[p̂N W (x; h)] =p (x) +
,
4
N1 Kh
∫1
∫1
v(x) = p(x)(1 − p(x)), R(ϕ) = −1 ϕ(z)2 dz and µr (ϕ) = −1 z r ϕ(z)dz. Using the
Taylor expansion, we obtain the MSE expression of the ideal estimator p∗ (x; h) as
)
(
h
6
∗
∗
MSE[p (x; h)] = AMSE[p (x; h)] + O
+h ,
K
where
v(x)R(ϕ)(3.5)
AMSE[p∗ (x; h)] = h4 µ2 (ϕ)2 f (x)2 +
N1 Kh
and
{
}
1
2(1 − 2p(x))p(1) (x)2(2)
f (x) =
p (x) −.
2v(x)Next, we focus on the proposed estimator p̂(x; h). We have the following.
Theorem 3.2. Under Assumptions 2–4, we have(

MSE[p̂(x; h)] = AMSE[p̂(x; h)] + O
AMSE[p̂(x; h)] = AMSE[p∗ (x; h)] −

)
h
6
+h ,
K

h2
(1 − 2p(x))µ2 f (x).
N1

KERNEL BINARY REGRESSION

41

The proof of Theorems 3.2 is also given in Okumura and Naito [13]. The termv(x)R(ϕ){N1 Kh}−1 appears in all MSE expressions above, which is the leading
term of the variances of p̂N W (x; h), p∗ (x; h) and p̂(x; h). Thus, we note that theessential diﬀerence in these estimators appears in the bias. We can immediately
obtain the following.
{}
Corollary 3.3. If (1−2p(x)) v(x)p(2) (x) − (1 − 2p(x))p(1) (x)2 ≥ 0 for x ∈ [0, 1],then
MSE[p∗ (x; h)] ≤ MSE[p̂N W (x; h)]by neglecting the terms of O(hK −1 + h6 ) and smaller.
The proposed WNWE p̂(x; h) is constructed by plugging the variance estimators
into p∗ (x; h). However, Theorem 3.2 reveals that there is a diﬀerence of O(h2 N1−1 )between the AMSE expressions of p̂(x; h) and p∗ (x; h). We also understand thatthe eﬀect of the plug-in variance estimators appears as O(h2 N1−1 ) in the sense of
MSE. The order of O(h2 N1−1 ) is important since it dominates the order in the MSEof p̂(x; h) under Assumption 2. Here, we discuss a manipulation that makes theMSE of p̂(x; h) close to that of p∗ (x; h). We have from the calculations presented
in the Appendix of Okumura and Naito [13] thatBias[p̂(x; h)] = h2 µ2 f (x) −1 − 2p(x)+ o(h2 ).
N1Hence, the bias-adjusted version of p̂(x; h) can be obtained as
(
)
2
1
+ 1−
p̂(x; h).
p̃(x; h) =N1
N1
Then, we have the following.
Theorem 3.4. Under Assumptions 2–4, we have∗

MSE[p̃(x; h)] = AMSE[p (x; h)] + O

(

)
h
6
+h .
K

Further, we have the following asymptotic normality of p̃(x; h):Theorem 3.5. Assume that Assumptions 2–4 hold and also that there exists a
constant ρ ≥ 0 such that N1 Kh5 → ρ2 . Then,√
N1 Kh {p̃(x; h) − p(x)} →D N (ρf (x), v(x)R(ϕ)).
It can be easily veriﬁed that p̂N W (x; h) also has asymptotic normality and that its
asymptotic variance equals v(x)R. The same argument for p̂LM L (x; h) can been
seen in Theorem 2 of Fan et al. [5]. These facts along with Theorem 3.5 revealthat p̃(x; h) has asymptotically the same precision for estimating p(x) as p̂N W (x; h)and p̂LM L (x; h), and we again note that p̂LM L (x; h) does not always exist. Hence,
the use of p̃(x; h) is justiﬁed since it does not have a fault on existence and it

42

H. OKUMURA

has the same asymptotic precision. Using Theorem 3.5, we obtain an asymptotic
100(1 − β)% conﬁdence interval for p(x) as
[
]
β
β
−1
−1
(3.6)
p̃(x; h) − Φ (1 − )V (x), p̃(x; h) + Φ (1 − )V (x) ,22√where V (x) = p̃(x; h)(1 − p̃(x; h))R(ϕ)/(N1 Kh) and Φ(z) is the standard normal
distribution function. The bias adjusted estimator p̃(x; h) has a nonzero asymptotic
bias ρf (x)µ2 in Theorem 3.5 if f (x) ̸= 0 and the bandwidth is optimally selected
as h ∼ (N1 K)−1/5 . However, this bias will be small for α = 0.5 since p(2) (x) is close
to 0, that is, f (x) is close to 0, if p(x) can be well approximated by symmetricand sigmoid models such as the probit models or the logistic models. On the other
hand, the bias may also varnish since ρ = 0 for a small h. Hence, we neglect thebias term when the conﬁdence interval is constituted as in Müller and Schmitt [10].
3.3. Quantal Bioassay. Okumura and Naito [12] proposed quantal bioassay usingp̃(x; h). We assume that p(x) is strictly monotone, and here discuss the estimationfor Θα = p−1 (α)(0 < α < 1). The estimator of Θα is similarly deﬁned as the wayof Müller and Schmitt [10]. Put Mα = {x ∈ [0, 1] : p̃(x; h) = α, p̃(1) (x; h) > 0},then the estimator of Θα is deﬁned by Θ̃α = (inf Mα + sup Mα )/2. Note that if α isincreasing, then Θ̃α is also increasing. The proposed estimator Θ̃α has consistent
property.
Theorem 3.6. Let p(x) be strictly monotone. Under Assumptions 2–4,supα:Θα ∈[0,1]|Θ̃α − Θα | → 0 in probability.

The asymptotic normality of Θ̃α is given as follows.
Theorem 3.7. Let p(x) be strictly monotone. Assume that there exists a constantρ ≥ 0 such that KN1 h5 → ρ2 . Then under Assumptions 2–4, for any Θα ∈[h, 1 − h],()√ρf (Θα )µ2 α(1 − α)R(ϕ)
KN1 h(Θ̃α − Θα ) −→ N,.
dp(1) (Θα )p(1) (Θα )2We also neglect the bias term from the same reasons in the construction of the conﬁdence interval of p(x) given by (3.6). Hence an asymptotic 100(1−β)% conﬁdence
interval for Θα can be obtained as


√
√
β
β
−1
−1
Φ (1 − 2 ) α(1 − α)R(ϕ)Φ (1 − 2 ) α(1 − α)R(ϕ)
Θ̃α −.
, Θ̃α +
KN1 h
KN1 hp̃(1) (Θ̃α )p̃(1) (Θ̃α )This conﬁdence interval is incalculable for all α such that Mα = ∅.

KERNEL BINARY REGRESSION

43

4. Bandwidth selection for WNWE
4.1. Bandwidth selection. From (3.5), the integrated mean squared error(IMSE) with a weight function κ(x) of the bias-adjusted estimator p̃(x; h) canbe approximated under Assumptions 1–3 with m ≥ 0 as follows:∫
IMSE[p̃(·; h)] =E[(p̃(x; h) − p(x))2 ]κ(x)dx
∫
∫
R(ϕ)
4
2
2
(4.1)
≈ h µ2 (ϕ)κ(x)f (x) dx +κ(x)v(x)dx.
N1 Kh
In Section 4, we employ κ(x) = 1[δ1 ,1−δ2 ] (x), where 1[δ1 ,1−δ2 ] (x) is the indicator function of the interval [δ1 , δ2 ] for some small positive constants δ1 and δ2 . Therefore,
the optimal bandwidth hopt that minimizes the right hand side of (4.1) is given as( )1/5θ2(4.2)
hopt = C(ϕ)(N1 K)−1/5 ,θ1∫∫where C(ϕ) = {R(ϕ)/(4µ2 (ϕ)2 )}1/5 , θ1 = κ(x)f (x)2 dx and θ2 = κ(x)v(x)dx.
The unknown amounts θ1 and θ2 that are functionals of p must be estimated. First,
in order to construct an estimator of θ1 , we adopt∑K(i)
i=1 Yi ϕg (xi − x)c(i)p (x; g) = ∑K
i=1 Ni ϕg (xi − x)
as a simple convenient estimator of p(i) (x) having the bandwidth g. The NWE for(0) (x; g). Then, a consistent estimator of f (x) is given by
B1 is pc
(2) (x; g)(0) (x; g))pc(1) (x; g)2pc(1 − 2pc
f¯(x; g) =−,2(0) (x; g))(1 − τ −r (pc(0) (x; g)))τK −r (pcK

where r is a positive integer andt ≤ c,
 c,t,
c < t < 1 − c,τc (t) =
1 − c, t ≥ 1 − c.
(0) (x; g))− pc(0) (x; g) = o (K −r ).
Note that it can be easily shown to hold that τK −r (pcpHence, our proposed estimator of θ1 is given by∫θ̄1 (g) =κ(x)f¯(x; g)2 dx.

To select the optimal bandwidth g, we make the following assumption instead ofAssumption 2.
Assumption 5 g → 0 as K → ∞ and Ni = N1 → ∞ for i = 2, . . . , K in such amanner that Kg 2m+2+ε = O(1) and N1 g 1−ε = O(1) for some 0 < ε < 1.
Then, we have the following theorem:

44

H. OKUMURA

Theorem 4.1. Under Assumptions 3–5 with m = 2, we have[]2
∆2 (p)∆3 (p)
2(1) (2) (3) (4)¯(4.3)
MSE[θ1 (g)] ≃ g ∆1 (p, p , p , p , p ) +
,
+
5
N1 Kg
N1 K 2 g 9
where∆1 = ∆1 (p, p(1) , p(2) , p(3) , p(4) )[∫µ4 (ϕ(2) )p(2) (x)p(4) (x) µ2 (ϕ)(1 − 2p(x))(1 − 2v(x))p(1) (x)4 p(2) (x)
=
κ(x)
−
48
v(x)3
2µ2 (ϕ(1) )(1 − 2p(x))2 p(1) (x)3 p(3) (x) µ2 (ϕ)(1 − 2v(x))p(1) (x)2 p(2) (x)2+3v(x)22v(x)2()](1 − 2p(x))p(1) (x) µ4 (ϕ(2) )p(1) (x)p(4) (x) µ3 (ϕ(1) )p(2) (x)p(3) (x)
−
−
dx,
v(x)
24
3
∫
R(ϕ(2) )∆2 = ∆2 (p) =κ(x)v(x)dx,
4∫
R(ϕ(2) ∗ ϕ(2) )κ(x)p(x)2 v(x)dx,∆3 = ∆3 (p) =4∫ϕ ∗ ϕ(x) = ϕ(t)ϕ(x − t)dt.
−

Futhermore, it holds that√
K N1 g 9 {θ¯1 (g) − θ1 − g 2 ∆1 − (N1 Kg 5 )−1 ∆2 } →d N (0, ∆3 ).
In the right hand side of (4.3), the order of the ﬁrst term, which is the asymptotic
squared bias (ASB) of θ¯1 (g), is larger than that of the second term, which is theasymptotic variance. Therefore, the optimal bandwidth gopt that minimizes the
ASB is given by
(
)1/7
∆2
gopt = gopt (∆1 , ∆2 ) = C1 (∆1 )(N1 K)−1/7 ,
|∆1 |
{

where
C1 (∆1 ) =We adopt
θ̄2 =1, ∆1 < 0,
1/7(5/2), ∆1 > 0.

1 ∑∗ Ni(Ȳi − Ȳi2 )K ∗ i Ni − 1

as an∑estimator of θ2 , where K ∗ is the number of xi s which belongs to [δ1 , 1 − δ2 ]and i ∗ designates the summation for those xi s. Then, under Assumptions 3 and4 with n ≥ 0, it can be shown that θ̄2 − θ2 = Op ((N1 K)−1/2 ). Hence, from (4.2), a
selection of hopt can be performed on the basis of the following:
(
)1/5
θ̄2
h̄PI = h̄PI (gopt ) = C(ϕ)(N1 K)−1/5 .
θ̄1 (gopt )KERNEL BINARY REGRESSION45

The relative error of h̄PI has the following property.
Theorem 4.2. Under Assumptions 3–5 with m = 2,
{
}
(
)
h̄PI
1
∆3
−1/7 5/14
−1 2
−5
N1 K
θ (c ∆1 + c ∆2 ) →D N 0,
,
−1+
hopt
5(N1 K)2/7 1
25c9 θ12
where c = (N1 K)1/7 gopt = C1 (∆1 )(∆2 /|∆1 |)1/7 .
From Theorem 4.2, we immediately have
(
)
1
h̄PI
− 1 = OP.
hopt
(N1 K)2/7
4.2. Implementation. The subsection brieﬂy illustrates the proposed practicalalgoritheorem to select the bandwidth from the data. Our devices in the algoritheorem include a scale adjustment pointed out by Yang and Tschernig [19] and the
GSK approach discussed in Grizzle et al. [6] for the ROT step.
4.2.1. Scale adjustment. As pointed out in Yang and Tschernig [19], the bandwidththat minimizes an ASB as described in Section 4.1 is smaller than the optimal
bandwidth that minimizes an asymptotic MSE (AMSE) in ﬁnite samples. Hence,
they recommended using a bandwidth slightly larger than the bandwidth thatminimizes the ASB.
†
We describe an algoritheorem to ditermine the constant ρ∗ > 1 such that gopt=†∗ρ gopt , where gopt is the optimal bandwidth minimizing the AMSE of θ̄1 (g) whichis the right hand side of (4.3). Denote the AMSE[θ̄1 (ρg)] as α(ρ : g, ∆1 , ∆2 , ∆3 ),which can be regarded as a function of ρ. We apply the Newton–Raphson method
to the derivation of the solution minimizing α(ρ : gopt , ∆1 , ∆2 , ∆3 ), where theexistence of the minimizer is guaranteed from the form of the function. The m-th
approximation of ρ is given byρ[m] = ρ[m−1] −α′ (ρ[m−1] : gopt , ∆1 , ∆2 , ∆3 ),α′′ (ρ[m−1] : gopt , ∆1 , ∆2 , ∆3 )where the derivatives are calculated with respect to ρ and ρ[0] = 1. From the fundamental result of the Newton–Raphson method, it follows that ρ∗ = limm→∞ ρ[m] .
†ˆ i (i = 1, 2, 3) are quantities= ρ̂∗ ĝROT , where ĝROT and ∆In practice, we adopt ĝopt
given by using a ROT method in the sequent discussion and ρ̂∗ is the minimizer ofˆ 1, ∆
ˆ 2, ∆
ˆ 3 ). Finally, we have the data-driven bandwidthα(ρ : ĝROT , ∆†h̃†PI = h̄PI (ĝopt).

4.2.2. GSK approach for ROT. The quantities gopt and ∆i (i = 1, 2, 3) are determined by the ROT method using certain parametric estimation. In the parametricestimation for p(x), it will be desirable that estimators of p(x) take values in (0, 1)
for any x and can be derived without iterative calculations. We present here amethod of ROT endowed with such desirable properties. The method we employ
is called the GSK approach given in Grizzle et al. [6], which is based on a generalized least squares method. The method is as follows. We consider the following

46

H. OKUMURA

polynomial logistic regression model: p(x; β) = G(zT β), where z = (1, x, . . . , xr )Tand unknown parameter β = (β0 , β1 , . . . , βr )T . The estimator β̂ is obtained by∑−1T2
minimizing the weighted sum of squared errors K
i=1 Ni v̂i (G (p̂i ) − zi β) , wherep̂i s are the Bayes estimators given in Section 3.1 and xi = (1, xi , . . . , xri )T . We
can explicitly express β̂ as β̂ = (X T Ω̂−1 X)X T Ω̂−1 η, where Z = (z1 , . . . , zK )T ,
Ω̂−1 = diag[N1 v̂1 , . . . , NK v̂K ] and η = (G−1 (p̂1 ), . . . , G−1 (p̂K ))T . Note that theempirical logit G−1 (Ȳi ) is asymptotically normally distributed with mean zTi β andvariance (Ni vi )−1 and is not well deﬁned if Ȳi is 0 or 1 in practice. The ROT estimator of p(x) and its ith derivative p(i) (x) are obtained by p(x; β̂) and p(i) (x; β̂),ˆ i (i = 1, 2, 3) of ∆i (i = 1, 2, 3), rerespectively. We can obtain the estimators ∆spectively, by substituting the ROT estimators deﬁned above to the correspondingˆ 1, ∆
ˆ 2 ).
parts in the deﬁnitions of ∆i (i = 1, 2, 3). Hence, we have ĝROT = gopt (∆4.2.3. Fully ROT method. As a simpler alternative method, the ROT method forhopt can be directly applied. This bandwidth selector is written as(∫)1/7
κ(x)v̂(x; β̂)dxhROT = C(ϕ) ∫(N1 K)−1/7 ,2ˆ
κ(x)f (x; β̂) dx
where v̂(x; β̂) and fˆ(x; β̂) are the ROT estimators of v(x) and f (x), respectively.
Note that hROT is not consistent in general.
5. Nonparametric binary regression
5.1. Kernel estimators. We consider the binary regression problem for the binary data B4 . Let us introduce the standard kernel-based estimators of the re−1
gression function p(x). Set wi = Kh (Xi − x) = Πdj=1 h−1
j ϕ(hj (Xij − xj )), whereh = (h1 , . . . , hd ). Then, the NWE can be expressed as∑n
Kh (Xi − x)Yi.
p̂N W (x; h) = ∑i=1n
i=1 Kh (Xi − x)
Table 1. Asymptotic conditional biases of three kernel-based estimators.
Estimator

Asymptotic bias{}d∑pj (x)φj (x)12p̂N W (x; h)µ2 (ϕ)
hj pjj (x) + 2
2
φ(x)
j=1
d
∑
1
µ2 (ϕ)
h2j pjj (x)p̆LLT (x; h)
2
j=1
{
}
d
∑
(1 − 2p(x))pj (x)2
1
2
p̂LLL (x; h)µ2 (ϕ)
hj pjj (x) −2p(x)(1 − p(x))
j=1

KERNEL BINARY REGRESSION

47

In order to derive local linear type estimators (LLTEs), we employ the following
criterion:
n
∑
(5.1)
qω (β; x) =ω(Xi )Kh (Xi − x)(Yi − β T Zi )2 ,i=1
T

where β = (β0 , β1 , . . . , βd ) and ω(t) is a weight function of t. Let LLTE be deﬁned
as β0 of β = (β0 , β1 , . . . , βd )T that minimizes qω (β; x). LLTE can be expressed as
the following matrix formula:p̄LLT (x; h) = eT1(d+1) (XTx Wx Xx )−1 XTx Wx Y,
where Wx = ΩKx , in which Ω = diag(ω(X1 ), . . . , ω(Xn )) and Kx = diag(Kh (X1 −x), . . . , Kh (Xn − x)). Note that p̄LLT (x; h) depends on ω. The truncated versionof LLTE proposed by Aragaki and Altman [1] can be expressed as p̆LLT (x; h) =
τ0 (p̄LLT (x; h)). In particular, OLLE can be deﬁned by p̆LLT (x; h) with ω(t) = 1
for any t, which is expressed as p̂OLL (x; h) = τ0 (p̄OLL (x; h)), where p̄OLL (x; h) =eT1(d+1) (XTx Kx Xx )−1 XTx Kx Y.
The LLLE of p(x) is deﬁned by p̂LLL (x; h) = G(eT1(d+1) β̂) if the following solution
exists uniquely:β̂ = arg min
β

n
∑

{
}2
Kh (Xi − x) Yi log(G(β T Zi )) + (1 − Yi ) log(1 − G(β T Zi )) .

i=1

The following asymptotic properties of p̂N W (x; h), p̆LLT (x; h) and p̂LLL (x; h) canbe obtained through the standard calculations used in kernel smoothing. The
asymptotic conditional variances of the three estimators given Ẋ = (X1 , . . . , Xn )are equal and are given as(5.2)

R(ϕ)d v(x),nh1 · · · hd φ(x)where v(x) = V [Yi |Xi = x] = p(x)(1 − p(x)). The asymptotic conditional biases of
the three estimators given Ẋ are shown in Table 1, where pj (x) and pjj (x) denotethe ﬁrst and second order partial derivatives of p(x) with respect to the jth variable
xj , respectively. The asymptotic conditional biases of LLTEs given Ẋ dependingon ω obtained using (5.1) are found to be the same. If φ(x) is not the uniform
density, then the bias of NWE only includes functionals of φ. The bias of LLLE
is more complex than that of LLTEs. If p(x) is a logistic linear model, then theasymptotic bias of LLLE becomes 0.
5.2. Weighted local linear estimators. We will focus on LLTE because it has
several advantages. In ordinary parametric linear regression with heteroskedasticity, the advantage of weighting by inverse variance is given by the Gauss-Markov
theorem. To improve the inverse variance weighting in local linear regression, we
consider p̆LLT (x; h) weighted through ω(x) = ωα (x), whereωα (t) = ηα (p(t)) = {p(t)α (1 − p(t))1−α }−248

H. OKUMURA

for some α(0 < α < 1), which we refer to as the weighted local linear estimator(WLLE) for 100α%. The corresponding estimator is given by p̆ωα (x; h). If α =0.5, then ω(t) = v(t)−1 . The diﬀerence between the variances of p̂OLL (x; h) andp̆ωα (x; h) is asymptotically given by the following theorem.
Theorem 5.1. Under Assumptions 6–8 in Appendix, we have()dν(ϕ) ∑ 2∥h∥2
V [p̂OLL (x; h)|Ẋ] = V [p̆ωα (x; h)|Ẋ] +h δαj (x) + op,
nh1 · · · hd j=1 j
nh1 · · · hd√
where ∥h∥ = hT h, ν(ϕ) = R(ϕ)d−1 {µ2 (ϕ2 ) − R(ϕ)µ2 (ϕ)} andδαj (x) =2
φ(x)2 p(x)2(1−α) (1 − p(x))2α× [(α − p(x))p(x)(1 − p(x)) {−2φj (x)pj (x) + φ(x)pjj (x)}{}]+ −α(1 − α) + (α − p(x))2 φ(x)pj (x)2 .

It can be shown that µ2 (ϕ2 ) − R(ϕ)µ2 (ϕ) < 0 when ϕ(x) = ϕB (x, r) r = 1, 2, ....
∩
Hence, ν(ϕ) < 0. Let Aα = dj=1 {x : δαj (x) ≤ 0}. For x ∈ Aα , it can be expectedthat the variance of p̆ωα (x; h) is smaller than that of p̂OLL (x; h). At least, for x
with α = p(x), δαj (x) ≤ 0 for any j.
Moreover, if φ is uniform on [0, 1]d and the value of p(x) is close to α and pjj (x)is rather small, that is, p(x) is rather ﬂat, p̆ωα (x; h) can be superior to p̂OLL (x; h)in terms of variance.
Furthermore, if ϕ(x) = ϕB (x; 0), then µ2 (ϕ2 )−R(ϕ)µ2 (ϕ) = 0. Note that ϕB (x; 0)
is a uniform kernel. In addition, if α = 0.5, then we obtain the following exact
result:V [p̂OLL (x; h)|Ẋ] ≥ V [p̆ω0.5 (x; h)|Ẋ].
In practice, p̆ωα (x; h) requires the estimation of Ω = Ωα , where Ωα = diag(ωα (X1 ), . . . , ωα (X1 )). Let p̂LLT (x; h) denote p̆LLT (x; h) in which Ω is replaced byΩ̂, where Ω̂ = diag(ω̂(X1 ), . . . , ω̂(Xn )) is an estimator of Ω. Then, we obtain theconditional AMSE of p̂LLT (x; h) given Ẋ in general as follows.
Theorem 5.2. Under Assumptions 6–9 in Appendix, if E[n−1 tr{Kx (Ω̂−Ω)2 }|Ẋ] =Op (||h||4 ), we have}2
{ d
R(ϕ)d v(x)µ2 (ϕ)2 ∑ 2.
AMSE[p̂LLT (x; h)|Ẋ] =hj pjj (x) +4
nh1 · · · hd φ(x)j=1
Let Ω̂α = diag(ω̂α (X1 ), . . . , ω̂α (Xn )), where ω̂α (Xi ) = ηα (τn−1 (p̂OLL (Xi ; h))).
Then, it follows that E[n−1 tr{Kx (Ω̂α − Ωα )2 }|Ẋ] = Op (||h||4 ). Practical WLLE
for 100α% can be expressed as p̂ωα (x; h) = τ0 (p̄ωα (x; h)), wherep̄ωα (x; h) = eT1 (XTx Ŵx Xx )−1 XTx Ŵx Y

KERNEL BINARY REGRESSION

49

and Ŵx = Ω̂α Kx . On the other hand, assume that p̃(x) is a parametric estimator of p(x) and p̃∗ (x) is the best parametric estimator of p̃(x). Let Ω̃α =diag(ω̃α (X1 ), . . . .ω̃α (Xn )) and Ω̃∗α = diag(ω̃α∗ (X1 ), . . . .ω̃α∗ (Xn )), where ω̃α (Xi ) =ηα (p̃(Xi )) and ω̃α∗ (Xi ) = ηα (p̃∗ ( Xi )). If the true function p(x) is not included in a
family of parametric models, then E[n−1 tr{Kx (Ω̃α − Ωα )2 }|Ẋ] does not converge
to 0 in general but it will follow that E[n−1 tr{Kx (Ω̃α − Ω̃∗α )2 }|Ẋ] = Op (n−1 ) =op (||h||4 ) under certain regularity conditions. Hence, if the parametric ﬁtting ispossible, the parametric estimator Ω̃α can be used as an alternative for Ω̂α . We
obtain an alternative estimator for p̂ωα (x; h) by replacing Ω̂α with Ω̃α , and denote
it as p̆ωα (x; h). Moreover, we adopt p̂ω (x; h) and p̃ω (x; h) as the estimators of p(x)for any x, where p̂ω (x; h) is p̂ωα (x; h) with α = τn−1 (p̂OLL (x; h)) and p̃ω (x; h) isp̃ωα (x; h) with α = p̃(x). It can also be shown that the condition of Theorem 5.2
is satisﬁed for p̂ω (x; h) and p̃ω (x; h).
6. Bandwidth selection for WLLEs
6.1. Preliminaries. Using a known weight function κ(x), the weighted IMSE ofp̂LLT (·; h) is deﬁned as∫
IMSE[p̂LLT (·; h)] = E[ {(p̂LLT (x; h) − p(x)}2 κ(x)φ(x)dx].
From the result of Theorem 5.2, the asymptotic IMSE of p̂LLT (·; h) is given asµ2 (ϕ)2 ∑ ∑ 2 2
R(ϕ)d
(6.1)
AIMSE[p̂LLT (·; h)] =hα hλ Cαλ +B,
4 α=1 λ=1nh1 · · · hd∫∫
where Cαλ = κ(x)pαα (x)pλλ (x)dx and B = κ(x)v(x)dx. Under some regularity
assumptions, given in the Appendix of Yang and Tschernig [19], it is conﬁrmedthat the optimal bandwidth exists uniquely. In such a case, the solution can becalculated numerically. We will deﬁne the optimal bandwidth as follows:ddhopt = arg min AIMSE[p̂LLT (·; h)].
h

We will present the ROT and PI methods for bandwidth selection using (6.1),
which are constructed as in Yang and Tschernig [19]. In practice, Cαλ and B are
to be estimated. For each α, λ = 1, . . . , d, let1∑= C̄αλ (pαα , pλλ ) =pαα (Xi )pλλ (Xi )κ(Xi ),
n i=1
n

C̄αλ

1 ∑ (Yi − p(Xi ))2 κ(Xi )n i=1
φ(Xi )nB̄ = B̄(p, φ) =
and

R(ϕ)d
µ2 (ϕ)2 ∑ ∑ 2 2B̄.
hα hλ C̄αλ +4 α=1 λ=1nh1 · · · hdd

AIM SE[h; {C̄αλ (pαα , pλλ )}, B̄(p, φ)] =d

50

H. OKUMURA

To derive the ROT and PI bandwidth selectors, C̄αλ and B̄ are estimated in parametric and nonparametric approaches, respectively.
6.2. ROT method. To estimate p(x) and its higher order partial derivatives,
Yang and Tschernig [19] have proposed blocked polynomial ﬁtting. A polynomial
model is ﬁtted to each block constructed by separating the covariate space. Mallows’s Cp is used to select the optimal separation for the various separated blocks.
Since p(x) is known to be bounded on (0, 1), we will carry out the blocked ﬁttingof the following logistic quadratic model:
(
)
d ∑
2
∑
p(x; γ) = G γ0 +γ2(j−1)+k xkj .
j=1 k=1

The maximum likelihood estimation of the model in each block is required. We
will select the optimal separation of blocks by the Bayesian information criterion.
The maximum likelihood estimator (MLE) of p(x) obtained by this approach isdenoted as p̃(x). As an estimator of φ(x), we adopt the uniform density function,
which is denoted as φ̇(x). Therefore, the ROT bandwidth of h is given byĥROT = arg min AIMSE[h; {C̄αλ (p̃αα , p̃λλ )}, B̄(p̃, φ̇)].
h6.3. PI method. We will use the estimator with a scalar bandwidth proposed by
Yang and Tschernig [19] as an estimator of pαα (x). The estimator of pαα (x) canbe expressed as p̂αα (x; g) = 2eT(α+1)(5d−1) (XTα ΞXα )−1 XTα ΞY, where g = (g, ..., g),Ξ = diag(Kg (X1 − x), . . . , Kg (Xn − x)), Xα is the n × (5d − 1) matrix deﬁned by
Xα = [1, {(Xij − xj )2 }1≤j≤d , {(Xij − xj )(Xiα − xα )}j̸=α , {(Xij − xj )}1≤j≤d ,{(Xij − xj )(Xiα − xα )2 }j̸=α , {(Xij − xj )2 (Xiα − xα )}j̸=α , (Xiα − xα )3 ]ni=1∑
and [(ai1 , . . . , aim )]ni=1 = (aij )n×m . Put Ĉαλ (gαλ ) = n−1 ni=1 p̂αα (Xi ; gαλ )p̂λλ (Xi ; gαλ ) for each α, λ = 1, . . . , d. Then, from the result of Yang and Tschernig [19], an expression of the asymptotic MSE of Ĉαλ (gαλ ) can be expressed as
}2
{
σ2
BVαλ
+ 2 αλ
,(6.2)
AMSE[Ĉαλ (gαλ )] = gαλ Dαλ + d+4
n g d+2ngαλin which∫Dαλ =4R(ϕ)d−2{δαλ µ4 (ϕ2 )R(ϕ) + (1 − δαλ )µ2 (ϕ2 )2(µ4 (ϕ) − µ2 (ϕ)2 )2−2µ2 (ϕ)µ2 (ϕ2 )R(ϕ) + µ2 (ϕ)2 R(ϕ)2 },∫= cαλ (ϕ) p(x)2 (1 − p(x))2 κ(x)2 dx,

Vαλ =

2
σαλ

{pαα (x)bλλ (x) + pλλ (x)bαα (x)}φ(x)κ(x)dx,

KERNEL BINARY REGRESSION

51

where δαλ = 1 if α = λ, δαλ = 0 if α ̸= λ,µ6 (ϕ) − µ2 (ϕ)µ4 (ϕ)µ2 (ϕ) ∑bαα (x) =p(x)+
pjjαα (x),αααα
12(µ4 (ϕ) − µ2 (ϕ)2 )2 j̸=α∫32
cαλ (ϕ) =Fαα (x)Fλλ (x)dx,(µ4 (ϕ) − µ2 (ϕ)2 )4Fαλ (x) = (1 − δαλ )(ϕ∗ ∗ ϕ)(xα )(ϕ∗ ∗ ϕ)(xλ )Πj̸=α,λ ϕ(2) (xj )+δαλ ϕ∗ (2) (xα )Πj̸=α ϕ(2) (xj ),∫(2)
ϕ (t) =
ϕ(u − t)ϕ(u)du and ϕ∗ (u) = ϕ(u)(u2 − µ2 (ϕ)).
Note that the ﬁrst and second terms in the right hand side of (6.2) are theleading terms of the ASB and the asymptotic variance of Ĉαλ (gαλ ), respectively.
We substituten1∑D̄αλ ={pαα (Xi )bλλ (Xi ) + pλλ (Xi )bαα (Xi )}κ(Xi )
n i=1
and

1 ∑ p(Xi )2 (1 − p(Xi ))2 κ(Xi )2
n i=1
φ(Xi )n2σ̄αλ=

for Dαλ and σαλ in (6.2), respectively. Instead of (6.2), we will consider
{
}2
2
B̄Vαλ
σ̄αλ
(6.3) AMSE[gαλ ; p, pαα , pλλ , {pjjαα , pjjλλ }, φ = gαλ D̄αλ + d+4+ 2 d+2 .
ngαλ
n gαλ
Yang and Tschernig [19] adopted for each α, λ = 1, . . . , d,ĝαλ = argmingαλ ASB[gαλ ; p̃, p̃αα , p̃λλ , {p̃jjαα , p̃jjλλ }, φ̇],
where

}2
{
B̄Vαλ
ASB[gαλ ; p, pαα , pλλ , {pjjαα , pjjλλ }, φ] = gαλ D̄αλ + d+4 .
ngαλ

We obtain Ĉαλ (ĝαλ ) as an estimator of C̄αλ . Therefore, the PI bandwidth selectorof h can be obtained as follows:ĥPI = arg min AIMSE[h; {Ĉαλ (ĝαλ )}, B(p̄, φ̄)],h∑where p̄(x) = p̂OLL (x; ĥROT ) and φ̄(x) = n−1 ni=1 KĥROT (Xi − x) + n−2 . SinceB(p̄, φ̄)/B = 1 + Op (n−2/(d+4) ), under the regularity assumptions given in Yang
and Tschernig [19], we haveĥPI = hopt {1 + Op (n−2/(d+6) )}.
The bandwidth that minimizes the AMSE given by (6.3) can be selected throughthe method proposed in Okumura and Naito [14]. The proposed bandwidth selectorof gαλ is given asĝαλ = arg min AMSE[gαλ ; p̃, p̃αα , p̃λλ , {p̃jjαα , p̃jjλλ }, φ̇],
gαλ

52

H. OKUMURA

which was used in the simulation study in the thesis.
In practice, the d × d matrix whose (α, λ)-element is Ĉαλ is not always a positive∑deﬁnite. Therefore, we replace Ĉαλ by n−1 ni=1 p̂αα (Xi ; gαα )p̂λλ (Xi ; gλλ ). The
ROT and PI bandwidth selectors for NWE and LLLE can also be derived similarly
based on each conditional bias in Table 1 and the same conditional variance (5.2).
7. Conclusion
We proposed the weighted kernel estimators of the regression function in thebinomial and binary regression problems.
For the binomial regression problem with a single covariate, WNWE was proposed in Section 3. This estimator is NWE with a kernel weighted by the inverse of
the variance estimator at each covariate. The bias-adjusted version of WNWE wasalso obtained through the asymptotic properties. The MSE of the bias-adjusted
WNWE is asymptotically smaller than that of NWE under some conditions. The
simulation study in Okumura and Naito [13] also demonstrated the good performance of the bias-adjusted WNWE. Two eﬃcient data-driven bandwidth selectors
were also proposed: the PI bandwidth selector and the ROT bandwidth selector.
Their eﬃciencies are described in Okumura and Naito [14].
For the binary regression problem with multiple covariates, two practical WLLEs,
which are LLTEs with a kernel weighted at each observed covariate, were proposed
in Section 5. The behavior of the estimators and the bandwidth selectors was
discussed. WLLE with weights constructed by using the true regression function
is asymptotically shown in Theorem 5.1 to have a variance smaller than that of
OLLE. With regard to practical applications, the WLLEs with weights constructedby using the parametric and nonparametric estimators were proposed. The practical WLLEs cannot necessarily guarantee the variance reduction asymptoticallyas given in Theorem 5.1. A simulation study to compare estimators described
was done in the thesis. From the simulation study, we inferred that the proposed
WLLEs reduce the variance and have good performance in terms of IMSE. Further,
it was shown that it is eﬀective to use the ROT bandwidth selector in comparatively small samples if the parametric model ﬁts well. The proposed WLLEs canbe recommended as a nonparametric estimator of the regression function, considering the disadvantages of the boundary bias problem inherent in NWE and thetroubles involved in the calculation of the LLLE. The WLLEs are also useful forthe estimation of eﬀective covariates corresponding to a response probability.
Appendix
The following assumptions are needed in Theorems 5.1 and 5.2.
Assumption 6. At x ∈ supp(φ), φ(x) is second continuously diﬀerentiable.
Assumption 7. At x ∈ supp(φ), p(x) is second continuously diﬀerentiable.
Assumption 8. As n → ∞, hj → 0(j = 1, . . . , d) in such a manner that hk /hj =O(1)(j ̸= k) and ||h||4 /(nh1 . . . hd ) = O(1).

KERNEL BINARY REGRESSION

53

Assumption 9. At x ∈ supp(φ), ω(x) is second continuously diﬀerentiable.
References[1] A. Aragaki and N. S. Altman, Local polynomial regression for binary response, Computer
Science and Statistics: Proceedings of the 29th Symposium on the Interface (1997), 467-472.
[2] J. B. Copas, Plotting p against x, Applied Statistics 32 (1983), 25-31.
[3] J. Fan, Local linear regression smoothers and their minimax eﬃciencies, Ann. Statist. 21(1993), 196-216.
[4] J. Fan and I. Gijbels, Local Polynomial Modeling and its Applications, Chapman and Hall,
London, 1996.
[5] J. Q. Fan, N. E. Heckman and M. P. Wand, Local polynomial kernel regression for generalized
linear models and quasi-likelihood functions, J. Amer. Statist. Assoc. 90 (1995), 141-150.
[6] J. E. Grizzel, C. F. Starmer and G. G. Koch, Analysis of growth and dose-response curves,
Biometrics 25 (1969), 489-504.
[7] M. L. Hazelton, Bias reduction in kernel binary regression, Comput. Stat. Data Anal. 51(2007), 4393-4402.
[8] C. J. Lloyd, Statistical Analysis of Categorical Data, John Wiley and Sons, New York, 1999.
[9] H.-G. Müller, Kernel estimators of zeros and of location and size of extrema of regressionfunctions, Scand. J. Statist. 12 (1985), 221-232.
[10] H.-G. Müller and T. Schmitt, Kernel and probit estimates in quantal bioassay, J. Amer.
Statist. Assoc. 83 (1988), 750-759.
[11] H. Okumura, A note on local likelihood regression for binary response data, Commun. Statist.Simul. Comp. 38 (2009), 1019-1025.
[12] H. Okumura and K. Naito, Kernel smoothing in quantal bioassay, Japanese Journal of Applied Statistics 32 (2003), 127-144.
[13] H. Okumura and K. Naito, The weighted estimators in nonparametric binomial regression,
J. Nonparam. Statist. 16 (2004), 39-62.
[14] H. Okumura and K. Naito, Bandwidth selection for kernel binomial regression, J. Nonparam.
Statist. 18 (2006), 343-356.
[15] D. Ruppert and M. P. Wand, Multivariate locally weighted least squares regression, Ann.
Statist. 22 (1994), 1346-1370.
[16] D. F. Signorini and M. C. Jones, Kernel estimators for univariate binary regression, J. Amer.
Statist. Assoc. 99 (2004), 119-126.
[17] J. G. Staniswalis and V. Cooper, Kernel estimates of dose response, Biometrics 44 (1988),
1103-1119.
[18] M. P. Wand and M. C. Jones, Kernel Smoothing Chapman and Hall, London, 1995.
[19] L. Yang and R. Tschernig, Multivariate bandwidth selection for local linear regression, J. R.
Statist. Soc. B 61 (1999), 793-815.
Hidenori Okumura: Department of Information Science and Business Management, Chugoku Junior College, Niwase, Okayama, 701-0197, Japan
E-mail address: okumura@cjc.ac.jp