Overview of High-Dimensional Measurement Error Regression Models

Abstract

High-dimensional measurement error data are becoming more prevalent across various fields. Research on measurement error regression models has gained momentum due to the risk of drawing inaccurate conclusions if measurement errors are ignored. When the dimension p is larger than the sample size n, it is challenging to develop statistical inference methods for high-dimensional measurement error regression models due to the existence of bias, nonconvexity of the objective function, high computational cost and many other difficulties. Over the past few years, some works have overcome the aforementioned difficulties and proposed several novel statistical inference methods. This paper mainly reviews the current development on estimation, hypothesis testing and variable screening methods for high-dimensional measurement error regression models and shows the theoretical results of these methods with some directions worthy of exploring in future research.

1. Introduction

Measurement error data inevitably exist in applications, causing significant concern in various fields including biology, medicine, epidemiology, economics, finance and remote sensing. So far, there have been a wealth of research achievements on classical low-dimensional measurement error regression models under various assumptions. Numerous studies focus on parameter estimation for low-dimensional measurement error regression models, with the primary techniques listed below: (1) Corrected regression estimation methods [1]; (2) Simulation–Extrapolation (SIMEX) estimation methods [2,3]; (3) Deconvolution methods [4]; (4) Corrected empirical likelihood methods [5,6]. For more detailed discussions on other estimation and hypothesis testing methods for classical low-dimensional measurement error models, please refer to the literature [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29], as well as the monographs [30,31,32,33,34,35].

As one of the most popular research fields in statistics, high-dimensional regression has been widely used in various fields including genetics, economics, medical imaging, meteorology and sensor networks. Over the past two decades, various high-dimensional regression methods have been widely proposed such as Lasso [36], smoothly clipped absolute deviation (SCAD) [37], Elastic Net [38], Adaptive Lasso [39], Dantzig Selector [40], smooth integration of counting and absolute deviation (SICA) [41], and minimax concave penalty (MCP) [42], among many others. These methods have been widely applied to estimate regression coefficients while also achieving the goal of variable selection by adding penalties to objective functions; please refer to the literature review [43,44,45,46] as well as the monographs [47,48,49].

For the variable screening methods of ultrahigh-dimensional regression models where dimension p and sample size n satisfy $logp=O\left(n^{\kappa }\right),\kappa >0$, Fan and Lv [50] proposed the sure independence screening (SIS) method, which is a pioneering method in this field. For the estimation and variable selection of ultrahigh-dimensional regression models, it is suggested to apply the SIS method for variable screening first. Then, based on the variables screened in the first step, we can utilize regularization methods with penalties to estimate the regression coefficients and identify the significant variables simultaneously. Due to the operability and effectiveness of the SIS method in applications, numerous works have extended the method; see [51,52,53,54,55,56,57,58,59,60].

However, most of the aforementioned theories and applications for high-dimensional regression models focused on clean data. In the era of big data, researchers frequently collect high-dimensional data with measurement errors. Typical instances include gene expression data [61] and sensor network data [62]. The imprecise measurements are the result of poorly managed and defective data collection processes as well as the imprecise measuring instruments. It is well known that ignoring the influence of measurement errors will result in biased estimators and erroneous conclusions. Therefore, developing statistical inference methods for high-dimensional measurement error regression models has drawn a lot of interest.

Based on the types of measurement errors, research on high-dimensional measurement error regression models can be divided into the following three categories: covariates containing measurement errors; response variables containing measurement errors; both covariates and response variables containing measurement errors. In this paper, we mainly focus on the category where covariates contain measurement errors. When the dimension p is larger than the sample size n, parameter estimation can be challenging due to the nonconvexity of the penalized objective function caused by correction for the bias. This further makes it impossible to obtain the optimal solution of optimization problem. We utilize the following linear regression model to illustrate this problem

$$\mathit{y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon },$$

(1)

where $\mathit{y}={(y_1,\dots ,y_n)}^{\mathrm{T}}\in {\mathbb{R}}^n$ is the $n\times 1$ response vector, $\mathbf{X}={({\mathit{X}}_1,\dots ,{\mathit{X}}_n)}^{\mathrm{T}}\in {\mathbb{R}}^{n\times p}$ is the $n\times p$ fixed design matrix with ${\mathit{X}}_i={(x_{i1},\dots ,x_{ip})}^{\mathrm{T}}$, $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^{\mathrm{T}}\in {\mathbb{R}}^p$ is the sparse regression coefficient vector with only s nonzero components, and assume that model error vector $\mathit{\epsilon }={({\epsilon }_1,\dots ,{\epsilon }_n)}^{\mathrm{T}}\in {\mathbb{R}}^n$ is independent of $\mathbf{X}$. In order to obtain a sparse estimator of the true regression coefficient vector ${\mathit{\beta }}_0={({\beta }_{01},\cdots ,{\beta }_{0p})}^{\mathrm{T}}\in {\mathbb{R}}^p$, we can minimize the following penalized least-square objective function

$$\frac{1}{2n}{\parallel \mathit{y}-\mathbf{X}\mathit{\beta }\parallel }_2^2+{\parallel p_{\lambda }\left(\mathit{\beta }\right)\parallel }_1,$$

(2)

which is equivalent to minimizing

$$\frac{1}{2}{\mathit{\beta }}^{\mathrm{T}}\Sigma \mathit{\beta }-{\mathit{\rho }}^{\mathrm{T}}\mathit{\beta }+{\parallel p_{\lambda }\left(\mathit{\beta }\right)\parallel }_1,$$

(3)

where $\Sigma =n^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{X}$, $\mathit{\rho }=n^{-1}{\mathbf{X}}^{\mathrm{T}}\mathit{y}$, $p_{\lambda }(·)$ is a penalty function with regularization parameter $\lambda \ge 0$. If the covariates matrix $\mathbf{X}$ can be precisely measured, the penalized objective functions (2) and (3) are convex. Thus, we can obtain a sparse estimator of ${\mathit{\beta }}_0$ by minimizing the penalized objective function (2) or (3).

However, it is common that the covariates matrix $\mathbf{X}$ cannot be accurately observed in practice. Let $\mathbf{W}={({\mathit{W}}_1,\dots ,{\mathit{W}}_n)}^{\mathrm{T}}={\left(w_{ij}\right)}_{n\times p}$ be the observed covariates matrix with additive measurement errors satisfying $\mathbf{W}=\mathbf{X}+\mathbf{U}$, where $\mathbf{U}={({\mathit{U}}_1,\dots ,{\mathit{U}}_n)}^{\mathrm{T}}$ is the matrix of measurement errors, ${\mathit{U}}_i={(u_{i1},\dots ,u_{ip})}^{\mathrm{T}}$ follows a sub-Gaussian distribution with mean zero and covariance matrix ${\Sigma }_u$, and it is assumed to be independent of $(\mathbf{X},\mathit{y})$. To reduce the influence of measurement errors, Loh and Wainwright [63] proposed replacing $\Sigma$ and $\mathit{\rho }$ in the penalized objective function (3) by their consistent estimators $\widehat{\Sigma }=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}-{\Sigma }_u$ and $\tilde{\mathit{\rho }}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathit{y}$, respectively. Then, we can obtain the sparse estimator of ${\mathit{\beta }}_0$ by minimizing the following penalized objective function

$$\frac{1}{2}{\mathit{\beta }}^{\mathrm{T}}\widehat{\Sigma }\mathit{\beta }-{\tilde{\mathit{\rho }}}^{\mathrm{T}}\mathit{\beta }+{\parallel p_{\lambda }\left(\mathit{\beta }\right)\parallel }_1.$$

(4)

Note that when the dimension p is fixed or smaller than the sample size n, it can be guaranteed that $\widehat{\Sigma }$ is a positive definite or semi positive-definite matrix. It further ensures that the penalized objective function (4) remains convex. Thus, the global optimal solution of $\mathit{\beta }$ can be obtained by minimizing the penalized objective function (4).

However, for high-dimensional or ultrahigh-dimensional regression models, i.e., $p>n$ or $p\gg n$, there are two key problems: (i) the penalized objective function (4) is no longer convex and unbounded from below because the corrected estimator $\widehat{\Sigma }$ of $\Sigma$ is no longer a semi-positive definite matrix. This further makes it impossible to obtain the estimator of ${\mathit{\beta }}_0$ by minimizing the penalty objective function (4). (ii) In order to construct an objective function similar to that of standard Lasso and solve the corresponding optimization problem using R package “glmnet” or “lars”, it is necessary to decompose $\widehat{\Sigma }$ by Cholesky’s decomposition method and obtain the substitution of response vector and covariates matrix. However, this process results in an error accumulation and makes it challenging to guarantee valid theoretical results; please see the detailed discussions in [64,65].

For problem (i), Loh and Wainwright [63] changed the unconstrained optimization problem into a constrained optimization problem by adding restrictions to $\mathit{\beta }$. They suggested applying the projected gradient descent algorithm to solve the restricted optimization problem and acquire the global optimal solution of true regression coefficient vector ${\mathit{\beta }}_0$. Nevertheless, the penalized objective function of the optimization problem is still nonconvex. To address this issue, Datta and Zou [64] suggested substituting $\widehat{\Sigma }$ by its semi-positive definite projection matrix $\tilde{\Sigma }$, and they proposed convex conditioned Lasso (CoCoLasso). Furthermore, Zheng et al. [65] introduced a balanced estimation that prevented overfitting while maintaining the estimation accuracy by combining $l_1$ and concave penalty. Zhang et al. [66] further proposed an estimation method based on $L_0$ regularization. Tao et al. [67] constructed a modified least-squares loss function using a semi-positive definite projection matrix for the estimated covariance matrix and proposed a calibrated zero-norm regularized least squares (CaZnRLS) estimation of regression coefficients. Rosenbaum and Tsybakov [68,69] proposed a matrix uncertainty (MU) selector and its improved version compensated MU selector for high-dimensional linear models with additive measurement errors in covariates. Sørensen et al. [70] extended the MU selector to generalized linear models and developed the generalized matrix uncertainty (GMU) selector. Sørensen et al. [71] showed the theoretical results of relevant variable selection methods. Based on the MU selector, Belloni et al. [72] introduced an estimator that can achieve the minimax efficiency bound. They proved that the corresponding optimization problem can be converted into a second-order cone programming problem, which can be solved in polynomial time. Romeo and Thoresen [73] evaluated the performance of the MU selector in [68], nonconvex Lasso in [63], and CoCoLasso in [64] using simulation studies. Brown et al. [74] proposed a path-following iterative algorithm called Measurement Error Boosting (MEBoost), which is a computationally effective method for variable selection in high-dimensional measurement error regression models. Nghiem and Potgieter [75] introduced a new estimation method called simulation–selection–extrapolation (SIMSELEX), which used Lasso in the simulation step and group Lasso in the selection step. Li and Wu [76] established minimax convergence rates for the estimation of regression coefficients under a more general situation. Bai et al. [77] proposed a variable selection method for ultrahigh-dimensional linear quantile regression models with measurement errors. Jiang and Ma [78] drew on the idea of nonconvex Lasso in [63] and proposed an estimator of the regression coefficients for high-dimensional Poisson models with measurement errors. Byrd and McGee [79] developed an iterative estimation method for high-dimensional generalized linear models with additive measurement errors based on the imputation-regularized optimization (IRO) algorithm in [80]. However, the error accumulation issue mentioned in problem (ii) has not been addressed in the literature.

The aforementioned works place more emphasis on estimation and variable selection problems rather than hypothesis testing. For high-dimensional regression models with clean data, research on hypothesis testing problems has made significant progress under various settings in [81,82,83,84,85,86,87,88]. For high-dimensional measurement error models, the hypothesis testing methods are equally crucial. However, the bias and instability caused by measurement errors make hypothesis testing extremely difficult. Recently, some progress has been achieved in statistical inference methods. Based on a multiplier bootstrap, Belloni [89] constructed simultaneous confidence intervals for the target parameters in high-dimensional linear measurement error models. Focused on the case where a fixed number of covariates contain measurement errors, Li et al. [90] proposed a corrected decorrelated score test for parameters corresponding to the error-prone covariates and created asymptotic confidence intervals for them. Huang et al. [91] proposed a new variable selection method based on debiased CoCoLasso and proved that it can achieve false discovery rate (FDR) control. Jiang et al. [92] developed Wald and score tests for high-dimensional Poisson measurement error models.

Compared to the above estimation and hypothesis testing methods, there are relatively few screening techniques for ultrahigh-dimensional measurement error models. Nghiem et al. [93] introduced two screening methods named corrected penalized marginal screening (PMSc) and corrected sure independence screening (SISc) for ultrahigh-dimensional linear measurement error models.

This paper gives an overview of the estimation and hypothesis testing methods for high-dimensional measurement error regression models as well as the variable screening methods for ultrahigh-dimensional measurement error models. The rest of this paper is organized as follows. In Section 2, we review some estimation methods for linear models. We survey the estimation methods for generalized linear models in Section 3. Section 4 presents the recent advances in hypothesis testing methods for high-dimensional measurement error models. Section 5 introduces the variable screening techniques for ultrahigh-dimensional linear measurement error models. We conclude the paper with some discussions in Section 6.

Notation 1. 

Let ${\mathbb{S}}^p$ be the set of all $p\times p$ real symmetric matrices and ${\mathbb{S}}_{+}^p$ be the subset of ${\mathbb{S}}^p$ containing all positive semi-definite matrix in ${\mathbb{S}}^p$. We use $\left|\mathcal{A}\right|$ to denote the cardinality of set $\mathcal{A}$. Let $\mathcal{S}=\{j:{\beta }_{0j}\ne 0,j=1,\dots ,p\}$ be the index set of nonzero parameters. For a vector $\mathit{a}=(a_1,\dots ,a_m)\in {\mathbb{R}}^m$, let ${\parallel \mathit{a}\parallel }_q=({\sum }_{\ell =1}^m|a_{\ell }{{|}^q)}^{1/q},1\le q<\infty$ denote its $l_q$ norm, and write ${\parallel \mathit{a}\parallel }_{\infty }={max}_{1\le \ell \le m}\left|a_{\ell }\right|$. Denote by ${\mathit{a}}_{\mathcal{A}}\in {\mathbb{R}}^{\left|\mathcal{A}\right|}$ the subvector of $\mathit{a}$ with index set $\mathcal{A}\subset \{1,\dots ,m\}$. Denote by $\mathit{e}$ the vector of all ones. For a matrix $\mathbf{B}=\left(b_{ij}\right)$, let ${\parallel \mathbf{B}\parallel }_1={max}_j{\sum }_i\left|b_{ij}\right|,{\parallel \mathbf{B}\parallel }_{max}={max}_{i,j}\left|b_{ij}\right|$ and ${\parallel \mathbf{B}\parallel }_{\infty }={max}_i{\sum }_j\left|b_{ij}\right|.$ For constants a and b, define $a\vee b=max\{a,b\}$. We use c and C to denote positive constants that may vary throughout the paper. Finally, let $\stackrel{d}{\to }$ denote convergence in distribution.

2. Estimation Methods for Linear Models

This section mainly focuses on the linear model (1) with high-dimensional settings where the dimension p is larger than the sample size n. When the data can be observed precisely, we can estimate the true regression coefficient vector ${\mathit{\beta }}_0$ by minimizing the penalized objective function (2) or (3). However, we frequently come across cases where the measured covariates contain measurement errors. There are various types of measurement error data, and we primarily focus on the two categories below.

(1) Covariates with additive errors. The observed error-prone covariate ${\mathit{W}}_i={\mathit{X}}_i+{\mathit{U}}_i$, where the measurement error ${\mathit{U}}_i$ is independent of ${\mathit{X}}_i$ and independently generated from a distribution with mean zero and known covariance matrix ${\Sigma }_u$.

(2) Covariates with multiplicative errors. The observed error-prone covariates ${\mathit{W}}_i={\mathit{X}}_i\odot {\mathit{M}}_i$, where ⊙ denotes the Hadamard product, and the measurement error ${\mathit{M}}_i$ is independent of ${\mathit{X}}_i$ and follows from a distribution with mean ${\mathit{\mu }}_M$ and known covariance matrix ${\Sigma }_M$.

Our main goal is to obtain the sparse estimator $\widehat{\mathit{\beta }}$ of true regression coefficient vector ${\mathit{\beta }}_0$ in the presence of measurement errors. As we introduced in Section 1, we will run into the issue of the penalized objective function being nonconvex and unbounded from below after correcting the bias caused by measurement errors. This prevents us solving the optimization problem. Several works focused on this issue and proposed some estimation methods.

2.1. Nonconvex Lasso

In order to resolve the issue of the objective function being unbounded from below and unsolvable in the presence of measurement errors, Loh and Wainwright [63] added restrictions to regression coefficients $\mathit{\beta }$ and adopted an $l_1$ penalty. Then, the estimator of ${\mathit{\beta }}_0$ can be obtained by the following $l_1$-constrained quadratic program

$${\widehat{\mathit{\beta }}}_{\mathrm{NCL}}\in arg\underset{{\parallel \mathit{\beta }\parallel }_1\le c_0\sqrt{s}}{min}\left\{\frac{1}{2}{\beta }^{\mathrm{T}}\widehat{\Sigma }\beta -{\tilde{\rho }}^{\mathrm{T}}\mathit{\beta }+\lambda {\parallel \mathit{\beta }\parallel }_1\right\}=:arg\underset{{\parallel \mathit{\beta }\parallel }_1\le c_0\sqrt{s}}{min}\left\{\mathcal{L}\left(\beta \right)+{\lambda \parallel \mathit{\beta }\parallel }_1\right\},$$

(5)

where $c_0>0$ is a constant, $s=\left|\mathcal{S}\right|$ denotes the number of nonzero components of ${\mathit{\beta }}_0$, $\mathcal{L}\left(\mathit{\beta }\right)=2^{-1}{\beta }^{\mathrm{T}}\widehat{\Sigma }\beta -{\tilde{\rho }}^{\mathrm{T}}\mathit{\beta }$ is the loss function, $\widehat{\Sigma }$ and $\tilde{\rho }$ are the consistent estimators of covariance matrix $\Sigma$ of ${\mathit{X}}_i$ and marginal correlation coefficient vector $\mathit{\rho }$ of $({\mathit{X}}_i,{\mathit{y}}_i)$, and they may differ in terms of various kinds of measurement error data. Under the additive error setting,

$${\widehat{\Sigma }}_{\mathrm{add}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}-{\Sigma }_u,{\tilde{\rho }}_{\mathrm{add}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathit{y}.$$

(6)

Under the multiplicative error setting,

$${\widehat{\Sigma }}_{\mathrm{mul}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}\oslash ({\Sigma }_m+{\mu }_m{\mu }_m^{\mathrm{T}}),{\tilde{\rho }}_{\mathrm{mul}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathit{y}\oslash {\mu }_m,$$

(7)

where ⊘ denotes the elementwise division operator, and let $\widehat{\Sigma }={\widehat{\Sigma }}_{\mathrm{add}}$ or ${\widehat{\Sigma }}_{\mathrm{mul}}$ throughout the sequel. The reason for using “∈” rather than “=” in (5) is that several local minima might exist in the objective function. Note that this method still relies on a nonconvex objective function to obtain the estimator of ${\mathit{\beta }}_0$. Thus, we refer to it as “nonconvex Lasso”. It can be implemented by the R package “hdme” [94] at https://cran.r-project.org/web/packages/hdme/vignettes/hdme.html (accessed on 13 June 2023).

The nonconvexity of the penalized objective function makes it challenging to obtain the global minimum of the optimization problem (5). To solve the optimization problem (5), Loh and Wainwright [63] applied the projected gradient descent algorithm and demonstrated that even if the penalized objective function is nonconvex, the solution produced by this algorithm can reach the global minimum with high probability. The algorithm finds the global minimum in an iterative way as follows. At the $(k+1)$th iteration,

$${\mathit{\beta }}_{\mathrm{NCL}}^{(k+1)}=arg\underset{{\parallel \mathit{\beta }\parallel }_1\le c_0\sqrt{s}}{min}\left\{\mathcal{L}\left({\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}\right)+\nabla \mathcal{L}{\left({\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}\right)}^{\mathrm{T}}(\mathit{\beta }-{\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)})+\frac{\eta }{2}\parallel \mathit{\beta }-{\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}{\parallel }_2^2+\lambda {\parallel \mathit{\beta }\parallel }_1\right\},$$

(8)

where $\nabla \mathcal{L}\left(\mathit{\beta }\right)=\widehat{\Sigma }\mathit{\beta }-\tilde{\mathit{\rho }}$ is the gradient of loss function $\mathcal{L}\left(\mathit{\beta }\right)$, $\eta >0$ denotes the step-size parameter. For details of this algorithm, please see [63,95,96,97]. Loh and Wainwright [63] proved that the solution obtained by iteration (8) is quite near to the global minimum in both $l_1$-norm and $l_2$-norm under some conditions. Specifically, for all $t\ge 0$,

$$\begin{array}{cc}\hfill \parallel {\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}-{\widehat{\mathit{\beta }}}_{\mathrm{NCL}}{\parallel }_2^2& \le {\gamma }^k\parallel {\mathit{\beta }}_{\mathrm{NCL}}^{\left(0\right)}-{\widehat{\mathit{\beta }}}_{\mathrm{NCL}}{\parallel }_2^2+C_1\frac{logp}{n}\parallel {\widehat{\mathit{\beta }}}_{\mathrm{NCL}}-{\mathit{\beta }}_0{\parallel }_1^2+C_2{\parallel {\widehat{\mathit{\beta }}}_{\mathrm{NCL}}-{\mathit{\beta }}_0\parallel }_2^2,\hfill \\ \hfill \parallel {\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}-{\widehat{\mathit{\beta }}}_{\mathrm{NCL}}{\parallel }_1& \le 2\sqrt{k}\parallel {\mathit{\beta }}_{\mathrm{NCL}}^{\left(k\right)}-{\widehat{\mathit{\beta }}}_{\mathrm{NCL}}{\parallel }_2+2\sqrt{k}\parallel {\widehat{\mathit{\beta }}}_{\mathrm{NCL}}-{\mathit{\beta }}_0{\parallel }_2+2{\parallel {\widehat{\mathit{\beta }}}_{\mathrm{NCL}}-{\mathit{\beta }}_0\parallel }_1,\hfill \end{array}$$

where $C_1$ and $C_2$ are positive constants, and $\gamma \in (0,1)$ is a contraction coefficient independent of $(n,p,k)$. For the estimator ${\widehat{\mathit{\beta }}}_{\mathrm{NCL}}$ of the true regression coefficient vector ${\mathit{\beta }}_0$, Loh and Wainwright [63] showed that with any $c_0\ge {\parallel {\mathit{\beta }}_0\parallel }_2$ and $\lambda =O\left(\sqrt{logp/n}\right)$, the $l_q$-estimation error of ${\widehat{\mathit{\beta }}}_{\mathrm{NCL}}$ satisfies the bounds

$$\parallel {\widehat{\mathit{\beta }}}_{\mathrm{NCL}}-{\mathit{\beta }}_0{\parallel }_q=O\left(s^{1/q}\sqrt{\frac{logp}{n}}\right),q=1,2.$$

When $q=1$, the $l_1$-estimation error can reach the convergence rate $s\sqrt{logp/n}$; when $q=2$, the $l_2$-estimation error can reach the convergence rate $\sqrt{slogp/n}$. However, Loh and Wainwright [63] did not establish the variable selection consistency and oracle inequality for the prediction error of the nonconvex Lasso estimator.

2.2. Convex Conditioned Lasso

The nonconvex Lasso [63] overcomes the problem of unsolvability caused by the nonconvex objective function in the presence of measurement errors. However, there are some drawbacks to this method. First, the nonconvex Lasso solves the problem by adding constraints to $\mathit{\beta }$, but the penalized objective function remains nonconvex. It is well recognized that the convexity of the penalized objective function will be incredibly useful for theoretical analysis and computation. Second, two important unknown parameters $c_0$ and s are included in the optimization problem (5). These two parameters have a direct impact on the estimation results, but we are not sure about their magnitudes in applications. Third, Loh and Wainwright [63] have not established the variable selection results of the nonconvex Lasso estimator. To remedy these issues, Datta and Zou [64] proposed Convex Conditioned Lasso (CoCoLasso) based on a convex objective function, which possesses computational and theoretical superiority brought by convexity.

In order to construct the convex objective function, Datta and Zou [64] introduced a nearest positive semi-definite matrix projection operator for the square matrix, which is defined as

$${\left(\mathbf{A}\right)}_{+}=arg\underset{{\mathbf{A}}_1\ge 0}{min}{\parallel \mathbf{A}-{\mathbf{A}}_1\parallel }_{max},$$

(9)

where $\mathbf{A}$ is a square matrix. Let $\tilde{\Sigma }={\left(\widehat{\Sigma }\right)}_{+}$, and the alternating direction method of multipliers (ADMM) algorithm [98] can be utilized to derive $\tilde{\Sigma }$ from $\widehat{\Sigma }$. Based on $\tilde{\Sigma }$, the following convex objective function can be constructed, and it yields the CoCoLasso estimator

$${\widehat{\mathit{\beta }}}_{\mathrm{coco}}=arg\underset{\mathit{\beta }}{min}\left\{\frac{1}{2}{\beta }^{\mathrm{T}}\tilde{\Sigma }\beta -{\tilde{\rho }}^{\mathrm{T}}\mathit{\beta }+\lambda {\parallel \mathit{\beta }\parallel }_1\right\}.$$

(10)

When the covariates contain additive measurement errors,

$${\tilde{\Sigma }}_{\mathrm{add}}={\left({\widehat{\Sigma }}_{\mathrm{add}}\right)}_{+},{\tilde{\rho }}_{\mathrm{add}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathit{y},{\widehat{\Sigma }}_{\mathrm{add}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}-{\Sigma }_u.$$

(11)

When the covariates contain multiplicative measurement errors,

$${\tilde{\Sigma }}_{\mathrm{mul}}={\left({\widehat{\Sigma }}_{\mathrm{mul}}\right)}_{+},{\tilde{\rho }}_{\mathrm{mul}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathit{y}\oslash {\mu }_m,{\widehat{\Sigma }}_{\mathrm{mul}}=n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}\oslash ({\Sigma }_m+{\mu }_m{\mu }_m^{\mathrm{T}}).$$

(12)

Note that $\tilde{\Sigma }$ not only contributes to the construction of the convex objective function but also possesses the same level of estimation accuracy as $\widehat{\Sigma }$ in [63]. It can be guaranteed by the following equation

$$\parallel \tilde{\Sigma }{-\Sigma \parallel }_{max}\le \parallel \tilde{\Sigma }-\widehat{\Sigma }{\parallel }_{max}+\parallel \widehat{\Sigma }{-\Sigma \parallel }_{max}\le 2{\parallel \widehat{\Sigma }-\Sigma \parallel }_{max}.$$

Since $\tilde{\Sigma }$ is semi-positive definite, we can perform Cholesky decomposition on $\tilde{\Sigma }$. Then, the Cholesky factor of $\tilde{\Sigma }$ can be used to simplify computations by rewriting (10) as

$${\widehat{\mathit{\beta }}}_{\mathrm{coco}}=arg\underset{\mathit{\beta }}{min}\frac{1}{2n}\parallel \tilde{\mathit{y}}-\tilde{\mathbf{W}}{\mathit{\beta }\parallel }_2^2+\lambda {\parallel \mathit{\beta }\parallel }_1,$$

(13)

where $\tilde{\mathbf{W}}$ denotes the Cholesky factor of $\tilde{\Sigma }$ satisfying $n^{-1}{\tilde{\mathbf{W}}}^{\mathrm{T}}\tilde{\mathbf{W}}=\tilde{\Sigma }$, and $\tilde{\mathit{y}}$ is the vector satisfying $n^{-1}{\tilde{\mathbf{W}}}^{\mathrm{T}}\tilde{\mathit{y}}=\tilde{\mathit{\rho }}$. The penalized objective function in (13) is similar to that of the standard Lasso. Thus, we can utilize the coordinate descent algorithm to obtain the CoCoLasso estimator; please see the details in [64,99,100]. Theoretically, Datta and Zou [64] established the $l_q$-estimation $(q=1,2)$ and prediction error bounds of the CoCoLasso estimator. Suppose that

$$\psi =\underset{\mathit{\delta }\ne 0,\parallel {\mathit{\delta }}_{{\mathcal{S}}^c}{\parallel }_1\le 3{\parallel {\mathit{\delta }}_{\mathcal{S}}\parallel }_1}{min}\frac{{\mathit{\delta }}^{\mathrm{T}}\Sigma \mathit{\delta }}{{\parallel \mathit{\delta }\parallel }_2^2}>0.$$

For $s\sqrt{\zeta logp/n}<\lambda \le min\{{ϵ}_0,12{ϵ}_0\parallel {\mathit{\beta }}_{0\mathcal{S}}{\parallel }_{\infty }\}$, where $\zeta =max\{{\sigma }_{\epsilon }^4,{\sigma }_U^4,1\}$, ${ϵ}_0={\sigma }_U^2$, ${\sigma }_{\epsilon }^2$ and ${\sigma }_U^2$ are sub-Gaussian parameters of model error and measurement error, respectively, the CoCoLasso estimator ${\widehat{\mathit{\beta }}}_{\mathrm{coco}}$ satisfies that with probability at least $1-Cexp(-clogp)$,

$$\begin{array}{cc}\hfill \parallel {\widehat{\mathit{\beta }}}_{\mathrm{coco}}-{\mathit{\beta }}_0{\parallel }_q=& O\left(\frac{\lambda s^{1/q}}{\psi }\right),q=1,2,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill n^{-1/2}{\parallel \mathit{X}({\widehat{\mathit{\beta }}}_{\mathrm{coco}}-{\mathit{\beta }}_0)\parallel }_2=& O\left(\lambda \sqrt{\frac{s}{\psi }}\right).\hfill \end{array}$$

(15)

The fomulas (14) and (15) show the oracle inequalities for the $l_q$-estimation error with $q=1,2$ and prediction error. Furthermore, Datta and Zou [64] established the sign consistency of the CoCoLasso estimator under an additional irrepresentable condition and minimum signal strength condition. Meanwhile, there was no variable selection result provided for the nonconvex Lasso estimator ${\widehat{\mathit{\beta }}}_{\mathrm{NCL}}$ in [63]. Thus, the CoCoLasso estimation method not only enjoys the computational convenience of convexity but also possesses excellent theoretical results. However, when the dimension of covariates p is large, the computation of $\tilde{\Sigma }$ is expensive. To improve the computational efficiency, Escribe et al. [101] applied a two-step block descent algorithm and proposed a block coordinate descent convex conditioned Lasso (BDCoCoLasso), which is designed for the case in which the covariate matrix is only partially corrupted. CoCoLasso and BDCoCoLasso are now available in R package “BDcocolasso” at https://github.com/celiaescribe/BDcocolasso (accessed on 13 June 2023).

2.3. Balanced Estimation

CoCoLasso is effective in the parameter estimation of high-dimensional measurement error models, but it suffers from overfitting. To overcome this drawback, Zheng et al. [65] replaced the Lasso penalty in CoCoLasso with the combined $l_1$ and concave penalty and developed the balanced estimator, which can be obtained by

$${\widehat{\mathit{\beta }}}_{\mathrm{bal}}=arg\underset{\mathit{\beta }}{min}\left\{\frac{1}{2}{\beta }^{\mathrm{T}}\tilde{\Sigma }\beta -{\tilde{\rho }}^{\mathrm{T}}\mathit{\beta }+{\lambda }_0{\parallel \mathit{\beta }\parallel }_1+{\parallel p_{\lambda }\left(\mathit{\beta }\right)\parallel }_1\right\},$$

(16)

where ${\lambda }_0=c_1\sqrt{logp/n}$ is the regularization parameter for the $l_1$ penalty with $c_1$ being a positive constant, $p_{\lambda }\left(\mathit{\beta }\right)=[p_{\lambda }\left(\right|{\beta }_1\left|\right),\dots ,p_{\lambda }\left(\right|{\beta }_p{\left|\right)]}^{\mathrm{T}}$, and $p_{\lambda }\left(u\right),u\in [0,+\infty )$ is a concave penalty function with the tuning parameter $\lambda \ge 0$. The definitions of $\tilde{\Sigma }$ and $\tilde{\mathit{\rho }}$ are the same as those in (11) and (12) with the two kinds of measurement error data. In contrast to the CoCoLasso estimator, the balanced estimator strikes a perfect balance between prediction and variable selection. Surprisingly, excellent variable selection results promote the estimation and prediction accuracy of the balanced estimator. The simulation studies in [65] demonstrate the estimation and prediction accuracy as well as the better variable selection results of the balanced estimator. As for the asymptotic properties of ${\widehat{\mathit{\beta }}}_{\mathrm{bal}}$, Zheng et al. [65] established the oracle inequalities for the $l_q$-estimation and prediction error,

$$\begin{array}{cc}\hfill \parallel {\widehat{\mathit{\beta }}}_{\mathrm{bal}}-{\mathit{\beta }}_0{\parallel }_q=& O_p\left(\frac{{\lambda }_0{s}^{1/q}}{{\varphi }^2}\right),q=1,2,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill n^{-1/2}{\parallel \mathit{X}({\widehat{\mathit{\beta }}}_{\mathrm{bal}}-{\mathit{\beta }}_0)\parallel }_2=& O_p\left(\frac{{\lambda }_0\sqrt{s}}{\varphi }\right),\hfill \end{array}$$

(18)

where

$$\varphi =\underset{\mathit{\delta }\ne 0,\parallel {\mathit{\delta }}_{{\mathcal{S}}^c}{\parallel }_1\le 7{\parallel {\mathit{\delta }}_{\mathcal{S}}\parallel }_1}{min}\frac{n^{-1/2}{\parallel \mathit{X}\mathit{\delta }\parallel }_2}{\parallel {\mathit{\delta }}_{\mathcal{S}}{\parallel }_2\vee {\parallel {\mathit{\delta }}_{{\mathcal{S}}^c}^{*}\parallel }_2}>0,$$

and ${\mathit{\delta }}_{{\mathcal{S}}^c}^{*}\in {\mathbb{R}}^s$ contains the s largest absolute vaules of ${\mathit{\delta }}_{{\mathcal{S}}^c}$. It can be seen from (17) and (18) that the bounds of $l_q$-estimation ($q=1,2$) and prediction error are free of regularization parameter $\lambda$ for the concave penalty. Also, the upper bound of falsely discovered signs is provided in [65]. Denote $\mathrm{FS}\left(\widehat{\beta }\right)=\left|\{1\le j\le p:sgn\left({\widehat{\beta }}_j\right)\ne sgn\left({\beta }_{0,j}\right)\}\right|$; then

$$\mathrm{FS}\left(\widehat{\beta }\right)=O_p\left(\frac{{\lambda }_0^{2}s}{{\lambda }^2{\varphi }^4}\right).$$

(19)

From (19), we can see that if ${min}_{j\in \mathcal{S}}\left|{\beta }_{0j}\right|\gg \sqrt{slogp/n}$ such that ${\lambda }^2\gg {\lambda }_0^{2}s$, a balanced estimator can achieve sign consistency, which is stronger than the variable selection consistency. Compared with the balanced estimator, the CoCoLasso estimator requires an additional irrepresentable condition to achieve this property.

2.4. Calibrated Zero-Norm Regularized Least Square Estimation

The nearest positive semi-definite matrix projection operator defined in [64] solves the problem that the penalized objective function is nonconvex in high-dimensional measurement error models. However, with the constraint of the positive semi-definite matrix, the computation cost of $\tilde{\Sigma }$ is high. Tao et al. [67] demonstrated that as the dimension p increases, the time required to calculate $\tilde{\Sigma }$ using the ADMM algorithm will increase significantly. Thus, Tao et al. [67] suggested replacing $\tilde{\Sigma }$ with an approximation of $\widehat{\Sigma }$ that is easy to obtain but less precise. To achieve this purpose, consider the eigendecomposition of $\widehat{\Sigma }$ as follows

$$\widehat{\Sigma }=\mathbf{V}\mathrm{diag}({\theta }_1,\dots ,{\theta }_p){\mathbf{V}}^{\mathrm{T}},$$

where $\mathrm{diag}({\theta }_1,\dots ,{\theta }_p)$ is a diagonal matrix containing the eigenvalues of $\widehat{\Sigma }$ with ${\theta }_1\ge {\theta }_2\ge \cdots \ge {\theta }_p$, and $\mathbf{V}\in {\mathbb{R}}^{p\times p}$ is an orthonormal matrix consisting of the corresponding eigenvectors. Then, Tao et al. [67] substituted the Frobenius norm for the elementwise maximum norm in (9) and obtained a positive definite approximation of $\widehat{\Sigma }$ as follows

$${\tilde{\Sigma }}_F=arg\underset{\mathbf{W}\ge \xi \mathit{I}}{min}{\parallel \widehat{\Sigma }-\mathbf{W}\parallel }_F\mathrm{for}\mathrm{some}\xi >0.$$

(20)

Note that the optimal solution of (20) is the same as that of the problem

$$\underset{\mathbf{W}\ge \xi \mathit{I}}{min}{\parallel \widehat{\Sigma }-\mathbf{W}\parallel }_F^2.$$

(21)

Thus, we have

$${\tilde{\Sigma }}_F=\xi \mathbf{I}+{\mathsf{\Pi }}_{{\mathbb{S}}_{+}^p}(\widehat{\Sigma }-\xi \mathbf{I})=\mathbf{V}\mathrm{diag}[max({\theta }_1,\xi ),\dots ,max({\theta }_p,\xi )]{\mathbf{V}}^{\mathrm{T}},$$

(22)

where ${\mathsf{\Pi }}_{{\mathbb{S}}_{+}^p}(·)$ denotes the projection of a matrix on ${\mathbb{S}}_{+}^p$. Similar to $\tilde{\Sigma }$, we have ${\tilde{\Sigma }}_F=n^{-1}{\tilde{\mathbf{W}}}_F^{\mathrm{T}}{\tilde{\mathbf{W}}}_F$, where $n^{-1/2}{\tilde{\mathbf{W}}}_F$ is the Cholesky factor of ${\tilde{\Sigma }}_F$. Let ${\tilde{\mathit{y}}}_F$ be the vector satisfying $n^{-1}{\tilde{\mathbf{W}}}_F^{\mathrm{T}}{\tilde{\mathit{y}}}_F=\tilde{\mathit{\rho }}$. By some simple calculation, we can obtain that

$$\left\{\begin{array}{cc}\hfill {\tilde{\mathbf{W}}}_F& =\sqrt{n}\mathbf{V}diag\left(\sqrt{max\left({\theta }_1,\xi \right)},\dots ,\sqrt{max\left({\theta }_p,\xi \right)}\right){\mathbf{V}}^{\mathrm{T}},\hfill \\ \hfill {\tilde{\mathit{y}}}_F& =\sqrt{n}\mathbf{V}diag\left(\frac{1}{\sqrt{max\left({\theta }_1,\xi \right)}},\dots ,\frac{1}{\sqrt{max\left({\theta }_p,\xi \right)}}\right){\mathbf{V}}^{\mathrm{T}}\tilde{\mathit{\rho }}.\hfill \end{array}\right.$$

(23)

Based on Equation (22), ${\tilde{\Sigma }}_F$ can be obtained easily. This implies that computing ${\tilde{\Sigma }}_F$ requires substantially less time than computing $\tilde{\Sigma }$. However, the approximation accuracy of ${\tilde{\Sigma }}_F$ to $\widehat{\Sigma }$ is not as good as that of $\tilde{\Sigma }$ because minimizing the Frobenius norm may yield larger components compared with the elementwise maximum norm. To obtain an excellent estimator of ${\mathit{\beta }}_0$, it is reasonable to find a more effective regression method to replace Lasso. Tao et al. [67] considered the zero norm penalty and defined the following calibrated zero-norm regularized least squares (CaZnRLS) estimator

$${\widehat{\mathit{\beta }}}_{\mathrm{zn}}\in arg\underset{\mathit{\beta }\in {\mathbb{R}}^p}{min}\left\{\frac{1}{2n\lambda }\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F{\parallel }_2^2+{\parallel \mathit{\beta }\parallel }_0\right\}.$$

(24)

However, it is difficult to solve (24) directly. Thus, to give an equivalent form for (24) that can be solved, Tao et al. [67] defined

$$\varphi \left(u\right):=\frac{a-1}{a+1}{u}^2+\frac{2}{a+1}u(a>1),u\in \mathbb{R}.$$

It is easy to verify that for any $\mathit{\beta }\in {\mathbb{R}}^p$,

$${\parallel \mathit{\beta }\parallel }_0=\underset{\mathit{w}\in {\mathbb{R}}^p}{min}\left\{\sum _{i=1}^p\varphi \left(w_i\right):{(\mathit{e}-\mathit{w})}^{\mathrm{T}}\left|\mathit{\beta }\right|=0,0\le \mathit{w}\le \mathit{e}\right\},$$

(25)

where $\left|\mathit{\beta }\right|=\left(\right|{\beta }_1|,\dots ,|{\beta }_p{\left|\right)}^{\mathrm{T}}$. The Formula (25) implies that the optimization problem (24) can be rewritten as the following mathematical program with equilibrium constraints (MPEC)

$$\underset{\mathit{\beta },\mathit{w}\in {\mathbb{R}}^p}{min}\left\{\frac{1}{2n\lambda }\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F{\parallel }_2^2+\sum _{i=1}^p\varphi \left(w_i\right):{(\mathit{e}-\mathit{w})}^{\mathrm{T}}\left|\mathit{\beta }\right|=0,0\le \mathit{w}\le \mathit{e}\right\}.$$

(26)

Note that if the optimal solution of optimization problem (24) is ${\widehat{\mathit{\beta }}}^{*}$, then the corresponding optimal solution of optimization problem (26) is $({\widehat{\mathit{\beta }}}^{*},\mathrm{sign}\left(\right|{\widehat{\mathit{\beta }}}^{*}\left|\right))$.

However, it can be seen that the annoying nonconvexity is introduced by the restriction ${(\mathit{e}-\mathit{w})}^{\mathrm{T}}\left|\mathit{\beta }\right|=0$ in (26), and it is the cause of the difficulty in obtaining the estimator ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}$. Accordingly, Tao et al. [67] considered the following penalized version of optimization problem (26)

$$\underset{\mathit{\beta },\mathit{w}\in {\mathbb{R}}^p}{min}\left\{\frac{1}{2n\lambda }\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F{\parallel }_2^2+\sum _{i=1}^p\varphi \left(w_i\right)+\rho {(\mathit{e}-\mathit{w})}^{\mathrm{T}}\left|\mathit{\beta }\right|,0\le \mathit{w}\le \mathit{e}\right\},$$

(27)

where $\rho >0$ is the penalty parameter. Tao et al. [67] proved that the global optimal solution of optimization problem (27) with $\rho \ge \overline{\rho }:=\left(4a{L}_f\right){\left[(a+1)\lambda \right]}^{-1}$ is the same as that of optimization problem (26), where $L_f$ is the Lipschitz constant of the function $f\left(\mathit{\beta }\right):={\left(2n\right)}^{-1}{\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F\parallel }_2^2$ on the ball $\{\mathit{\beta }\in {\mathbb{R}}^p:\parallel \mathit{\beta }{\parallel }_2\le R\}$, and R is a constant. Thus, ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}$ can be obtained by solving the following optimization problem with $\rho \ge \overline{\rho }$

$${\widehat{\mathit{\beta }}}_{\mathrm{zn}}\in arg\underset{\mathit{\beta }\in {\mathbb{R}}^p,\mathit{w}\in [\mathit{0},\mathit{e}]}{min}\left\{\frac{1}{2n}{\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F\parallel }_2^2+\sum _{i=1}^p\lambda \left[\varphi \left(w_i\right)+\rho (1-w_i)\left|{\beta }_i\right|\right]\right\}.$$

(28)

Tao et al. [67] recommended using the multi-stage convex relaxation approach (GEP–MSCRA) to obtain ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}$. This approach solves (28) in an iterative way with the main steps summarized as follows.

Step 1. Initialize the algorithm with ${\mathit{w}}^{\left(0\right)}\in [\mathit{0},2^{-1}\mathit{e}]$, ${\rho }^{\left(0\right)}=1$, $\lambda >0$, $k=1$.

Step 2. Solve the following optimization problem and obtain ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}$

$${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}=arg\underset{\mathit{\beta }\in {\mathbb{R}}^p}{min}\left\{\frac{1}{2n}\parallel {\tilde{\mathbf{W}}}_F\mathit{\beta }-{\tilde{\mathit{y}}}_F{\parallel }_2^2+\lambda \sum _{i=1}^p(1-w_i^{(k-1)})\left|{\beta }_i\right|\right\}.$$

Step 3. If $k=1$, choose an appropriate ${\rho }^{\left(1\right)}>{\rho }^{\left(0\right)}$ using the information from $\parallel {\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(1\right)}{\parallel }_{\infty }$; if $1<k\le 3$, choose ${\rho }^{\left(k\right)}$ satisfying ${\rho }^{\left(k\right)}>{\rho }^{(k-1)}$; if $k>3$, let ${\rho }^{\left(k\right)}={\rho }^{(k-1)}$.

Step 4. Obtain $w_i^{\left(k\right)}(i=1,\dots ,p)$ through the following optimization problem

$$w_i^{\left(k\right)}=arg\underset{0\le w_i\le 1}{min}\left\{\varphi \left(w_i\right)-{\rho }^{\left(k\right)}{w}_i\left|{\widehat{\beta }}_{\mathrm{zn},i}^{\left(k\right)}\right|\right\}.$$

Step 5. Let $k\leftarrow k+1$ and repeat Steps 2–4 until the stopping conditions are satisfied.

Note that the initial ${\mathit{w}}^{\left(0\right)}$ in Step 1 is an arbitrary vector from the interval $[\mathit{0},2^{-1}\mathit{e}]$ rather than the feasible set $[\mathit{0},\mathit{e}]$ in (28). The reason is to obtain a better initial estimator ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(1\right)}$. In addition, $w_i^{\left(k\right)}$ in Step 4 has the following closed form based on the convexity of $\varphi$

$$w_i^{\left(k\right)}=min\left[1,max\left(\frac{(a+1){\rho }^{\left(k\right)}\left|{\beta }_i^{\left(k\right)}\right|-2}{2(a+1)},0\right)\right],i=1,\dots ,p.$$

Consequently, the primary calculation in each iteration is to solve a weighted $l_1$-norm regularized least square problem. Under some regularity conditions, ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}$ satisfies

$$\parallel {\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}-{\mathit{\beta }}_0{\parallel }_2=O_p\left(\lambda \sqrt{s}\right)\forall k\in {\mathbb{N}}^{+}.$$

(29)

It can be seen from (29) that the $l_2$-estimation error bound of CaZnRLS estimator possesses the same order as those of nonconvex Lasso and CoCoLasso estimators. Tao et al. [67] further showed that the error bound of ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{(k+1)}$ is better than that of ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}$ for all $k\in {\mathbb{N}}^{+}$. Furthermore, Tao et al. [67] demonstrated that GEP-MSCRA will produce a ${\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}$ such that $supp\left({\widehat{\mathit{\beta }}}_{\mathrm{zn}}^{\left(k\right)}\right)=supp\left({\mathit{\beta }}_0\right)$ in a finite number of iterations if the minimum nonzero value of the smallest nonzero entries of ${\mathit{\beta }}_0$ is not too small.

2.5. Linear and Conic Programming Estimation

In addition to the approaches mentioned above, another class of methods is based on the idea of the Dantzig selector to acquire an estimator of true regression coefficients ${\mathit{\beta }}_0$. Rosenbaum and Tsybakov [68] proposed the following matrix uncertainty (MU) selector

$${\widehat{\mathit{\beta }}}_{\mathrm{MU}}=arg\underset{\mathit{\beta }}{min}\left\{{\parallel \mathit{\beta }\parallel }_1:\parallel n^{-1}{\mathbf{W}}^{\mathrm{T}}{(\mathit{y}-\mathbf{W}\mathit{\beta })\parallel }_{\infty }\le \delta {\parallel \mathit{\beta }\parallel }_1+\lambda \right\},$$

(30)

where $\delta \ge 0$ and $\lambda \ge 0$ are tuning parameters depending on the level of measurement error U and model error $\epsilon$, respectively. The MU selector is available in R package “hdme” [94].

However, the $n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}$ is included in (30) rather than $n^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{X}$ due to the unobservability of $\mathbf{X}$. Obviously, the matrix $n^{-1}{\mathbf{W}}^{\mathrm{T}}\mathbf{W}$ contains bias caused by measurement errors. To address this issue, Rosenbaum and Tsybakov [69] proposed an improved version of the MU selector called the compensated MU selector. It is applicable to the case that the entries of measurement error ${\mathit{U}}_i$ are independent such that ${\sigma }_{U,j}^2=n^{-1}{\sum }_{i=1}^n\mathbb{E}\left(U_{ij}^2\right)$ is finite for $j=1,\dots ,p$. The compensated MU selector is defined as

$${\widehat{\mathit{\beta }}}_{\mathrm{CMU}}=arg\underset{\mathit{\beta }}{min}\left\{{\parallel \mathit{\beta }\parallel }_1:\parallel n^{-1}{\mathbf{W}}^{\mathrm{T}}(\mathit{y}-\mathbf{W}\mathit{\beta })+\widehat{\mathbf{D}}{\mathit{\beta }\parallel }_{\infty }\le \delta {\parallel \mathit{\beta }\parallel }_1+\lambda \right\},$$

(31)

where $\widehat{\mathbf{D}}$ is a diagonal matrix consisting of ${\widehat{\sigma }}_{U,j}^2,j=1,\dots ,p$, and constants $\delta$ and $\lambda$ are the same as those in (30). Rosenbaum and Tsybakov [69] showed that the $l_q$-estimation error of the estimator ${\widehat{\mathit{\beta }}}_{\mathrm{CMU}}$ satisfies

$$\parallel {\widehat{\mathit{\beta }}}_{\mathrm{CMU}}-{\mathit{\beta }}_0{\parallel }_q=O_p\left(s^{1/q}(\parallel {\mathit{\beta }}_0{\parallel }_1+1)\sqrt{\frac{logp}{n}}\right),1\le q\le \infty .$$

The MU selector and compensated MU selector provide two alternative estimation methods for high-dimensional measurement error models, but there remains an issue. The optimization problem in (31) may be nonconvex, and Rosenbaum and Tsybakov [69] did not offer a suitable algorithm to the general case. To remedy this issue, Belloni et al. [72] proposed the conic-programming-based estimator ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$. Consider the following optimization problem

$$\begin{array}{cc}\hfill & \underset{\mathit{\beta },t}{min}\left\{{\parallel \mathit{\beta }\parallel }_1+\kappa t\right\},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.\parallel n^{-1}{\mathbf{W}}^{\mathrm{T}}(\mathit{y}-\mathbf{W}\mathit{\beta })& +\widehat{\mathbf{D}}{\mathit{\beta }\parallel }_{\infty }\le \delta t+\lambda ,{\parallel \mathit{\beta }\parallel }_2\le t,t\in {\mathbb{R}}^{+},\hfill \end{array}$$

(32)

where $\kappa ,\delta$ and $\lambda$ are positive tuning parameters. Suppose that the solution of (32) is $({\widehat{\mathit{\beta }}}_{\mathrm{cp}};\widehat{t})$, then, ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$ is defined as the conic-programming-based estimator of true regression coefficients ${\mathit{\beta }}_0$. It is easy to verify that the optimization problem (32) can be solved efficiently in a polynomial time as it is a second-order cone programming problem. To analyze the asymptotic properties of ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$, assume that $\kappa \in [2^{-1},2]$, $\delta =O\left(\sqrt{logp/n}\right)$, and $\lambda =O\left(\sqrt{logp/n}\right)$. Then, Belloni et al. [72] showed that the $l_q$-estimation $(1\le q\le \infty )$ and prediction error of ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$ satisfy

$$\begin{array}{cc}\hfill \parallel {\widehat{\mathit{\beta }}}_{\mathrm{cp}}-{\mathit{\beta }}_0{\parallel }_q& =O_p\left(s^{1/q}(\parallel {\mathit{\beta }}_0{\parallel }_2+1)\sqrt{\frac{logp}{n}}\right),1\le q\le \infty ,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill n^{1/2}{\parallel \mathit{X}({\widehat{\mathit{\beta }}}_{\mathrm{cp}}-{\mathit{\beta }}_0)\parallel }_2& =O_p\left(s^{1/2}(\parallel {\mathit{\beta }}_0{\parallel }_2+1)\sqrt{\frac{logp}{n}}\right).\hfill \end{array}$$

(34)

In contrast to nonconvex Lasso in [63], the conic-programming-based estimator ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$ can achieve the convergence rate in (33) and (34) without any information of the parameters $\parallel {\mathit{\beta }}_0{\parallel }_1$, $\parallel {\mathit{\beta }}_0{\parallel }_2$ or s. Compared with the compensated MU selector in [69], the conic-programming based estimator ${\widehat{\mathit{\beta }}}_{\mathrm{cp}}$ can be obtained in the general case without the computational difficulty of nonconvexity.

3. Estimation Methods for Generalized Linear Models

The above methods are mainly for linear models. This section introduces the estimation methods for high-dimensional generalized linear models with measurement errors.

3.1. Estimation Method for Poisson Models

Count data are commonly encountered in various fields including finance, economics and social sciences. In order to analyze count data, Poisson regression models are a popular choice in practice. Jiang and Ma [78] studied the high-dimensional Poisson regression models with additive measurement errors and proposed a novel optimization algorithm to obtain the estimator of true regression coefficient vector ${\mathit{\beta }}_0$. Suppose that $Y_i$ is the response variable following a Poisson distribution satisfying $\mathbb{E}\left(Y_i\right|{\mathit{X}}_i)=exp\left({\mathit{X}}_i^{\mathrm{T}}\beta \right)$, where ${\mathit{X}}_i\in {\mathbb{R}}^p$ is an unobservable covariate. Its error-prone surrogate ${\mathit{W}}_i={\mathit{X}}_i+{\mathit{U}}_i$, and the measurement error ${\mathit{U}}_i$ follows from a sub-Gaussian distribution with known covariance matrix ${\Sigma }_u$. It is easy to verify that

$$\mathbb{E}\left\{Y_i{\mathit{W}}_i^{\mathrm{T}}\beta -exp\left({\beta }^{\mathrm{T}}{\mathit{W}}_i-{\beta }^{\mathrm{T}}{\Sigma }_u\beta /2\right)\mid {\mathit{X}}_i,Y_i\right\}=Y_i{\mathit{X}}_i^{\mathrm{T}}\beta -exp\left({\beta }^{\mathrm{T}}{\mathit{X}}_i\right).$$

(35)

From (35), Jiang and Ma [78] imposed a restriction on $\mathit{\beta }$ similar to it in [63] and estimated $\mathit{\beta }$ by solving the following optimization problem

$${\widehat{\mathit{\beta }}}_{\mathrm{p}}=arg\underset{{\parallel \mathit{\beta }\parallel }_1\le c_p\sqrt{s},{\parallel \mathit{\beta }\parallel }_2\le c_p}{min}\left\{\mathcal{L}\left(\mathit{\beta }\right)+{\lambda \parallel \mathit{\beta }\parallel }_1\right\},$$

(36)

where

$$\mathcal{L}\left(\beta \right)=-\frac{1}{n}\sum _{i=1}^n\left\{Y_i{\mathit{W}}_i^{\mathrm{T}}\beta -exp\left({\beta }^{\mathrm{T}}{\mathit{W}}_i-{\beta }^{\mathrm{T}}{\Sigma }_u\beta /2\right)\right\}.$$

(37)

The estimator ${\widehat{\mathit{\beta }}}_{\mathrm{p}}$ can be obtained by the composite gradient descent algorithm. Specifically, at the $(k+1)$th iteration, first solve the following optimization problem without any restrictions on $\mathit{\beta }$

$${\tilde{\beta }}_{\mathrm{p}}^{(k+1)}=arg\underset{\beta }{min}\left\{\partial \mathcal{L}\left({\beta }_{\mathrm{p}}^{\left(k\right)}\right)/\partial {\beta }^{\mathrm{T}}(\beta -{\beta }^{\left(k\right)})+\eta /2\parallel \beta -{\beta }^{\left(k\right)}{\parallel }_2^2+\lambda {\parallel \beta \parallel }_1\right\},$$

where $\eta >0$ is a stepsize parameter. Next, apply the projection method in [95] to project ${\tilde{\beta }}_{\mathrm{p}}^{(k+1)}$ onto the $l_1$ ball with radius $c_p\sqrt{s}$ and produce ${\breve{\beta }}_{\mathrm{p}}^{(k+1)}$. If $\parallel {\breve{\beta }}_{\mathrm{p}}^{(k+1)}{\parallel }_2>c_p$, let ${\widehat{\mathit{\beta }}}_{\mathrm{p}}^{(k+1)}={\breve{\beta }}_{\mathrm{p}}^{(k+1)}{c}_p/{\parallel {\breve{\beta }}_{\mathrm{p}}^{(k+1)}\parallel }_2$; otherwise, let ${\widehat{\mathit{\beta }}}_{\mathrm{p}}^{(k+1)}={\breve{\beta }}_{\mathrm{p}}^{(k+1)}$. Repeat the above steps until the stopping condition is satisfied. Jiang and Ma [78] proved the convergence of this algorithm. Under some regularity conditions, they further showed that the global minimum ${\widehat{\mathit{\beta }}}_{\mathrm{p}}$ of (36) satisfies

$$\parallel {\widehat{\beta }}_{\mathrm{p}}-{\beta }_0{\parallel }_q=O\left(s^{1/q}\lambda \right).$$

(38)

There is an usual requirement that $\lambda ⩾2{∥\partial \mathcal{L}\left(\beta \right)/\partial \beta ∥}_{\infty }$ in Poisson models, and the term ${∥\partial \mathcal{L}\left(\beta \right)/\partial \beta ∥}_{\infty }=O\left(\sqrt{n/logp}\right)$. Thus, the convergence rate of ${\widehat{\mathit{\beta }}}_{\mathrm{p}}$ is slower than those of nonconvex Lasso, CoCoLasso and balanced estimators in linear models.

3.2. Generalized Matrix Uncertainty Selector

The method in [78] is only designed for high-dimensional Poisson models with measurement errors. To develop a method that is applicable to generalized linear models, Sørensen et al. [70] drew on the idea of the MU selector and proposed the generalized matrix uncertainty (GMU) selector for high-dimensional generalized linear models with additive measurement errors.

Consider a generalized linear model with response variable Y distributed according to

$$f_Y(y;\theta ,\varphi )=exp\left\{\frac{y\theta -b\left(\theta \right)}{a\left(\varphi \right)}+c(y,\varphi )\right\},$$

where $\theta ={\mathit{X}}^{\mathrm{T}}{\mathit{\beta }}_0$, $\mathit{X}\in {\mathbb{R}}^p$ are the covariates. The expected response is given by the mean function $\mu \left(\theta \right)=b^{\prime }\left(\theta \right)$, and the Taylor expansion of the mean function $\mu \left({\mathit{X}}_i^{\mathrm{T}}{\mathit{\beta }}_0\right)$ at point ${\mathit{W}}_i^{\mathrm{T}}{\mathit{\beta }}_0$ is

$$\mu \left({\mathbf{X}}_i^{\mathrm{T}}{\beta }_0\right)=\sum _{\ell =0}^{\infty }\frac{{\mu }^{(\ell )}\left({\mathit{W}}_i^{\mathrm{T}}{\beta }_0\right)}{\ell !}{(-{\mathit{U}}_i^{\mathrm{T}}{\beta }_0)}^{\ell },$$

(39)

where ${\mu }^{(\ell )}(·)$ is the ℓth derivative of function $\mu (·)$. With the Taylor expansion (39) of the mean function, the generalized matrix uncertainty selector can be defined as

$$\begin{array}{cc}\hfill & {\widehat{\mathit{\beta }}}_{\mathrm{GMU}}^L=arg\underset{\mathit{\beta }}{min}\left\{{\parallel \mathit{\beta }\parallel }_1:\mathit{\beta }\in {\Theta }^L\right\},\hfill \\ \hfill {\Theta }^L=\left\{\mathit{\beta }\in {\mathbb{R}}^p:\underset{1\le j\le p}{max}\right.& \left.\left|\frac{1}{n}{w}_{ij}[Y_i-\mu \left({\mathit{W}}_i^{\mathrm{T}}\mathit{\beta }\right)]\right|\le \lambda +\sum _{\ell =1}^L\frac{{\delta }^{\ell }}{\ell !\sqrt{n}}{\parallel \mathit{\beta }\parallel }_1^{\ell }{\parallel {\mathit{\mu }}^{(\ell )}\left(\mathit{W}\mathit{\beta }\right)\parallel }_2\right\},\hfill \end{array}$$

(40)

where $\delta$ is the positive parameter satisfying

$${\parallel \mathit{U}\parallel }_{\infty }\le \delta ,{\mathit{\mu }}^{(\ell )}\left(\mathit{W}\mathit{\beta }\right)={[{\mu }^{(\ell )}\left({\mathit{W}}_1^{\mathrm{T}}\mathit{\beta }\right),\dots ,{\mu }^{(\ell )}\left({\mathit{W}}_n^{\mathrm{T}}\mathit{\beta }\right)]}^{\mathrm{T}}.$$

In practice, Sørensen et al. [70] recommended using $L=1$ for computational convenience and demonstrated that the first-order approximation produces satisfactory results.

To solve the optimation problem (40) and obtain the estimator ${\widehat{\mathit{\beta }}}_{\mathrm{GMU}}^L$, we can utilize the iterative reweighing algorithm. The main iteration step of the algorithm is stated as follows

$${\widehat{\beta }}_{\mathrm{GMU}}^{(k+1)}=arg\underset{\beta }{min}\left\{{\parallel \beta \parallel }_1:\frac{1}{n}{∥{\tilde{\mathbf{W}}}_{\mathrm{g}}^{\left(k\right)\mathrm{T}}({\tilde{\mathit{z}}}^{\left(k\right)}-{\tilde{\mathbf{W}}}_{\mathrm{g}}^{\left(k\right)}\beta )∥}_{\infty }\le \lambda +\sum _{\ell =1}^L\frac{{\delta }^{\ell }}{\ell !\sqrt{n}}{\parallel \beta \parallel }_1^{\ell }{\parallel {\mathit{V}}^{(\ell ,k)}\parallel }_2\right\},$$

(41)

where ${\tilde{\mathbf{W}}}_{\mathrm{g}}\in {\mathbb{R}}^{n\times p}$ is a matrix of the weighted error-prone surrogate of covariates with elements ${\tilde{w}}_{\mathrm{g},ij}^{\left(k\right)}=w_{ij}\sqrt{V_i^{(1,k)}}$, ${\tilde{\mathit{z}}}^{\left(k\right)}\in {\mathbb{R}}^n$ is a vector with the elements ${\tilde{z}}_i^{\left(k\right)}=z_i^{\left(k\right)}\sqrt{V_i^{(1,k)}}$,

$$z_i^{\left(k\right)}={\mathit{W}}_i^{\mathrm{T}}{\widehat{\beta }}_{\mathrm{GMU}}^{\left(k\right)}+\left[Y_i-\mu \left\{{\mathit{W}}_i^{\mathrm{T}}{\widehat{\beta }}_{\mathrm{GMU}}^{\left(k\right)}\right\}\right]{\mu }^{\prime }{\left\{{\mathit{W}}_i^{\mathrm{T}}{\widehat{\beta }}_{\mathrm{GMU}}^{\left(k\right)}\right\}}^{-1},i=1,\dots ,n,$$

and

$${\mathit{V}}^{(\ell ,k)}={\left[{\mu }^{(\ell )}\left\{{\mathit{W}}_1^{\mathrm{T}}{\widehat{\beta }}_{\mathrm{GMU}}^{\left(k\right)}\right\},\dots ,{\mu }^{(\ell )}\left\{{\mathit{W}}_n^{\mathrm{T}}{\widehat{\beta }}_{\mathrm{GMU}}^{\left(k\right)}\right\}\right]}^{\mathrm{T}}={\left(V_1^{(\ell ,k)},\dots ,V_n^{(\ell ,k)}\right)}^{\mathrm{T}},\ell =1,\dots ,L$$

is the weight vector in Taylor expansion with L terms. When $L=1$ is applied, it is easy to verify that (41) is a linear program. For more details about the algorithm, please see [70,102]. The GMU selector can be implemented by R package “hdme” [94]. However, Sørensen et al. [70] did not establish any asymptotic properties of the GMU selector.

4. Hypothesis Testing Methods

The aforementioned works on high-dimensional measurement error models mainly investigate estimation problems and numerical algorithms of optimization problems as well as the theoretical properties of estimators. Recently, some works have studied the hypothesis testing problems for high-dimensional measurement error regression models, which will be introduced in this section.

4.1. Corrected Decorrelated Score Test

The above methods are proposed under the setting that all covariates are corrupted. In practice, it is common that not all covariates are measured with errors. Thus, Li et al. [90] investigated high-dimensional measurement error models where a fixed number of covariates contain measurement errors and proposed statistical inference methods for the regression coefficients corresponding to these covariates.

Consider the following high-dimensional linear model with one of the covariates containing additive errors

$$\left\{\begin{array}{c}{y}_i={\beta }_0{X}_i+{\mathit{\gamma }}_0^{\mathrm{T}}{\mathit{Z}}_i+{\epsilon }_i,\hfill \\ W_i=X_i+U_i,i=1,\dots ,n,\hfill \end{array}\right.$$

where $X_i\in \mathbb{R}$ is an unobservable covariate, $W_i$ is its error-prone surrogate, and ${\mathit{Z}}_i\in {\mathbb{R}}^{p-1}$ is an observed covariate precisely. The measurement error $U_i$ follows from sub-Gaussian distribution with mean zero and variance ${\sigma }_U^2$, and $U_i$ is independent of $(X_i,{\mathit{Z}}_i,{\epsilon }_i)$. Denote $\mathit{y}={(y_1,\dots ,y_n)}^{\mathrm{T}},\mathit{X}={(X_1,\dots ,X_n)}^{\mathrm{T}},\mathit{W}={(W_1,\dots ,W_n)}^{\mathrm{T}}$ and $\mathit{Z}={({\mathit{Z}}_1,\dots ,{\mathit{Z}}_n)}^{\mathrm{T}}.$ This subsection aims to test the hypothesis:

$$H_0:{\beta }_0={\beta }^{*}⟷H_1:{\beta }_0\ne {\beta }^{*}({\beta }^{*}\in \mathbb{R}),$$

and construct a confidence interval for ${\beta }_0$ under high-dimensional settings.

Since we are only concerned with the inference of the parameter $\beta$, then the parameter $\mathit{\gamma }$ is regarded as a nuisance. Following the idea in [85], Li et al. [90] defined the corrected score function as

$${\mathit{S}}_{\mathit{\theta }}\left(\mathit{\theta }\right)=\widehat{\Sigma }\mathit{\theta }-\widehat{\mathit{\rho }}=\frac{1}{n}\sum _{i=1}^n{\mathit{S}}_{i\mathit{\theta }}\left(\mathit{\theta }\right)=\left(\begin{array}{c}{S}_{\beta }(\beta ,\mathit{\gamma })\hfill \\ {\mathit{S}}_{\mathit{\gamma }}(\beta ,\mathit{\gamma })\hfill \end{array}\right)=\left(\begin{array}{c}{\widehat{\Sigma }}_{11}\beta +{\widehat{\Sigma }}_{12}\mathit{\gamma }-{\widehat{\rho }}_1\hfill \\ {\widehat{\Sigma }}_{21}\beta +{\widehat{\Sigma }}_{22}\mathit{\gamma }-{\widehat{\mathit{\rho }}}_2\hfill \end{array}\right),$$

where $\mathit{\theta }={(\beta ,{\mathit{\gamma }}^{\mathrm{T}})}^{\mathrm{T}}$,

$$\widehat{\Sigma }=\left(\begin{array}{cc}{\widehat{\Sigma }}_{11}& {\widehat{\Sigma }}_{12}\\ {\widehat{\Sigma }}_{21}& {\widehat{\Sigma }}_{22}\end{array}\right)=\left(\begin{array}{cc}{\mathit{W}}^{\mathrm{T}}\mathit{W}/n-{\sigma }_U^2& {\mathit{W}}^{\mathrm{T}}\mathit{Z}/n\\ {\mathit{Z}}^{\mathrm{T}}\mathit{W}/n& {\mathit{Z}}^{\mathrm{T}}\mathit{Z}/n\end{array}\right)\mathrm{and}\widehat{\mathit{\rho }}=\left(\begin{array}{c}{\widehat{\rho }}_1\\ {\widehat{\mathit{\rho }}}_2\end{array}\right)=\left(\begin{array}{c}{\mathit{W}}^{\mathrm{T}}\mathit{y}/n\\ {\mathit{Z}}^{\mathrm{T}}\mathit{y}/n\end{array}\right)$$

are consistent estimators of $\Sigma ={(\mathit{X},\mathit{Z})}^{\mathrm{T}}(\mathit{X},\mathit{Z})/n$ and $\mathit{\rho }={(\mathit{X},\mathit{Z})}^{\mathrm{T}}\mathit{y}/n$, respectively. The corrected score covariance matrix is defined as

$$\mathbf{I}\left(\theta \right)=\mathbb{E}\left\{{\mathit{S}}_{i\mathit{\theta }}\left(\mathit{\theta }\right){\mathit{S}}_{i\mathit{\theta }}{\left(\mathit{\theta }\right)}^{\mathrm{T}}\right\}=\left(\begin{array}{cc}{I}_{\beta \beta }\hfill & {\mathit{I}}_{\beta \mathit{\gamma }}\hfill \\ {\mathit{I}}_{\mathit{\gamma }\beta }\hfill & {\mathbf{I}}_{\mathit{\gamma }\mathit{\gamma }}\hfill \end{array}\right).$$

To conduct statistical inference on the target parameter $\beta$, it is crucial to eliminate the influence of nuisance parameter $\mathit{\gamma }$. Thus, Li et al. [90] developed the corrected decorrelated score function for the target parameter $\beta$ as

$$S(\beta ,\mathit{\gamma })=S_{\beta }(\beta ,\mathit{\gamma })-{\mathit{\omega }}^{\mathrm{T}}{\mathit{S}}_{\mathit{\gamma }}(\beta ,\mathit{\gamma }),$$

where ${\mathit{\omega }}^{\mathrm{T}}={\mathit{I}}_{\beta \mathit{\gamma }}{\mathbf{I}}_{\mathit{\gamma }\mathit{\gamma }}^{-1}=\mathbb{E}\left(X_i{\mathit{Z}}_i^{\mathrm{T}}\right)\mathbb{E}{\left({\mathit{Z}}_i{\mathit{Z}}_i^{\mathrm{T}}\right)}^{-1}$. It easy to verify that $\mathbb{E}\left[S({\beta }_0,{\mathit{\gamma }}_0){\mathit{S}}_{\mathit{\gamma }}({\beta }_0,{\mathit{\gamma }}_0)\right]=0$, which indicates that $S(\beta ,\mathit{\gamma })$ and nuisance score function ${\mathit{S}}_{\mathit{\gamma }}(\beta ,\mathit{\gamma })$ are uncorrelated. Obviously, we can obtain that $Var\left[S(\beta ,\mathit{\gamma })\right]=I_{\beta \beta }-{\mathit{I}}_{\beta \mathit{\gamma }}{\mathbf{I}}_{\mathit{\gamma }\mathit{\gamma }}^{-1}{\mathit{I}}_{\mathit{\gamma }\beta }=:{\sigma }_{\beta \mid \mathit{\gamma }}^2$. Then, Li et al. [90] constructed the test statistic and the confidence interval for ${\beta }_0$ based on the estimated decorrelated score function. This statistical inference procedure is summarized as follows.

Step 1. Apply the CoCoLasso estimation method in [64] to calculate initial estimator $\tilde{\mathit{\theta }}={(\tilde{\beta },{\tilde{\mathit{\gamma }}}^{\mathrm{T}})}^{\mathrm{T}}$, and utilize the following Dantzig-type estimator to estimate $\mathit{\omega }$

$$\widehat{\mathit{\omega }}=arg\underset{\mathit{\omega }}{min}{\parallel \mathit{\omega }\parallel }_1,\mathrm{s}.\mathrm{t}.{\parallel {\widehat{\Sigma }}_{12}-{\mathit{\omega }}^{\mathrm{T}}{\widehat{\Sigma }}_{22}\parallel }_{\infty }\le {\lambda }^{\prime },$$

where ${\lambda }^{\prime }=O\left(\sqrt{logp/n}\right)$.

Step 2. Estimate the decorrelated score function by

$$\widehat{S}(\beta ,\tilde{\mathit{\gamma }})=S_{\beta }(\beta ,\tilde{\mathit{\gamma }})-{\widehat{\mathit{\omega }}}^{\mathrm{T}}{\mathit{S}}_{\mathit{\gamma }}(\beta ,\tilde{\mathit{\gamma }}),$$

and calculate the test statistic $\widehat{T}=\sqrt{n}\widehat{S}\left({\beta }^{*},\tilde{\mathit{\gamma }}\right){\left({\widehat{\sigma }}_{\beta \mid \mathit{\gamma },H_0}^2\right)}^{-1/2}$, where

$$\begin{array}{cc}\hfill {\widehat{\sigma }}_{\beta \mid \gamma ,H_0}^2& ={\left.\left\{{\widehat{I}}_{\beta \beta }-{\widehat{\mathit{\omega }}}^{\mathrm{T}}{\widehat{\mathit{I}}}_{\mathit{\gamma }\beta }\right\}\right|}_{\beta ={\beta }^{*}}\hfill \\ & =({\widehat{\sigma }}_{\epsilon ,H_0}^2+{\beta }^{*2}{\sigma }_U^2)(1-{\widehat{\mathit{\omega }}}^{\mathrm{T}}{\widehat{\Sigma }}_{21})+{\beta }^{*2}\mathbb{E}\left(U_i^4\right)+{\widehat{\sigma }}_{\epsilon ,H_0}^2{\sigma }_U^2-{\beta }^{*2}{\sigma }_U^4.\hfill \end{array}$$

Step 3. Estimate $\beta$ as

$$\widehat{\beta }=\tilde{\beta }-\widehat{S}\left(\tilde{\theta }\right)/({\widehat{\Sigma }}_{11}-{\widehat{\omega }}^{\mathrm{T}}{\widehat{\Sigma }}_{21}),$$

and construct the $(1-\alpha )100\%$ confidence interval for ${\beta }_0$ as

$$\left[\widehat{\beta }-u_{1-\alpha /2}\sqrt{{\widehat{\sigma }}_{\beta }^2/n},\widehat{\beta }+u_{1-\alpha /2}\sqrt{{\widehat{\sigma }}_{\beta }^2/n}\right],$$

where $u_{1-\alpha /2}$ is the $(1-\alpha /2)$ quantile of standard normal distribution,

$${\widehat{\sigma }}_{\beta }^2={(1-{\widehat{\omega }}^{\mathrm{T}}{\widehat{\Sigma }}_{21})}^{-2}\left\{({\widehat{\sigma }}_{\epsilon }^2+{\widehat{\beta }}^2{\sigma }_U^2)(1-{\widehat{\omega }}^{\mathrm{T}}{\widehat{\Sigma }}_{21})+{\widehat{\beta }}^2\mathbb{E}\left(U_i^4\right)+{\widehat{\sigma }}_{\epsilon }^2{\sigma }_U^2-{\widehat{\beta }}^2{\sigma }_U^4\right\}$$

is the estimator of the asymptotic variance ${\sigma }_{\beta }^2$ of $\widehat{\beta }$, and ${\widehat{\sigma }}_{\epsilon }^2=n^{-1}{\sum }_{i=1}^n{(y_i-\widehat{\beta }{W}_i-{\tilde{\mathit{\gamma }}}^{\mathrm{T}}{\mathit{Z}}_i)}^2-{\widehat{\beta }}^2{\sigma }_U^2$ is the estimator of the variance ${\sigma }_{\epsilon }^2$ of ${\epsilon }_i$.

Note that the methods used to estimate $\mathit{\theta }$ and $\mathit{\omega }$ in Step 1 can be varying, as long as the corresponding estimators are consistent; please see more discussions in [90]. Li et al. [90] showed that under some regularity conditions,

$$\sqrt{n}\widehat{S}\left({\beta }^{*},\tilde{\mathit{\gamma }}\right){\left({\widehat{\sigma }}_{\beta \mid \mathit{\gamma },H_0}^2\right)}^{-1/2}\stackrel{d}{\to }\mathcal{N}(0,1)\mathrm{as}n\to \infty .$$

Furthermore, the asymptotic normality of the test statistic ${\widehat{T}}_n$ at local alternatives was also established in [90] without any additional condition. Li et al. [90] also constructed the asymptotic confidence interval for target parameter $\beta$ in Step 3 based on the asymptotic normality of $\widehat{\beta }$, which is given as follows

$$\sqrt{n}(\widehat{\beta }-{\beta }_0)=-{\left[\mathbb{E}\left\{{\left.\frac{\partial S\left(\beta ,{\gamma }_0\right)}{\partial \beta }\right|}_{\beta ={\beta }_0}\right\}\right]}^{-1}\sqrt{n}S\left({\beta }_0,{\gamma }_0\right)+o_P\left(1\right)\stackrel{d}{\to }\mathcal{N}(0,{\sigma }_{\beta }^2)\mathrm{as}n\to \infty ,$$

where ${\sigma }_{\beta }^2={\left\{\mathbb{E}\left(X_i^2\right)-{\mathit{\omega }}^{\mathrm{T}}\mathbb{E}\left(X_i{\mathit{Z}}_i\right)\right\}}^{-2}{\sigma }_{\beta \mid \mathit{\gamma },0}^2$, and

$${\sigma }_{\beta \mid \mathit{\gamma },0}^2=({\sigma }_{\epsilon }^2+{\beta }_0^2{\sigma }_U^2)\left\{1-{\omega }^{\mathrm{T}}\mathbb{E}\left(X_i{\mathit{Z}}_i\right)\right\}+{\beta }_0^2\mathbb{E}\left(U_i^4\right)+{\sigma }_{\epsilon }^2{\sigma }_U^2-{\beta }_0^2{\sigma }_U^4.$$

4.2. Wald and Score Tests for Poisson Models

In addition to linear models, researchers have made some progress on hypothesis-testing problems for Poisson models. Jiang et al. [92] studied hypothesis-testing problems for high-dimensional Poisson measurement error models, and they proposed Wald and score tests for the linear function of regression coefficients.

Consider the following hypothesis test

$$H_0:\mathbf{C}{\mathit{\beta }}_{0\mathcal{M}}=\mathit{b}⟷H_1:\mathbf{C}{\mathit{\beta }}_{0\mathcal{M}}=\mathit{b}+{\mathit{h}}_n\mathrm{for}\mathrm{some}{\mathit{h}}_n\in {\mathbb{R}}^r,$$

where $\mathbf{C}\in {\mathbb{R}}^{r\times m}$ is a matrix with $r\le m$, and ${\mathit{\beta }}_{0\mathcal{M}}\in {\mathbb{R}}^{m\times 1}$ is a subvector of the true regression coefficient vector ${\mathit{\beta }}_0={({\beta }_{01},\dots ,{\beta }_{0p})}^{\mathrm{T}}$ formed by ${\beta }_{0j}(j\in \mathcal{M})$. To construct a valid test statistic, Jiang et al. [92] drew on the idea of the estimation method in [78] and suggested estimating regression coefficients under the null hypothesis by

$${\widehat{\beta }}_{\mathrm{p}\mathrm{n}}=arg\underset{{\parallel \beta \parallel }_1⩽R_1,{\parallel \beta \parallel }_2⩽R_2}{min}\left\{\mathcal{L}\left(\beta \right)+p_{\lambda }\left({\beta }_{{\mathcal{M}}^c}\right)\right\},\mathrm{s}.\mathrm{t}.\mathbf{C}{\beta }_{\mathcal{M}}=\mathit{b},$$

(42)

where $p_{\lambda }(·)$ is a penalty function, and $\mathcal{L}\left(\mathit{\beta }\right)$ is defined in (37). Similarly, the following estimator of ${\mathit{\beta }}_0$ can be considered without assuming the null hypothesis

$${\widehat{\beta }}_{\mathrm{p}\mathrm{w}}=arg\underset{{\parallel \beta \parallel }_1⩽R_1,{\parallel \beta \parallel }_2⩽R_2}{min}\left\{\mathcal{L}\left(\beta \right)+p_{\lambda }\left({\beta }_{{\mathcal{M}}^c}\right)\right\}.$$

(43)

The estimators ${\widehat{\mathit{\beta }}}_{\mathrm{p}\mathrm{n}}$ and ${\widehat{\mathit{\beta }}}_{\mathrm{p}\mathrm{w}}$ can be obtained by the ADMM algorithm; for more details, please see [92]. It can be seen that optimization problems (42) and (43) can be distinguished from the method in (36) because we do not impose penalties on the components of the target parameter ${\mathit{\beta }}_{\mathcal{M}}$ to avoid forcing them to be zeros. Then, based on the above estimators of ${\mathit{\beta }}_0$, Jiang et al. [92] proposed the following score statistic and Wald statistic to test whether $\mathbf{C}{\beta }_{0\mathcal{M}}=\mathit{b}$ or not

$$\begin{array}{cc}\hfill T_S=& n{\left\{\frac{\partial \mathcal{L}\left(\widehat{\beta }\right)}{\partial {\beta }^{\mathrm{T}}}\right\}}_{\mathcal{M}\cup \mathcal{S}}{\mathbf{A}}^{\mathrm{T}}{\Psi }^{-1}({\widehat{\Sigma }}^r,\widehat{\mathbf{Q}},\widehat{\beta })\mathbf{A}{\left\{\frac{\partial \mathcal{L}\left(\widehat{\beta }\right)}{\partial \beta }\right\}}_{\mathcal{M}\cup \mathcal{S}},\hfill \\ \hfill T_W=& n{(\mathbf{C}{\widehat{\beta }}_{\mathrm{pw},\mathcal{M}}-\mathit{b})}^{\mathrm{T}}\Psi {({\widehat{\Sigma }}^r,\widehat{\mathbf{Q}},{\widehat{\beta }}_{\mathrm{pw}})}^{-1}(\mathbf{C}{\widehat{\beta }}_{\mathrm{pw},\mathcal{M}}-\mathit{b}),\hfill \end{array}$$

where $\mathbf{A}=\mathbf{C}\left[{\mathit{I}}_{m\times m},{\mathit{0}}_{m\times k}\right]{\widehat{\mathbf{Q}}}_{\mathcal{M}\cup \mathcal{S},\mathcal{M}\cup \mathcal{S}}^{-1}\left(\widehat{\beta }\right)$,

$$\Psi (\Sigma ,\mathbf{Q},\beta )\equiv \mathbf{C}\left[{\mathit{I}}_{m\times m},{\mathit{0}}_{m\times k}\right]{\mathbf{Q}}_{\mathcal{M}\cup \mathcal{S},\mathcal{M}\cup \mathcal{S}}^{-1}\left(\beta \right){\Sigma }_{\mathcal{M}\cup \mathcal{S},\mathcal{M}\cup \mathcal{S}}\left(\beta \right){\mathbf{Q}}_{\mathcal{M}\cup \mathcal{S},\mathcal{M}\cup \mathcal{S}}^{-1}\left(\beta \right){\left[{\mathit{I}}_{m\times m},{\mathit{0}}_{m\times k}\right]}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}},$$

${\widehat{\Sigma }}^r\left(\beta \right)$ and $\widehat{\mathbf{Q}}\left(\beta \right)$ are estimators of ${\Sigma }^r\left(\beta \right)$ and $\mathbf{Q}\left(\beta \right)=\mathbb{E}\left\{exp\left({\beta }^{\mathrm{T}}\mathit{X}\right)\mathit{X}{\mathit{X}}^{\mathrm{T}}\right\}$, respectively, and

$${\Sigma }^r\left(\beta \right)=\mathbb{E}\left[{\left\{Y_i{\mathbf{W}}_i-exp({\beta }^{\mathrm{T}}{\mathbf{W}}_i-{\beta }^{\mathrm{T}}{\Sigma }_u\beta /2)\left({\mathbf{W}}_i-{\Sigma }_u\beta \right)\right\}}^{\otimes 2}\right]$$

is the covariance of the residuals.

Jiang et al. [92] established the consistency of ${\widehat{\mathit{\beta }}}_{\mathrm{p}\mathrm{n}}$ and ${\widehat{\mathit{\beta }}}_{\mathrm{p}\mathrm{w}}$ with $\lambda$ larger than $O\left({\{logp/n\}}^{1/4}\right)$, $m=o\left({\{logp/n\}}^{1/2}\right)$ and $s=o\left({\{logp/n\}}^{1/2}\right)$. Furthermore, the asymptotic distributions of the two test statistics are established; specifically, as $n\to \infty$, we have

$$T_S\stackrel{d}{\to }{\chi }^2\left(r,n{\mathbf{h}}_n^{\mathrm{T}}{\Psi }^{-1}\left(\Sigma ,\mathit{Q},{\beta }_t\right){\mathbf{h}}_n\right),T_W\stackrel{d}{\to }{\chi }^2\left(r,n{\mathbf{h}}_n^{\mathrm{T}}{\Psi }^{-1}\left(\Sigma ,\mathit{Q},{\beta }_t\right){\mathbf{h}}_n\right).$$

Thus, we reject the null hypothesis if $T_S>{\chi }_{1-\alpha }^2\left(r\right)$ for the score test with the nominal significance level $\alpha >0$, and we reject the null hypothesis if $T_W>{\chi }_{1-\alpha }^2\left(r\right)$ for the Wald test, where ${\chi }_{1-\alpha }^2\left(r\right)$ is the $(1-\alpha )$ quantile of the chi-square distribution ${\chi }^2\left(r\right)$.

5. Screening Methods

As the dimensions of data become higher and higher, we often encounter ultrahigh-dimensional data. For the ultrahigh-dimensional models, we frequently reduce the dimensions using variable screening techniques and then apply other estimation or hypothesis-testing methods. The variable screening technique SIS [50] designed for ultrahigh-dimensional clean data has achieved great success and has been extended to various settings. SIS screens the variables according to the magnitudes of their marginal correlations with the response variable. Nghiem et al. [93] drew inspiration from the ideas of SIS in [50] and marginal bridge estimation in [103], and they proposed the corrected sure independence screening (SISc) method and corrected penalized marginal screening method (PMSc). Consider the following optimization problem

$$\begin{array}{cc}\hfill {\tilde{\mathit{\beta }}}_{\mathrm{sc}}=& arg\underset{\mathit{\beta }}{min}L\left(\mathit{\beta }\right)=arg\underset{\mathit{\beta }}{min}\left\{\sum _{j=1}^p{L}_j\left({\beta }_j\right)\right\}\hfill \\ \hfill =& arg\underset{\mathit{\beta }}{min}\left\{\frac{1}{n}\sum _{j=1}^p\left[\sum _{i=1}^n{\left(y_i-w_{ij}{\beta }_j\right)}^2-{\sigma }_j^2{\beta }_j^2+p_{\lambda }\left(\left|{\beta }_j\right|\right)\right]\right\},\hfill \end{array}$$

(44)

where $p_{\lambda }(·)$ is a penalty function, and the bridge penalty is adopted in [93]. Based on (44), Nghiem et al. [93] proposed PMSc and SISc methods. For the PMSc method, it suggested taking the selected submodel as

$${\widehat{\mathcal{S}}}_{\mathrm{PMSc}}=\left\{j:{\tilde{\beta }}_{\mathrm{sc},j}\ne 0\right\}.$$

Under some regularity conditions, Nghiem et al. [93] showed that $\mathbb{P}(\mathcal{S}\subset {\widehat{\mathcal{S}}}_{\mathrm{PMSc}})\to 1$. Furthermore, when $\lambda =0$, we can obtain that

$${\tilde{\beta }}_{\mathrm{sc},j}=\frac{{\sum }_{i=1}^n{w}_{ij}{y}_i}{{\sum }_{i=1}^n{w}_{ij}^2-n{\sigma }_{u,j}^2},j=1,\dots ,p,$$

which measures the marginal correlation between the jth variable and the response variable. The SISc selects the variable according to the magnitude of ${\tilde{\beta }}_{\mathrm{sc},j}$. The corresponding selected set is

$${\widehat{\mathcal{S}}}_{\mathrm{SISc}}=\left\{1\le j\le p:|{\tilde{\beta }}_j|\mathrm{is}\mathrm{among}\mathrm{the}d\mathrm{largest}\mathrm{of}\mathrm{all}\right\}.$$

Nghiem et al. [93] proved that $\mathbb{P}(\mathcal{S}\subset {\widehat{\mathcal{S}}}_{\mathrm{SISc}})=1-O\{pexp(-Cn)\}$ for some constant $C>0$ under some regularity conditions.

6. Conclusions

With the advent of the big data era, high-dimensional measurement error data have proliferated in various fields. Over the past few years, many statistical inference methods for high-dimensional measurement error regression models have been developed to overcome the difficulties in scientific research and provide effective approaches for tackling problems in applications. This paper reviews the research advances in estimation and hypothesis testing methods for high-dimensional measurement error models as well as variable screening methods for ultrahigh-dimensional measurement error models. The aforementioned estimation methods can be classified into the following three categories: (i) methods based on a nonconvex objective function with restrictions on the regression coefficients, such as the nonconvex Lasso and the estimation method for Poisson models in [78]; (ii) methods with a convex objective function including CoCoLasso, the balanced estimation method and the CaZnRLS estimation method; (iii) methods that draw on the idea of a Dantzig selector, such as the MU selector, compensated MU selector, GMU selector, and conic-programming-based estimation method. Many methods are now available in R packages “hdme” and “BDcocolasso”. Thus, we can apply these methods to analyze high-dimensional measurement error data. For the use of estimation methods, it is recommended to use CoCoLasso and balanced estimation methods due to their operability. If a higher computational efficiency is required, the CaZnRLS estimation method can be considered. If covariates are only partially corrupted by measurement errors, it is better to apply BDCoCoLasso.

Due to the prevalence of high-dimensional measurement error data in daily life and the growing demand for the statistical inference methods of measurement error regression models in applications, the related research is still one of the crucial aspects in statistical research. At present, the statistical inference methods and the theoretical system of high-dimensional measurement error models are far from complete. To the best of our knowledge, the study of high-dimensional measurement error regression models is currently limited to linear models and generalized linear models. However, it is common that covariates and response variables show a complicated relationship rather than a simple linear relationship in practice. Therefore, in order to meet the urgent needs of applications, it is necessary to develop more general statistical inference methods for high-dimensional nonlinear measurement error models. Further research in this area includes the following aspects.

• Existing estimation methods for high-dimensional measurement error regression models are mainly for linear or generalized linear models. Therefore, it is urgent to develop estimation methods for nonlinear models with high-dimensional measurement error data such as nonparametric and semiparametric models.
• Existing works mainly focus on independent and identically distributed data. It is worthwhile to extend the estimation and hypothesis-testing methods to measurement error models with complex data such as panel data and functional data.
• In most studies of high-dimensional measurement error models, it is assumed that the covariance structure of the measurement errors is specific or the covariance matrix of measurement errors is known. Thus, it is a challenging problem to develop estimation and hypothesis-testing methods in the case that the covariance matrix of measurement errors is completely unknown.