Semiparametric model averaging for high-dimensional quantile regression with nonignorable nonresponse

TL;DR

This paper introduces a novel semiparametric model averaging approach for high-dimensional quantile regression with nonignorable missing data, improving model selection and estimation accuracy in incomplete datasets.

Contribution

It develops a semiparametric model averaging method that handles nonignorable missingness and asymptotically concentrates on the correct model, advancing quantile regression techniques.

Findings

The proposed SMA estimator achieves asymptotic optimality.

The method effectively handles complex missingness mechanisms.

Empirical results demonstrate improved efficiency over existing methods.

Abstract

Model averaging has demonstrated superior performance for ensemble forecasting in high-dimensional framework, its extension to incomplete datasets remains a critical but underexplored challenge. Moreover, identifying the parsimonious model through averaging procedure in quantile regression demands urgent methodological innovation. In this paper, we propose a novel model averaging method for high-dimensional quantile regression with nonignorable missingness. The idea is to relax the parametric constraint on the conditional distribution of respondents, which is constructed through the two-phase scheme: (i) a semiparametric likelihood-based estimation for the missing mechanism, and (ii) a semiparametric weighting procedure to combine candidate models. One of pivotal advantages is our SMA estimator can asymptotically concentrate on the optimally correct model when the candidate set involves at least one correct model. Theoretical results show that the estimator achieves asymptotic optimality even under complex missingness conditions. Empirical conclusions illustrate the efficiency of the method.

Tables

Table 1. Table 1: Simulation results: MS, FP, FN, FPR and its standard deviation (in parentheses) for different methods.

$(\tau ,\rho )$	Missingness	Method	MS1	$\mathrm{FP}$2	$\mathrm{FN}$3	FPR(SD) ($\times {10}^2$)
(0.05, 0)	37.87$\%$	SMA	14	9.860	0.860	78.84(18.47)
		PS-SMA	8.215	4.475	1.260	78.97(19.36)
		GPS-SMA	10.210	6.110	0.900	76.86(19.43)
(0.05, 0.5)	39.82$\%$	SMA	14	9.010	0.010	71.70(22.42)
		PS-SMA	7.045	2.130	0.100	61.86(18.31)
		GPS-SMA	5.200	0.445	0.245	59.71(18.16)
(0.5, 0)	38.54$\%$	SMA	14	10.43	1.430	64.00(10.40)
		PS-SMA	7.860	4.665	1.805	62.72(10.88)
		GPS-SMA	10.385	6.870	1.485	62.45(10.94)
(0.5, 0.5)	37.67$\%$	SMA	14	9.270	0.270	58.33(11.33)
		PS-SMA	6.250	2.250	1.015	55.72(10.79)
		GPS-SMA	6.345	2.265	0.920	54.17(10.89)

Equations

$$Q_{\tau }(y_i|{\mathbf{x}}_i)={\mathbf{\beta }}^{\top }(\tau ){\mathbf{x}}_i\text{triangleeq}{\mu }_i(\tau ),$$

$$y_i={\mu }_i(\tau )+{\epsilon }_i(\tau ),$$

$$\mathrm{Pr}(r_i=1|{\mathbf{x}}_i,y_i)=\frac{1}{1+\exp \left({\phi }_0+{\sum }_{j=1}^p{\phi }_j{x}_{ij}+\gamma y_i\right)}\text{triangleeq}{\pi }_i(\mathbf{\theta }),$$

$$\frac{1-\mathrm{Pr}(r_{i}=1|\bm{x}_{i};\bm{\theta})}{\mathrm{Pr}(r_{i}=1|\bm{x}_{i};\bm{\theta})}=E\left(\frac{1-\pi_{i}(\bm{\theta})}{\pi_{i}(\bm{\theta})}\bigg{|}\bm{x}_{i},r_{i}=1\right)=\exp\left\{\varphi_{0}+\sum^{p}_{j=1}\varphi_{j}x_{ij}+m_{1}(\bm{x}_{i};\gamma)\right\},$$

$$m_1(\mathbf{x};\gamma )=\log E\{r\exp (\gamma y)|\mathbf{x}\}-\log \mathrm{Pr}(r=1).$$

$${\hat{m}}_1({\mathbf{x}}_i;\gamma )=\log \frac{{\sum }_{j=1}^n{r}_j\exp (\gamma y_j)K_h({\mathbf{x}}_i-{\mathbf{x}}_j)}{{\sum }_{j=1}^n{K}_h({\mathbf{x}}_i-{\mathbf{x}}_j)}-\log \frac{{\sum }_{j=1}^n{r}_j}{n}.$$

$$-\mathrm{logit}\left\{{\mathrm{Pr}}^{(k)}(r=1|\mathbf{x};{\mathbf{\theta }}^{(k)})\right\}={\mathbf{\phi }}^{{(k)}^{\top }}{\mathbf{x}}^{(k)}+m_1(\mathbf{x};\gamma )=\sum _{j=1}^{|{\mathcal{M}}_k|}{\phi }_j{x}_j+m_1(\mathbf{x};\gamma ),$$

$${\hat{\mathbf{\theta }}}_n^{(k)}={\left({\hat{\mathbf{\phi }}}_n^{{(k)}^{\top }},{\hat{\gamma }}_n^{(k)}\right)}^{\top }={\mathrm{argmin}}_{{\mathbf{\theta }}^{(k)}\in {\Phi }_k\times \Gamma }{L}_n^{(k)}\left({\mathbf{\theta }}^{(k)}\right),$$

$$L_n^{(k)}({\mathbf{\theta }}^{(k)})=n^{-1}\sum _{i=1}^n\left[(r_i-1)\left({\mathbf{\phi }}^{{(k)}^{\top }}{\mathbf{x}}_i^{(k)}+{\hat{m}}_1({\mathbf{x}}_i;\gamma )\right)+\log \left\{1+\exp \left({\mathbf{\phi }}^{{(k)}^{\top }}{\mathbf{x}}_i^{(k)}+{\hat{m}}_1({\mathbf{x}}_i;\gamma )\right)\right\}\right]$$

$${\hat{\eta }}_i(\mathbf{\omega })=\sum _{k=1}^{S_1}{\omega }_k{\hat{\mathbf{\phi }}}_n^{{(k)}^{\top }}{\mathbf{x}}_i^{(k)}+{\hat{m}}_1\left({\mathbf{x}}_i;\hat{\gamma }(\mathbf{\omega })\right),$$

$$\mathrm{SL}(\hat{\eta }(\mathbf{\omega }))=\sum _{i=1}^n\left[(1-r_i){\hat{\eta }}_i(\mathbf{\omega })-\log \left\{1+\exp \left({\hat{\eta }}_i(\mathbf{\omega })\right)\right\}\right]-{\lambda }_n\sum _{k=1}^{S_1}{\omega }_k|{\mathcal{M}}_k|,$$

$$\hat{\mathbf{\theta }}(\hat{\mathbf{\omega }})={\left({\hat{\mathbf{\phi }}}^{\top }(\hat{\mathbf{\omega }}),\hat{\gamma }(\hat{\mathbf{\omega }})\right)}^{\top }=\sum _{k=1}^{S_1}{\hat{\omega }}_k{\left({\hat{\mathbf{\phi }}}_n^{{(k)}^{\top }},0_{S_1-k}^{\top },{\hat{\gamma }}_n^{(k)}\right)}^{\top },$$

$$\mathrm{KL}^{(k)}(\bm{\theta})=E\left\{r\log\frac{\mathrm{Pr}(r=1|\bm{x};\bm{\theta}_{0})}{\mathrm{Pr}^{(k)}(r=1|\bm{x};\bm{\theta}^{(k)})}+(1-r)\log\frac{\mathrm{Pr}(r=0|\bm{x};\bm{\theta}_{0})}{\mathrm{Pr}^{(k)}(r=0|\bm{x};\bm{\theta}^{(k)})}\bigg{|}\bm{x}\right\}$$

$${\mathbf{\theta }}_{\ast }^{(k)}={\mathrm{argmin}}_{{\mathbf{\theta }}^{(k)}\in {\Phi }_k\times \Gamma }E\left({\mathrm{KL}}^{(k)}(\mathbf{\theta })\right).$$

$${\hat{\mathbf{\theta }}}_n^{(k)}\stackrel{p}{\to }{\mathbf{\theta }}_{\ast }^{(k)}\text{as}n\to \infty .$$

$$\mathrm{KL}({\eta }_{\ast }(\mathbf{\omega }))=\sum _{i=1}^n\left\{E(r|{\mathbf{x}}_i)\log \frac{\mathrm{Pr}(r=1|{\mathbf{x}}_i;{\mathbf{\theta }}_0)}{\mathrm{Pr}(r=1|{\mathbf{x}}_i;{\mathbf{\theta }}_{\ast }(\mathbf{\omega }))}+E(1-r|{\mathbf{x}}_i)\log \frac{\mathrm{Pr}(r=0|{\mathbf{x}}_i;{\mathbf{\theta }}_0)}{\mathrm{Pr}(r=0|{\mathbf{x}}_i;{\mathbf{\theta }}_{\ast }(\mathbf{\omega }))}\right\}.$$

$$\mathrm{KL}({\eta }_{\ast }({\mathbf{\omega }}_u))=\sum _{i=1}^n\left\{E(r|{\mathbf{x}}_i)\log \frac{\mathrm{Pr}(r=1|{\mathbf{x}}_i;{\mathbf{\theta }}_0)}{\mathrm{Pr}(r=1|{\mathbf{x}}_i;{\mathbf{\theta }}_{\ast }({\mathbf{\omega }}_u))}+E(1-r|{\mathbf{x}}_i)\log \frac{\mathrm{Pr}(r=0|{\mathbf{x}}_i;{\mathbf{\theta }}_0)}{\mathrm{Pr}(r=0|{\mathbf{x}}_i;{\mathbf{\theta }}_{\ast }({\mathbf{\omega }}_u))}\right\}.$$

$${\mu }_i^{(m)}={\mathbf{\beta }}^{{(m)}^{\top }}{\mathbf{x}}_i^{(m)}=\sum _{j\in {\mathcal{I}}_m}{\beta }_j{x}_{ij},$$

$$E{\rho }_{\tau }\left(\epsilon \right)=E\left\{r{\rho }_{\tau }\left(\epsilon \right)+(1-r)M_0(\mathbf{x};\gamma )\right\},$$

$$M_0(\mathbf{x};\gamma )=\frac{E_1\left\{\exp (\gamma y){\rho }_{\tau }\left(\epsilon \right)|\mathbf{x}\right\}}{E_1\left\{\exp (\gamma y)|\mathbf{x}\right\}},$$

$${\hat{\mathbf{\beta }}}_n^{(m)}={\mathrm{argmin}}_{{\mathbf{\beta }}^{(m)}\in {\mathcal{B}}_m}{n}^{-1}\sum _{i=1}^n\left\{r_i{\rho }_{\tau }\left(y_i-{\mathbf{\beta }}^{{(m)}^{\top }}{\mathbf{x}}_i^{(m)}\right)+(1-r_i){\hat{M}}_0\left({\mathbf{x}}_i;\hat{\gamma }(\hat{\mathbf{\omega }}),{\mathbf{\beta }}^{(m)}\right)\right\},$$

$${\hat{M}}_0\left({\mathbf{x}}_i;\hat{\gamma }(\hat{\mathbf{\omega }}),{\mathbf{\beta }}^{(m)}\right)=\frac{{\sum }_{j=1}^n{r}_j{K}_h\left({\mathbf{x}}_j^{(m)}-{\mathbf{x}}_i^{(m)}\right)\exp (\hat{\gamma }(\hat{\mathbf{\omega }})y_j){\rho }_{\tau }\left(y_j-{\mathbf{\beta }}^{{(m)}^{\top }}{\mathbf{x}}_i^{(m)}\right)}{{\sum }_{j=1}^n{r}_j{K}_h\left({\mathbf{x}}_j^{(m)}-{\mathbf{x}}_i^{(m)}\right)\exp (\hat{\gamma }(\hat{\mathbf{\omega }})y_j)}.$$

$${\widehat{\mathrm{QL}}}_n(\mathbf{w})=n^{-1}\sum _{i=1}^n\left\{r_i{\rho }_{\tau }\left(y_i-\sum _{m=1}^{S_2}{w}_m{\hat{\mu }}_i^{(m)}\right)+(1-r_i){\hat{M}}_0\left({\mathbf{x}}_i;\hat{\gamma }(\hat{\mathbf{\omega }}),\mathbf{w}\right)\right\}+\frac{{\varphi }_n}{n}\sum _{m=1}^{S_2}{w}_m\left|{\mathcal{I}}_m\right|,$$

$${\hat{M}}_0\left({\mathbf{x}}_i;\hat{\gamma }(\mathbf{\omega }),\mathbf{w}\right)$$

$$=\frac{{\sum }_{j=1}^n{r}_l{K}_h\left({\mathbf{x}}_j^{(m)}-{\mathbf{x}}_i^{(m)}\right)\exp (\hat{\gamma }(\hat{\mathbf{\omega }})y_j){\rho }_{\tau }\left(y_j-{\sum }_{m=1}^{S_2}{w}_m{\hat{\mu }}_i^{(m)}\right)}{{\sum }_{j=1}^n{r}_j{K}_h\left({\mathbf{x}}_j^{(m)}-{\mathbf{x}}_i^{(m)}\right)\exp (\hat{\gamma }(\hat{\mathbf{\omega }})y_j)}.$$

$$\hat{\mathbf{w}}={\mathrm{argmin}}_{\mathbf{w}\in {\mathcal{H}}_2}{\widehat{\mathrm{QL}}}_n\left(\mathbf{w}\right),$$

$${\hat{Q}}_{\tau }(y|\mathbf{x})=\sum _{m=1}^{S_2}{\hat{w}}_m{\hat{\mathbf{\beta }}}_n^{{(m)}^{\top }}{\mathbf{x}}^{(m)}.$$

$$\mathrm{Pr}\left(\underset{1\le m\le S_2}{\max }\Vert {\hat{\mathbf{\beta }}}_n^{(m)}-{\mathbf{\beta }}_{\ast }^{(m)}\Vert \ge C{\left|\bar{\mathcal{I}}\right|}^{{\kappa }_1}{n}^{-{\kappa }_2}\right)\le 2S_2\exp \left(-C_4{\underline{c}}^2{n}^{1-4{\kappa }_2}{\left|\bar{\mathcal{I}}\right|}^{4{\kappa }_1-2}\right)+S_2\exp \left(-C_{6}n\right).$$

$$\mathrm{FPR}_{n,\tau}\left(\bm{w}\right)=E\left\{\rho_{\tau}\left(y-\sum^{S_{2}}_{m=1}w_{m}\hat{\bm{\beta}}^{{(m)}^{\top}}_{n}\bm{x}^{(m)}\right)\bigg{|}\mathcal{D}_{n}\right\}.$$

$$\frac{{\mathrm{FPR}}_{n,\tau }\left(\hat{\mathbf{w}}\right)}{{\inf }_{\mathbf{w}\in {\mathcal{H}}_2}{\mathrm{FPR}}_{n,\tau }\left(\mathbf{w}\right)}=1+O_p\left(R_n+n^{-1}{\varphi }_n|\bar{\mathcal{I}}|\right)=1+o_p(1),$$

$${\mathrm{EFPR}}_{n,\tau }(\mathbf{w})={\mathrm{FPR}}_{n,\tau }(\mathbf{w})-E{\rho }_{\tau }(\epsilon ).$$

Full text

Keywords: Asymptotic optimality, Missing not at random (MNAR), Frequentist model averaging, Quantile regression, Weight consistency.

1 Introduction

Statistical inference of high-dimensional data has become increasingly important across many scientific fields. In the past two decades, there has been a number of contributive literature such as Fan and Lv (2008), Wang (2011), Li et al. (2012). Compared to least squares model, quantile regression (Koenker 2005) offers a robust alternative to check the effect of high-dimensional covariates on the entire response distribution. To efficiently select active predictors, He et al. (2013) introduced a criterion of variable screening by marginal quantile utility. More recently, Kong et al. (2019) and Jiang et al. (2024) modified the method and used it to refine variable selection in ultra-high dimensional quantile model, respectively.

One of central challenges in high-dimensional regression analysis is prediction. Model averaging emerges as a compelling alternative to classical selection procedures (such as AIC and BIC), effectively addressing the inherent uncertainty associated with selecting a single model. Following a series of pioneering contributions by Yang (2001), Hjort and Claeskens (2003), and Hansen (2007), theoretical development of frequentist model averaging has been rapidly advancing. For quantile regression, Shan and Yang (2009) applied the adaptive aggregating algorithm for mixing different specifications, Lu and Su (2015) proposed the jackknife model averaging (JMA) procedure for the divergent number of predictors, Wang et al. (2023) combined JMA with variable screening for the ultra-high dimensional case, Xu et al. (2023) studied the application of cross-validation model averaging (CVMA) in functional linear model, and among others.

Abovementioned methodologies are established with fully available dataset. In practice, nonresponse is frequently encountered during data collecting. Naive exclusion of missingness may causes biased results when the response mechanism depends on unobserved parts — the scenario defined as missing not at random (MNAR) (Little and Rubin 2002; Kim and Shao 2021). Statistical inference of MNAR data has been greatly developed, including dimension reduction and estimation (e.g., Kim and Yu 2011, Miao and Tchetgen 2016, Shao and Wang 2016, Morikawa and Kim 2021, Zhang et al. 2020, Zhao et al. 2021). In recent years, researchers begin to take focus on model averaging with missing data, and most of studies are regarded with ignorable mechanisms. For example, Wei et al. (2021), Xie et al. (2021) used frequentist model averaging and inverse probability weighting (IPW) to predict the conditional mean function, respectively. However, for MNAR, the literature currently lacks model averaging procedures for either the response mechanism or the conditional quantile function (CQF). We wish to address the void in this paper.

The objective of the paper is to contribute a novel semiparametric model averaging (SMA) procedure for predicting CQF. To the best of our knowledge, most of methodologies can hardly skirt the estimation of MNAR mechanism. As pointed out by Li et al. (2023), instability will be occurred through the moment estimating equation. They suggested a likelihood estimation under the parametric assumption of the conditional distribution of respondents. In the high-dimensional case, we derive a weighted estimator that combines semiparametric likelihood and model averaging. The advantage is to avoid potential misspecification of the conditional distribution, and the computational burden from nonparametric parts can be practically reduced by a sampling technique. We verify that the estimator of nuisance parameter is consistent when the correct model are contained in the candidate set.

Asymptotic theory of the weight estimator is another interesting topic in model averaging literature, while related studies are almost focused on least squares regression (Zhang and Liu 2019; Zhang et al. 2020; Fang et al. 2023). Recent advance in Xu et al. (2023) addressed the gap by establishing the over-consistency (the concept refers to Chen et al. 2023) of weights of CVMA, while the behavior of the weight assigned to the optimally correct model remains to be developed. The current paper enriches the theoretical results of proposed procedure in two scenarios: (i) when the candidate set contains at least one correct model, SMA can asymptotically put the weight one to the optimally correct model, i.e., parsimony; (ii) when all candidates are misspecified, SMA preserves the same asymptotic optimality derived from full data in the sense that its prediction risk via quantile loss is asymptotically identified to the best–but–infeasible averaging estimator. Note that above conclusions are investigated on the continuous weight support with “nonnegative and sum-up-to-one”, it is regarded as a supplement beyond conventional least squares framework.

Numerically we propose a sampling approximation for SMA procedure to alleviate computational burden induced by nonparametric regression. Besides, to address the curse of dimensionality and efficiently construct candidate sets, we also develop a greedy post-selection algorithm based on grouped variable importance. We empirically compare the proposed estimation to existing methods of model selection and model averaging. Numerical results demonstrate that SMA can effectively identify the optimal candidate and has satisfactory predictive performance. Further, our method is more robust than others to against misspecification of the response mechanism.

The reminder of this paper is organized as follows. Section 2 describes the model framework, the SMA procedure for MNAR mechanism, and investigates the consistency of nuisance parameter. Section 3 proposes the SMA criterion to predict CQF. Section 4 establishes the theoretical properties of the CQF estimator. Section 5 introduces an efficient computational sampling method to evaluate conditional expectations of respondents. Numerical studies are presented in Section 6. A post-selection algorithm for preparing candidate sets is described in Appendix. Supplementary simulation results and technical details of proofs are attached in online supplementary material, respectively.

2 Model and estimation for MNAR mechanism

2.1 Model specification

Let $\mathcal{D}_{n}=\{(\bm{x}^{\top}_{i},y_{i})^{\top}:i=1,...,n\}$ be the set of independent and identically distributed (i.i.d.) copies of $(\bm{x}^{\top},y)^{\top}$, where $\bm{x}$ is a $p$-dimensional covariate vector that can be fully observed, and $y$ is the outcome with some of observations are nonresponse. Given a probability level $\tau\in(0,1)$, the CQF of the response takes the form as

[TABLE]

where $\bm{\beta}(\tau)$ is the parameter vector that may depends on $\tau$. Obviously (2.1) is equivalent to the following quantile regressive model:

[TABLE]

where $\varepsilon_{i}(\tau)$’s are i.i.d. random errors such that $\mathrm{Pr}(\varepsilon_{i}(\tau)\leq 0|\bm{x}_{i})=\tau$. Hereafter we omit the notation “$(\tau)$” for simplicity. Since some of $y_{i}$’s are subject to missingness, we define $r_{i}$ as the binary response indicator such that $y_{i}$ is observed when $r_{i}=1$ and $y_{i}$ is missing when $r_{i}=0$. The response mechanism is assumed as a logistic model that

[TABLE]

where $\bm{\theta}=(\varphi_{0},...,\varphi_{p},\gamma)^{\top}$ is the unknown parameter vector.

Through the paper, $p$ is allowed to be divergent with respect to $n$, and the main interest is to effectively predict (2.1) by the pairs $(\bm{x}^{\top}_{i},y_{i})^{\top}$ among $r_{i}=1$. Due to MNAR, we firstly need to evaluate the nuisance parameter in (2.2).

2.2 Semiparametric model averaging for response mechanism

As we can see, evaluating $\bm{\theta}$ is challenging because (i) the nonresponse variable is involved in the mechanism model and (ii) the distribution of $y|\bm{x},r=1$ may be inconsistent with that of $y|\bm{x}$, one cannot directly construct the likelihood function for (2.2). Motivated by Kim and Shao (2021), we can find

[TABLE]

where $m_{1}(\bm{x}_{i};\gamma)=\log E(\exp(\gamma y)|\bm{x}_{i},r_{i}=1)$ replaces the nonresponse part. To avoid the potential misspecification, we allow a nonparametric form on $y|\bm{x},r=1$, hence the question is transformed into estimating parameter in $\mathrm{Pr}(r_{i}=1|\bm{x}_{i};\bm{\theta})$ by a semiparametric method.

To describe the SMA procedure, we first use nonparametric regression for $m_{1}(\bm{x}_{i};\gamma)$. Note that

[TABLE]

For a kernel function $K_{h}(\cdot)=K(\cdot/h)$ and an approximated bandwidth $h=h_{n}$, the corresponding estimator is

[TABLE]

Having parametrized $\mathrm{Pr}(r_{i}=1|\bm{x}_{i};\bm{\theta})$, we can decompose it into multiple candidates and aggregate them through model averaging.

Let $\bm{\theta}^{(k)}=(\bm{\varphi}^{{(k)}^{\top}},\gamma)^{\top}$, and $\mathcal{S}_{1}=\{\mathrm{Pr}^{(k)}(r=1|\bm{x};\bm{\theta}^{(k)}):k=1,...,S_{1}\}$ be the candidate model set of $\mathrm{Pr}(r=1|\bm{x};\bm{\theta})$. Note that the estimation of $\bm{\theta}$ will participate in predicting CQF, it imposes the requirement on $\mathcal{S}_{1}$ that: there exists at least one candidate model that sufficiently covers active predictors, i.e., it is preferable to contain correct models. Of this view we make a nested assumption on $\mathcal{S}_{1}$ for both theoretical and practical convenience. Suppose that all linear predictors are decreasingly ranked by their importance, the $k$-th candidate is formed as

[TABLE]

where $\mathcal{M}_{1}\subset...\subset\mathcal{M}_{S_{1}}\subset\{1,...,p\}$ are index sets of linear predictors and $|\mathcal{M}_{k}|$ is the cardinality. The parameter can be estimated as

[TABLE]

where

[TABLE]

is the negative log-likelihood function of the $k$-th model and $\hat{m}_{1}(\bm{x}_{i};\gamma)$ is given in (2.3).

The remainder is to match weights over candidates. Let $\bm{\omega}=(\omega_{1},...,\omega_{S_{1}})^{\top}\in\mathcal{H}_{1}$ be a vector where the $k$-th component $\omega_{k}$ represents the weight of $\hat{\bm{\theta}}^{(k)}_{n}$, and $\mathcal{H}_{1}=\{\bm{\omega}\in[0,1]^{S_{1}}:\sum^{S_{1}}_{k=1}\omega_{k}=1\}$. Define $\hat{\gamma}(\bm{\omega})=\sum^{S_{1}}_{k=1}\omega_{k}\hat{\gamma}^{(k)}_{n}$ and

[TABLE]

the weighting criterion of SMA is proposed as

[TABLE]

where $\lambda_{n}$ is a penalty scale depending on $n$. The regularized term is imposed to shrink the weight corresponding to overfitted models. Consequently, the SMA estimator of $\bm{\theta}$ is

[TABLE]

where $\bm{0}_{S_{1}-k}$ is a $(|\mathcal{M}_{S_{1}}|-|\mathcal{M}_{k}|)$-dimensional vector of [math], and $\hat{\bm{\omega}}=\mathop{\arg\max}_{\bm{w}\in\mathcal{H}_{1}}\mathrm{SL}(\hat{\eta}(\bm{\omega}))$.

2.3 Asymptotic properties of the estimator

This section devotes to study the theoretical performance of (2.4). Under the linear specification, it is apparent to conform the identification of MNAR mechanism through distinguishing instrumental variable (Morikawa and Kim 2021, Li et al. 2023). In the following discussion we default the estimator is identified. Another fundamental permission is that both $S_{1}$ and $|\mathcal{M}_{k}|$’s can diverge with $n$ due to the high-dimensional setting (e.g., $|\mathcal{M}_{S_{1}}|$ is allowed to be upper-bounded by $O(p)$). Since (2.2) performs as an agent of quantile regression, we particularly pay attention to the convergency of nuisance parameter. Denote $\bm{\theta}_{0}$ as the “true value” of $\bm{\theta}$, and

[TABLE]

as the Kullback–Leibler (KL) divergence between $\mathrm{Pr}^{(k)}(r=1|\bm{x};\bm{\theta}^{(k)})$ and $\mathrm{Pr}(r=1|\bm{x};\bm{\theta}_{0})$. Commonly conditions (see in, White 1982) imply the existence of a unique “pseudo-true” $\bm{\theta}^{(k)}_{*}$, which satisfies

[TABLE]

To prevent confusion we state some notations: (i) $\|\cdot\|$ denotes the Euclidean norm (i.e., for a matrix $A=(a_{ij})$, $\|A\|=(\sum_{i,j}a_{ij})^{1/2}$). (ii) $\mathcal{X}$ denotes the support of $\bm{x}$. (iii) $a\lesssim(\gtrsim)b$ denotes $a\leq(\geq)C\cdot b$ with a positive constant $C$. (iv) $f_{r}(\cdot)$ and $E_{r}(\cdot)$ denote conditional density function and conditional expectation among $r\in\{0,1\}$, respectively. Firstly, we investigate consistency of candidate estimators, which requires the following regularity conditions:

1. (A1)

$\Phi_{k}$’s and $\Gamma$ are compact over $k=1,...,S_{1}$. 2. (A2)

$f_{1}(\bm{x})$ is bounded away from [math] and $\infty$ uniformly over $\bm{x}\in\mathcal{X}$. 3. (A3)

The kernel $K(u)$ has a compact support and bounded derivatives until order $k$ ($k\geq 2$), $\int K(u)\mathrm{d}u=1$ and has zero moments of order up to $m-1$ and nonzero $m$-th order moment. 4. (A4)

There exists $v\geq 4$ such that $E(y^{v})$ is finite. On the other hand, $E\{\exp(v\gamma y)\}$ and $E_{1}\{\exp(v\gamma y)|\bm{x}\}f_{1}(\bm{x})$ are bounded away from [math] and infinity on $\mathcal{X}$. 5. (A5)

$0<\pi_{i}(\bm{\theta})<1$ uniformly over $\bm{\theta}\in\Phi\times\Gamma$. 6. (A6)

As $n\to\infty$, the bandwidth satisfies: (i) $h\to 0$, (ii) $n^{1-2/v}h^{p}/\log(n)\to\infty$, (iii) $\sqrt{n}h^{p+2k}/\log(n)\to\infty$ and (iv) $\sqrt{n}h^{2m}\to 0$. 7. (A7)

$\max_{1\leq j\leq p}Ex^{4}_{j}=O(1)$.

Conditions (A1)–(A6) are widely employed in the literature of semiparametric inference for missing data (Kim and Yu 2011; Shao and Wang 2016; Morikawa and Kim 2021). In addition to establish the convergence of nonparametric estimator, they also imply that $Em_{1}(\bm{x};\gamma)$ has finite derivatives to order $2$. Condition (A7) is a mild restriction of the moment in high-dimensional data analysis (Wang et al. 2023), and is alternative for deriving the asymptotic properties of SMA estimator.

Lemma 1**.**

Under conditions (A1)–(A7), if $\max_{1\leq k\leq S_{1}}|\mathcal{M}_{k}|\lesssim n^{1/2}$, we have for $k=1,...,S_{1}$,

[TABLE]

From Lemma 1 it is easy to see that the estimator of a correct specification (in $\mathcal{S}_{1}$) is consistent to $(\bm{\theta}^{\top}_{0},0,...,0)^{\top}$. This result will play as a crucial premise to establish the consistency of SMA estimator. To the best of our knowledge, a well-working model averaging procedure should perform the largest propensity to the “optimally correct” (abbreviated OC) candidate. Let $\mathcal{M}_{T}$ be the index set of significant linear predictors in $\mathrm{Pr}(r=1|\bm{x};\bm{\theta})$, the $k$-th nested candidate is OC if $\mathcal{M}_{T}\subset\mathcal{M}_{k}$ and $\mathcal{M}_{k-1}\bigcap\mathcal{M}_{T}\subsetneq\mathcal{M}_{T}$ (i.e., the first $k-1$ candidates are underfitted and the last $S_{1}-k$ candidates are overfitted). We remain to make structural regularities on the candidate set to guarantee the reliability of subsequent inference.

We consider the scenario that $\mathcal{S}_{1}$ contains the OC model. Without loss of generality, the index of the OC candidate is specified as $1\leq s_{*}\leq S_{1}$. Similar to model averaging literature for generalized linear model, we define the infeasible KL loss among candidates as $\xi_{n}=\inf_{\bm{\omega}\in\mathcal{H}_{1}}\mathrm{KL}(\eta_{*}(\bm{\omega}))$, where

[TABLE]

Rewriting $\bm{\omega}=(\bm{\omega}^{\top}_{u},\omega_{s_{*}},\bm{\omega}^{\top}_{o})^{\top}$ and $\bar{\mathcal{H}}_{1}=\{\bm{\omega}_{u}\in[0,1]^{s_{*}-1}:0<\sum^{s_{*}-1}_{k=1}\omega_{k}\leq 1\}$, we further present a weighted infeasible KL loss among underfitted models as $\bar{\xi}_{n}=\inf_{\bm{\omega}_{u}\in\bar{\mathcal{H}}_{1}}(\sum^{s_{*}}_{k=1}\omega_{k})^{-1}\mathrm{KL}(\eta_{*}(\bm{\omega}_{u}))$, where

[TABLE]

The following conditions are emerged for deriving the theorem:

1. (A8)

$n^{1/2}\xi^{-1}_{n}=o_{p}(S_{1}^{-1/2}|\mathcal{M}_{S_{1}}|^{-1})$. 2. (A9)

The penalty factor satisfies: (i) $\frac{n^{1/2}S^{1/2}_{1}|\mathcal{M}_{S_{1}}|}{\lambda_{n}(|\mathcal{M}_{s_{*}+1}|-|\mathcal{M}_{s_{*}}|)}=o(1)$ and (ii) $\bar{\xi}^{-1}_{n}\lambda_{n}|\mathcal{M}_{S_{1}}|=o_{p}(1)$.

Condition (A8) reflects that orders of the model size and the number of candidates are both related to KL divergence of the best-but-infeasible weighted model. It is regarded as an extension in high-dimensional case. When $S_{1}$ and $|\mathcal{M}_{S_{1}}|$ are finite, this condition degenerates to $n\xi_{n}^{-2}=o_{p}(1)$, a classical constraint on model averaging inference (e.g., Ando and Li 2014; Zhang et al. 2016; Zhang et al. 2020; Fang et al. 2022). Condition (A9) ensures that the OC model can be consistently selected by SMA. Specifically, the first part in (A9) illustrates that $\lambda_{n}$ should not be a constant, while the second part shows the order is dominated by $\bar{\xi}_{n}$, the latter is separated from KL divergence for demonstrating the “over-consistency” of weights (see condition (11) in Chen et al. 2023).

Theorem 1**.**

Under conditions (A1)–(A9), when $\mathcal{S}_{1}$ contains the OC model, we have for $n\to\infty$, (i) $\hat{\omega}_{s_{*}}\stackrel{{\scriptstyle p}}{{\to}}1$; (ii) $\hat{\bm{\theta}}(\hat{\bm{\omega}})\stackrel{{\scriptstyle p}}{{\to}}\bm{\theta}_{0}$.

Theorem 1 guarantees the consistency of nuisance parameter estimator. Notably, this result relies on the structural condition for candidates and will raise two potential problems: (1) how to construct $\mathcal{S}_{1}$ to cover the OC model, and (2) whether the estimator is valid for subsequent inference when the true mechanism deviates from (2.2). The first will be addressed in Section 5, while the second will be examined through numerical studies to assess prediction robustness.

3 SMA for conditional quantile function

In this section we intend to estimate (2.1) by SMA. Given $0<\tau<1$, we still need to divide covariates to construct submodels. Denote $\mathcal{S}_{2}=\{\mu^{(m)}_{i}:m=1,...,S_{2}\}$ as the set of total $S_{2}$ candidates, where the $m$-th model is written as

[TABLE]

and $\mathcal{I}_{m}\subset\{1,...,p\}$ is the index set of predictors. As mentioned previously, the SMA procedure is implemented by matching weights for candidate estimators. In fact, (2.2) implies that

[TABLE]

where

[TABLE]

and $\rho_{\tau}(\varepsilon)=\varepsilon(\tau-I(\varepsilon\leq 0))$ is named as quantile loss or check loss. Detailed derivation of (3.1) is attached in Part C of online supplementary material and, under MNAR, the parameter in $\mu^{(m)}_{i}$ hence can be consistently estimated as

[TABLE]

where $\hat{\gamma}(\hat{\bm{\omega}})$ is given in (2.4), and $\hat{M}_{0}\left(\cdot\right)$ is obtained by nonparametric regression such that

[TABLE]

Let $\bm{w}=(w_{1},...,w_{S_{2}})^{\top}\in\mathcal{H}_{2}$ and $\mathcal{H}_{2}=\{\bm{w}\in[0,1]^{S_{2}}:\sum^{S_{2}}_{m=1}w_{k}=1\}$. Similar to the procedure in Section 2.2, the SMA estimator of the $\tau$-th CQF has the form $\sum^{S_{2}}_{m=1}w_{m}\hat{\mu}^{{(m)}}$. Note that the complexity of (3.2) is almost higher than $O(pn^{2})$, adopting multifolds cross-validation will sensibly deteriorate computing efficiency (e.g., the complexity of JMA is approximating $O(pn^{3})$). To suppress computational costs and overfitting, we utilize regularization to optimize $\bm{w}$, of which the criterion is modified by

[TABLE]

here $\phi_{n}$ is a penalty factor that diverges to infinity, and

[TABLE]

Therefore, the weight is chosen by

[TABLE]

and the SMA estimator of (2.1) is

[TABLE]

4 Asymptotic behaviors

In this section, we demonstrate the asymptotic properties of CQF estimator. We firstly establish the consistency of the candidate estimator and then study the predicting efficiency. The latter is investigated by the asymptotic performance of selected weights. Let $\bm{\beta}^{(m)}_{*}\in\mathcal{B}_{m}$ be the unique minimizer of $E\rho_{\tau}(y-\bm{\beta}^{{(m)}^{\top}}\bm{x}^{(m)})$. The following conditions are required for theoretical inference:

1. (B1)

$\mathcal{B}_{m}$’s are compact over $m=1,...,S_{2}$. 2. (B2)

$\underline{c}=\min_{1\leq m\leq S_{2}}\underline{c}_{(m)}>0$ and is nonincreasing as $n\to\infty$, where $\underline{c}_{(m)}$ is the minimum eigenvalue of $E\bm{x}^{(m)}\bm{x}^{{(m)}^{\top}}$. 3. (B3)

(i) Uniformly for $\bm{x}\in\mathcal{X}$, $f_{y|x}(\cdot|\bm{x})$ is bounded away from [math] and $\infty$ uniformly on its support; (ii) the marginal densities of $x_{j}$’s ($j=1,...,p$) are uniformly bounded away from [math] and $\infty$. 4. (B4)

Let $|\bar{\mathcal{I}}|=\max_{1\leq m\leq S_{2}}|\mathcal{I}_{m}|$, (i) $|\bar{\mathcal{I}}|^{1/2}(\hat{\gamma}(\hat{\bm{\omega}})-\gamma_{0})=o_{p}(1)$; (ii) $|\bar{\mathcal{I}}|^{\kappa_{1}}n^{-\kappa_{2}}=o(1)$, where $\kappa_{1}$, $\kappa_{2}$ are both positive.

Condition (B1) is similar to (A1), and (B2) is common in high-dimensional random matrix theory. It is easy to see that $\underline{c}$ is inversely proportional to $n$ while the order is refined to promote the consistency. Condition (B3) is sufficient for the concentration of estimators. Condition (B4) restricts the order of candidate model size.

Lemma 2**.**

Under conditions (A1)–(A9), (B1)–(B4), when $\mathcal{S}_{1}$ contains the OC model, then for any $C>0$, there exist positive constants $C_{4}$ and $C_{6}$ such that

[TABLE]

Above lemma shows the tail behavior of estimators uniformly over candidate models. One benefit of this result is to determine the optimal order of $\left|\bar{\mathcal{I}}\right|$ and $S_{2}$ separately. Note that (4.1) does not directly reveal the consistency of parameter estimators unless further constraints are imposed. For instance, if we regard $\underline{c}=O(1)$ and $\left|\bar{\mathcal{I}}\right|=n^{\alpha}$ for $\alpha>0$, to promote the convergence of the exponential tail, one can get a necessary condition that $\alpha\in\left(0,1/2\right)$ (detailed discussion is found in the proof of Lemma 3). Of such case considered, $\log(S_{2})=o\left(n+n^{1-4\kappa_{2}}\left|\bar{\mathcal{I}}\right|^{4\kappa_{1}-2}\right)$, which means the number of candidates can be grown by an exponential order as well.

Next, we establish theoretical properties of $\hat{\bm{w}}$. Denote the out-of-sample finial prediction risk (FPR) of weighted CQF as

[TABLE]

Following conventional frameworks of frequentist model averaging, we assume additional conditions.

1. (B5)

$E\mu^{4}=O(1)$. 2. (B6)

(i) $n^{-1}\left|\bar{\mathcal{I}}\right|^{2}=o(\underline{c}^{2}\left|\bar{\mathcal{I}}\right|^{4\kappa_{1}}n^{-4\kappa_{2}})$; (ii) $S_{2}\lesssim n^{2\kappa_{2}}\left|\bar{\mathcal{I}}\right|^{-2\kappa_{1}-1}$; (iii) $\phi_{n}\left|\bar{\mathcal{I}}\right|=o(n)$.

Condition (B5) is same as assumption 4 (iii) in Wang et al. (2023). Condtion (B6) (i) guarantees an exponentially decaying tails in above lemma, and (iii) theoretically prevents the regularization term from dominating the weighting process. Condition (B6) (ii) determines the maximum dimensionality of candidates to establish the asymptotic optimality. Specially, for finite-dimensional candidate models it demonstrates that $\mathcal{S}_{2}$ can accommodate up to $O(n^{1/2})$ estimators with exponential tail. In high-dimensional settings, the dimensionality bounds can be determined through the following lemma:

Lemma 3**.**

Suppose the size of the largest candidate of CQF has the order $n^{\alpha}$, if Lemma 2 holds with condition (B6), then $\alpha\in(0,1/4)$.

Lemma 3 verifies the divergence of $|\bar{\mathcal{I}}|$ under complexity order assumptions, and the upper-bound is consistent with that in most high-dimensional inference. Crucially, (B6) (ii) is satisfied with $\phi_{n}\lesssim n^{1-\alpha}$, which takes an analogous penalty restriction in model selection (Shao 1997). Consequently, we establish the following property:

Theorem 2** (Asymptotic optimality).**

Given $0<\tau<1$, conditions (A1)–(A9) and (B1)–(B6), when $\mathcal{S}_{1}$ contains the OC model, (3.3) is asymptotically optimal in the sense that

[TABLE]

where $R_{n}=\sum^{5}_{l=1}r_{ln}$ and $r_{ln}$’s are given in the proof.

Theorem 2 reveals the primal result in model averaging theory. When all candiates in $\mathcal{S}_{2}$ are misspecified, it investigates that the SMA procedure has the optimal predicting accuracy in the sense that its FPR is asymptotically identical to that of the best-but-infeasible averaging estimator. Although weights are selected from respondents, FPR is calculated with fully out-of-sample observations. The property is asymptotically equivalent to that obtained in the full case. In fact, above theorem implies a principle that $|E\rho_{\tau}(\varepsilon)|$ is bounded away from [math] and infinity, further requires $E\varepsilon\neq 0$ for $\tau\in(0,1)$. To cover this situation we reform asymptotic optimality via “excess final prediction risk” (EFPR) defined as

[TABLE]

As discussed in Xu et al. (2023), $\mathrm{EFPR}_{n,\tau}(\bm{w})\geq 0$. The advantage of EFPR is given as an approximating unbiasedness in quantile framework, while it induces an extra order of $\inf_{\bm{w}\in\mathcal{H}_{2}}\rm{EFPR}_{n,\tau}(\bm{w})$. Based on Theorem 2, it can be naturally promoted into our SMA estimator by the following corollary:

Corollary 1**.**

Let

[TABLE]

Under same conditions in Theorem 2, if $\zeta^{-1}_{n}(R_{n}+n^{-1}\phi_{n}|\bar{\mathcal{I}}|)=o_{p}(1)$, (3.3) is asymptotically optimal in the sense that

[TABLE]

Remark 1**.**

It is worthy noting that asymptotic optimalities are established by the premise of the existence of the OC candidate in $\mathcal{S}_{1}$. This is because (3.1) yields an unbiased estimating equation iff the response probability is correctly specified. However, above properties may persist even under model misspecification in $\mathcal{S}_{1}$. In such scenario, $\hat{\bm{\beta}}^{(m)}_{n}$’s are no longer consistent to the pseudo-true parameters associated with the full case, and further propagates to FPR, EFPR.

Hereafter we focus on the property of $\hat{\bm{w}}$ in the sense that $\mathcal{S}_{2}$ contains correct models. We anticipate the estimator in (3.3) achieves the same convergency to least squares model averaging weights (see in, Zhang et al. 2020). Fortunately, the equivalence emerges through an algebraic transformation between quantile loss and mean squared error, and enables linearization of underfitted components. By Taylor’s expansion, we have

[TABLE]

where $\bar{\mu}(\bm{w})$ lies between $\mu$ and $\sum^{S_{2}}_{m=1}w_{m}\hat{\mu}^{(m)}$, and $\mathcal{I}_{\mathrm{cor}}\subset\{1,...,S_{2}\}$ is the index set of correct candidates in $\mathcal{S}_{2}$. Moreover, denote $\bar{\mathcal{H}}_{2}=\{\bm{w}\in[0,1]^{S_{2}}:\sum^{S_{2}}_{m=1}w_{m}=1,~\sum_{m\in\mathcal{I}_{\mathrm{cor}}}w_{m}\neq 1\}$ and

[TABLE]

Similar with (A9) (ii), we impose the following condition:

1. (B7)

$\max\left\{\left|\bar{\mathcal{I}}\right|^{\kappa_{1}+1/2}n^{-\kappa_{2}},n^{-1}\phi_{n}\left|\bar{\mathcal{I}}\right|,R_{n}\right\}=o(\bar{\zeta}_{n})$.

Lemma 4** (Over-consistency).**

Given $0<\tau<1$, conditions (A1)–(A9) and (B1)–(B7), SMA weight estimator is over-consistent, i.e.,

[TABLE]

The lemma implies that weight estimators asymptotically exclude incorrect models. However, it is not enough to investigate the consistency that assigns with the optimally correct model if $|\mathcal{I}_{\mathrm{cor}}|>1$. Let $s_{*}$ be the index of the OC model, for convenience we suppose that the first $s_{*}-1$ candidates in $\mathcal{S}_{1}$ are incorrect and the last $S_{2}-s_{*}$ are correct but $|\mathcal{I}_{s_{*}}|<|\mathcal{I}_{m}|$ for $m=s_{*}+1,...,S_{2}$. We have the following result:

Theorem 3** (Parsimony).**

Under conditions in Lemma 4, if $\phi_{n}$ further satisfies

[TABLE]

(3.3)* asymptotically concentrates on the OC candidate, i.e.,*

[TABLE]

Theorem 3 presents that SMA can control the number of nonzero coefficient estimator to mitigate overfitting in high-dimensional settings, contributing it as a parsimonious model averaging approach. The divergence constraint for $\phi_{n}$ highlights the necessity of incorporating the regularization term, which is by no means contradictory to (B6) (iii) because $|\bar{\mathcal{I}}|^{\kappa_{1}+1/2}n^{-\kappa_{2}}+R_{n}=o(1)$. Adopting theoretical bounds of $\phi_{n}$ requires very large $n$ and $|\bar{\mathcal{I}}|$, it is able to empirically select the order to balance computational flexibility (e.g., $\phi_{n}=O(\log(n))$).

5 Implement

In previous sections we illustrate the SMA procedure and its propertie, where the estimation of $E(\cdot|\bm{x},r=1)$ is obtained by nonparametric regression. The above method is advantageous for theoretically demonstrating the order of model size, however, it increases the complexity of objective functions and further leads to numerical confusions. In this section we consider an alternative semiparametric procedure, which is implemented by the idea of nonparametric sampling for conditional expectation.

Observed that for a measurable $g(\bm{x},y)$ and a large scalar $L$,

[TABLE]

where $\{\hat{y}_{l}:l=1,...,L\}$ is a set of random draws from the distribution of $y|\bm{x},r=1$. Based on the fact that ${f}_{1}(y|\bm{x})\propto f(y,\bm{x},r=1)$, the latter can be definitely estimated by kernel density technique. For a smoothed multivariate kernel function $K_{h}$, we have

[TABLE]

where $\bm{z}=(\bm{x}^{\top},y)^{\top}$. In this sense, $(\ref{sec2.eq3})$ is reformed as

[TABLE]

and conditional expectations in (3.2), (3.3) are replaced by

[TABLE]

The SMA procedure hence can be numerically executed through the following steps:

1. Step 1.

Draw $\{\hat{y}_{l}:l=1,...,L\}$ from (5.1). 2. Step 2.

Implement SMA procedure to estimate $\bm{\theta}$, where $\hat{m}_{1}(\bm{x}_{i};\gamma)$ is given in (5.2). 3. Step 3.

Implement SMA procedure for CQF, where $\hat{M}_{0}(\bm{x}_{i};\gamma)$ and $\hat{M}_{0}\left(\bm{x}_{i};\hat{\gamma}(\bm{\omega}),\bm{w}\right)$ are given in (5.3).

In numerical studies we employ Markov Chain Monte-Carlo (MCMC) sampling in step 1. However, the proposed methods with high-dimensional covariates inherently face the curse of dimensionality. Firstly, nonparametric components suffer from an exponentially deteriorating efficiency as dimension increases. To address this we use sufficient dimension reduction to prepare a $p\times d$ matrix $B_{y|x,r=1}$ satisfying $y\perp\bm{x}|B^{\top}_{y|x,r=1}\bm{x},r=1$, hence the conditional expectation in Sections 2–3 is replaced by $E_{1}(\cdot|B^{\top}_{y|x,r=1}\bm{x})$. Specifically we can apply dr package in R to determine the central subspace.

Secondly, one need to constrain the cardinality of the nested candidate set $\mathcal{S}_{1}$ to relieve computational burden. We implement this via feature screening, which is efficiently executed by existing programs (e.g., screenIID function in R). Since the screened index set is prone to severe overfitting, we wish to develop a shrinkage protocol to refine survived predictors. An enhanced algorithm is proposed in Appendix.

6 Numerical studies

In this section, we conduct numerical experiments and empirical analysis to compare the proposed estimation against existing model selection and model averaging methods. We implement MCMC sampling with $500$ total iterations and take the last $50$ draws as $\{\hat{y}_{l}:l=1,...,50\}$. We denote the estimator from Sections 2–3 as SMA, and the estimator from Section 5 as SMA-s. The competitive methods include: (i) AIC; (ii) BIC; (iii) smoothed AIC (denoted by SAIC); (iv) smoothed BIC (denoted by SBIC); (v) averaging estimator via quantile regression information criterion proposed by Lu and Su (2015) (denoted by QRIC); (vi) averaging estimator via $5$-folds cross-validation without regularization (denoted by CVMA); (vii) averaging estimator by inverse probability weighting (IPW) criterion (denoted by IPW-MNAR); (viii) averaging estimator by IPW criterion under MAR specification $\mathrm{logit}\{\pi_{i}(\bm{\theta})\}=\bm{\theta}^{\top}\bm{x}_{i}$ (denoted by IPW-MAR). For (iii) and (iv), corresponding weight estimators are assigned as

[TABLE]

For (v) and (vi), the dominant term in the criterion is same as that in $\widehat{\mathrm{QL}}_{n}(\bm{w})$, while the penalty in (vi) is adopted as in (3.9) of Lu and Su (2015). For (vii), the weight is optimized as

[TABLE]

where $\hat{\bm{\theta}}_{n}$ is obtained by SMA and each $\hat{\mu}^{(m)}_{i}$ is obtained by IPW. For (viii), the response probability is estimated by group-lasso in R package grplasso. The kernel function is chosen by multivariate Gaussian with $h=1.5n^{-1/3}\hat{\sigma}_{z}$, where $\hat{\sigma}_{z}$ is the estimated standard deviation of observations. Besides, penalty terms in (ii), (iv), (vii) and (viii) are equally set as $\phi_{n}=\log(n)$.

It should be pointed out that for SMA, SMA-s and (i)–(vi), the estimation procedure of $\gamma$ is same as that of QR model. For example, $\hat{Q}_{\tau}(y|\bm{x})$ and $\hat{\gamma}_{n}$ are both estimated by AIC in (i). To save the computational cost, all competitors are implemented with MCMC sampling estimation for $E_{1}(\cdot|\bm{x})$.

6.1 High-dimensional simulation design

We consider the following data generating process:

[TABLE]

Here the sample size $n$ is set as either $80$ or $150$, each of $\bm{x}_{i}$’s are i.i.d. from $N(\bm{0}_{p},\bm{I}_{p})$ and $p=5n$. The error term are generated by the following scenarios:

1. Scenario I (heteroskedasticity): $\varepsilon_{i}=e_{i}\sum^{p}_{j=p-9}x^{2}_{ij}/8$, where $e_{i}\sim N(0,1)$;
2. Scenario II (homoscedasticity): $\varepsilon_{i}=e_{i}$, where $e_{i}\sim N(0,1)$.

The active linear components are indexed as $j=1,...,s_{0}$ with $s_{0}=p/10$, $\beta_{j}=(-1)^{j}j^{-1}$ for $j=1,...,s_{0}$ and $\beta_{j}=0$ for $j=s_{0}+1,...,p$. As designed in many model averaging literature, $c$ is supplied to control the population $R^{2}=\{\mathrm{Var}(y_{i})-\mathrm{Var}(e_{i})\}/\mathrm{Var}(y_{i})$ such that $R^{2}$ varies across $\{0.1,0.3,0.5,0.7,0.9\}$.

To produce nonignorable nonresponse, we use $2$ response mechanisms:

1. Case (a) (logistic): $-\mathrm{logit}\{\pi_{i}(\bm{\theta})\}=\theta_{1}x_{1i}+\theta_{2}x_{2i}+\theta_{3}+\theta_{4}y_{i}$;
2. Case (b) (log-log): $\log\{-\log(1-\pi_{i}(\bm{\theta}))\}=\theta_{1}x_{1i}+\theta_{2}x_{2i}+\theta_{3}+\theta_{4}x_{1i}\sin(y_{i})$.

In each case, the true value of the parameter is set by: (a) $\bm{\theta}=(1,-1,-1,0.7)^{\top}$ in scenario I; $\bm{\theta}=(1,-1,-1,1)^{\top}$ in scenario II, and (b) $\bm{\theta}=(1,-1,0.3,1)^{\top}$ in scenario I; $\bm{\theta}=(1,-1,0.3,0.8)^{\top}$ in scenario II. We uniformly apply the logistic specification for MNAR to evaluate $\bm{\theta}$ during estimations (i)–(vii), SMA and SMA-s. It can be seen that the mechanism is correctly specified in case (a), while it is structurally misspecified in case (b). On the other hand, there is always a misspecification in (viii). To check the robustness for case (b), we additionally design a SMA procedure based on the true response mechanism, denoted as SMA-T. The variation of the missing rate with $R^{2}$ are visualized in Part A of online supplementary material, which floats around $31\%-37\%$.

To make preparations for constructing candidate models, we employ the idea of Zhang et al. (2020) to screen covariates. Specifically, we implement DC-SIS method (Li et al. 2012) to estimate $\mathcal{A}_{ry|x}\triangleeq\{j:ry\text{ functional depends on }x_{j}\}$, which retains the top $d_{n}=\lfloor 2n/(3\log(n))\rfloor$ corresponding to a decreasing order of distance correlation. Set $\hat{\mathcal{A}}_{ry|x}\triangleeq\{j_{1},...,j_{d_{n}}\}$, the linear components in the $m$-th nested models of $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are identically formed as: $\{x_{j_{1}},...,x_{j_{2m}}\}$ for $1\leq m<\lceil d_{n}/2\rceil$ and $\{x_{j_{1}},...,x_{j_{d_{n}}}\}$ for $m=\lceil d_{n}/2\rceil$. Besides, we use sliced inverse regression to compute $\hat{B}_{y|x,r=1}$.

Based on $200$ independent replications, we assess the normalized final prediction risk to quantify the performance of CQF prediction, which is given as their FPR divided by the optimal full-case single-estimator (Hansen 2007, Lu and Su 2015). FPR under the $r$-th replication is calculated by 100 pairs of out-of-sample observations $\{(\dot{\bm{x}}^{\top}_{s,r},\dot{y}_{s,r})^{\top}:s=1,...,100\}$ and is further averaged as

[TABLE]

To save space, we only report results for scenario I in Figures 1–4 (the other is found in supplementary material), which are summarized as: (1) When $\pi_{i}(\bm{\theta})$ is correct, two SMA estimators have outperformance than model selections, SAIC, SBIC, QRIC, CVMA and IPW-MAR, particularly in $n=80$. For $(\tau,n,R^{2})=(0.05,150,0.9)$, SMAs performs close to CVMA and, for $(\tau,n,R^{2})=(0.5,150,0.9)$ SMAs performs close to BIC, SBIC and CVMA. (2) When $\pi_{i}(\bm{\theta})$ is misspecified, two SMA estimators are significantly superior to others except IPW-MNAR, such advantage increases as $n$ decreasing. Particularly, SMA is more close to SMA-T for most of $R^{2}$. (3) Compared to IPW-MNAR, SMA generally performs better for $R^{2}<0.9$.

Above phenomenons reflect that SMA combines parsimony of BIC with the asymptotic optimality of CVMA. The outperformance of SMA over CVMA can be attributed to two key factors: firstly, multi-folds CVMA requires a data splitting process, which may compromise numerical precision with small sample size; secondly, CVMA incurs significantly higher computational costs. Note that IPW-MAR consistently exhibits the poorest performance, indicating that the misjudgement of the type of missingness can substantially impair prediction accuracy. Furthermore, the comparison between SMA and IPW-MNAR reveals that the latter suffers from prediction instability, particularly when the fitted model has lower $R^{2}$ or the sample size is small. This instability arises from extreme estimators of propensity score. To further illustrate the robustness under misspecified mechanism, Figure 5 shows boxplots for the case that $(n,R^{2})=(150,0.5)$. It is clearly shown that SMA has smaller estimation variance than IPW-MNAR. This advantage arises because SMA compensates for bias induced by leveraging distributional information from the observed respondents.

6.2 Performance of weight consistency

In this section, we design a simulation process to illustrate the consistency of SMA weight estimators. For the sake of space, the detail discussion is presented in Part A.2 of online supplementary material.

6.3 Empirical application

In this section we apply the SMA method for a micro-array genetic dataset. The detail discussion is presented in Part B of online supplementary material.

Conclusion

In the background of nonignorable nonresponse, we propose a semiparametric model averaging criterion for high-dimensional quantile regression. We demonstrate that the weight estimator is parsimonious when the optimal correct model is contained in the candidate set. The asymptotic optimality in the sense of FPR and EFPR are established as well. Furthermore, we develop two nonparametric estimators for conditional expectation. Numerical results confirmed the promise of the proposed methodology.

It should be pointed out that the asymptotic optimality of quantile model averaging theory is not strictly unbiased (even for EFPR, it has been verified to be weighted-unbiased with the weight corresponding to the conditional distribution of response variable). While Lu and Su (2015) established a asymptotic $C_{p}$-type criterion for quantile regression with homoscedastic errors, the extension to heteroscedastic models remains an open research problem. Another promising direction for future studies involves applying model averaging to obtain efficient propensity score estimation across different type of response mechanisms in nonresponse analysis.

Acknowledgements

The authors gratefully acknowledge Prof. Xinyu Zhang for his valuable suggestions on model averaging theory. We also thank the co-editor and anonymous referees for their insightful and constructive comments that significantly improved this work.

Funding

Xiong’s work is supported by National Natural Science Foundation of China (No.12401352), Postdoctoral Fellowship Program of CPSF (No.GZC20231022). Deng’s work is supported by Natural Sciences and Engineering Research Council of Canada (NSERC). Wang’s work is supported by National Natural Science Foundation of China (No.12271231, 12001229).

Statements and declarations

The authors report there are no competing interests to declare.

Appendix. A discussion for selecting nested models

Through the discussion of asymptotic properties for SMA, an important sufficient condition is the existence of the correct candidate in $\mathcal{S}_{1}$, that is, all active predictors are embraced in the largest nested model. For sparse linear logistic regression, one can employ the regularition path method to insert significant components. However, the nonparametric form in $\mathrm{Pr}(r=1|\bm{x})$ deteriorates both computational burden and numerical implementation difficulty. On the other hand, applying model-free methods to screen $x_{j}$’s dependence on $r$ is sensitive to the failure of efficiency under MNAR, which is mainly because active parts in $y|\bm{x},r=1$ may be incorrectly selected while they are inactive in $r=1|\bm{x},y$. In this section we discuss a fast and effect algorithm to create sub-models.

Let $\mathcal{A}_{r|x,y}=\{k:r\text{~depends on~}x_{k}\text{~among~}y\}$ and $\mathcal{A}_{y|x}=\{k:y\text{~depends on~}x_{k}\}$. Despite the considerable difficulty in selecting $\mathcal{A}_{r|x,y}$, we instead identify covariates that depend on $ry$ through feature screening, denoted as $\mathcal{A}_{ry|x}$. Existing studies have demonstrated the relationship $\mathcal{A}_{r|x,y}\bigcup\mathcal{A}_{y|x}\subset\mathcal{A}_{ry|x}$ (see in, Zhang et al. 2020, Zhao et al. 2021), which reduces our original problem to the evaluation of $\mathcal{A}_{ry|x}$. Different from variable screening process, it makes more sense to quantify the discrepancy between successive nested candidates through grouped variable importance (GVI). Formally, the GVI between models of size $m$ and $m-1$ is

[TABLE]

where $D_{m}$ measures the dependence between $\bm{x}^{(m)}$ and response in the $m$-th model (e.g., distance correlation in DC-SIS of Li et al. 2012, and weighted Cramér-von Mises distance in MV-SIS of Cui et al. 2015). Intuitively, GVI dynamics in the nested set reflects the incremental effect of covariates: $\mathrm{GVI}(m)>0$ if newly added covariates are influential for predicting, whereas $\mathrm{GVI}(m)\leq 0$ if relatively insignificant covariates are incorporated. This phenomenon mirrors the trade-off of the goodness-of-fit by information criteria. Another contribution for estimand $\mathcal{A}_{ry|x}$ is to assess the nested candidate for QR as well. Motivated by these properties, we propose a hybrid method integrating GVI with feature screening to determine the maximum size of nested models.

We firstly summarize the naive implementation in algorithm 1, where the sample distance correlation on $\mathbb{R}^{p}$ is typically adopted by DC-SIS (referred as $\widehat{\mathrm{dc}}$), which can be well-calculated in (2.10) of Sźekely et al. (2007).

As we can see, aforementioned algorithm is implemented through GVI between adjacent variables. However, it may inherently ignore variables that exhibit weak individual correlations but demonstrate strong joint correlations. Conversely, a substantial number of predictors with negligible contributions to inter-group associations may be redundantly selected. To address the limitation, we incorporate forward selection (denoted as F-step) to supplement the search for eliminated-but-active variables in $\hat{\mathcal{F}}_{n}\backslash\hat{\mathcal{A}}_{S}$, and use endogenous competition (denoted as E-step) to remove those with poor grouped importance out of $\hat{\mathcal{A}}_{S}$. The enhanced algorithm is formulated as a greedy method, which can be explicitly described as follows:

We employ the following DGP to simulate the proposed algorithm:

[TABLE]

for $i=1,...,n$, where $\beta_{1}=\beta_{3}=\beta_{5}=\sqrt{8}|\tau-0.5|$, $\beta_{2}=\beta_{4}=1$, $\beta_{6}=\beta_{7}=2(\tau-0.05)$ and $\beta_{j}=0$ with other $j$’s. $\bm{x}_{i}$ is drawn from $N(\bm{0}_{p},\bm{\Sigma})$ with $\Sigma=(\rho)_{p\times p}$ and $p=1000$. The error term $\varepsilon_{i}$ is generated from $N(0,1)$, and the response mechanism is specified as

[TABLE]

The estimation of CQF is implemented by SMA with different variable screening procedures. We set thresholds $d_{n}=\lfloor 2n/(3\log(n))\rfloor$, $s_{n}=\min_{j=1,...,d_{n}}\mathrm{GVI}(j)$, and the construction of $\mathcal{S}_{1}$, $\mathcal{S}_{2}$ are same as those in Section 5.1.

We consider (a) SMA by feature-screening (denoted as SMA) and (b) SMA by algorithm 1 (denoted as PS-SMA) as competitors, carry the simulation in terms of $\rho=\{0,0.5\}$ and $\tau=\{0.05,0.5\}$, respectively, and present the results averaged through $200$ replications. From Table 1 we can find that all scenarios exhibit nonresponse rates exceeding 35$\%$. SMA was implemented with a fixed size $|\hat{\mathcal{F}}_{n}|=14$ and PS-SMA, GPS-SMA both selected predictors from $\hat{\mathcal{F}}_{n}$. When covariates are uncorrelated, GPS-SMA conducted a very close FN to SMA’s, but PS-SMA often left out some significant parts. Conversely, when covariates are moderately correlated, GPS-SMA performed a very similar FP with PS-SMA but a smaller FN under $\tau=0.5$. In all cases GPS-SMA has a smallest FPR than others. These results imply that our greedy algorithm can efficiently enhance the post-selection of candidates, further improve the prediction of model averaging. Although the FN of GPS-SMA is unsatisfactory under $\tau=0.05$ and $\rho=0.5$, it did not affect the prediction of CQF. This mainly arises from the over-exclusion of significant variables under high missingness, however, the effect of excluded parts and response is inherently captured by other strongly correlated participants.