Inference for L-estimators of location using a bootstrap warping approach

Abstract

In this note we propose a new semi-parametric bootstrap procedure for hypothesis tests about a statistical function and termed bootstrap warping. This procedure was motivated by empirical likelihood and bootstrap tilting techniques. The procedure is computationally efficient and has a fixed number of parameters. We show that the warping procedure has good type I error control and has monotone power as a function of sample size and shift alternatives.

Keywords:

• L-estimators
• L-moments
• LQ-moments
• Order statistics

1. Introduction

Let $X_1,X_2,\dots ,X_n$ denote an i.i.d. sample drawn from an absolutely continuous population with cumulative distribution function (c.d.f.) denoted as F and corresponding quantile function denoted as $Q(u)=F^{-1}(u).$ For the application described in this note we are interested in making inferences about a one-dimensional parameter of the form $\theta =T(F),$ where in our methodology T(F) denotes a specific smooth statistical functional for measuring expectation and having the form (1.1) $$T\left(F\right)=E\left[g\right(X\left)j\right(F\left(X\right)\left)\right]=\int g\left(X\right)j\left(F\right(X\left)\right)\mathit{dF}={\int }_0^{1}g\left(Q\right(u\left)\right)j\left(u\right)\mathit{du}.$$(1.1)

In our applications we restrict $j(F(X))$ to be a smooth absolutely continuous weighting function such that $\int j(F(x))\mathit{dF}={\int }_0^{1}j(u)\mathit{du}=1$ is essentially a p.d.f. for $j(u)>0,\forall u\in (0,1)$ and we assume $\int g(x)\mathit{dF}$ is bounded, e.g. suppose the parameter of interest is the population mean, then the statistical functional has the well-known form $\theta =T(F)=E(g(X)j(F(X)))=\int xdF={\int }_0^{1}Q(u)\mathit{du},$ with $j(F(x))=1$ or alternatively j(u) = 1 and g(x) = x. The classic “bootstrap” estimator of T(F) is given by replacing the c.d.f. F with its empirical counterpart, $\widehat{F}(x)={\sum }_{i=1}^n{I}_{(x_i\le x)}/n,$ in (1.1)(1.1) $$T\left(F\right)=E\left[g\right(X\left)j\right(F\left(X\right)\left)\right]=\int g\left(X\right)j\left(F\right(X\left)\right)\mathit{dF}={\int }_0^{1}g\left(Q\right(u\left)\right)j\left(u\right)\mathit{du}.$$(1.1) or alternatively replacing Q with its empirical counterpart, $\widehat{Q}(u)=x_{[\mathit{nu}]+1:n},$ where $[·]$ denotes the floor function and $X_{i:n}$ denotes the ith order statistic. Substituting $\widehat{F}$ into (1.1)(1.1) $$T\left(F\right)=E\left[g\right(X\left)j\right(F\left(X\right)\left)\right]=\int g\left(X\right)j\left(F\right(X\left)\right)\mathit{dF}={\int }_0^{1}g\left(Q\right(u\left)\right)j\left(u\right)\mathit{du}.$$(1.1) for F yields the empirical estimator of T(F), which has the well-known form (1.2) $$T\left(\widehat{F}\right)=\sum _{i=1}^{n}g\left(X_{i:n}\right){\int }_{(i-1)/n}^{i/n}j\left(u\right)\mathit{du}.$$(1.2)

Some classic examples for $T(\widehat{F})$ include, kernel density and quantile estimators, sample moment estimators and L-estimators, e.g. see Serfling (1980) for a technical overview of estimators having this form relative to their asymptotic properties.

Now, suppose we are interested in testing a hypothesis about the given statistical functional T(F) having the form $H_0:T(F_0)={\theta }_0$ versus $H_1:T(F_0)>{\theta }_0,$ or without loss of generality $H_1:T(F_0)<{\theta }_0,$ in a semiparametric fashion. Note that we will touch on two-sided tests later in this note. Popular nonparametric approaches for testing hypotheses of this type are given by the well-known empirical likelihood (EL) method due to Owen (1988) and bootstrap tilting methods such as exponential tilting or other multinomial based resampling schemes, e.g. see Davison and Hinkley (1997).

The nonparametric EL and bootstrap approaches provide the motivation for our new semiparametric testing methodology. The key idea behind the EL and bootstrap tilting approaches is to find the nonparametric maximum likelihood estimator for the probability density function (p.d.f.) f₀ and c.d.f F₀ given the constraint $T(F_0)={\theta }_0$ prescribed under H₀, as estimated by its empirical counterpart $T({\widehat{F}}_0)={\theta }_0,$ where ${\widehat{F}}_0={\sum }_{i=1}^n{v}_i{I}_{(x_i\le x)}$ and the v_i parameters corresponding to f₀ sum to 1 and are bounded between 0 and 1. The common definition for the v_i parameters in the continuous case for the discretized model is given as $v_i=F(x_i)-F(x_i-),$ where $F(x)=P(X\le x)$ and $F(x-)=P(X<x),$ respectively, e.g. see Owen (1988) for the technical argument pertaining to this definition.

In the most common scenario the likelihood under the “unconstrained” alternative hypothesis yields the classic estimates of ${\widehat{v}}_i=1/n$ for simple statistics such as the sample mean. Other weights may occur for functionals corresponding to trimmed estimators, e.g. see Qin and Taso (2002) with respect to the weights for trimmed mean. The weights, v_i, under the null hypothesis are generally determined by minimizing a given distance measure such as the Kullback-Leibler distance $D(\mathbf{v},{\mathbf{v}}_0)={\sum }_{i=1}^n{v}_i\hspace{0.17em}\text{log}\hspace{0.17em}(n{v}_i),$ where the $1\times n$ vectors $\mathbf{v}=(1/n,1/n,\dots ,1/n)$ and ${\mathbf{v}}_0=(v_1,v_2,\dots ,v_n).$ Alternatively, one may use constrained maximum likelihood approaches for determining the vector ${\mathbf{v}}_0=(v_1,v_2,\dots ,v_n)$ under H₀, e.g. see Vexler and Gurevich (2010) a typical model scenario. We use this idea of a discretized model as a starting point for developing an alternative based on smooth statistical functional inferential procedure using what we term statistical warping as defined in the next section.

In this Section 2 we outline the bootstrap warping procedure and follow this with simulation study in Section 3

2. Bootstrap warping and hypothesis testing

The key features of bootstrap warping, as contrasted with the EL and bootstrap tilting approaches, is that it is a semi-parametric approach and that the number of parameters in the model is reduced from n – 1 to 2, conditional on the observed data, i.e. we “warp” the observed e.d.f. $\widehat{F}(x)={\sum }_{i=1}^n{I}_{(x_i\le x)}/n$ versus treating each discretized segment on the continuum as a parameter. In addition, our resampling scheme follows the classic bootstrap multinomial resampling scheme with cell probabilities of $1/n$ versus bootstrap tilting, which requires the weights to be determined conditional on the dataset under investigation, thus adding a layer of complexity to the computational components to these problems. The direct benefit of this parameterization is computational ease without suffering the “curse of dimensionality” associated with big data scenarios. This will be described in detail below. Additionally, in terms of future work, covariate adjustments may be made through the warping model parameters, thus extending the utility of this approach to more complex settings.

In terms of testing $H_0:T(F_0)={\theta }_0$ versus $H_1:T(F_0)\ne {\theta }_0$ we need to first define $T({\widehat{F}}_0)$ relative to obtaining an empirical version of the constraint $T(F_0)={\theta }_0.$ Towards this end we define the warped empirical estimator of $T(F_0)$ based on formulation (1.1) as (2.1) $$T({\widehat{F}}_0;\alpha ,\beta )=\sum _{i=1}^{n}g\left(X_{i:n}\right)\left\{K_{\alpha ,\beta }\left[J\left(\frac{i}{n}\right)\right]-K_{\alpha ,\beta }\left[J\left(\frac{i-1}{n}\right)\right]\right\},$$(2.1) (2.2) $$=\sum _{i=1}^{n}g\left(X_{i:n}\right)f\left(X_{i:n}\right),$$(2.2) where $J(v)={\int }_0^{v}j(u)\mathit{du},f(X_{i:n})$ is a weighting function defined more formally below at (2.5)(2.5) $$f_{\alpha ,\beta }\left(X_{i:n}\right)=F_{\alpha ,\beta }\left(X_{i:n}\right)-F_{\alpha ,\beta }(X_{i:n}-)=K_{\alpha ,\beta }\left[J\left(\frac{i}{n}\right)\right]-K_{\alpha ,\beta }\left[J\left(\frac{i-1}{n}\right)\right],$$(2.5) , $j(·)$ is defined at (1.1) and we define $K_{\alpha ,\beta }$ such that $K_{\alpha ,\beta }(t)=t$ under H₀. For example, if we were interested in testing about $E(g(X))$ then (2.3) $$T({\widehat{F}}_0;\alpha ,\beta )=\sum _{i=1}^{n}g\left(X_{i:n}\right)\left[K_{\alpha ,\beta }\left(\frac{i}{n}\right)-K_{\alpha ,\beta }\left(\frac{i-1}{n}\right)\right]$$(2.3) and under H₀ true we would have $T({\widehat{F}}_0)={\sum }_{i=1}^{n}g(X_{i:n})/n.$

Comment. Note that $T({\widehat{F}}_0;\alpha ,\beta )$ at (2.3)(2.3) $$T({\widehat{F}}_0;\alpha ,\beta )=\sum _{i=1}^{n}g\left(X_{i:n}\right)\left[K_{\alpha ,\beta }\left(\frac{i}{n}\right)-K_{\alpha ,\beta }\left(\frac{i-1}{n}\right)\right]$$(2.3) is used specifically to generate the null distribution for the estimator $T(\widehat{F})$ at (1.2)(1.2) $$T\left(\widehat{F}\right)=\sum _{i=1}^{n}g\left(X_{i:n}\right){\int }_{(i-1)/n}^{i/n}j\left(u\right)\mathit{du}.$$(1.2) and is not meant to be alternative estimator for θ.

The components of the weighting function in (2.2)(2.2) $$=\sum _{i=1}^{n}g\left(X_{i:n}\right)f\left(X_{i:n}\right),$$(2.2) , denoted as $K_{\alpha ,\beta }(·),$ are defined as the c.d.f. of a Kumaraswamy distribution and given as (2.4) $$K_{\alpha ,\beta }\left(u\right)=1-{(1-u^{\alpha })}^{\beta },$$(2.4) where $0<u<1,\alpha >0$ and $\beta >0.$ The choice of the Kumaraswamy distribution in terms of a weighting function is due to its numerical tractability and flexibility in terms of the relative shapes it contains, i.e. our test will be sensitive to a number of alternatives given H₁ via the choice of this weighting function.

Our semi-parametric density utilized within (2.2)(2.2) $$=\sum _{i=1}^{n}g\left(X_{i:n}\right)f\left(X_{i:n}\right),$$(2.2) is now defined as discretized type model similar to what is used in the EL methodology and bootstrap tilting and given as (2.5) $$f_{\alpha ,\beta }\left(X_{i:n}\right)=F_{\alpha ,\beta }\left(X_{i:n}\right)-F_{\alpha ,\beta }(X_{i:n}-)=K_{\alpha ,\beta }\left[J\left(\frac{i}{n}\right)\right]-K_{\alpha ,\beta }\left[J\left(\frac{i-1}{n}\right)\right],$$(2.5) where $F(x)=P(X\le x)$ and $F(x-)=P(X<x),$ and $f(X_{i:n})$ represents a point mass corresponding to the ith order statistic. The Kumaraswamy distribution was chosen over other candidate distributions, e.g. the beta distribution, due to its well-behaved numerical properties and relatively straightforward parameterization. See Jones (2009) for a detailed description of the Kumaraswamy distribution and a description of its close relationship to the beta distribution. In essence $f(X_{i:n})$ serves as a standard weighting function such that when α = 1 and β = 1 then $T(\widehat{F};\alpha ,\beta )$ equates to $T({\widehat{F}}_0;\alpha ,\beta ).$

The test of interest in this note is given as $H_0:T(F_0)={\theta }_0$ versus $H_1:T(F_0)>{\theta }_0.$ As in EL methods and bootstrap tilting the first step is to maximize the constrained pseudo-likelihood (2.6) $$L_{\alpha ,\beta }=\prod _{i=1}^n{f}_{\alpha ,\beta }\left(X_{i:n}\right),$$(2.6) with respect to α and β and under the constraint $H_0:T(F_0)={\theta }_0,$ where $f(X_{i:n})$ is defined at (2.5)(2.5) $$f_{\alpha ,\beta }\left(X_{i:n}\right)=F_{\alpha ,\beta }\left(X_{i:n}\right)-F_{\alpha ,\beta }(X_{i:n}-)=K_{\alpha ,\beta }\left[J\left(\frac{i}{n}\right)\right]-K_{\alpha ,\beta }\left[J\left(\frac{i-1}{n}\right)\right],$$(2.5) . Clearly, α = 1 and β = 1 given H₀ is true.

The bootstrap resampling scheme for our inferential method then is as follows:

1. Calculate the observed test statistic $T(\widehat{F}).$

2. Obtain $\widehat{\alpha }$ and $\widehat{\beta }$ from (2.6)(2.6) $$L_{\alpha ,\beta }=\prod _{i=1}^n{f}_{\alpha ,\beta }\left(X_{i:n}\right),$$(2.6) .

3. Generate B nonparametric bootstrap samples of size n, i.e. generate n uniform (0,1) random variables and apply $\widehat{Q}(u)=x_{[\mathit{nu}]+1:n}$ to those randomly generated uniform variates.

4. Calculate $T_i^{*}({\widehat{F}}_0|\widehat{\alpha },\widehat{\beta })$ from (2.2)(2.2) $$=\sum _{i=1}^{n}g\left(X_{i:n}\right)f\left(X_{i:n}\right),$$(2.2) , replacing α with $\widehat{\alpha }$ and β with $\widehat{\beta }$ in step (2), for $i=1,2,\dots ,B.$

5. Calculate the approximate one-sided bootstrap p-value $p_{\mathit{boot}}={\sum }_{j=1}^{B}I(T_i^{*}({\widehat{F}}_0|\widehat{\alpha },\widehat{\beta })>T(\widehat{F}))/B,$ where $T(\widehat{F})$ is the observed estimator defined at (1.2)(1.2) $$T\left(\widehat{F}\right)=\sum _{i=1}^{n}g\left(X_{i:n}\right){\int }_{(i-1)/n}^{i/n}j\left(u\right)\mathit{du}.$$(1.2) .

For a test of $H_0:T(F_0)={\theta }_0$ versus $H_1:T(F_0)<{\theta }_0$ simply reverse the inequality in step (5) above. For the test $H_0:T(F_0)={\theta }_0$ versus $H_1:T(F_0)\ne {\theta }_0$ there is an added assumption of the symmetry of the distribution of $T_i^{*}({\widehat{F}}_0|\widehat{\alpha },\widehat{\beta })$ under H₀. Under this assumption the two-sided p-value is given as $p_{\mathit{boot}}={\sum }_{j=1}^{B}I(|T_i^{*}({\widehat{F}}_0|\widehat{\alpha },\widehat{\beta })|>|T(\widehat{F})|)/B.$ In general, for most of the tests of interest the test statistic will have an asymptotic normal distribution thus most two-sided tests should satisfy approximately a symmetry assumption, i.e. the statistics are based on smooth functions, which then in turn lend themselves to well-behaved and symmetric bootstrap resampling distributions.

As with similar bootstrapping methodologies for inference. e.g. see Davison and Hinkley (1997), the key is that $T({\widehat{F}}_0;\alpha ,\beta )$ is a consistent estimator of $T(F_0)$ under H₀, which in the methodology presented above holds given α = 1 and β = 1 under H₀, $\widehat{F}$ converges to F and $T(\widehat{F};\alpha ,\beta )$ converges to T(F) given the smoothness conditions outlined earlier, which is by definition the statistical function of interest, e.g. see van der Vaart (1998). The large sample proof of this concept is given by the following theorem:

Theorem 1.

Under $H_0:T(F_0)={\theta }_0$ and as $n\to \infty ,$ $\sqrt{n}(\widehat{\alpha }-1,\widehat{\beta }-1)$ has a centered bivariate normal distribution with variance-covariance matrix $B^{-1}\Sigma B^{-1}$, where B is the standard maximum likelihood based information matrix associated with the Kumaraswamy density, $k_{\alpha ,\beta }=d{K}_{\alpha ,\beta }$, and Σ is the variance-covariance matrix of a 2-dimensional random vector whose components are given by (2.7) $$\partial \hspace{0.17em}\text{log}\hspace{0.17em}{f}_{\alpha ,\beta }\left(X_{i:n}\right)/\partial \alpha +W_{\alpha }\left(X_{i:n}\right),$$(2.7) (2.8) (2.9) $$\partial \hspace{0.17em}\text{log}\hspace{0.17em}{f}_{\alpha ,\beta }\left(X_{i:n}\right)/\partial \beta +W_{\beta }\left(X_{i:n}\right),$$(2.8) (2.9) where (2.10) $$W_{\alpha }\left(X_{i:n}\right)=\int I_{\left(F_{\alpha ,\beta }\right(X_{i:n})<u)}\frac{{\partial }^2}{\partial \alpha \partial u}\text{log}\hspace{0.17em}\left(f_{\alpha ,\beta }\right(u\left)\right)d{F}_{\alpha ,\beta }\left(u\right),$$(2.10) (2.11) $$W_{\beta }\left(X_{i:n}\right)=\int I_{\left(F_{\alpha ,\beta }\right(X_{i:n})<u)}\frac{{\partial }^2}{\partial \beta \partial u}\text{log}\hspace{0.17em}\left(f_{\alpha ,\beta }\right(u\left)\right)d{F}_{\alpha ,\beta }\left(u\right).$$(2.11)

Proof.

The technical details have been worked out in an elegant fashion for the case of a semi-parametric copula model with marginal distribution functions estimated by the empirical distribution function estimator. The result in Theorem 1 follows directly from the theoretical developments used in the copula approach in Section 4 of the copula paper Genest, Ghoudi, and Rivest (1995) by simply replacing the multivariate copula function with the univariate beta density, which is essentially a special case of the higher dimension copula model. Estimates of the variance-covariance matrix $B^{-1}\Sigma B^{-1}$ are not as straightforward to obtain and we recommend bootstrap resampling for this purpose.

3. Simulation results

For our simulation study we focused on the trimmed mean with known statistical functional given as $T_{\gamma }(F)=\frac{1}{1-2\gamma }{\int }_{Q(\gamma )}^{Q(1-\gamma )}xdF=\frac{1}{1-2\gamma }{\int }_{\gamma }^{1-\gamma }Q(u)\mathit{du},$ $0<\gamma <1/2.$ Similar results in terms of behavior hold for moment estimators and kernel estimators and are not presented here. We centered our simulation study for the trimmed mean on the hypothesis test $H_0:T(F_0)=0$ versus $H_1:T(F_0)>0$ at type I error rate 0.05 for symmetric distributions with trimming proportions $\alpha =0,0.1,0.2$ and samples of size n = 10, 20, 50. For each simulation result we utilized 1000 Monte Carlo resamples with the number of bootstrap resamples set to B = 250. For the exponential distribution we tested $H_0:T(F_0)=T_{\gamma }(F)$ versus $H_1:T(F_0)>T_{\gamma }(F)$ with shifted exponential alternatives. For the γ trimmed mean we can simplify (2.3) such that (3.1) $$T({\widehat{F}}_0;\alpha ,\beta )=\sum _{i=k+1}^{n-k}g\left(X_{i:n}\right)K_{\alpha ,\beta }\left(\frac{i}{n-2k}\right)-K_{\alpha ,\beta }\left(\frac{i-1}{n-2k}\right),$$(3.1) where $k=[n\tau ].$ For power examinations we used shift alternatives of $\delta =0.5$ and δ = 1 different from H₀. The Type I error was estimated from the 1000 Monte Carlo resamples as the number of times p_boot was less than the Type error of 0.05.

The results of our simulation study are presented in Table 1. We see that the type I error is controlled at the nominal level and that fluctuations about that level are primarily due to simulation error. The power is monotone in increase δ and n. As compared with an optimal scenario such as a t-test under normality with $\delta =0.5$ and $\alpha =0.0$ the one-sample t-test has power of 0.427, 0.695 and 0.967 as compared to the warping powers of 0.369, 0.656 and 0.966 for samples of size n of 10, 20 and 50, which yields relative efficiencies of 86.4%, 94.4% and 99.8%,respectively. By comparison the empirical likelihood approach yields powers of 0.414, 0.640 and 0.960 for samples of size n of 10, 20 and 50 under the same scenarios. However, it should be noted that the Type I error control for the empirical likelihood approach was inflated at 0.084 at $\delta 0$ and n = 10, hence in turn the power value of 0.414 is inflated at this same sample size due to a much higher than desired Type I error level as compared to the warping approach.

Table 1. Type 1 error set to 0.05(δ = 0), power ($\delta >0$) at trimming proportions $\alpha =0,0.1$ and 0.2.

	Normal			Logistic			Exponential
n	10	20	50	10	20	50	10	20	50
δ = 0
$\alpha =0.0$	0.060	0.048	0.059	0.056	0.065	0.051	0.049	0.049	0.049
$\alpha =0.1$	0.043	0.045	0.051	0.038	0.053	0.056	0.042	0.056	0.040
$\alpha =0.2$	0.037	0.061	0.048	0.039	0.045	0.049	0.031	0.049	0.042
$\delta =0.5$
$\alpha =0.0$	0.369	0.656	0.966	0.210	0.319	0.637	0.641	0.872	0.999
$\alpha =0.1$	0.314	0.624	0.946	0.199	0.324	0.637	0.695	0.961	0.999
$\alpha =0.2$	0.0167	0.576	0.949	0.105	0.286	0.625	0.349	0.949	0.999
δ = 1
$\alpha =0.0$	0.826	0.997	0.999	0.419	0.714	0.943	0.999	0.999	0.999
$\alpha =0.1$	0.720	0.987	0.999	0.448	0.727	0.985	0.947	0.999	0.999
$\alpha =0.2$	0.235	0.941	0.999	0.202	0.694	0.981	0.182	0.970	0.999

Acknowledgments

This work was supported by Roswell Park Cancer Institute and National Cancer Institute (NCI) grant P30CA016056, NRG Oncology Statistical and Data Management Center grant U10CA180822 and IOTN Moonshot grant U24CA232979-01. We wish to thank the reviewers for their thoughtful comments, which led to an improved version of this work.