Sampling Importance Resampling Algorithm with Nonignorable Missing Response Variable Based on Smoothed Quantile Regression

: The presence of nonignorable missing response variables often leads to complex conditional distribution patterns that cannot be effectively captured through mean regression. In contrast, quantile regression offers valuable insights into the conditional distribution. Consequently, this article places emphasis on the quantile regression approach to address nonrandom missing data. Taking inspiration from fractional imputation, this paper proposes a novel smoothed quantile regression estimation equation based on a sampling importance resampling (SIR) algorithm instead of nonparametric kernel regression methods. Additionally, we present an augmented inverse probability weighting (AIPW) smoothed quantile regression estimation equation to reduce the inﬂuence of potential misspeciﬁcation in a working model. The consistency and asymptotic normality of the empirical likelihood estimators corresponding to the above estimating equations are proven under the assumption of a correctly speciﬁed parameter working model. Furthermore, we demonstrate that the AIPW estimation equation converges to an IPW estimation equation when a parameter working model is misspeciﬁed, thus illustrating the robustness of the AIPW estimation approach. Through numerical simulations, we examine the ﬁnite sample properties of the proposed method when the working models are both correctly speciﬁed and misspeciﬁed. Furthermore, we apply the proposed method to analyze HIV— CD4 data, thereby exploring variations in treatment effects and the inﬂuence of other covariates across different quantiles.


Introduction
Missing data analysis has gained significant attention in recent years.To analyze missing data, it is crucial to understand the response mechanism that leads to missing data.If the missingness of the variable of interest is conditionally independent of that variable, the response mechanism is considered to be random or ignorable.Otherwise, the response mechanism is considered to be nonrandom or nonignorable.Dealing with nonrandom missing data presents greater challenges, which are evident in two aspects: Firstly, the assumed response model cannot be validated solely based on observed data; secondly, the model parameters may be unidentifiable.
To obtain meaningful inferences from incomplete data with nonrandom missingness, it is necessary to satisfy a set of identifying conditions [1,2].Moreover, the accuracy of the methods based on parameter models is greatly influenced by the correct specification of the assumed parameter model [3].Consequently, researchers aim to impose weaker model assumptions on the response mechanism to achieve robust results.The semiparametric response model was initially considered by Kim and Yu [4], but their proposed method necessitated a validation sample to estimate the model parameters.Similarly, Shao and Wang [5] examined the same semiparametric exponential tilting model and proposed a parameter estimation approach based on calibration estimation equations.A comprehensive review of parameter estimation methods for nonrandom missing data is provided by Kim and Shao [6].
Quantile regression, introduced by Koenker and Bassett [7], has become a widely used statistical analysis tool.It offers more adaptability and flexibility compared to mean regression.Notably, quantile regression does not require the assumption of error term distribution and demonstrates robustness against heavy-tailed errors and outliers.Furthermore, by considering regressions at different quantiles of the response variable, quantile regression enables the assessment of covariate effects at various quantiles and yields a more comprehensive understanding of the conditional distribution.However, there is a scarcity of literature on quantile regression for nonrandom missing data.
The nonsmoothness of the check function for standard quantile estimators makes it impossible to directly estimate the asymptotic covariance matrix [8].As a result, the existing theoretical results for nonrandom missing mean regression cannot be directly extended to quantile regression.
The idea of smoothing nondifferentiable objective functions can be traced back to Horwitz [9], while Whang [10] introduced the smoothed empirical likelihood approach for quantile regression.Luo et al. [11] extended the aforementioned method to analyze data with random missingness; Zhang and Wang [12] further expanded it to handle cases of nonignorable missingness.
However, on the one hand, this method relies on the assumption of a parametric propensity missingness model, which introduces the risk of model misspecification.On the other hand, this method corrects estimation biases caused by missing data through inverse probability weighting but may not fully utilize the information from incomplete observations.Regarding nonrandom missing data, previous studies have addressed the issue in different regression settings.Specifically, Niu et al. [13] and Bindele and Zhao [14] focused on estimation equation imputation in linear regression and rank regression, respectively.In the context of quantile regression, Chen et al. [15] introduced three missing quantile regression estimation equations: inverse probability weighting, estimation equation imputation, and an enhanced approach combining both methods.It is important to note that these studies assume a response mechanism with random missingness.
Moreover, the existing literature commonly utilizes kernel estimation methods [16] to estimate the conditional means involved in the imputation estimation equation.However, when the dimension of the covariates is high, the kernel estimation results can become unstable.To overcome the curse of dimensionality associated with multivariate nonparametric kernel estimation, Kim [17] proposed a parametric fractional imputation method for handling missing data.Additionally, Riddles et al. [18] extended this method to address the scenario of nonignorable missing data.They developed an EM algorithm based on a parameter working model derived from observed data and incorporated the parametric fractional imputation (FI) method.Nevertheless, these approaches heavily rely on parameter-based response models, which renders them sensitive to model misspecification.Furthermore, the likelihood-based EM algorithm is not directly applicable to quantile regression.
Utilizing estimation equations, Paik and Larsen [19] incorporated a working model for observed data and employed a sampling importance resampling (SIR) algorithm to estimate the missing data and corresponding estimation equations.Building upon this, Wang et al. [20] and Song et al. [21] extended the logistic response model used in the aforementioned approaches to a develop a semiparametric exponential tilt model.However, in the absence of knowledge about the tilting parameter, these methods relied on validation samples.
In this paper, we propose a smoothed empirical likelihood approach for imputing quantile regression estimation equations with nonignorable missing data based on a semiparametric response model.The novel estimation equation guarantees the second-order differentiability of the objective function with respect to the parameter vector.The imputed values for the missing data were derived from a parameter working model and obtained using sampling importance resampling.
Although imputation estimation equations applying information from missing data compared to IPW estimation equations can enhance estimation efficiency, both theoretical and numerical experiments have shown that imputation estimation equations are sensitive to misspecification of the parameter working model.Therefore, to mitigate the impact of misspecification in the working model, this paper further proposes the AIPW smoothed quantile regression estimation equation.It is demonstrated that, when the working model is correctly specified, the asymptotic variance of the AIPW estimation equation shares the same form as the asymptotic variance of the nonparametric model estimator.Furthermore, even when the working model is misspecified, the estimates remain consistent.
The remaining sections of this paper are organized as follows.Section 2 establishes the semi-parametric response model and the AIPW quantile regression estimation equation, along with the algorithmic procedure for estimating the skewness parameter and quantile regression coefficients using importance resampling.Section 3 presents the large sample properties of the parameter estimators.Section 4 demonstrates the finite sample properties of the estimators through numerical simulations.Section 5 applies the proposed methodology to analyze the HIV-CD4 dataset.

Proposed Method
Consider a linear quantile regression model as follows: where Y i is the response variable, Z i is a fully observed q-dimensional covariate vector, θ τ represents the unknown regression coefficient vector, i denotes the random error term satisfying P( i ≤ 0|Z i ) = τ, τ ∈ (0, 1), and the i values are mutually independent.In the subsequent discussion, we will abbreviate θ τ as θ.
If the response variables Y i , i = 1, . . ., n are fully observed, the quantile regression estimator of θ is obtained by minimizing the following equation: where ρ τ (u) = u(τ − I(u < 0)) is the check function, and I(•) is the indicator function.For a given τ, θ satisfies the following estimation equation: where In the scenario where the missingness of response variable Y i is nonignorable, let δ i denote the missing indicator.If Y i is missing, δ i = 1; otherwise, δ i = 0. (Z i , Y i , δ i ), i = 1, . . ., n represents an independent and identically distributed sample from (Z, Y, δ).We establish a semiparametric exponential tilting model for missing propensity as follows: where g(•) is an unspecified function, X ⊂ Z is a d-dimensional vector, and there exists an instrumental variable V = Z/X that is unrelated to δ given (X, Y).
Let f κ (Y, Z|X), where κ ∈ 0, 1 denotes the conditional density of (Z, Y) of X when δ = κ.Specifically, we have the following: where O(X, Y) = P(δ=0|X,Y) P(δ=1|X,Y) = exp(−g(X) + γY).For the quantile estimation equation ψ(Y, Z; θ), let it can be easily shown that where The nonparametric kernel estimate of Equation ( 6) is given by m0 (X; where Due to the instability of the nonparametric multivariate kernel estimation of the aforementioned conditional expectation, this paper adopts Monte Carlo methods to estimate m 0 ψ (X; θ, γ).For simplicity of discussion, we consider the parameter assumption of the conditional distribution f (Y|Z, δ = 1; β) of the observed response.This assumption can be verified easily using fully observed samples.Consequently, the conditional distribution of the response with nonrandom missingness satisfies Let Y (j) i , j = 1, . . ., M be independent and identically distributed samples from f (Y|Z i , δ = 0; β, γ).According to the law of large numbers, as M → ∞, we have To obtain a set of random realizations from f 0 (Y|Z; β, γ), the SIR algorithm [19] can be employed based on the parametric representation in (7) for a given (β, γ): (1) Random samples (2) Calculate the adjustment weights for each sample point in S as (3) Resample from S i according to the probabilities ω i1 (γ), . . ., ω iM 2 (γ) to obtain Y i , . . ., Y (M) i .To ensure the convergence of the aforementioned process, it is crucial to have M 2 → ∞ and M/M 2 → 0.
Due to the nonsmoothness of the aforementioned estimation equation, obtaining the sandwich estimator of the asymptotic covariance matrix directly is not feasible.Therefore, this paper proposes using a smooth function G(Z i θ − Y i ) as a substitute for the indicator function I(Y i − Z i θ < 0) in the quantile estimation equation, thus resulting in a smooth approximation of ψ τ (Y i − Z i θ): For nonignorable missing data, we have the following representation for the smoothed SIR-based quantile regression estimation equation: i (β, γ); θ).
The estimation equation based on imputation is susceptible to the misspecification of f (Y|Z, δ = 1; β).Due to the relative robustness of the semiparametric response model, we consider the AIPW (augmented inverse probability weighting) estimation equation: where it can be proven that the AIPW estimation equation is consistent in the case of the misspecification of the parameter model f (Y|Z, δ = 1; β).In practice, (β, γ) are often unknown and need to be estimated separately.The maximum likelihood estimation of β, denoted as β, is the solution to the following score function: Then, we consider the estimation of the tilting parameter γ.The semiparametric missing propensity model is analyzed by considering two estimation approaches for the tilting parameter γ: the profile two-step generalized method of the moments estimation and the kernel regression estimation for the nonparametric component g(•).To estimate the skewness parameter γ, we define the profile estimation equation as follows: where h(V ) is an arbitrary specified function of the instrumental variable V , and g γ (•) satisfies the following: Under the assumption of a correctly specified missingness mechanism, it holds that E{ξ i (g γ , γ)} = 0, and the vector ξ i (g γ , γ) is overidentified with respect to γ.The profile two-stage generalized method of the moments estimation for γ is given by the following: . The estimator ĝγ (•) represents the kernel regression estimate of g(•) and satisfies the following equation: , where K h (u 1 , . . ., u d ) represents the d-dimensional kernel function with a bandwidth h.Define ψ(l) hi (θ, β, γ) and l = 1, 2, which satisfy the following: where m0 Let p i represent the probability mass of ψ(l) hi (θ, β, γ), where i = 1, 2, • • • , n.The empirical log-likelihood ratio function with respect to θ is defined as follows: Using the method of Lagrange multipliers, it can be shown that R(l) (θ) can be expressed as follows: R(l where λ satisfies the following: The empirical likelihood estimators of the quantile regression coefficients based on the two proposed estimation equations in this paper, denoted as θ(l) , l = 1, 2, are given by the following: θ(l) = arg min θ R(l) (θ).

Theoretical Analysis
To elucidate the theoretical properties of the proposed estimators in this paper, we first define the matrix as follows: The density of Z is bounded and has continuous and bounded second-order derivatives; (b) the density of Z and the propensity π(Z, Y, γ 0 ) are bounded away from 0 ; and (c) and K denote a generic notation for a d-dimensional kernel function, and the value of d is determined by the context of use.K is a bounded, uniformly continuous, symmetric function of the m th order satisfying the following conditions: K(s)ds = 1, s = (s 1 , . . . ,s d ), s t l K(s)ds = 0, and s m l K(s)ds = 0 for any l = 1 . . ., /∂γ exists at γ 0 , E sup γ ξ(g γ , γ) < ∞, γ 0 is the unique solution to E{ξ(g γ , γ)} = 0, and W(γ 0 ) is positive definite. (C5) {(Z i , Y i , δ i ) : i = 1, . . ., n} are independent and identically distributed random vectors.
The support of θ denoted by B is a compact set in R q , and θ 0 ∈ B is the unique solution to 3 are bounded by an integrable function H(x, y) within a neighborhood of θ 0 .(C6) For all in a neighborhood of zero and for almost every Z, F( | Z), f ( | Z) and f ( | Z, δ = 0) to exist, they are bounded away from zero and are r times continuously differentiable with r ≥ 2. There exists a function is either strictly positive or strictly negative within (a l−1 , a l ) for l = 1, . . ., L + 1. (C8) The positive bandwidth parameter h satisfies nh 2r → 0 and nh/ log(n) → ∞ as n → ∞. (C9) Z has a bounded support, E Z 4 < ∞ and the matrices Γ and T l ; additionally, l = 1, 2 are nonsingular.
(C10) Under complete observation of (Z i , Y i ) for i = 1, . . ., n, the unique solution β to the score equation in (3) satisfies for some Σ for sufficiently large n.
To ensure the requirements of Lemma 8.11 by Newey and McFadden [22] and Theorem 6.18 by Van [23], the conditions (C1)-(C4), which are commonly found in the literature on missing data and nonparametric method [24,25], are primarily employed.These conditions encompass the following: (1) random equivalence and continuity conditions; (2) linearity conditions on the objective function with respect to nonparametric components and convergence rate conditions for nonparametric estimators; and (3) the differential continuity condition of the estimating equations with respect to the parameter of interest.Conditions (C5)-(C9) ensure the consistency and asymptotic normality of the empirical likelihood estimator for quantile regression smoothing [10].To simplify the discussion on the asymptotic properties of maximum likelihood estimation in the working model, we introduce condition (C10).
Under the fulfillment of the assumed conditions, we define the following: .
In addition, we have the following lemma, whose proof is given in the Appendix A: Lemma 2. Under the assumption conditions (C1)-(C10), with the notation from Section 3, the following results hold as n → ∞: (1) (3) Theorem 1.Under conditions (C1)-(C10), if the parameter working model is correctly specified, as n → ∞ for l = 1, 2, we have If there is no missing data, we have P(δ|Z, Y) = 0, which implies that H 1 , H 2 , and H 3 are all zero.Additionally, we have The above results are consistent with the asymptotic normality conclusion of classical quantile regression.
The different forms of Λ 1 and Λ 2 indicate that if the parameter working model is no longer consistent, while θ(2) remains a consistent estimator of θ 0 .The following procedure demonstrates the double robustness of the AIPW estimation equation.For misspecified f (Y|Z, δ = 1) values, there exists C(Z; θ, γ) such that m0 ψ (Z i ; θ, β, γ) = C(Z i ; θ, γ) + o p (1).This illustrates the double robustness property of the AIPW estimation equation.
It can be shown that If π(X i , Y i ) is correctly specified, the IPW estimation equation is consistent, which implies that the AIPW quantile regression estimation equation remains consistent in this case.
Theorem 2. Under the conditions of Theorem 1, if the parameter working model is correctly specified, for l = 1, 2 and as n → ∞, where r First, if there is no data missingness, we have q , and the Wilks' Theorem holds.Furthermore, it is worth noting that if β and γ are known, we still have T l = A l , and in this case, the Wilks' Theorem still holds.The above conclusion is consistent with Zhao [26].

Simulation Study
To investigate the finite-sample properties of the proposed method, this study conducted numerical simulations under both correctly specified and misspecified working model scenarios.

Simulation 1: Correctly Specified Working Model
In the numerical simulation, we generated a random vector (x, y, δ), where x is the independent variable, y is the response variable of interest, and δ is the indicator variable for the observation of y.When δ = 1, y has observed values; otherwise, the observation of y is missing.Let x ∼ N(0, 0.5), and generate the observed response variable y according to the following equation: where µ(x) = −1 + x.We considered two different distributions for the random error term e: (a) N(0, 0.9) and (b) N(0, 0.49(1 + x 2 )).For the working model f (y|x, δ = 1; β), the former follows a homoscedastic structure N(µ(x), σ 2 1 ), while the latter follows a heteroscedastic structure N(µ(x), σ 2 2 (x)).
The indicator variable δ follows a Bernoulli distribution with parameter p, i.e., δ ∼ Bernoulli(1, p).The conditional probability of δ given (x, y) is defined as follows: where (φ 0 , φ 1 ) = (0.8, −0.2).In this case, the missingness mechanism for the response variable y is nonrandom, with x serving as a missingness instrument.The average observed rate in the sample was approximately 73%.
We establish a quantile regression model of the response variable y on x as follows: where τ represents the quantile of interest, specifically τ = (0.25, 0.5, 0.75).
Consider the following five quantile regression estimation equations: (1) Full Estimation Equation: To generate a sample of size n = 500 that meets the requirements of the simulation, we can use the law of total probability and express f (y|x) as follows: .
Under the specified nonrandom missingness mechanism, we have O(x, y) = exp(−φ 0 − φ 1 y) for the homoscedastic case of f (y|x, δ = 1) = N(µ, σ 2 1 ).In this case, we can express the ratio of the conditional probabilities as follows: where Since x is completely observed, we can draw a sample of size n = 500 from the mixed distribution of f (y|x).A similar approach can be applied under the heteroscedastic assumption.
It should be noted that the response variable y originates from a distribution with a complex, mixed form.As a result, discussing the true values of the parameters (θ 0 , θ 1 ) poses a formidable challenge.This complexity renders it difficult to assess the performance of the estimation methods using conventional measures such as bias or the root mean squared error (RMSE).Consequently, we introduce the following approxmate relative evaluation metrics: where ARMSE(Method * ) = SD(Method * ) 2 + (Mean(Method * ) − Mean(Full)) 2 .Tables 1 and 2 summarize the mean and variance of the five coefficient estimates at different quantiles based on 1000 Monte Carlo simulations under the homoscedastic case (a) and heteroscedastic case (b) of f (y|x, δ = 1).From the estimation results, it can be observed that the coefficient estimates based on complete observations have larger bias compared to the other estimation methods.When the working model f (y|x, δ = 1; β) was correctly specified, the proposed imputation estimates yielded smaller variances compared to the IPW estimates.In this case, the performance of the AIPW estimates was similar to the IPW estimates.Comparing the results at different quantiles, it can be seen that the variances of the five estimation methods at the 0.5 quantile were smaller than those at the 0.25 and 0.75 quantiles, which is due to the larger sample size at the central quantile compared to the tails.Under the homoscedastic assumption, the variances of the estimates at the 0.25 and 0.75 quantiles were similar.Under the existing missing mechanism, as the value of the response variable y increased, the missing propensity also increased, thereby indicating higher missing rates at the upper quantiles.Consequently, the estimation variances of the IPW and AIPW estimates were higher at the high quantile of τ = 0.75 compared to the low quantile of τ = 0.25.However, proper imputation could greatly improve the estimation efficiency at the high quantile of τ = 0.75.This improvement was more pronounced under the heteroscedastic model.These results demonstrate that the imputation estimates are nearly unbiased when the working model is correctly specified and have higher estimation efficiency compared to the IPW and AIPW estimates.

Simulation 2: Misspecification of the Working Model
In practical situations, the true data generation mechanism is unknown, and it is challenging to accurately specify the working model f (y|x, δ = 1; β) for the observed data.In this study, we investigated the finite sample properties of the proposed imputation estimator and calibration estimator under the misspecification of the working model.The simulation model includes two covariates: X 1 ∼ N(0, 1) and X 2 ∼ Exp(0.2).Given the covariates, the response variable Y is generated as follows: where ε ∼ N(0, 1), and X 1 , X 2 , and ε are mutually independent.
In this simulation setup, the error term distribution of Y is heteroscedastic.The missing data mechanism for Y is nonrandom and follows The average observed rate in the model was approximately 73%, and X 2 served as an instrumental variable.We generated a random sample of size n = 500 denoted as (X i , Y i , δ i ) : i = 1, . . ., n.For the aforementioned simulation model, we consider the following quantile regression model: where θ 0 (τ) = (1 + 0.5Q τ (ε), 1, 1 + 0.25Q τ (ε)) .Under the aforementioned data generating mechanism, obtaining an explicit expression for f (Y|X, δ = 1; β) is challenging and requires specifying the working model based on the observed data.In this simulation model, we consider two possible working models: (1) N( μ(X), σ2 ) and ( 2) N( μ(X), (0.5 + 0.25X 2 ) 2 ).The estimation results of the five types of quantile regression estimates obtained from 1000 random simulations at different quantiles are summarized in Tables 3-5.These tables include two types of imputation estimates based on the parameter working models ( 1) and ( 2), as well as the corresponding AIPW estimates based on the parameter working models (1) and ( 2), and the combined estimation equations.The results show that the imputation estimates based on the erroneously specified working model (1) exhibited significant estimation bias.On the other hand, although the working model (2) was also misspecified, it took into account the heteroscedasticity in the conditional distribution of the response variable, thus resulting in smaller estimation bias compared to model (1) and better estimation performance.These findings highlight the sensitivity of imputation methods to misspecified working models.Across the three quantiles, the IPW estimates performed well, thus indicating the robustness of the semiparametric response assumption.Even in the presence of misspecified parameter working models, both of the AIPW estimates had similar median absolute deviations to IPW, which were significantly smaller than the misspecified imputation estimates, thereby demonstrating the robustness of the AIPW estimation.Comparing the two AIPW estimates, it is observed that the estimate based on the correctly specified parameter working model had smaller estimation bias and higher estimation efficiency.

Real Data Application
We applied our proposed method to the data of 2139 HIV-infected patients enrolled in the ACTG175 study [27].The ACTG175 study evaluated the efficacy of monotherapy or combination therapy in HIV-infected patients with CD4 cell counts between 200 and 500 cells/mm³.Following the studies by Davidian et al. [28] and Zhang et al. [29], we categorized all the treatment regimens into two groups.The first group consisted of the standard zidovudine (ZDV) monotherapy arm, while the second group included three newer treatment arms: ZDV and dual nucleoside analogue (ddl), ZDV and zalcitabine (ddC), and ddl monotherapy.The first group comprised 532 subjects, while the second group comprised 1697 subjects.We investigated the effect of the treatment arm (trt, 0 = ZDV monotherapy only) on the τ quantile of the CD4 cell count (CD4 96 ) measured at baseline and adjusted for the baseline CD4 cell count (CD4 0 ) and other baseline covariates, including age, weight, race (0 = Caucasian), gender (0 = female), history of reverse transcriptase inhibitor use (0 = no), and whether the subject discontinued treatment before 96 weeks (offtrt, 0 = no).
Consider fitting a linear quantile regression model as follows: The dataset used in this study is sourced from the R package "speff2trial".The study population consists of 1522 Caucasian individuals and 617 non-Caucasian individuals, with 1171 males and 368 females.The average age of the participants is 35 years, with a standard deviation of 8.7 years.Among the participants, 1253 individuals had a history of antiretroviral therapy, and 776 individuals discontinued treatment before the 96th week.
Due to attrition during the study period, approximately 37% of the participants have missing values for the variable CD496.Although complete measurements of other variables related to CD496, such as baseline CD4 and CD8 cell counts CD40 and CD80, as well as CD4 and CD8 cell counts at 20 ± 5 weeks CD420 and CD820, were obtained at baseline and follow-up visits, these variables may not fully explain the propensity for participants to drop out.In other words, we cannot assume that the missingness of CD4 96 is random.Therefore, in our analysis, we consider a more comprehensive semiparametric nonrandom missingness mechanism: where S represents the set of variables associated with attrition, and g(S) is a function capturing the relationship between these variables and the missingness indicator R. Figure 3 displays the histograms of the observed CD4 96 and its logarithm.From the figure, it can be observed that the conditional distribution f (y|X, R = 1) of observed CD4 96 is right-skewed.However, the logarithmic transformation did not result in improved symmetry, thus indicating that the normality assumption did not hold.In our analysis, we can assume that CD4 96 follows a truncated normal distribution with left truncation at 0, where its mean is primarily determined by the influence of eight covariates and three auxiliary variables.
The parameters β in the working model f (y|X, R = 1; β) are estimated using the truncation regression model in R package "truncreg".The parameter γ in π(S, Y) is estimated using the method of the profile generalized method of moments (GMMs).
Figures 4 and 5 illustrate the normality properties of the residuals from the truncated regression working model.Visually, the distribution of residuals appears to be symmetric.The calculated sample skewness is 0.05, thus indicating a slight deviation from perfect symmetry.The Q-Q plot reveals that the distribution of residuals has a kurtosis less than 3. Further computation reveals a kurtosis of 2.11, thus indicating that the residual distribution is flatter than a standard normal distribution.Table 6 presents the estimates of the quantile regression coefficients and corresponding 95% confidence intervals at the τ = 0.25, 0.5, 0.75 quantile levels.The four estimation methods considered include complete case (CC) estimation, inverse probability weighting (IPW) estimation, multiple imputation (MI) estimation, and augmented inverse probability weighting (AIPW) estimation.The MI estimation is based on averaging over L = 20 randomly generated imputations.Confidence intervals for the coefficient estimates were obtained using the bootstrap method with B = 200 resampling iterations.
From Table 6, it can be observed that for the three given quantile levels and four estimation methods, patients receiving the three new combined treatment methods had significantly higher CD4 cell counts at 96 ± 5 weeks compared to the traditional treatment method.In other words, the new treatment methods had significantly slowed down the progression of AIDS compared to the traditional method.Comparing the four estimation methods, it is evident that the complete case estimation overestimated the performance of the treatment group.The results of the IPW estimation and AIPW estimation were similar and higher than the imputation estimation.When comparing the treatment effects at different quantile levels, both the IPW estimation and imputation estimation reflected a decreasing trend in treatment effect from the 0.25th to the 0.75th quantile.Although the AIPW estimation and complete case estimation did not show a similar trend, the coefficient estimates of the AIPW estimation also indicate a more significant improvement in treatment effect for patients at lower quantiles.
Upon examining the effects of the other covariates, it is found that for all four estimation methods, the baseline CD4 level CD4 0 had a positive impact on the CD4 cell count at 96 ± 6 weeks, while patients with a history of antiretroviral therapy or early treatment discontinuation exhibited poorer CD4 cell levels at 96 ± 5 weeks.In comparison to the covariates directly related to the disease progression mentioned above, the effects of age, weight, race, gender, and other covariates on the CD4 cell count at 96 ± 6 weeks were minimal.The impact directions and significance obtained from different methods were also not consistent.Therefore, although these variables needed to be considered in the modeling process, conclusions regarding their effects should be drawn with caution.

Discussion
In this study, we address the bias in quantile regression estimates by constructing imputation and AIPW estimation equations, with both involving the estimation of conditional means under nonrandom missingness.Many existing methods rely on kernel regression to estimate conditional means.However, nonparametric estimation methods may suffer from the curse of dimensionality when the dimension of the covariates is high.Paik and Larsen [19] proposed using importance resampling to obtain Monte Carlo estimates of conditional means, and Song et al. [21] further applied this method to estimation equations.In this study, we extend these methods to quantile regression and overcome the theoretical and computational challenges caused by the nonsmoothness of the checking function in classical quantile regression by employing convolution smoothing.
Common parameter working models are based on linear regression for observed data.Song et al.'s [21] simulation results showed that model misspecification does not lead to estimation bias.However, their simulation study was based on a regression model that satisfied the Gauss-Markov assumption, with missing response variables following a normal distribution with homoscedasticity concerning the covariates.Misspecification was reflected in the estimation of the mean or location variables.However, the advantages of quantile regression are more evident in situations involving skewness, heavy tails, and heteroscedasticity.In this study, our simulation results under heteroscedasticity showed that imputation estimation based on the assumption of a linear regression working model leads to significant estimation bias, while the AIPW estimation equation can mitigate the impact of model misspecification.We also provide theoretical proof of the consistency of AIPW estimation.
Our simulation results demonstrate that, under the ideal scenario of correctly specified parameter working models, the imputation estimator is more efficient than the IPW and AIPW estimators.The AIPW estimator based on the correctly specified model was found to be more efficient than that based on the misspecified model.Therefore, in practical applications, it is crucial to appropriately specify the parameter working model based on the observed data.Fortunately, the effectiveness of the model specification can be assessed using various methods such as Q-Q plots and histograms.For the observed response conditional distributions that do not conform to the linear regression assumption, a Box-Cox transformation can be applied to approximate a normal parameter working model.If such a parameter working model is difficult to obtain, the AIPW estimator proposed in this study can still provide relatively reliable estimates.This is because the proposed response mechanism model is semiparametric and offers certain flexibility.However, the response model constructed in this study does not consider the interaction effects between covariates X and the response variable Y or the potential nonlinear effects of the response variable Y on the missingness propensity.

Similarly, we have
Based on the assumptions and the Taylor expansion, we have Notice that Under the assumption conditions, we have where lim .
By the law of large numbers, as n → ∞, we have 1 n ∑ n i=1 ∂ ψ(l) hi (θ 0 , β, γ) ∂θ p → Γ.Finally, we demonstrate (4).From n −1 (max i ψ(l) hi (θ 0 , β, γ) ) 2 n −1 ∑ n i=1 ψ(l) hi (θ 0 , β, γ) ⊗2 → 0, it can be easily shown that, for l = 1, 2, max i ψ(l) hi (θ 0 , β, γ) = o p (n 1/2 ).Proof of Theorem 1.By applying the Lagrange multiplier method, we obtain the empirical log-likelihood ratio function with respect to the parameter vector θ: 2n ( θ, 0) = 0.Under the assumption conditions, according to Lemma A.1 in Newey and Smith [30] and Theorem 1(a) in Leng and Tang Leng [31], it can be shown that θ is a consistent estimator of θ0.By Taylor expanding T r, for almost all Z and in a neighborhood of zero, and E C(X) X 2 < ∞. (C7) The kernel function K(•) is a probability density function such that (a) it is bounded and has a compact support; (b) K(•) is an rth order kernel, i.e., K(•) satisfies u Figures 1 and 2 illustrate that the residual distribution of the working model (1) exhibited peakedness, thus violating the normality assumption and indicating model misspecification.In contrast, working model (2) took into account the correct specification of the variance.

Figure 1 .
Figure 1.Histogram of the residual distribution for the parameter working model N( μ(X), σ2 ).

Figure 2 .
Figure 2. QQ plot of the residual distribution for the parameter working model N( μ(X), σ2 ).

Figure 4 .
Figure 4. Histogram of residuals from the parameterized working model.

Figure 5 .
Figure 5. Q-Q plot of residuals from the parameterized working model.

Table 1 .
Monte Carlo mean, standard deviation (SD), and approximate relative performance (ARE) of the five methods for error term (a).

Table 2 .
Monte Carlo mean, standard deviation (SD), and approximate relative performance (ARE) of the five methods for error term (b).

Table 3 .
The bias (Bias), standard deviation (SD), and median absolute deviation (MAD) of the five types of quantile regression coefficient estimates at τ = 0.25.

Table 4 .
The bias (Bias), standard deviation (SD), and median absolute deviation (MAD) of the five types of quantile regression coefficient estimates at τ = 0.5.

Table 5 .
The bias (Bias), standard deviation (SD), and median absolute deviation (MAD) of the five types of quantile regression coefficient estimates at τ = 0.75.