The Lagrange approach in the monotone single index model

TL;DR

This paper introduces a Lagrange-based method for estimating parameters in the monotone single index model, avoiding reparametrization, and compares various estimators including score and least squares methods.

Contribution

It proposes a novel Lagrange-type approach to estimate parameters without reparametrization and applies it to several existing estimators in the monotone single index model.

Findings

The Lagrange approach simplifies the estimation process.

The estimators are compared with existing methods like LSE, MRE, and EDR.

Effects of random starting values on the algorithms are analyzed.

Abstract

The finite-dimensional parameters of the monotone single index model are often estimated by minimization of a least squares criterion and reparametrization to deal with the non-unicity. We avoid the reparametrization by using a Lagrange-type method and replace the minimization over the finite-dimensional parameter alpha by a `crossing of zero' criterion at the derivative level. In particular, we consider a simple score estimator (SSE), an efficient score estimator (ESE), and a penalized least squares estimator (PLSE) for which we can apply this method. The SSE and ESE were discussed in Balabdaoui, Groeneboom and Hendrickx (2018}, but the proofs still used reparametrization. Another version of the PLSE was discussed in Kuchibhotla and Patra (2017), where also reparametrization was used. The estimators are compared with the profile least squares estimator (LSE), Han's maximum rank estimator (MRE), the effective dimension reduction estimator (EDR) and a linear least squares estimator, which can be used if the covariates have an elliptically symmetric distribution. We also investigate the effects of random starting values in the search algorithms.

Figures

Tables

Table 1. Table 1 : Simulation, model ( X i ∼ N ​ ( 0 , 1 ) , d = 3 formulae-sequence similar-to subscript 𝑋 𝑖 𝑁 0 1 𝑑 3 X_{i}\sim N(0,1),d=3 ): The mean value ( μ ^ i subscript ^ 𝜇 𝑖 \hat{\mu}_{i} = mean( α ^ i ​ n ) , i = 1 , 2 , 3 \hat{\alpha}_{in}),i=1,2,3 ) and n 𝑛 n times the variance-covariance ( σ ^ i ​ j = n ⋅ \hat{\sigma}_{ij}=n\cdot cov ( α ^ i ​ n , α ^ j ​ n ) subscript ^ 𝛼 𝑖 𝑛 subscript ^ 𝛼 𝑗 𝑛 (\hat{\alpha}_{in},\hat{\alpha}_{jn}) , i , j = 1 , 2 , 3 formulae-sequence 𝑖 𝑗 1 2 3 i,j=1,2,3 ) of the simple score estimate (SSE), the efficient score estimate (ESE), the least squares estimate (LSE), the maximum rank correlation estimate (MRE), the linear estimate (this is H-LFLSE in [ 2 ] ), the Penalized Least Squares Estimate (PLSE), the Penalized Least Squares Estimate as computed by the R package simest (PLSE2) and the Efficient Dimension Reduction estimate as computed by the R package edr (EDR), for different sample sizes n 𝑛 n with X i ∼ N ​ ( 0 , 1 ) similar-to subscript 𝑋 𝑖 𝑁 0 1 X_{i}\sim N(0,1) . The line, preceded by ∞ \infty , gives the asymptotic values. The values are based on 1000 1000 1000 replications.

Method	$n$	${\hat{\mu }}_1$	${\hat{\mu }}_2$	${\hat{\mu }}_3$	${\hat{\sigma }}_{11}$	${\hat{\sigma }}_{22}$	${\hat{\sigma }}_{33}$	${\hat{\sigma }}_{12}$	${\hat{\sigma }}_{13}$	${\hat{\sigma }}_{23}$
SSE	100	0.5735	0.5769	0.5752	0.2434	0.2482	0.2399	-0.1234	-0.1165	-0.1229
	500	0.5771	0.5768	0.5774	0.1493	0.1426	0.1470	-0.0721	-0.0771	-0.0700
	1000	0.5772	0.5776	0.5770	0.1230	0.1150	0.1268	-0.0558	-0.0671	-0.0594
	2000	0.5770	0.5775	0.5775	0.1118	0.1091	0.1004	-0.0600	-0.0516	-0.0491
	$\infty$	0.5774	0.5774	0.5774	0.0741	0.0741	0.0741	-0.0370	-0.0370	-0.0370
ESE	100	0.5756	0.5779	0.5758	0.1061	0.1058	0.1105	-0.0510	-0.0549	-0.0548
	500	0.5773	0.5770	0.5776	0.0491	0.0489	0.0518	-0.0231	-0.0261	-0.0258
	1000	0.5775	0.5773	0.5771	0.0438	0.0406	0.0420	-0.0212	-0.0226	-0.0193
	2000	0.5773	0.5773	0.5774	0.0362	0.0367	0.0327	-0.0201	-0.0161	-0.0166
	$\infty$	0.5774	0.5774	0.5774	0.0247	0.0247	0.0247	-0.0123	-0.0123	-0.0123
LSE	100	0.5751	0.5766	0.5759	0.1627	0.1729	0.1753	-0.0799	-0.0820	-0.0923
	500	0.5773	0.5770	0.5771	0.1048	0.1100	0.1044	-0.0550	-0.0497	-0.0548
	1000	0.5773	0.5776	0.5761	0.0985	0.0884	0.0995	-0.0437	-0.0550	-0.0445
	2000	0.5775	0.5772	0.5772	0.0921	0.0968	0.0930	-0.0474	-0.0436	-0.0494
	$\infty$	0.5774	0.5774	0.5774	?	?	?	?	?	?
MRE	100	0.5683	0.5731	0.5736	0.6614	0.6764	0.6181	-0.3258	-0.2999	-0.3235
	500	0.5760	0.5770	0.5763	0.5208	0.4969	0.5125	-0.2521	-0.2670	-0.2437
	1000	0.5771	0.5776	0.5761	0.4815	0.4512	0.4842	-0.2259	-0.2555	-0.2262
	2000	0.5772	0.5770	0.5773	0.4704	0.4502	0.4170	-0.2506	-0.2186	-0.1996
	$\infty$	0.5774	0.5774	0.5774	0.3583	0.3583	0.3583	-0.1791	-0.1791	-0.1791
linear	100	0.5746	0.5736	0.5724	0.4675	0.4588	0.4644	-0.2309	-0.2323	-0.2243
	500	0.5780	0.5751	0.5762	0.5072	0.5057	0.4966	-0.2635	-0.2733	-0.2445
	1000	0.5771	0.5768	0.5768	0.5182	0.5099	0.5085	-0.2600	-0.2583	-0.2491
	2000	0.5772	0.5770	0.5710	0.5129	0.5288	0.5067	-0.2673	-0.2453	-0.2769
	$\infty$	0.5774	0.5774	0.5774	0.5185	0.5185	0.5185	-0.2593	-0.2593	-0.2593
PLSE	100	0.5765	0.5772	0.5766	0.0657	0.0682	0.0629	-0.0352	-0.0303	-0.0326
	500	0.5775	0.5770	0.5774	0.0368	0.0362	0.0384	-0.0173	-0.0195	-0.0189
	1000	0.5775	0.5774	0.5771	0.0329	0.0314	0.0305	-0.0169	-0.0160	-0.0145
	2000	0.5773	0.5775	0.5774	0.0301	0.0313	0.0282	-0.0166	-0.0135	-0.0147
	$\infty$	0.5774	0.5774	0.5773	0.0276	0.0278	0.0282	-0.0137	-0.0140	-0.0142
PLSE2	100	0.5767	0.5771	0.5766	0.0642	0.0661	0.0614	–0.0342	-0.0301	-0.0314
	500	0.5777	0.5768	0.5773	0.0496	0.0485	0.0505	-0.2635	-0.2733	-0.2445
	1000	0.5776	0.5774	0.5771	0.0609	0.0542	0.0565	-0.0169	-0.0317	-0.0248
	2000	0.5771	0.5777	0.5772	0.0724	0.0713	0.0705	-0.0365	-0.0356	-0.0350
	$\infty$	0.5774	0.5774	0.5774	0.0247	0.0247	0.0247	-0.0123	-0.0123	-0.0123
EDR	100	0.5735	0.5755	0.5783	0.1771	0.1827	0.1836	-0.0815	-0.0882	-0.0993
	500	0.5772	0.5770	0.5775	0.0609	0.0598	0.0642	-0.0293	-0.0328	-0.0315
	1000	0.5773	0.5776	0.5769	0.0515	0.0486	0.0521	-0.0241	-0.0275	-0.0245
	2000	0.5773	0.5774	0.5773	0.0450	0.0439	0.0417	-0.0236	-0.0214	-0.0417
	$\infty$	0.5774	0.5774	0.5774	?	?	?	?	?	?

Equations

$$Y={\psi }_0({\mathbf{\alpha }}_0^T\mathbf{X})+\epsilon ,$$

$$\psi \mapsto n^{-1}\sum _{i=1}^n{\left\{Y_i-\psi ({\mathbf{\alpha }}^T{\mathbf{X}}_i)\right\}}^2$$

$$\mathbf{\alpha }\mapsto n^{-1}\sum _{i=1}^n{\left\{Y_i-{\hat{\psi }}_{n,\mathbf{\alpha }}({\mathbf{\alpha }}^T{\mathbf{X}}_i)\right\}}^2$$

$$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X})-Y_{i}\bigr{\}}\left\{\bm{X}_{i}\frac{d}{du}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}+\frac{\partial}{\partial\bm{\alpha}}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}\right\}+\lambda\bm{\alpha}=\bm{0},$$

$$\displaystyle\lambda=-\frac{1}{n}\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X}_{i})-Y_{i}\bigr{\}}\bm{\alpha}^{T}\left\{\bm{X}_{i}\frac{d}{du}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}+\frac{\partial}{\partial\bm{\alpha}}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}\right\},$$

$$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X})-Y_{i}\bigr{\}}^{2}+\lambda\left\{\|\bm{\alpha}\|^{2}_{2}-1\right\},$$

$$\displaystyle\frac{1}{n}\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right)\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X})-Y_{i}\bigr{\}}\left\{\bm{X}_{i}\frac{d}{du}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}+\frac{\partial}{\partial\bm{\alpha}}\hat{\psi}_{n,\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}\right\}=\bm{0}.$$

$$\displaystyle\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right){\mathbb{E}}\left[\bigl{\{}\psi_{\bm{\alpha}}(\bm{\alpha}^{T}\bm{X})-Y\bigr{\}}\left\{\bm{X}\frac{d}{du}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}+\frac{\partial}{\partial\bm{\alpha}}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}\right\}\right]$$

$$\displaystyle=\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right){\mathbb{E}}\bigl{\{}\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})|\bm{\alpha}^{T}X)-\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})\bigr{\}}\bm{X}\frac{d}{du}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}=\bm{0},$$

$$\psi \mapsto \int {\left\{{\psi }_0({\mathbf{\alpha }}_0^T\mathbf{x})-\psi ({\mathbf{\alpha }}^T\mathbf{x})\right\}}^{2}dG(\mathbf{x}),$$

$$\displaystyle{\mathbb{E}}\bigl{\{}\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})|\bm{\alpha}^{T}X)-\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})\bigr{\}}\frac{\partial}{\partial\bm{\alpha}}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}$$

$$\displaystyle={\mathbb{E}}\left[\frac{\partial}{\partial\bm{\alpha}}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}{\mathbb{E}}\left(\left\{\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})|\bm{\alpha}^{T}X)-\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})\right\}\Bigm{|}\bm{\alpha}^{T}X\right)\right]$$

$$=0.$$

$$\displaystyle\frac{1}{n}\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right)\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X}_{i})-Y_{i}\bigr{\}}\bm{X}_{i}=\bm{0},$$

$${\tilde{\psi }}_{n,h,\mathbf{\alpha }}^{\prime }(u)=\frac{1}{h}\int K\left(\frac{u-x}{h}\right)d{\hat{\psi }}_{n,\mathbf{\alpha }}(x),$$

$$\displaystyle\frac{1}{n}\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right)\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{X}_{i})-Y_{i}\bigr{\}}\bm{X}_{i}\tilde{\psi}_{n,h,\bm{\alpha}}^{\prime}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}=\bm{0},$$

$$\sum _{i=1}^n{\left\{f({\mathbf{\alpha }}^T{\mathbf{X}}_i)-Y_i\right\}}^2+\mu {\int }_a^b{f}^{\prime \prime }(x{)}^{2}dx,$$

$$a=\underset{i}{\min }{\mathbf{\alpha }}^T{\mathbf{X}}_i,b=\underset{i}{\max }{\mathbf{\alpha }}^T{\mathbf{X}}_i,$$

$$\mathbf{\beta }\mapsto \sum _{i=1}^n\{Y_i-{\hat{\psi }}_{S(\mathbf{\beta }),\lambda }(S(\mathbf{\beta }{)}^T{\mathbf{X}}_i){\}}^2,$$

$$\displaystyle\frac{1}{n}\left(\bm{I}-\bm{\alpha}\bm{\alpha}^{T}\right)\sum_{i=1}^{n}\bigl{\{}\hat{\psi}_{n,\bm{\alpha},\lambda}(\bm{\alpha}^{T}\bm{X}_{i})-Y_{i}\bigr{\}}\bm{X}_{i}\tilde{\psi}_{n,\bm{\alpha},\lambda}^{\prime}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}_{i}}=\bm{0},$$

$$Y={\psi }_0({\mathbf{\alpha }}_0^T\mathbf{X})+\epsilon ,{\psi }_0(\mathbf{x})=x^3,{\alpha }_{01}={\alpha }_{02}={\alpha }_{03}=1/\sqrt{3},X_1,X_2,X_3\stackrel{i.i.d}{\sim }N[0,1],\epsilon \sim N(0,1),$$

$${\psi }_{\mathbf{\alpha }}(u)=\mathbb{E}\{{\mathbf{\alpha }}_0^T\mathbf{X}|{\mathbf{\alpha }}^T\mathbf{X}=u\}=\frac{\left({\sum }_i{\alpha }_i\right)u\left\{6\left({\sum }_i{\alpha }_i^2\right)\left({\sum }_i{\alpha }_i^2-{\sum }_{i<j}{\alpha }_i{\alpha }_j\right)+{\left({\sum }_i{\alpha }_i\right)}^2{u}^2\right\}}{3\sqrt{3}{\left({\alpha }_1^2+{\alpha }_2^2+{\alpha }_3^2\right)}^3},$$

$$n^{-1}\sum _{i=1}^n{\left\{Y_i-{\hat{\psi }}_{n,\mathbf{\alpha }}({\mathbf{\alpha }}^T{\mathbf{X}}_i)\right\}}^2,$$

$${\hat{\psi }}_{n,\mathbf{\alpha }}({\mathbf{\alpha }}^T\mathbf{x})={\hat{\psi }}_{n,c\mathbf{\alpha }}(c{\mathbf{\alpha }}^T\mathbf{x}),$$

$${\mathbf{S}}_n\mathbf{\alpha }=n^{-1}\sum _{i=1}^n\left({\mathbf{X}}_i-{\bar{\mathbf{X}}}_n\right)Y_i$$

$${\mathbf{S}}_n=n^{-1}\sum _{i=1}^n\left({\mathbf{X}}_i-{\bar{\mathbf{X}}}_n\right){\left({\mathbf{X}}_i-{\bar{\mathbf{X}}}_n\right)}^T,{\bar{\mathbf{X}}}_n=n^{-1}\sum _{i=1}^n{\mathbf{X}}_i,$$

$$n^{-1}\sum _{i=1}^n{\left\{Y_i-{\mathbf{\alpha }}^T({\mathbf{X}}_i-{\bar{\mathbf{X}}}_n)\right\}}^2$$

$$\sum _{i=1}^n{\left\{Y_i-{\psi }_0({\mathbf{\alpha }}^T{\mathbf{X}}_i)\right\}}^2.$$

$$\mathbf{\Sigma }\tilde{\mathbf{\alpha }}=\text{cov}(\mathbf{X},Y),$$

$$\text{cov}(\mathbf{X},Y)$$

$$=\mathbb{E}\left\{{\psi }_0({\mathbf{\alpha }}_0^T\mathbf{X})\mathbb{E}\{\mathbf{X}-\mathbb{E}\mathbf{X}|{\mathbf{\alpha }}_0^T\mathbf{X}\}\right\}.$$

$$\mathbb{E}\{(\mathbf{X}-\mathbb{E}\mathbf{X})|{\mathbf{\alpha }}_0^T\mathbf{X}\}={\mathbf{\alpha }}_0^T(\mathbf{X}-\mathbb{E}\mathbf{X})\mathbf{\beta },$$

Full text

62G05,

62N01,

62-04,

single index model, natural cubic smoothing splines, LSE, Lagrange,

journalname \endlocaldefs

1 Introduction

The monotone single index model tries to predict a response from the linear combination of a finite number of parameters and a function linking this linear combination of parameters to the response via a monotone link function $\psi_{0}$ which is unknown. So, more formally, we have the model

[TABLE]

where $Y$ is a one-dimensional random variable, $\bm{X}=(X_{1},\ldots,X_{d})^{T}$ is a $d$-dimensional random vector with distribution function $G$, $\psi_{0}$ is monotone and $\varepsilon$ is a one-dimensional random variable such that ${\mathbb{E}}[\varepsilon|\bm{X}]=0$ $G$-almost surely. The regression parameter $\bm{\alpha}_{0}$ is a vector of norm $\|\bm{\alpha}_{0}\|_{2}=1$, where $\|\,\cdot\|_{2}$ denotes the Euclidean norm in $\mathbb{R}^{d}$.

The profile least squares estimate of $\bm{\alpha}_{0}$ is an $M$ estimate in two senses: for fixed $\bm{\alpha}$ the least squares criterion

[TABLE]

is minimized for all monotone functions $\psi$ (either decreasing or increasing) which gives an $\bm{\alpha}$ dependent function $\hat{\psi}_{n,\bm{\alpha}}$, and the function

[TABLE]

is minimized over $\bm{\alpha}$. This gives the profile least squares estimator $\hat{\bm{\alpha}}_{n}$ of $\bm{\alpha}_{0}$, which we will call LSE in the sequel. Although this estimate of $\bm{\alpha}_{0}$ has been known now for a very long time (more than 30 years probably), it is not known whether it is $\sqrt{n}$ convergent (under appropriate regularity conditions), let alone that its asymptotic distribution is known. Also, simulation studies are rather inconclusive. For example, it is conjectured in [9] on the basis of simulations that the rate of convergence of $\hat{\bm{\alpha}}_{n}$ is $n^{9/20}$. Other simulation studies were presented in [2] which suggest that even if the LSE is $\sqrt{n}$ convergent under some conditions (for example elliptic symmetry of the distribution of the covariate $\bm{X}$), there are other estimates which have a better performance in these circumstances.

In order to get more insight in the matter, a combination of an $M$ estimate and an $L$ estimate was studied in [2]. For fixed $\bm{\alpha}$ we use the least squares estimator $\hat{\psi}_{n,\bm{\alpha}}$ of the link function, minimizing (1.1) (the $M$ part), just as is done for the LSE, but we try to find the estimate of $\bm{\alpha}_{0}$ by solving an equation in $\bm{\alpha}$ for these least squares functions $\hat{\psi}_{n,\bm{\alpha}}$ (the $L$ part) instead of minimizing again over $\bm{\alpha}$. This means, using the Lagrange approach for the side condition that $\|\bm{\alpha}\|_{2}=1$, that in principle we want to solve an equation of the type

[TABLE]

(differentiating both w.r.t. the argument of the function $x\mapsto\hat{\psi}_{n,\bm{\alpha}}(\bm{\alpha}^{T}\bm{x})$ and w.r.t. its parameter $\bm{\alpha}$), where the Lagrange multiplier $\lambda$ satisfies

[TABLE]

using the condition $\bm{\alpha}^{T}\bm{\alpha}=1$. Note that the left-hand side of (1.3) is $1/2$ times the derivative w.r.t. $\bm{\alpha}$ of the least squares criterion

[TABLE]

where one would try to choose the Lagrange multiplier $\lambda$ in such a way that $\|\bm{\alpha}\|_{2}^{2}=1$, if $\bm{\alpha}$ is the minimizer of (1.5).

Eliminating $\lambda$ using (1.4), this would give the following equation in $\bm{\alpha}$:

[TABLE]

The derivatives in the last factor of the left-hand side of this expression do not exist for the direct least squares estimate $\hat{\psi}_{n,\bm{\alpha}}$. But the equation corresponds in the underlying model to the equation

[TABLE]

where $u\mapsto\psi_{\bm{\alpha}}(u)=E\{\psi_{0}(\bm{\alpha}_{0}^{T}\bm{x})|\bm{\alpha}^{T}\bm{x}=u)$ minimizes

[TABLE]

over nondecreasing functions $\psi$ on $\mathbb{R}$ (see Proposition 3 in [2]), if $G$ is the distribution function of the covariate $\bm{X}$, where we can assume that the derivatives in (1) exist. The first equality in (1) follows from

[TABLE]

and expresses that $\psi_{\bm{\alpha}}$ is the minimizer of (1.8) in the model, for fixed $\bm{\alpha}$.

One can show that the sample equivalent of equation (1), but without the derivative, yields a consistent and $\sqrt{n}$ convergent estimate of the regression parameter $\bm{\alpha}_{0}$. So the Simple Score Estimator (SSE) is defined (approximately) as the value of $\bm{\alpha}$ such that

[TABLE]

where we add “approximately”, since we cannot hope to achieve exact equality to zero because of the discrete character of the left-hand side of the equation. Instead we look for a “crossing of zero” of the left-hand side, as explained in [2]. In practice, this means that we minimize the norm of the left-hand side as a function of $\bm{\alpha}$ in the computer program for computing the SSE, see [4]. Note that we do not need any tuning parameter for this estimator.

The method of computing $\hat{\bm{\alpha}}_{n}$ via (1.9) was introduced in [2], but consistency and asymptotic normality of the estimate defined in this way were not proved. Instead, the proofs of these facts were based on a lower dimensional parametrization and computation of the estimate via the lower dimensional parametrization, involving the Jacobian of the transformation from $\mathbb{R}^{d-1}$ to $\mathbb{R}^{d}$, and not using the Lagrange approach. But in fact, many of the computations in the simulations were done using the Lagrange approach after we noticed empirically that these gave the same results in our simulations.

It is also possible to use the least squares estimator $\hat{\psi}_{n,\bm{\alpha}}$ in the sample equivalent of (1), and use an estimate of the derivative \frac{d}{du}\psi_{\bm{\alpha}}(u)\bigr{|}_{u=\bm{\alpha}^{T}\bm{X}}, constructed from $\hat{\psi}_{\bm{\alpha}}$. This estimate is defined by:

[TABLE]

where $K$ is one of the usual kernels, symmetric around zero and with support $[-1,1]$, and where $h$ is a bandwidth of order $n^{-1/7}$ for sample size $n$. Note that this estimate is rather different from the derivative of a Nadaraya-Watson estimate which is also used in this context and which is in fact the derivative of a ratio of two kernel estimates. Replacing (1.9) by the equation

[TABLE]

then yields, under the appropriate regularity conditions, the so-called efficient score estimator (ESE), which is not only consistent and square root $n$ consistent, but also asymptotically efficient (see [2]).

So, surprisingly, the replacement of the minimization (1.2) in the profile least squares estimator by the search of the solution of (1.9) or (1.11) in $\bm{\alpha}$ produces a $\sqrt{n}$ convergent estimator $\hat{\bm{\alpha}}_{n}$, which is even efficient if we use (1.11).

We can also use as our estimator of the link function a cubic spline $\hat{\psi}_{n,\bm{\alpha}}$, which is defined as the function minimizing

[TABLE]

over the class of functions ${\cal S}_{2}[a,b]$ of differentiable functions $f$ with an absolutely continuous first derivative, where

[TABLE]

see [3], pp. 18 and 19, where $\mu>0$ is the penalty parameter. Using these estimators of the link function, the estimate $\hat{\bm{\alpha}}_{n}$ of $\bm{\alpha}_{0}$ is then found in [8] by using a $(d-1)$-dimensional parametrization $\beta$ and a transformation $S:\bm{\beta}\mapsto S(\bm{\beta})=\bm{\alpha}$, where $S(\bm{\beta})$ belongs to the surface of the unit sphere in $\mathbb{R}^{d}$, and minimizing the criterion

[TABLE]

over $\bm{\beta}$, where $\hat{\psi}_{S(\bm{\beta}),\lambda}$ minimizes (1.12) for fixed $\bm{\alpha}=S(\bm{\beta})$.

Analogously to our approach above we can skip the reparametrizion, and, moreover, can try to solve directly the equation

[TABLE]

where $\tilde{\psi}_{n,\bm{\alpha},\lambda}$ minimizes (1.12) for fixed $\bm{\alpha}$ and $\tilde{\psi}^{\prime}_{n,\bm{\alpha},\lambda}$ is its derivative. The method can also be used for not necessarily monotone functions $\psi_{0}$. We call this estimator the Penalized Least Squares Estimator (PLSE).

Another important point not addressed in [2] is the starting value of the iterations in the computer program. It is rather essential to have a starting point for the search of the regression parameter which is sufficiently close to the real underlying parameter $\bm{\alpha}_{0}$, in particular for the efficient score estimator and spline estimator. The simulations in [2] were therefore done by starting at $\bm{\alpha}_{0}$, something one can of course not do in practice. We therefore investigate how the methods perform if one uses random starting values. To this end we first compute the (cube root $n$ consistent, see [1]) profile least squares estimate of $\bm{\alpha}_{0}$, and next use this as a starting value for the algorithms for the SSE, ESE and PLSE.

It is the purpose of the present paper to provide the changes in the proofs of the consistency and asymptotic distribution of the SSE, ESE and PLSE, computed via the Lagrange-type method given above, but without going through all technicalities again, which are given in [2] and [8]. So we will give the essential changes in the proofs and refer for the needed entropy machinery to these papers. For the simple score estimator (SSE) we go through all the needed steps, though, as an example of the skeleton of the proof for the other estimators. We further illustrate the validity of the Lagrange approach by simulations and also show that one can start the algorithms at a random starting point.

2 An illustrative example

We consider the following model, that was also studied in [2]:

[TABLE]

where $\varepsilon$ is independent of the covariate vector $\bm{X}=(X_{1},X_{2},X_{3})^{T}$.

If we would not make the effort to estimate the link function, but just assume that this function is linear, we get, for a sample of size $n=200$, Figure 1(a), where the points represent the values $(\hat{\bm{\alpha}}_{n}^{T}\bm{X}_{i},Y_{i})$. Although it might seem somewhat unlikely from Figure 1, this linear regression method produces in fact a $\sqrt{n}$ convergent estimate of $\bm{\alpha}_{0}$ in the present model.

On the other hand, if we estimate the function $\psi_{0}$ by the profile Least Squares Estimator (LSE), shown in Figure 1(b), where one uses for fixed $\bm{\alpha}$ the least squares estimate under the restriction that the estimate is nondecreasing in the ordered variables $\bm{\alpha}^{T}\bm{X}_{i}$, and minimize over all such $\bm{\alpha}$, one would hope that the estimate of $\bm{\alpha}_{0}$ is better because here the link function is estimated rather well. However, as mentioned in Section 1, its limiting distribution is unknown, and it is even not known whether it attains the $\sqrt{n}$ rate. This was the motivation to introduce other estimates that are still built on the same pointwise least squares estimate $\hat{\psi}_{n,\bm{\alpha}}$ for fixed $\bm{\alpha}$, but for which we can prove $\sqrt{n}$ convergence and asymptotic normality.

The crucial model function $u\mapsto\psi_{\bm{\alpha}}(u)={\mathbb{E}}\{\bm{\alpha}_{0}^{T}\bm{X}|\bm{\alpha}^{T}\bm{X}=u\}$ is in this model given by:

[TABLE]

which is clearly differentiable w.r.t. $\bm{\alpha}$ and $u$.

In Section 1 we defined the Simple Score Estimator (SSE) and the Efficient Score Estimator (ESE), for which pictures of the corresponding functions $\hat{\psi}_{n,\hat{\bm{\alpha}}}$ are shown in Figure 2. For the ESE also the estimate of the derivative, defined by (1.10), is shown in Figure 2 (c). A picture of the Penalized Least Squares Estimate (PLSE) and its derivative, as defined in Section 1, is given in Figure 3, where a penalty parameter $\mu=0.1$ is used.

We compare the behavior of these estimates in the boxplot Figure 4 with the behavior of three other estimates, the so-called Effective Dimension Reduction (EDR) method, the linear estimate and Han’s maximum correlation estimate. The R script for this comparison can be found on [4].

The simple score estimator (SSE) was computed by minimizing the norm of the left-hand side of (1.9), the efficient score estimator (ESE) by minimizing the norm of the left-hand side of (1.11), and the penalized least squares estimator (PLSE) by minimizing the norm of the left-hand side of (1.13). For the ESE the derivative is estimated by (1.10), for the penalized least squares estimator (PLSE) by the derivative of the cubic spline, which can be computed explicitly, see [3], Section 2.5.1. For the PLSE we took the penalty parameter $\mu=0.1$. In all cases we get an approximate solution of either equation (1.9) (for the SSE) or equation (1.11) (ESE) or equation (1.13) (for the PLSE). The LSE, SSE and ESE estimate $\hat{\psi}_{n,\bm{\alpha}}$ is, for fixed $\bm{\alpha}$, the isotonic least squares estimate, as defined in Section 3 of [2]).

After some experimentation for the best behavior, we used the Hooke-Jeeves method for the PLSE and the Nelder-Mead algorithm for minimizing over $\bm{\alpha}$ in all other cases. The implementation of each method can again be found on [4] in the form of an R script and C++ code.

The LSE is computed in the following way. We generate a number of standard normal $d$-dimensional vectors and normalize these to have norm equal to 1. These are taken as initial values of $\bm{\alpha}$ in the search for the minimizing $\hat{\bm{\alpha}}$. In the example we generated $20$ of such normal vectors for each sample. These were taken as starting values for the search of the minimum of the criterion

[TABLE]

over $\bm{\alpha}$. We do not try to use the constraint $\|\bm{\alpha}\|_{2}=1$, but just let the Nelder-Mead algorithm look for an $\bm{\alpha}$ minimizing (2.2), where we use the following property of the least squares estimate $\hat{\psi}_{n,\bm{\alpha}}$ for fixed $\bm{\alpha}$:

[TABLE]

for each $c>0$ (only the ordering of the variables $\bm{\alpha}^{T}\bm{X}_{i}$ counts, which does not change if we multiply $\bm{\alpha}$ with $c>0$). So the criterion (2.2) also does not change by the scale change (2.3). We normalize the $\hat{\alpha}$ found for each starting value to have norm $1$.

Note that the SSE and the ESE use the same estimators as the LSE; they both minimize the criterion (2.2) for fixed $\bm{\alpha}$, just as the LSE is doing, yielding the same $\hat{\psi}_{n,\bm{\alpha}}$. Only the subsequent search for the estimate of $\bm{\alpha}_{0}$ is done differently: the LSE tries to minimize again, whereas the SSE and ESE try to solve an equation in $\bm{\alpha}$.

The LSE estimate of $\bm{\alpha}_{0}$ was used as starting value for the search procedures for the SSE, ESE, PLSE and Han’s MRE (see [6]). The EDR and the linear method do not need a starting value from outside. This is clear for the linear method. For the EDR the starting value is generated within the algorithm by a simple (not $\sqrt{n}$ convergent) average derivative (w.r.t $\bm{X}$) estimator, which is iteratively improved upon in the algorithm.

The results of this procedure are shown in Figure 4. The EDR was computed using the R package edr, the PLSE by our own R script on [4], and the PLSE2 is computed using the R package simest, where we also used the LSE as starting value. It is seen that Han’s Maximum Rank Correlation method and the linear regression method have in this case a behavior which is clearly inferior to the other methods. If we add dependence to the distribution of the covariates the linear estimate becomes even worse and the other estimates show roughly the same configuration (not shown here). The profile least squares estimate LSE has a remarkably good behavior, although for this estimator not even the $\sqrt{n}$ convergence has been proved, as in fact has been proved for Han’s MRE and the linear regression estimate. The linear regression method attains the $\sqrt{n}$ convergence because the covariates have an elliptically symmetric distribution in the present model, see Section 3.

3 The linear estimator and elliptically symmetric distributions

It is proved in [1] that in the case that the distribution of $\bm{X}$ is elliptically symmetric, ordinary least squares estimates of $\bm{\alpha}_{0}$ are $\sqrt{n}$ convergent, under some (moment) regularity condition on $\bm{X}$ and $Y$. Normal distributions are a particular case of elliptically symmetric distributions. So, no matter how far off the linear fit in Figure 1 (a) seems, the resulting estimate of $\bm{\alpha}_{0}$ is $\sqrt{n}$ convergent. We briefly explain this fact, which seems surprising at first sight.

In computing the linear estimate, we first do not attempt to solve for a vector $\bm{\alpha}$ with norm $1$. So we do not use the Lagrange approach, which we used for the SSE, ESE and PLSE, and simply solve the equation

[TABLE]

in $\bm{\alpha}$, where

[TABLE]

So the solution in fact minimizes

[TABLE]

instead of

[TABLE]

Next we claim that the solution $\hat{\bm{\alpha}}_{n}$ of the linear problem has the property that $\hat{\bm{\alpha}}_{n}/\|\hat{\bm{\alpha}}_{n}\|$ is a $\sqrt{n}$ convergent estimate of $\bm{\alpha}_{0}$.

Let us first investigate why $\hat{\bm{\alpha}}_{n}/\|\hat{\bm{\alpha}}_{n}\|$ would be consistent. The vector $\hat{\bm{\alpha}}_{n}$ converges in probability to a vector $\tilde{\bm{\alpha}}$, satisfying

[TABLE]

where $\bm{\Sigma}$ is the covariance matrix of $\bm{X}$. Furthermore,

[TABLE]

Using the elliptic symmetry, we get:

[TABLE]

for a vector $\bm{\beta}$ such that

[TABLE]

see [1]. Hence:

[TABLE]

and we obtain:

[TABLE]

So equation (3.3) turns into

[TABLE]

where

[TABLE]

showing that $\bm{\Sigma}\tilde{\bm{\alpha}}$ is a multiple of $\bm{\Sigma}\bm{\alpha}_{0}$. Assuming $\bm{\Sigma}$ to be non-singular, we therefore get that $\tilde{\bm{\alpha}}$ is a multiple of $\bm{\alpha}_{0}$. This means that $\hat{\bm{\alpha}}_{n}/\|\hat{\bm{\alpha}}_{n}\|$ converge in probability to $\bm{\alpha}_{0}$, since we may assume in the present case that $\text{cov}\left(\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X}),\bm{\alpha}_{0}^{T}\bm{X}\right)>0$, using the fact that $\psi_{0}$ is a nondecreasing function.

To prove the $\sqrt{n}$ convergence and asymptotic normality, we first derive the asymptotic distribution of $\tilde{\bm{\alpha}}_{n}$, solving (3.1). We have:

[TABLE]

A simple application of the central limit theorem, together with Slutsky’s theorem, give that

[TABLE]

where $c$ is defined by (3.4) and $\bm{A}$ is defined by

[TABLE]

where we follow the notation in [1]. Another application of Slutsky’s theorem now yields that

[TABLE]

where

[TABLE]

Details on the Jacobian for the last application of Slutsky’s theorem are given in [1].

The theorem in [1] runs as follows.

Theorem 3.1**.**

Let $(\bm{X}_{1},Y_{1}),\dots,(\bm{X}_{n},Y_{n})$ be an i.i.d. sample from $(\bm{X},Y)$ such that $E(Y|\bm{X})=\psi_{0}(\bm{\alpha}_{0}\bm{X})$ almost surely, where $\psi_{0}$ is non-decreasing and $\bm{\alpha}_{0}^{T}\bm{\alpha}_{0}=1$. Suppose that $\bm{X}$ has an elliptically symmetric distribution with finite mean $\bm{\mu}\in\mathbb{R}^{d}$ and a positive definite covariance matrix $\bm{\Sigma}$. Assume, moreover, that ${\mathbb{E}}\|Y\bm{X}\|<\infty$ and that there exists a nonempty interval $[a,b]$ on which $\psi_{0}$ is strictly increasing. Then, as $n\to\infty$, the estimator $\tilde{\bm{\alpha}}_{n}=\hat{\bm{\alpha}}_{n}/\|\hat{\bm{\alpha}}_{n}\|$, where $\hat{\bm{\alpha}}_{n}$ is defined by

[TABLE]

converges in probability to $\bm{\alpha}_{0}$. If, moreover, ${\mathbb{E}}Y^{2}\|\bm{X}\|^{2}<\infty$, $\sqrt{n}\{\tilde{\bm{\alpha}}_{n}-\bm{\alpha}_{0}\}$ converges in distribution to a normal distribution with mean $\bm{0}$ and covariance matrix

[TABLE]

where $\bm{\Gamma}$ is defined by (3.5).

4 The (profile) LSE

The consistency of the profile LSE can be derived in the following way. The function $\bm{x}\mapsto\hat{\psi}_{n,\hat{\bm{\alpha}}_{n}}(\hat{\bm{\alpha}}_{n}^{T}\bm{x})$ minimizes the expression:

[TABLE]

under the restriction $\|\bm{\alpha}\|_{2}=1$, where we subtract in the first line a term that does not influence the minimization. Since the function $\bm{x}\mapsto\hat{\psi}_{n,\hat{\bm{\alpha}}_{n}}(\hat{\bm{\alpha}}_{n}^{T}\bm{x})$ minimizes this expression, it gives in particular a smaller value of the expression than the function $\bm{x}\mapsto\psi_{0}(\bm{\alpha}_{0}^{T}\bm{x})$. This implies:

[TABLE]

and hence

[TABLE]

So we find:

[TABLE]

almost surely, which implies

[TABLE]

almost surely. From this and the fact that $\bm{\alpha}_{n}$ has norm 1, the consistency of $\hat{\bm{\alpha}}_{n}$ for $\bm{\alpha}_{0}$ can easily be deduced.

The convergence rate and asymptotic distribution of the LSE are unknown. In [1] a cube root convergence rate for both $\hat{\psi}_{n,\hat{\bm{\alpha}}_{n}}$ and $\hat{\bm{\alpha}}_{n}$ is derived.

5 Consistency and asymptotic normality of the SSE

We prove consistency under the conditions of [2]. Let $(\hat{\bm{\alpha}}_{n_{k}})$ be a subsequence converging to a vector $\tilde{\bm{\alpha}}$, satisfying $\|\tilde{\bm{\alpha}}\|_{2}=1$. The vector $\tilde{\bm{\alpha}}$ has to satisfy:

[TABLE]

where

[TABLE]

Since ${\mathbb{E}}\left\{Y|\bm{X}\right\}=\psi_{0}(\bm{\alpha}_{0}^{T}\bm{X})$, we can write this in the form

[TABLE]

Premultiplying this by $(\bm{\alpha}_{0}-\tilde{\bm{\alpha}})^{T}$ we get:

[TABLE]

So either $\tilde{\bm{\alpha}}=\pm\bm{\alpha}_{0}$ (using $\|\tilde{\alpha}\|_{2}=1$), or $\tilde{\bm{\alpha}}_{0}$ has a component in the space perpendicular to $\bm{\alpha}$, in which case we get that $(\bm{\alpha}_{0}-\tilde{\bm{\alpha}})^{T}\bm{X}$ is, conditionally on $\tilde{\bm{\alpha}}^{T}\bm{X}$, nondegenerate and that

[TABLE]

contradicting (5). So we obtain

[TABLE]

Since a sufficiently small neighborhood around $\bm{\alpha}_{0}$ excludes $-\bm{\alpha}_{0}$, the consistency follows if we only look for a solution in this neighborhood.

We now give the argument for the following theorem on the asymptotic normality of $\sqrt{n}\{\hat{\bm{\alpha}}_{n}-\bm{\alpha}_{0}\}$.

Theorem 5.1**.**

Let Assumptions A1-A6 of [2]) be satisfied and assume that $\|\bm{\alpha}_{0}\|_{2}=1$ and that

A7. There exists a $\delta_{0}>0$ such that for all $\bm{\alpha}\neq\bm{\alpha}_{0}$, satisfying $\|\bm{\alpha}-\bm{\alpha}_{0}\|<\delta_{0}$ and $\|\bm{\alpha}\|_{2}=1$,

[TABLE]

A9. {\mathbb{E}}\Bigl{[}\psi_{0}^{\prime}(\bm{\alpha}_{0}^{T}\bm{X})\,\text{\rm var}(\bm{X}|\bm{\alpha}_{0}^{T}\bm{X})\Bigr{]} has rank $d-1$.

Let $\hat{\bm{\alpha}}_{n}$ be defined as a crossing of zero of the left-hand side of (1.9). Let the matrices $\bm{A}$ and $\bm{\Sigma}$ be defined by:

[TABLE]

and

[TABLE]

Then

[TABLE]

where $\bm{A}^{-}$ is the Moore-Penrose inverse of $\bm{A}$ and $N_{d}$ denotes the degenerate $d$-dimensional normal distribution with mean zero and covariance matrix $\bm{A}^{-}\bm{\Sigma}\bm{A}^{-}$.

Remark 5.1**.**

This is a modification of the asymptotic normality part of Theorem 5 in [2], avoiding the manipulations with the Jacobians, which were needed there because of the lower dimensional parametrization. **

Before going into the proof, we highlight the essential elements. We first note that we may assume, as in [2], that the estimator $\hat{\bm{\alpha}}_{n}$ satisfies the equation:

[TABLE]

Using the characterization of the least squares estimator $\hat{\psi}_{n,\bm{\alpha}}$ for fixed $\bm{\alpha}$, we show that

[TABLE]

see (5) below. Next we show that this, together with (5.5), leads to the asymptotic equation

[TABLE]

see (5). The dominating terms on both sides are perpendicular to $\bm{\alpha}_{0}$. From this, the consistency and the relation $\|\bm{\alpha}\|_{2}=\|\bm{\alpha}_{0}\|_{2}=1$, the result follows.

Proof.

We define the piecewise constant function $\bar{E}_{n,\bm{\alpha}}$

[TABLE]

where the $\tau_{i}$ are the ordered points of jump of the least squares estimate $\hat{\psi}_{n\bm{\alpha}}$, and where $\psi_{\bm{\alpha}}(u)={\mathbb{E}}\left\{\psi_{0}(\bm{\alpha}_{0}^{T}\bm{x})|\bm{\alpha}^{T}\bm{x}=u\right\}$. The origin of this method is the treatment of smooth functionals in the interval censoring model, see, e.g., p. 290 of [5].

Then, since $\hat{\psi}_{n,\bm{\alpha}}$ is the least squares estimator for fixed $\bm{\alpha}$ and since the function $\bar{E}_{n,\bm{\alpha}}$ is constant between jumps:

[TABLE]

implying that we can write (5.5) in the form:

[TABLE]

Note that $\bar{E}_{n,\bm{\alpha}}$ cannot computed in practice, since we cannot compute the function $\psi_{\bm{\alpha}}$ using the data only.

We have:

[TABLE]

The last term can be written:

[TABLE]

Since, by the choice of $\bar{E}_{n,\hat{\bm{\alpha}}_{n}}$,

[TABLE]

(see for this type of argument for example the top of p. 308 in [5]), we find that the first integral is $o_{p}(n^{-1/2})$. Using (5.10) again, together with the fact that $|\psi_{0}(\bm{\alpha}_{0}^{T}\bm{x})-\psi_{\hat{\bm{\alpha}}_{n}}(\hat{\bm{\alpha}}_{n}^{T}\bm{x})|$ is bounded by a constant times $\|\hat{\bm{\alpha}}_{n}-\bm{\alpha}_{0}\|$, we find that the second term is of smaller order than $\|\hat{\bm{\alpha}}_{n}-\bm{\alpha}_{0}\|$. Finally, again using (5.10) we find from the Cauchy-Schwarz inequality that the third term is bounded above by a constant times

[TABLE]

So we get the relation

[TABLE]

Furthermore,

[TABLE]

where the last equality follows from:

[TABLE]

So, using

[TABLE]

we end up with the equation

[TABLE]

where we use that premultiplying with $\hat{\bm{\alpha}}_{n}\hat{\bm{\alpha}}_{n}^{T}$ yields zero for the dominating terms.

Using ${\mathbb{E}}(\bm{X}|\hat{\bm{\alpha}}_{n}^{T}\bm{X})-{\mathbb{E}}\left(\bm{X}|\bm{\alpha}_{0}^{T}\bm{X}\right)=O_{p}(\hat{\bm{\alpha}}_{n}-\bm{\alpha}_{0})$, it is easily seen that

[TABLE]

Likewise,

[TABLE]

where the last dominating term satisfies

[TABLE]

So we arrive at the following asymptotic equation:

[TABLE]

Defining

[TABLE]

we get:

[TABLE]

for some $k\in\mathbb{R}$, where $\bm{B}^{-}$ is the Moore-Penrose inverse of $\bm{B}$.

Define

[TABLE]

Then, since $\bm{\beta}_{n}$ is perpendicular to $\bm{\alpha}_{0}$, we get:

[TABLE]

implying

[TABLE]

and hence

[TABLE]

The second solution would lead to

[TABLE]

which is excluded by $\hat{\bm{\alpha}}_{n}-\bm{\alpha}_{0}=o_{p}(1)$ (the consistency of $\hat{\bm{\alpha}}_{n}$). Using $\|\bm{\beta}_{n}\|_{2}^{2}=O_{p}(n^{-1})$ we get:

[TABLE]

and hence

[TABLE]

∎

6 Consistency and asymptotic normality of the Efficient Score Estimate (ESE) and the Penalized Least Squares Estimate (PLSE) as defined by (1.13)

The proofs of the consistency and asymptotic normality of the ESE and PLSE are highly similar to the proofs of these facts for the SSE in the preceding section. The only extra ingredient is occurrence of the estimate of the derivative of the link function. We only discuss the asymptotic normality.

An essential step is again to show that

[TABLE]

For the ESE this is done by defining the piecewise constant function $\bar{\rho}_{n,\bm{\alpha}}$ for $u$ in the interval between successive jumps $\tau_{i}$ and $\tau_{i+1})$ of $\hat{\psi}_{n\bm{\alpha}}$ by:

[TABLE]

where $\bar{\rho}_{n,\bm{\alpha}}$ replaces $\bar{E}_{n,\bm{\alpha}}$ in (5.9), see Appendix E in the supplement of [2]. The remaining part of the proof runs along the same lines as the proof for the SSE. For additional details, see Appendix E in the supplement of [2].

The corresponding step in the proof for the PLSE is given by the following lemma.

Lemma 6.1**.**

Let the conditions of Theorem 5 in [8] be satisfied. In particular, let the penalty parameter $\mu_{n}$ satisfy $\mu_{n}=o_{p}(n^{-1/2})$. Then we have for all $\bm{\alpha}$ in a neighborhood of $\bm{\alpha}_{0}$ and for the corresponding penalized least squares estimate (natural cubic spline) $\hat{\psi}_{n\bm{\alpha}}$:

[TABLE]

Remark 6.1**.**

The result shows that we have as our basic equation in $\bm{\alpha}$:

[TABLE]

Proof.

Fix $\bm{\alpha}$ and let $\hat{\psi}_{n,\bm{\alpha}}$ be the natural cubic spline, minimizing

[TABLE]

over functions $g\in{\cal S}_{2}[a,b]$, where the $t_{i}$ are the the ordered values $\bm{\alpha}^{T}\bm{X}_{i}$, and where $\mu_{n}=o_{p}(n^{-1/2})$, and where $a=\min_{i}\bm{\alpha}^{T}\bm{X}_{i}$ and $b=\max_{i}\bm{\alpha}^{T}\bm{X}_{i}$. We can write the minimum in the following form:

[TABLE]

We extend the natural cubic spline $\hat{\psi}_{n,\bm{\alpha}}$ linearly to a function on $\mathbb{R}$, and define the function

[TABLE]

We have:

[TABLE]

This implies, assuming the boundedness of the derivative of the function $t\mapsto\frac{\partial}{\partial t}{\mathbb{E}}(\bm{X}|\bm{\alpha}^{T}\bm{X}=t)$ for $t\in[a,b]$,

[TABLE]

Also assuming that

[TABLE]

(see Theorem 2 in [8]), we obtain from (6.2):

[TABLE]

since the function $\phi$ attains its minimum at $\bm{0}$ by the definition of the (natural) cubic spline as a least squares estimate. It follows that

[TABLE]

∎

The remaining part of the proof of the asymptotic normality can either run along the same lines as the proof for the corresponding fact for the SSE, using the function $u\mapsto\psi_{\bm{\alpha}}(u)={\mathbb{E}}\{\psi_{0}(\bm{\alpha}^{T}\bm{x})|\bm{\alpha}^{T}\bm{X}=u\}$, or directly use the convergence of $\hat{\psi}_{n\hat{\bm{\alpha}}_{n}}$ to $\psi_{0}$ and of $\hat{\psi}^{\prime}_{n\hat{\bm{\alpha}}_{n}}$ to $\psi_{0}^{\prime}$ (see Theorem 3 in [8]). For the SSE and ESE we were forced to introduce the intermediate function $\psi_{\bm{\alpha}}$ to get to the derivatives, because there the derivative of $\hat{\psi}_{n\hat{\bm{\alpha}}_{n}}$ did not exist.

We get the following result.

Theorem 6.1**.**

Let either $\hat{\bm{\alpha}}_{n}$ be the ESE of $\bm{\alpha}_{0}$ and let Assumptions A1-A8 of [2] be satisfied and let $\psi_{\bm{\alpha}}$ be twice continuously differentiable or let $\hat{\bm{\alpha}}_{n}$ be the PLSE of $\bm{\alpha}_{0}$ and let Assumptions (A0) to (A6) and (B1) to (B3) of [8]) be satisfied. Moreover, let the bandwidth $h\asymp n^{-1/7}$ in the estimate of the derivative of $\psi_{\bm{\alpha}}$ for the ESE. Define the matrices,

[TABLE]

and

[TABLE]

Then

[TABLE]

where $\tilde{\bm{A}}^{-}$ is the Moore-Penrose inverse of $\tilde{\bm{A}}$.

This corresponds to Theorem 6 in [2] and Theorem 5 in [8]), but note that the formulation of Theorem 5 in [8] still contains the Jacobian connected with the lower dimensional parametrization.

7 Simulations

In this section we give results similar to those in Section 7 of [2]. First we want to show that the Lagrange approach and the random starting values give results similar to those in Table 2 of [2] for the model introduced in Section 2 of the present paper. This can be seen from Table 1, where the results for the SSE, ESE, LSE, MRE and linear estimate can be compared with the corresponding results inTable 2 of [2]. The MRE is called MRCE (Maximum Rank Correlation Estimator) and the linear estimator is called H-LFLSE (“Hybrid Link-Free Least Squares Estimator”) in Table 2 of [2]. Moreover, we have added the results for the penalized least squares estimator, as computed via the algorithm in [4] (PLSE) and via the R package simest (PLSE2), respectively.

It is clear that there is good correspondence between the results. It is also of interest to see (again) that the linear estimator and Han’s maximumum rank correlation estimator have a behavior that is inferior to that of the other estimators. The good behavior of the penalized least squares estimators is also of interest. The penalty parameter $\mu$ was taken $0.1$ for sample sizes $100,500$ and $1000$ and taken $0.05$ for $n=2000$, but the results do not seem very sensitive to this choice. For exampe, if we take $\mu=0.1$ also for sample size $n=2000$, the results are almost the same. It is also seen that the behavior of the PLSE as computed via the algorithm in [4] has a slightly better behavior for the larger sample sizes, possibly because the solution via equation (1.13), motivated by the Lagrange-type approach, gives a more stable behavior than the application of the argmax approach for $\bm{\alpha}$. The splines $\hat{\psi}_{n,\bm{\alpha}}$ are the same in the two methods, but the search for the best $\bm{\alpha}$ is different, also because different search methods are used; [4] uses the Hooke-Jeeves algorithm and simest uses Broyden’s algorithm.

Just as in [2], Section 7.2, we investigate the effect of an increasing dimension for this model. The result is given in Figure 5. The inferiority of the linear estimator and of Han’s maximum rank correlation estimator is again clear. The estimator EDR is deteriorating in comparison with the estimators SSE, ESE, LSE, PLSE and PLSE2 with increasing dimension and, moreover, its computing time increases steeply with increasing dimension (see also for this Section 7.2 of [2]).

There is again some evidence that the PLSE, computed by the algorithm in [4] has a slightly better behavior than the PLSE2, computed by simest. Both estimates seem to have the best behavior among the estimates studied here, though, with the ESE as a close second.

8 Conclusion

A natural estimate the finite dimensional parameter in the monotone single index is the profile least squares estimator which we called the LSE. But its asymptotic distribution is unknown and it is even not known whether this estimate is $\sqrt{n}$ consistent. The LSE is an $M$ estimate in two senses: for fixed $\bm{\alpha}$ the least squares criterion (1.1) is minimized over all monotone functions $\psi$ and the function (1.2) is minimized for the resulting $\hat{\psi}_{n,\bm{\alpha}}$ over $\bm{\alpha}$.

We proposed to replace the second minimization by the solution of an equation arising from a Lagrange formulation of the minimization problem. This method yields estimators for which one can prove $\sqrt{n}$ convergence and asymptotic normality. Introducing the estimate of the derivative of the link function (1.10), one can even construct an asymptotically efficient estimate in this way. The method avoids reparametrization, which is the usual way to go in this model, and continues work in [2].

It is shown that one can also apply this approach for a method which uses cubic splines in the first minimization step, which then becomes a penalized minimization problem. Our method avoids again reparametrization and continues work in this direction in [8], which was implemented in the R package simest.

The methods are illustrated by simulations, and the implementation can be found on [4] in the form of R scripts and C++ code. We also address the important issue of the starting points of the algorithms and show that one can start by using the result of the LSE as starting value. The LSE itself is computed from a number of random starting values and the estimate resulting from the starting value which gives a minimum of the criterion function is used as starting value for the other algorithms.

The estimates are compared with a number of estimates in the literature: the estimate in the R package simest, Han’s maximum rank correlation estimate (see [6]), a linear estimate, applicable if the distribution of the covariate is elliptically symmetric, and the effective dimension reduction estimate, introduced in [7], and available via the R package edr. In our comparison the penalized least squares estimates come out best, with the efficient score estimate ESE as a close second. The linear estimate and Han’s maximum rank correlation estimate seem clearly inferior to the other estimators.

In the simulation examples we considered covariates with a normal distribution, as a particular example of elliptically symmetric distributions, since in that case the conditions for the application of the linear method are satisfied. We also took normal errors, since this was a condition for the application of edr with which we also make a comparison. Other comparisons with edr and other packages can be found in [8] and [2]. A general finding is that the linear method completely breaks down if the condition of elliptic symmetry of the covariates is not fulfillled.

Further improvement of our Lagrange-type method is to be expected if we refine the search algorithms. Although most of the programming was done from first principles, in particular the programming of the spline method which used the Reinsch method, as exposed in [3], to make effective use of the band matrix structure, the Hooke-Jeeves and Nelder-Mead C++ algorithms were “taken off the (internet) shelf” (see [4]), and certainly can be improved upon. Because the criterion function is not continuous, at least for the SSE, ESE and LSE, doing the search for $\bm{\alpha}$ in the right way is rather non-trivial. Changes of $\bm{\alpha}$ lead to changes of the order of the values $\bm{\alpha}^{T}\bm{X}_{i}$ and this leads to discontinuities of the criterion function. Note that these searches are “grid-free”, looking for a solution on a grid is a much too crude method in these problems.