{'fraction_non_alphanumeric': 0.039125983748160474, 'fraction_numerical': 0.020026873120481157, 'mean_word_length': 4.576984834968778, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Semiparametric Approach for Regression with Covariate Subject to Limit of Detection

We consider generalized linear regression analysis with left-censored covariate due to the lower limit of detection. Complete case analysis by eliminating observations with values below limit of detection yields valid estimates for regression coefficients, but loses efficiency; substitution methods are biased; maximum likelihood method relies on parametric models for the unobservable tail probability distribution of such covariate, thus may suffer from model misspecification. To obtain robust and more efficient results, we propose a semiparametric likelihood-based approach for the estimation of regression parameters using an accelerated failure time model for the covariate subject to limit of detection. A two-stage estimation procedure is considered, where the conditional distribution of the covariate with limit of detection given other variables is estimated prior to maximizing the likelihood function. The proposed method outperforms the complete case analysis and the substitution methods as well in simulation studies. Technical conditions for desirable asymptotic properties are provided.


Introduction
Detection limit is a threshold below which measured values are not considered significantly different from background noise (Helsel, 2005). Hence, values measured below this threshold are unreliable. In environmental epidemiology, particularly exposure analysis, when exposure levels are low, measurement of chemicals has a large percentage falling below the limit of detection due to inadequate instrument sensitivity. This is a common issue in the National Health and Nutrition Examination Survey (Crainiceanu et al., 2008), where many exposure variables have large proportions of measurements below their limits of detection. For example, 27.8% of the urine arsenobetaine measures are below its limit of detection at 0.4 µg/l (Caldwell et al., 2009). In the Diabetes Prevention Program, 66 of the 301 eligible participants had their testosterone levels below the detection limit of 8.0 ng/dl (Kim et al., 2012). In an analysis for the Michigan Bone Health and Metabolism Study, up to 66% of the 50 study participants had anti-Mullerian hormone below the limit of detection at 0.05ng/ml (Sowers et al., 2008). For illustrative purpose in this article, we consider an analysis for the National Health and Nutrition Examination Survey, which examines the relationship between arsenic exposure and the prevalence of type 2 diabetes .
A variable with limit of detection can be either a response variable or a covariate in regression analysis. We focus on the latter in this article. Although many ad hoc methods have been implemented in practice, more appropriate statistical methods for regression models with a covariate subject to limit of detection are yet to be thoroughly studied (Schisterman and Little, 2010). The complete case analysis, simply eliminating observations with values below limit of detection, yields consistent estimates of the regression coefficients (Nie et al., 2010;Little and Rubin, 2002), but loses efficiency. Substitution methods are frequently used in epidemiology studies, where the values of covariate Z below limit of detection, denoted by L, are substituted by L, or L/ √ 2, or zero, see for example, Schisterman et al. (2006), Hornung and Reed (1990), Moulton et al. (2002), Koru-Sengul et al. (2011), Kroger et al. (2009, Boomsma et al. (2009), Bloom et al. (2008), and Gollenberg et al. (2010) among many others. These methods are easily implementable, but found to be inappropriate and can yield large biases (Helsel, 2006;Nie et al., 2010). Richardson and Ciampi (2003) proposed to replace the values below limit of detection L with E(Z|Z < L), which is obtained from an assumed known distribution of Z. There are two issues with this approach, however, (i) the distributional assumption is not verifiable; (ii) even if E(Z|Z < L) is correctly specified, the method only leads to consistent estimates in linear regression when the covariate subject to limit of detection Z is independent of all the other covariates.
Another widely used method is the maximum likelihood estimation based on a parametric distributional assumption to the unobservable tail probability of the covariate that is subject to limit of detection. For examples, Nie et al. (2010) and Arunajadai and Rauh (2012) considered the linear regression based on a normal and a generalized gamma distribution for the covariate subject to limit of detection, respectively; Cole et al. (2009), Wu et al. (2012 and Albert et al. (2010) considered parametric maximum likelihood approach under the generalized linear regression; D' Angelo and Weissfeld (2008) applied this approach to the Cox regression. In practice, however, the underlying covariate distribution is unknown. The test of the parametric assumption to the unobservable tail probability of the covariate that is subject to limit of detection is usually unavailable because there is no information/observation below the detection limit. Both Lynn (2001) and Nie et al. (2010) noted that a parametric assumption can yield large bias if misspecified and argued that such an approach should not be attempted. Nie et al. (2010) recommended the complete case analysis despite the fact that simply dropping data below the limit of detection can lose a significant amount of information.
To obtain more efficient and yet robust results, we propose a semiparamet-ric likelihood-based approach to fit generalized linear models with covariate subject to limit of detection. The tail distribution of the covariate beyond its limit of detection is estimated from a semiparametric accelerated failure model, conditional on all the fully observed covariates. Model checking can be done using martingale residuals for semiparametric accelerated failure time models. The proposed method is shown to be consistent and asymptotically normal, and outperforms existing methods in simulations. The proof of the asymptotic properties relies heavily on empirical process theory, which is provided in the online Supplementary Material.

A semiparametric approach
For a single observation, denote the response variable by Y , the covariate subject to limit of detection by Z, and the fully observed covariates by X = (X 1 , . . . , X p ) ′ , where p is the number of fully observed covariates. For simplicity, we only consider one covariate that is subject to limit of detection. Consider a generalized linear model with where g is the link function, D ′ θ is the linear predictor with D = (1, X ′ , Z) ′ and θ = (β ′ , γ) ′ , here β is a (p + 1)-dimensional vector and γ is a scalar. The variance of Y , typically a function of the mean, is denoted by We consider the exponential dispersion family in the natural form (Agresti, 2002;McCullagh and Nelder , 1989) given (Z, X) where φ is the dispersion parameter and ̟ is the natural parameter. We have µ = E(Y ) =ḃ(̟), and var(Y ) =b(̟)a(φ), whereḃ is the first derivative of b andb is the second derivative of b.
The actual value of Z is not observable when Z < L, where the constant L denotes the limit of detection, which is an example of left-censoring. In practice Z is a concentration measure of certain substance and thus non-negative. Consider a monotone decreasing transformation h that yields , then T is observed; otherwise T is right-censored by C. We denote the observed value by V = min(T, C) and the censoring indicator by ∆ = I(T ≤ C).
The proposed methodology works for a broad family of link functions defined by the regularity conditions given in the Appendix. For notational simplicity, we present the main material using canonical link function g, where g = (ḃ) −1 . Then, when T is observed, model (2) becomes Denote the conditional cumulative distribution function of T given X by F 1 (t|X) with density f 1 (t|X). The likelihood function for the observed data (V, ∆, Y, X) can be factorized into where f denotes the joint density of (V, ∆, Y, X), f 2 denotes conditional density of (V, ∆) given (Y, X), f 3 denotes conditional density of Y given X, and f 4 denotes marginal density of X. Going through conditional arguments using the Bayes' rule and dropping f 4 (X), we obtain the likelihood function where only f θ,φ contains the parameter of interest θ, whereas f 1 is a nuisance parameter in addition to φ.
There are two parts in (4) 1−∆ for subject with covariate below limit of detection. Complete case analysis is only based on the first part and, although it yields a consistent estimate of θ, it clearly loses efficiency. We see from the second part of (4) that the efficiency gain comparing to the complete case analysis depends on how well we can recover the right tail of the conditional distribution F 1 (t|X) beyond C. Parametric models for F 1 (t|X) are often considered in the literature, see Nie et al. (2010), but it may suffer from model misspecification. The nonparametric method degenerates to the complete case analysis because there is no actual observation beyond censoring time C. We consider a semiparametric approach that allows reliable extrapolation beyond C and is robust against any parametric assumption. Among all the commonly used semiparametric models for right-censored data, only the accelerated failure time model allows extrapolation beyond C, and model checking can be done by visualizing the cumulative sums of the martingale-based residuals (Lin et al., 1993(Lin et al., , 1996Peng and Fine, 2006). We hence propose a semiparametric accelerated failure time model for the transformed covariate subject to limit of detection given by where ς follows some unknown distribution, denoted by η, and is independent of X. We only consider a fixed h for T in this article. More flexible transformation, for example, the Box-Cox transformation (Box and Cox, 1964;Foster et al., 2001;Cai et al., 2005), is worth further investigation. Note that X appears in both models (1) and (5), but it may refer to different forms of covariates in these models. For example, X 1 is a covariate in (1) whereas X 2 1 is a covariate in (5). We use the same X to denote all fully observed covariates for notational simplicity. The log-likelihood function then becomes where τ is a truncation time at the residual scale defined in Condition 4 in the Appendix.

The pseudo-likelihood method
The log likelihood function (6) involves an unknown distribution function η and its derivative, hence a semiparametric maximum likelihood estimation, if it exists, can be complicated. We propose a tractable two-stage pseudo-likelihood approach in which the nuisance parameters (φ, α, η) are estimated in stage 1, and the parameter of interest θ is then estimated by maximizing the data version of (6) in stage 2 with nuisance parameters replaced by their estimators obtained in stage 1. Details are given below: Stage 1. Nuisance parameter estimation. Dispersion parameter φ is estimated by the complete case analysis of the generalized linear model (2); the accelerated failure time model regression coefficient α is estimated by either the rank based methods, see Wei et al. (1990), Jin et al. (2003, Nan et al. (2009) or the sieve maximum likelihood method, see ; and the accelerated failure time model error distribution η is estimated by the Kaplan-Meier estimator from the censored residuals.
The complete case analysis can be done by any standard statistical package for generalized linear models. The rank based estimates for the accelerated failure time model usually are obtained by using linear programming. The R package "rankreg" (Zhou, 2006), now archived by CRAN, can be implemented for small to moderate sample sizes because it solves a linear programming problem of size n 2 . This is the method we implemented in simulations and the arsenic exposure data example. An alternative approach is to modify the Newton algorithm for solving the discrete rank based estimating equation (Yu and Nan, 2006). Standard Newton-Raphson algorithm can be implemented to obtain the sieve maximum likelihood estimates  for the accelerated failure time model when the sample size is large.
Stage 2. Pseudo-likelihood estimation of θ. Replacing (φ, α, η) by their Stage 1 estimates (φ n ,α n ,η n,αn ) in the log likelihood function yields the following log pseudo-likelihood function for a random sample of n observations: Note that the term ∆ logη(T ) in (6) is dropped because it does not involve θ.
We maximize (7) by setting its derivative to zero and then solving the equation to obtain the pseudo-likelihood estimatorθ n . This is implemented by using standard Newton-Raphson algorithm with the initial value obtained from the complete case analysis in Stage 1.
Sinceθ n is obtained by solving an estimating equation, its asymptotic properties can be obtained from Z-estimation theory. It can be shown that all the estimates obtained in Stage 1 have desirable statistical properties for Stage 2 estimation. In particular,φ n obtained from the complete case analysis is n 1/2consistent by Little and Rubin (2002);α n is n 1/2 -consistent by Nan et al. (2009) or ; andη n,αn is also n 1/2 -consistent in a finite interval, and its proof is provided in the online Supplementary Material.
A set of regularity conditions is introduced in the Appendix. Some conditions are commonly assumed for the accelerated failure time models, and other conditions are for the generalized linear models, which are easily verifiable for linear, logistic and Poisson regression models. We then have the following asymptotic results forθ n . Then for the two-stage pseudo-likelihood estimatorθ n satisfying Ψθ n,n (φ n ,α n ,η n,αn ) = 0, we have: (i)θ n converges in outer probability to θ 0 , and (ii) n 1/2 (θ n −θ 0 ) converges weakly to a mean zero normal random variable with variance A −1 BA −1 , where A and B are provided in the online Supplementary Material.
Because the asymptotic variance ofθ n has a very complicated expression that prohibits the direct calculation of its estimate from observed data, we recommend using the bootstrap variance estimator.
The proof of Theorem 4.1 is based on the general Z-estimation theory of . Define a deterministic function and denote the true values of (φ, α, η) by (φ 0 , α 0 , η 0 ). We can show that Ψ θ,n (φ n ,α n ,η n,αn ) converges uniformly to Ψ θ (φ 0 , α 0 , η 0 ) as n → ∞. Then the consistency is achieved given that θ 0 is the unique solution of Ψ θ (φ 0 , α 0 , η 0 ) = 0. The asymptotic normality is derived by showing the asymptotic linear representation of n 1/2 (θ n − θ 0 ). The detailed proofs rely heavily on empirical process theory and can be found in the online Supplementary Material, where we only provide the analytic form of the asymptotic variance for the Gehan weighted estimate of α. The analytic forms of the asymptotic variance for other rank based estimates and the sieve maximum likelihood estimates can be obtained similarly.

Simulations
We conduct simulations to investigate the finite sample performance of the proposed method. Simulation data sets are generated from the generalized linear model where β 0 = −1, β 1 = 0.5, β 2 = −1, γ = 2, and g is chosen to be the canonical link function for normal, bernoulli and poisson distributions, respectively. The normal error variance is chosen to be 1 for the linear regression model. The three covariates are: X 1 ∼ Bernoulli(0.5), X 2 is normal with mean 1 and standard deviation 1 truncated at ±3, and Z = exp (−T ) is generated from the following linear model where α 0 = 0.25, α 1 = 0.25, α 2 = −0.5, ς ∼ 0.5N (0, 1/8 2 ) + 0.5N (0.5, 1/10 2 ). The limit of detection L for covariate Z is chosen to yield 30% right censoring for T . We simulate 1000 replications for each scenario and compare the proposed method with full data analysis, complete case analysis, and four different substitution methods. The full data analysis represents the case of no limit of detection, which serves as a benchmark. We conduct simulations with two different sample sizes: 200 and 400. The four substitution methods for Z < L are: (i) replacing Z by L, (ii) replacing Z by L/ √ 2, (iii) replacing Z by zero, and (iv) replacing Z by E(Z|Z < L). For the proposed two-stage method, we report the 90% and 95% coverage proportions for which the variances are obtained from 200 bootstrap samples. Empirical variances obtained from 1000 independent data sets are provided only for the full data analysis and the two valid methods for the limit of detection problem. The results are presented in Tables 1-3.
The results suggest that all the substitution methods yield biased estimates, including substituting Z by the true E(Z|Z < L). The biases for the proposed two-stage method are minimal, which are comparable to both the full data analysis and the complete case analysis. Clearly, the proposed method is much more efficient than the complete case analysis, and the bootstrap method performs well in estimating the variance, which yields reasonable coverage rate of the confidence intervals for all considered sample sizes.
Different censoring rates varying from 10% to 60% were considered in additional simulations. Both the proposed two-stage method and the complete case analysis yield unbiased estimates in all simulation scenarios. The efficiency gain of the two-stage method comparing to complete case analysis increases as the percentage of censoring increases. All the substitution methods yield biased estimates and the biases increase as the censoring rate increases as well. The detailed simulation results are not provided here.
The consistency of the estimates from the two-stage method depends on correctly specifying the semiparametric accelerated failure time model (5). Although the assumption for model (5) is much less restrictive than any parametric model, model checking should be done before applying the proposed two-stage method. Though less severe than ad hoc substitution methods, misspecification of the accelerated failure time model also yields biased results, and the bias increases as the severity of misspecification of (5) grows. Again, the detailed simulation results are not provided here.

The National Health and Nutrition Examination Survey
We consider the National Health and Nutrition Examination Survey 2003-2004 as an illustrative example for the regression with a covariate subject to limit of detection. In particular, we focus on the effect of left-censored urine arsenobetaine on the prevalence of type 2 diabetes, see Navas-Acien et al. (2008) National Health and Nutrition Examination Survey, conducted by the US National Center for Health Statistics, used a complex multistage sampling design to obtain a representative sample of the civilian noninstitutionalized individuals within the US population. The data set contains a subsample of 1542 study participants with arsenic measurements of the National Health and Nutrition Examination Survey 2003-2004. For each participant in this subsample, total urine arsenic and arsenic species including arsenobetaine were collected for      arsenic analysis. The limit of detection for total urine arsenic and urine arsenobetaine were 0.6 and 0.4 µg/l, respectively. The percentage of study participants with levels below the limit of detection were 1.3% for total urine arsenic and 27.8% for urine arsenobetaine (Caldwell et al., 2009). Navas-Acien et al. (2008) found that total urine arsenic was associated with increased prevalence of type 2 diabetes, hence it was adjusted in the analysis, and 19 participants with total urine arsenic below limit of detection were dropped from the study. The urine creatinine level that was used to account for urine dilution in spot urine samples was fully observed and was adjusted in the analysis as well. We further excluded 23 participants with missing values in other variables of interest.
For illustration purpose, we focus on the male participants, where 24.1% of 730 subjects had urine arsenobetaine below limit of detection. Age, race/ethnicity, body mass index, the logarithm of total urine arsenic and the logarithm of urine creatinine level are used as covariates to fit the accelerated failure time model for − log(arsenobetaine). All of them are significant with p-values < 0.0001 except one dummy variable for race. Dotted lines in Figure 1 are the plots of 50 realizations from the distributions of the score processes. The observed score processes are presented with solid lines which randomly fluctuate around zero. From Figure 1 we see that the considered accelerated failure time model for urine arsenobetaine fits the data reasonably well, with respective goodness-of-fit empirical p-values of 0.686, 0.104, 0.706, 0.782, 0.68, 0.794, 0.646 and 0.834 for age, log total urine arsenic, log urine creatinine, body mass index and race that contains five ethnic groups based on 500 simulated martingale residual score processes. The response of interest is the status of type 2 diabetes. Thus, a logistic regression is considered with age, body mass index, log total urine arsenic, log urine creatinine and log urine arsenobetaine as the covariates of interest, whereas race is not significantly associated. Table  4 shows the regression coefficient estimates, where we see that the proposed two-stage method yields similar point estimates with smaller variances/p-values comparing to the complete case analysis, indicating the efficiency gain of the proposed method. Substitution methods with values below limit of detection replaced by L or L/ √ 2 yield quite different point estimates comparing to the complete case analysis and the two-stage method, thus put the significant results of log(arsenobetaine) in question. The effect of log(arsenobetaine) is clearly biased when the urine arsenobetaine values below limit of detection are replaced by zero.

Disscussion
The estimates from the proposed method are consistent and asymptotically normal under much less restrictive assumptions than parametric approaches. Another advantage of the two-stage method over the parametric method is that the model checking tools are available. All the substitution methods could yield large bias including replacing the values below the limit of detection by the true E(Z|Z < L). The consistency and efficiency gain of the proposed two-  stage method rely on the correctly specified accelerated failure time model. We suggest to fit the accelerated failure time model for the covariate subject to limit of detection before applying the two-stage method. If the data fits the model reasonably well and there are at least some fully observed covariates significantly associated with the covariate subject to limit of detection, then the proposed method is recommended. Otherwise, we suggest to use the complete case analysis. The amount of efficiency gain of the proposed two-stage method depends on how far we can estimate F (t|X) reasonably well beyond the limit of detection. We assume some finite value τ for residuals in this article. In practice, the upper limit of the integral in the pseudo-likelihood function can go as far as the largest observed residual in the fitted accelerated failure time model, which is max i (T i − X ′ iα n ); and theoretically, this upper limit is ∞ when the support of X ′ α 0 is unbounded. In the latter case, it can be shown thatF n converges to F on the entire real line with a polynomial rate at n −1/8 , see Lai and Ying (1991) and Ding and Nan (2014), and we may still obtain consistent estimates for the parameters of interest. The asymptotic normality, however, will be difficult to show.
We only consider the case with one covariate subject to limit of detection in this article for simplicity. Regression with multiple covariates subject to limits of detection may occur in practice. Parametric models have been considered for such problems (May et al., 2011;D'Angelo and Weissfeld, 2008). To achieve robust results, the proposed semiparametric approach can be generalized to tackle the problem with multiple covariates subject to limits of detection. The critical step is to provide an valid nonparametric estimate for the multivariate survival function, for which available methods include Dabrowska (1988), Prentice and Cai (1992), van der Laan (1996), and Prentice and Moodie (2004). The constant limit of detection assumption considered in this article, though commonly seen in practice, is also for notational simplicity, and can be relaxed to cases with random limit of detection.
Limit of detection issue is a special missing data problem. Multiple imputation (Little and Rubin, 2002) may be considered as an alternative method if the tail distribution of the covariate subject to limit of detection conditional on all other variables, including the response variable, can be estimated reasonably well. This research is ongoing and will be presented elsewhere.
When the original lab measurements are available, which may or may not be below limit of detection, Murphy et al. (2010) and Buck Louis et al. (2012) directed used the machine-read lab values in their analysis to avoid the potential biases caused by substitution. More appropriate analysis should be treating these lab data as error prone measurements (Guo and Little, 2010), hence methods dealing with measurement errors would apply.
The proposed two-stage method is ready to be generalized to other regression models that have a likelihood function to work with, for example, the Cox regression model and mixed effects model. Extra care will be needed due to special features of these models. For example, handling the nonparametric baseline function in the Cox model in the context of limit of detection can be delicate. These are interesting topics worth further investigation.

Supplementary material
The online Supplementary Material contains the detailed proof of Theorem 1.

Acknowledgement
We thank Professor Sung Kyun Park for the National Health and Nutrition Examination Survey arsenic exposure data. We are grateful to the editor, associate editor and a referee for helpful and constructive comments. The research was partially supported by several grants from the U.S. National Science Foundation and the National Institute of Health.

Appendix: regularity conditions
Denote the sample space of response variable Y by Y, the sample space of covariate X by X , the parameter space of θ by Θ, the parameter space of α by A, and the parameter space of η by H. In addition to the assumptions of bounded support for (X, Z) and compact parameter spaces Θ and A, we provide a set of regularity conditions for Theorem 4.1 in the following.

Condition 7.b(·) is a bounded Lipschitz function.
Condition 8. There exist constants C i , i = 1, . . . , 5, such that for any constant U < ∞, sup y∈Y,θ∈Θ,|1/a(φ)|≤l,x∈X ,t∈ [C,U] Condition 9. There exist constants δ 1 > 0 and δ 2 > 0, such that Remark: Condition 1 is for the consistency, which may be unnecessarily strong for the proposed two-stage method. Direct calculation yieldṡ , which is negative definite. ThusΨ θ , a continuous matrix with θ, is also negative definite in a neighborhood of θ 0 , which guarantees that θ 0 is the unique solution of Ψ θ (φ 0 , α 0 , η 0 ) = 0 in a neighborhood of θ 0 . The initial value we use in the Newton-Raphson algorithm for solving Ψ θ,n (φ n ,α n ,η n,αn ) = 0 is obtained from the complete case analysis, which is consistent, thus the solution of the proposed two-stage method should also be consistent. Condition 2 holds for many commonly used transformations, for example, h(t) = exp(−t) and polynomial functions. Condition 3 and 4 are usual assumptions for accelerated failure time models (Tsiatis, 1990;Nan et al., 2009). Conditions 5-8 automatically hold for common generalized linear models, for example, linear, logistic or poisson regression.
Condition 9 is mainly for technical convenience. One way to obtain Condition 9 might be to truncate response variable Y such that |Y | ≤ M < ∞ for a large constant M and to further truncate the residual in the accelerated failure time model with some constant τ ′ < τ . In our simulations, however, we do not implement such truncations but still obtain satisfactory results.  1 General Z-estimation theory The proof of Theorem 1 in the main text is based on the general Z-estimation theory of , which is provided in the following Lemmas 1.1 and 1.2 for our problem setting. Detailed discussion and proofs of these two lemmas can be found in . Let | · | be the Euclidian norm We use P * to denote outer probability, which is defined as P * (A) = inf{pr(B) : B ⊃ A, B ∈ B} for any subset A of Ω in a probability space (Ω, B, P ).
Thenθ n is n 1/2 -consistent and further we have (3)

Technical lemmas
Now we provide technical preparations for the proof of Theorem 1, some of which are from Ying Ding's 2010 University of Michigan Ph.D. thesis. We adopt the empirical process notation of van der vaart and Wellner (1996).
where ⌈x⌉ is the smallest integer that is greater than or equal to x. Furthermore, denote t 0 = 0 and t ⌈1/ε⌉ = +∞. For i = 1, · · · , ⌈1/ε⌉, define brackets where K 2 = 2ξ. From Theorem 2.14.9 in van der vaart and Wellner (1996), we have for some constant D 1 depends on K 1 and constant D 2 depends on K 2 . We now where F X (x) is the distribution function of X, from Condition 3 we have Since it can be shown that m 3 (α 0 , ǫ 0 ; X, ∆) has finite second moment, we have is Donsker.
Lemma 2.3. Suppose X and A be the bounded covariate and parameter spaces. Let H be a collection of distribution functions satisfying Condition 3. We have F = {η(t − x ′ α), t ∈ T ⊂ R, x ∈ X , α ∈ A, η ∈ H} is Donsker.
Proof. From Condition 9, we have { τ C−x ′ α f θ,φ (y | t + x ′ α, x)dη(t)} bounded away from zero. From Section 2.10.2 of van der vaart and Wellner (1996), we only need to show that both the numerator and denominator in (14) belong to Donsker classes. By integration by parts, we have In the above,ḣ(·) is Lipschitz by Condition 2 and f θ,φ (y | t + x ′ α, x) is Lipschitz function for θ, φ and α by Conditions 2, 5 and 8, thus both belong to Donsker classes by Theorem 2.10.6 in van der vaart and Wellner (1996). By Lemma 2.3 we know that {η(C − x ′ α)} is Donsker. Since the class of indicator functions of half spaces is a VC-class, see e.g. Exercise 9 on page 151 and Exercise 14 on page 152 in van der vaart and Wellner (1996), thus the set of functions {1(t ≥ C − x ′ α)} is a Donsker class. By Theorem 2.10.3 in van der vaart and Wellner (1996), the permanence of the Donsker property for the closure of the convex hull, we have Hence the denominator of (14) belongs to a Donsker class.
Proof. Again, the proof is based on direct calculation.
Note that for other rank based estimates of α, m 3 in B is the corresponding influence function with different forms; For the sieve maximum likelihood estimates , m 3 is the efficient influence function (Ritov and Wellner, 1988). It is clearly seen that the analytic form of the asymptotic variance is too complicated to be useful for the asymptotic variance estimation. Hence in our numerical studies we use bootstrap to obtain the variance estimator.