Robust estimation for varying-coefficient partially linear measurement error model with auxiliary instrumental variables

1.   Introduction

Varying-coefficient partially linear model is one of the most popular semiparametric regression models. It takes the form as

where $Y$ is the response variable, $\mathit{\boldsymbol{X}}$, $\mathit{\boldsymbol{Z}}$ and $U$ are the associated covariates, ${\bf{\Theta}}$ is a $p$-dimensional vector of unknown parameters, $\boldsymbol{\alpha}(.)$ is a $q$-dimensional vector of nonparametric function, $\varepsilon$ is model error with $E(\varepsilon|\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U) = 0$ and $\text{Cov}(\varepsilon|\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U) = \sigma^2$. By combining the flexibility of nonparametric model and the adaptability of parametric model, model (1.1) is quite general which covers many important models as special cases, such as liner model, varying-coefficient model, partially linear model and others. However, due to the low accuracy of instrument equipment, the imperfection of measurement technique, and the high price or time cost, it is often impossible to observe accurate data in the experiments, but the data with errors-in-variables(EV) ^[1]. If one simply ignores the measurement error, known as naive method, it will result in unreliable statistical inference, such as biased estimators, loss of efficiency, decreasing power of hypothesis testing. According to the diversity of errors in variables, the treatment of measurement errors differs for various types.

Early research focused on the study of simple additive measurement errors in the linear component. Instead of true variable $\mathit{\boldsymbol{X}}$, we can observe $\mathit{\boldsymbol{\eta}}$ with $\mathit{\boldsymbol{\eta}} = {X\mathit{\boldsymbol{}}}+{e\mathit{\boldsymbol{}}}$ and ${e\mathit{\boldsymbol{}}}$ is the measurement error, which is independent of $(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, \varepsilon)$. You and Chen ^[2] proposed a modified profile least squares procedure and Wang et al. ^[3] applied empirical likelihood inference to study model (1.1) with measurement errors in the linear part. On the other hand, measurement error in the nonparametric part in model (1.1) may occur in reality. The surrogate variable $\mathit{\boldsymbol{\eta}}$ is observed by $\mathit{\boldsymbol{\eta}} = {Z\mathit{\boldsymbol{}}}+{e\mathit{\boldsymbol{}}}$ with the measurement error ${e\mathit{\boldsymbol{}}}$. Feng and Xue ^[4] proposed a local bias-corrected restricted profile least squares estimators and Fan et al. ^[5] applied some auxiliary information to construct empirical log-likelihood ratios for both parameter and nonparametric functions. It is noted that in the above literature, we assumed that the relationship between the unobservable variable and the observed EV data is simple additive. Furthermore, it is necessary that in the statistical inference, the variance of measurement error is known before. Even if the variance of error is unavailable in practice, we can use the estimators to replace it by repeatedly measuring actually observation data $\mathit{\boldsymbol{\eta}}$.

However, the additive measurement error structure may not describe the complexity of experimental data and sometimes there are no repeated observations to estimate the variance of measurement error. In this paper, we do not specify any model structure of the measurement error, and do not require the information of the variance of measurement error. Furthermore, repeated measurement data are not necessary. We assume that some variables are unobserved directly, but auxiliary information is available to remit the true variable. Specifically, the unobserved variable and the observed surrogate variable with measurement error is connected by a nonparametric structure with an instrument variable. Let ${X\mathit{\boldsymbol{}}} = (\mathit{\boldsymbol{\xi}}^T, \mathit{\boldsymbol{W}}^T)^T$, where $\mathit{\boldsymbol{\xi}}$ is a $p_1\times 1$ unobserved directly vector and $\mathit{\boldsymbol{W}}$ is a vector of the remaining exactly observed components. We assume that the true variable $\mathit{\boldsymbol{\xi}}$ is related to the observed measurement variable $\mathit{\boldsymbol{\eta}}$ and auxiliary variable $V$ through the nonparametric structure

From assumption (1.2), model (1.1) with error-prone linear covariates implies that

with ${\bf{\Theta}} = (\mathit{\boldsymbol{\beta}}^T, \mathit{\boldsymbol{\theta}}^T)^T$, $\mathit{\boldsymbol{X}} = (\mathit{\boldsymbol{\xi}}^T, \mathit{\boldsymbol{W}}^T)^T$ and $\mathit{\boldsymbol{e}}$ is the measurement error with $E(\mathit{\boldsymbol{e}}|V) = \mathit{\boldsymbol{0}}$ and the covariance matrix $\Sigma_{e}$. Therefore, the error model structures of $\mathit{\boldsymbol{\xi}}$ is a special cases of the additive error models. For the sake of simplicity, we assume that the auxiliary variables $V$ is scalar. This structure of measurement error is more general, introduced by Zhou and Liang ^[6], which includes denoise linear model and rational expectation model and others. Zhou and Liang ^[6] proposed a regression calibration technique to reduce the bias due to mismeasurement and then developed the profile least-square based estimation procedure for parametric and nonparametric components.

Meanwhile, this structure of error-prone covariates has been extensively studied by researchers for other semiparametric models. For the varying-coefficient models with error-prone covariates, Zhao and Xue ^[7] constructed confidence interval for the varying-coefficient functions using instrumental variable-based empirical likelihood inference and Xu et al. ^[8] presented an instrumental variable type local polynomial estimator. Zhang et al. ^[9] proposed the estimation and variable selection for partial linear single index model when some linear covariates are not observed, but their ancillary variables are available. Huang and Ding ^[10] studied a new partially linear error-in-variable models with error-prone covariates in the parametric part and mismeasured covariates in the nonparametric part, simultaneously. Sun et al. ^[11] developed an estimation procedure for semiparametric varying-coefficient partially linear model when both response and part of covariates are measured with error by functions of some auxiliary variables.

In this paper, we propose two step estimation procedure for model (1.3). In the first step, we give a bias correction for the error-prone covariate by using ancillary information and applying locally linear technique. Then the covariates are corrected based on instrumental variables. In the second stage, instead of profile least squares estimation, we propose a robust estimation procedure. It is remarkable that the least squares technique may be sensitive to outliers and inefficient for many non-normal error distributions. Modal regression, introduced by Yao et al. ^[12], is a powerful tool for revealing the underlying relationship between the response and corresponding covariates in the presence of outlier or non-normal error distribution, and it has been successfully applied to various semiparametric models (Zhao et al. ^[13], Yang and Yang ^[14] and Yang et al. ^[15], Lv et al. ^[16], Yu et al. ^[17] and others). This technique has some good properties, such as easy to implement, robust to outliers, full asymptotic efficiency under the non-normal error distribution. Motivated by these, we extend modal regression in the second stage to estimate the parameter and nonparametric function.

The rest of this paper is organized as follows. Section 2 introduces a two-stage robust instrumental variables-based modal regression estimator and Section 3 establishes its theoretical properties. In section 4, we discuss the selection of bandwidth and the specific estimation algorithm. Simulation studies are conducted to evaluate the performances of the proposed estimation procedure in Section 5. A real data analysis is illustrated in section 6 to show the effective of the proposed estimation procedures. We make our concluding remarks in Section 7 and leave the proofs of the main Theorems in Section 8.

2.   Robust estimation procedure

In this section, we propose a two-stage robust estimation procedure for model (1.3). Firstly, the local linear smoothing is adapted to calibrate the variable $\mathit{\boldsymbol{\xi}}$ by using auxiliary instrument variables $\mathit{\boldsymbol{\eta}}$ and $V$. Secondly, based on B-spline basis functions and modal regression, the final estimators of the parametric and nonparametric components are obtained.

2.1.   Covariate calibration

Suppose that $\{(Y_{i}, \mathit{\boldsymbol{\xi}}_{i}, \mathit{\boldsymbol{W}}_{i}, \mathit{\boldsymbol{Z}}_{i}, U_{i}, \mathit{\boldsymbol{\eta}}_{i}, V_{i})\}_{i = 1}^n$ is the independent and identically distributed sample from $({Y}, \mathit{\boldsymbol{\xi}}, \mathit{\boldsymbol{W}}, \mathit{\boldsymbol{Z}}, {U}, \mathit{\boldsymbol{\eta}}, {V})$. The varying-coefficient partially linear model with error-prone in linear covariates has the following form,

Firstly, we need to calibrate the covariate $\mathit{\boldsymbol{\xi}}_{i}$, which is not observed in model (2.1). Let $\eta_{ik}$ be the $k$th entry of vector $\mathit{\boldsymbol{\eta}}_{i}, i = 1, \ldots, n$. To estimate $\xi_{k}(v)$, the $k$th component of $\mathit{\boldsymbol{\xi}}(v)$, we employ the local linear smoothing technique ^[18] to minimize

where $L_b(.) = L(./b)/b$ is the kernel function, and $b = b_k$ for $k = 1, \ldots, p_1$ is the bandwidth for the $k$th component of $\mathit{\boldsymbol{\xi}}$. The local linear estimator of ${\xi}_k(v)$, $k = 1, \ldots, p_1$, denoted by $\hat{{\xi}}_k(v)$,

where $\mathit{\boldsymbol{\eta}}^{(k)} = (\eta_{1k}, \eta_{2k}, \ldots, \eta_{nk})^T$, $W_{v} = \mathrm{diag}(L_{b_k}(V_{1}-v), L_{b_k}(V_{2}-v), \ldots, L_{b_k}(V_{n}-v))$, and $V_{v}$ is defined in the following,

Therefore, unobserved covariate $\mathit{\boldsymbol{\xi}}_{i}$ can be replaced by its local linear estimator $\hat{\mathit{\boldsymbol{\xi}}}(V_{i})$, abbreviated as $\hat{\mathit{\boldsymbol{\xi}}}_i$. Then the robust estimators of parameters and nonparametric function can be constructed in the following procedure.

2.2.   Modal regression estimator

Following Schumaker ^[19], we can replace the nonparametric functions ${\alpha}_k(u), k = 1, 2, \ldots, q$ through their basis function approximations. Let $\mathit{\boldsymbol{B}}(u) = (B_1(u), \ldots, B_L(u))^T$ be the B-spline basis function with the order of $M$, where $L = K+M$ and $K$ is the number of interior knots. Then, the nonparametric function ${\alpha}_{k}(u)$ can be approximated by

where $\mathit{\boldsymbol{\gamma}}_{k} = ({\gamma_{k1}}, {\gamma_{k2}}, \ldots, {\gamma_{kL}})^T$. Then model (2.4) is substituted into the model (2.1), we can obtain

where $\mathit{\boldsymbol{\Pi}}_{i} = I_{q}\otimes\mathit{\boldsymbol{B}}(U_{i})\cdot \mathit{\boldsymbol{Z}}_{i}$ and $\mathit{\boldsymbol{\gamma}} = (\gamma_{1}^T, \gamma_{2}^T, \ldots, \gamma_{q}^T)^T$. Model (2.5) is a standard linear regression model, and it contains three parametric vectors $\mathit{\boldsymbol{\beta}}$, $\mathit{\boldsymbol{\theta}}$ and $\mathit{\boldsymbol{\gamma}}$. Each function $\alpha_k(.)$ in model (2.1) is characterized by $\mathit{\boldsymbol{\gamma}}_k$ in model (2.4).

Based on the thought of modal regression ^[12], the estimators of the parameters $\mathit{\boldsymbol{\beta}}$, $\mathit{\boldsymbol{\theta}}$ and $\mathit{\boldsymbol{\gamma}}$ are obtained through maximizing the following objective function,

where $\phi_{h}(.) = \phi({.}/{h})/h$, $\phi(.)$ is the kernel function and $h$ is a bandwidth. The value of bandwidth $h$ determines the level of the robustness of the estimator, then the specific selection of optimal bandwidth $h$ is described in the next section.

However, we can not maximize (2.6) directly due to the unobserved variable $\mathit{\boldsymbol{\xi}}_{i}$. Instead, we substitute the local linear estimator $\hat{\mathit{\boldsymbol{\xi}}}_i$ with the true variable $\mathit{\boldsymbol{\xi}}_{i}$. Further, the calibrated expression is rewrite as

Therefore, we get the spline modal regression estimators $\hat{\mathit{\boldsymbol{\beta}}}$, $\hat{\mathit{\boldsymbol{\theta}}}$ and $\hat{\mathit{\boldsymbol{\gamma}}}$ by maximizing the objective function (2.7). Then, the corresponding estimator of coefficient function ${\alpha}_k(u)$ is obtained by $\hat{{\alpha}}_{k}(u) = \mathit{\boldsymbol{B}}(u)^T \hat{\mathit{\boldsymbol{\gamma}}_{k}}, k = 1, 2, \ldots, q$.

3.   Theoretical properties

In this section, we discuss the asymptotic properties of the resulting estimators. Denote ${\bf{\Theta}}_0 = (\mathit{\boldsymbol{\beta}}_0^T, \mathit{\boldsymbol{\theta}}_0^T)^T$ and $\mathit{\boldsymbol{\alpha}}_0(.)$ to be the true value of ${\bf{\Theta}} = (\mathit{\boldsymbol{\beta}}^T, \mathit{\boldsymbol{\theta}}^T)^T$ and $\mathit{\boldsymbol{\alpha}}(.)$ in model (1.1). Correspondingly, $\mathit{\boldsymbol{\gamma}}_{0j}$ is the best approximation coefficient of $\alpha_{0j}(u)$ in the B-spline space. Let ${F(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u, h)} = \mathrm{E}(\phi_h''(\varepsilon)|{\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{x}}, \mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{z}}}, {U} = {u})$ and ${G(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u, h)} = \mathrm{E}(\phi_h'(\varepsilon)^2|{\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{x}}, \mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{z}}}, {U} = {u})$.

We begin with the following assumption conditions required to derive the main results. These conditions are quite mild and can be easily satisfied.

A1: The random variable $U$ has a bounded support ${\mathcal{U}}$. Its density function $f_u(.)$ is Lipschitz continuous and bounded away from 0 and infinite on its support. The density function of random function $V$, $f_v(.)$ is continuously differentiable and bounded away from 0 and infinite on its support $\mathcal{V}$.

A2: The nonparametric function $\alpha_k(u), k = 1, \ldots, q$ and $\xi_k(v), k = 1, \ldots, p_1$ have a continuous second derivative.

A3: The matrix $E(\mathit{\boldsymbol{XX}}^T|U)$ is nonsingular for each $U\in \mathcal{U}$. All elements of the matrices $E(\mathit{\boldsymbol{XX}}^T|U)$, $E(\mathit{\boldsymbol{XX}}^T|U)^{-1}$ and $E(\mathit{\boldsymbol{XZ}}^T|U)^{-1}$ are Lipschitz continuous.

A4: The kernel functions $L(.)$ and $\phi(.)$ are density functions with compact supports.

A5: There exists an $s > 0$ such that $\mathrm{E}||\textbf{X}||^{2s} < \infty$ and for some $r < 2-s^{-1}$ such that $n^{2r-1}h\rightarrow \infty$, $n^{2r-1}b_k\rightarrow \infty$ and $nhb_k^4\rightarrow 0$ as $n\rightarrow \infty$, $k = 1, 2, \ldots, p_1$.

A6: $nh^8\rightarrow 0$ and $nh^2/(\mathrm{log} n)^2\rightarrow \infty$.

A7: Let $c_1, \ldots, c_K$ be the interior knots on $[0, 1]$. Set $c_0 = 0, \; c_{K+1} = 1$ and $h_i = c_i-c_{i-1}$. Then, there exists a constant $C_0$ such that

A8: $F(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u, h)$ and $G(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u, h)$ are continuous with respect to $(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u)$. In addition, $F(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u, h) < 0$ for any $h > 0$.

A9: $\mathrm{E}(\phi_h'(\varepsilon)|\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u) = 0$ and $\mathrm{E}(\phi_h''(\varepsilon)^2|\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u)$, $\mathrm{E}(\phi_h'(\varepsilon)^3|\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u)$ and $\mathrm{E}(\phi_h'''(\varepsilon)|\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{z}}, u)$ are continuous with respect to $\mathit{\boldsymbol{x}}$, $\mathit{\boldsymbol{z}}$ and $u$.

The following Theorem 1 gives the consistency of spline modal estimators $\mathit{\boldsymbol{\hat\Theta}} = (\mathit{\boldsymbol{\hat\beta}}^T, \mathit{\boldsymbol{\hat\theta}}^T)^T$ and $\mathit{\boldsymbol{\hat\alpha}}_{k}(u), k = 1, \ldots, q$.

Theorem 1. Suppose that the assumptions A1-A9 hold, and the number of knots $K = O(n^{1/(2r+1)})$. Then we have

The following Theorem 2 shows the asymptotic normality for the parameter estimator $\mathit{\boldsymbol{\hat\Theta}}$.

Let $\Phi(U) = \mathrm{E}(F({X\mathit{\boldsymbol{}}}, {Z\mathit{\boldsymbol{}}}, U, h){\bf \Pi}{\bf X}^T|U)$, and $\Gamma(U) = \mathrm{E}(F({X\mathit{\boldsymbol{}}}, {Z\mathit{\boldsymbol{}}}, U, h){\bf \Pi}{\bf \Pi}^T|U)$. Denote ${\bf \tilde{X}}({Z\mathit{\boldsymbol{}}}, {X\mathit{\boldsymbol{}}}, U) = {\bf X}-\Phi(U)^T\Gamma^{-1}(U){\bf \Pi}$.

Theorem 2. Suppose that the assumptions A1-A9 hold, and the number of knots $K = O(n^{1/(2r+1)})$, then we have

where $\mathit{\boldsymbol{\Sigma}} = \mathrm{E}(F(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h){\bf \tilde{X}}{\bf \tilde{X}}^T)$ and $\mathit{\boldsymbol{\Omega}} = \mathrm{E}(G({X\mathit{\boldsymbol{}}}, {Z\mathit{\boldsymbol{}}}, U, h){\bf \tilde{X}}{\bf \tilde{X}}^T)+\mathit{\boldsymbol{\beta}}_0^T\Sigma_e\mathit{\boldsymbol{\beta}}_0\mathrm{E}(F(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h)^2\mathrm{E}(\tilde{\mathit{\boldsymbol{X}}}\tilde{\mathit{\boldsymbol{X}}}^T))$.

4.   Bandwidth selection and estimation algorithm

In this section, we discuss the selection of bandwidth and estimation procedure based on the modal expectation maximization(MEM) algorithm proposed by Li et al. ^[20].

4.1.   Bandwidth selection in practice

In order to obtain the spline modal estimators $\hat{\mathit{\boldsymbol{\beta}}}$, $\hat{\mathit{\boldsymbol{\theta}}}$ and $\hat{\mathit{\boldsymbol{\gamma}}}$ by maximizing function (2.7), we need to select an appropriate bandwidth $h$ due to the fact that its value affects directly the robustness of the estimators. According to the thought of Yao et al. ^[12], it is easy to show that the ratio of the asymptotic variance of our proposed estimator to that of the least square B-spline estimator is given by

where $\hat F(h) = \frac{1}{n}\sum_{i = 1}^n \phi''_h(\hat\varepsilon_i)$, $\hat G(h) = \frac{1}{n}\sum_{i = 1}^n \phi'_h(\hat\varepsilon_i)^2$, $\hat\varepsilon_i = Y_{i}-\hat{\mathit{\boldsymbol{\xi}}}_i^{T} \mathit{\boldsymbol{\beta}}^{(0)}-\mathit{\boldsymbol{W}}_{i}^T \mathit{\boldsymbol{\theta}}^{(0)}-\mathit{\boldsymbol{\Pi}}_{i}^T {\mathit{\boldsymbol{\gamma}}}^{(0)}$, the parameters $\mathit{\boldsymbol{\beta}}^{(0)}$, $\mathit{\boldsymbol{\theta}}^{(0)}$ and $\mathit{\boldsymbol{\gamma}}^{(0)}$ are the initial estimators of the parametric vectors $\mathit{\boldsymbol{\beta}}$, $\mathit{\boldsymbol{\theta}}$ and $\mathit{\boldsymbol{\gamma}}$ via the least square estimation procedure, and $\hat\sigma^2 = \frac{1}{n}\sum_{i = 1}^n \hat\varepsilon_i^2.$

The optimal bandwidth $h$ is calculated by minimizing (4.1). In the actual calculation, we use the grid search method to choose $h$, and the possible grids points for $h$ can be $h = 0.5\hat\sigma\times 1.02^j, j = 0, 1, \ldots, 100$ (See Yao et al. ^[12] and Lv et al. ^[16]).

4.2.   The MEM algorithm for parameters

In this subsection, the following modal expectation maximization(MEM) algorithm by Li et al. ^[20] is adopted to obtain the estimators of the parameters ${\bf{\Theta}}$ and $\mathit{\boldsymbol{\gamma}}$. Let $\mathit{\boldsymbol{Q}}_{i} = (\hat{\mathit{\boldsymbol{\xi}}}_i^{T}, \mathit{\boldsymbol{W}}_{i}^T, \mathit{\boldsymbol{\Pi}}_{i}^T)^{T}$, $\mathit{\boldsymbol{\Xi}} = ({\bf{\Theta}}^{T}, \mathit{\boldsymbol{\gamma}}^{T})^T$. Set $\mathit{\boldsymbol{\Xi}}^{(0)} = ({\bf{\Theta}}^{(0)T}, \mathit{\boldsymbol{\gamma}}^{(0)T})^{T}$, where ${\bf{\Theta}}^{(0)}$ and $\mathit{\boldsymbol{\gamma}}^{(0)}$ are initial estimators of ${\bf{\Theta}}$ and $\mathit{\boldsymbol{\gamma}}$, respectively with $m = 0$.

${\bf Step 1.(E-step)}$: We update $\pi(i\left|\mathit{\boldsymbol{\Xi}}^{(m)}\right)$ by

${\bf Step 2.(M-step)}$: Update $\mathit{\boldsymbol{\Xi}}^{(m+1)}$ by

where ${\bf Q} = (\mathit{\boldsymbol{Q}}_{1}, \mathit{\boldsymbol{Q}}_{2}, \ldots, \mathit{\boldsymbol{Q}}_{n})^{T}$, $\mathit{\boldsymbol{Y}} = (Y_{1}, Y_{2}, \ldots, Y_{n})^{T}$, $\mathit{\boldsymbol{\omega}}$ is $n\times n$ diagonal matrix with $\mathit{\boldsymbol{\omega}} = diag(\pi(1|\mathit{\boldsymbol{\Xi}}^{(m)}), \pi(2|\mathit{\boldsymbol{\Xi}}^{(m)}), \ldots, \pi(n|\mathit{\boldsymbol{\Xi}}^{(m)}))$.

${\bf Step 3}$: Iterate the E-step and M-step until the algorithm converges. We denote the final estimator of $\mathit{\boldsymbol{\Xi}}$ as $\hat{\mathit{\boldsymbol{\Xi}}}$. Then, the parameters for ${\bf{\Theta}}$ and $\mathit{\boldsymbol{\gamma}}$ based on modal regression are $\hat{{\bf{\Theta}}} = ({\bf I}_{p\times p}, {\bf 0}_{p\times qL})\hat{\mathit{\boldsymbol{\Xi}}}$ and $\hat{\mathit{\boldsymbol{\gamma}}} = ({\bf 0}_{qL\times p}, {\bf I}_{qL\times qL})\hat{\mathit{\boldsymbol{\Xi}}}$. Then, the corresponding estimator of $\alpha_k(.)$ is $\hat{{\alpha}}_{k}(u) = \mathit{\boldsymbol{B}}(u)^T \hat{\mathit{\boldsymbol{\gamma}}_{k}}, k = 1, 2, \ldots, q$.

5.   Simulation studies

In this section, we conduct the following simulation to assess the finite sample performance of the proposed estimation procedure. For this purpose, we consider the following varying-coefficient partially linear model with error-prone linear covariates,

where the parametric vector $\mathit{\boldsymbol{\beta}} = (\beta_{1}, \beta_{2}, \beta_{3})^{T} = (3, 1.8, 0.8)^{T}$, $\mathit{\boldsymbol{\theta}} = (\theta_{1}, \theta_{2}, \theta_{3})^{T} = (2, 0.5, -1.5)^{T}$, the nonparametric function $\mathit{\boldsymbol{\alpha}}(U) = (\alpha_{1}(U), \alpha_{2}(U))^{T}$ with $\alpha_{1}(U) = 3.5(\mathrm{exp}(-(4U-1)^2)+\mathrm{exp}(-(4U-3)^2))-1.5$, $\alpha_{2}(U) = \mathrm{cos}(2\pi U)$ and $U_{i}\sim U(0, 1)$, the covariates $\mathit{\boldsymbol{W}}_{i}$ are multivariate normal distribution with mean ${\bf 0}$ and covariance matrix $\Sigma_{1} = (\sigma_{ij})_{1\leq i, j\leq3}$ with $\sigma_{ij} = 0.5^{|i-j|}$ for $i, j = 1, 2, 3$ and the covariates $\mathit{\boldsymbol{Z}}_{i}$ are multivariate normal distribution with mean ${\bf 0}$ and covariance matrix $\Sigma_{2} = (\sigma_{ij})_{1\leq i, j\leq2}$ with $\sigma_{ij} = 0.5^{|i-j|}$ for $i, j = 1, 2$. The following three different model errors: (1) the normal distribution: $\varepsilon_{i}\sim N(0, 1)$; (2) the $t$ distribution: $\varepsilon_{i}\sim t(3)$; (3) the mixed normal distribution(MN): $\varepsilon_{i}\sim 0.9N(0, 1)+0.1N(0, 9^{2})$ are considered.

The variables $\mathit{\boldsymbol{\xi}}_{i}$ in model (5.1) are the function of the auxiliary instrumental variables $V_{i}$. It is noted that $\mathit{\boldsymbol{\xi}}_{i}$ are not observed, and the observed variables are $\mathit{\boldsymbol{\eta}}_{i}$ with $\mathit{\boldsymbol{\eta}}_{i} = \mathit{\boldsymbol{\xi}}(V_{i})+\mathit{\boldsymbol{e}}_{i}$, and $\xi_{1}(V_{i}) = V_{i}(1-V_{i})$, $\xi_{2}(V_{i}) = 2\mathrm{sin}(V_{i})$, $\xi_{3}(V_{i}) = 3V_{i}-2\mathrm{cos}(V_{i})$ with $V_{i}\sim N(2, 1.5)$. Furthermore, the measurement error $\mathit{\boldsymbol{e}}_{i}$ is independent with $V_{i}$ and $(\mathit{\boldsymbol{W}}_{i}, \mathit{\boldsymbol{Z}}_{i}, \varepsilon_{i})$ and $\mathit{\boldsymbol{e}}_{i}\sim N(\textbf{0}, I_3\sigma_{e}^{2})$, where $\sigma_{e}^2 = 0.4^2$ and $\sigma_{e}^2 = 0.6^2$ represent the different level of the measurement error.

In the simulation experiments, we choose the Epanechnikov kernel function $L(t) = \frac{3}{4}(1-t^{2})_{+}$, the bandwidth $b = \hat{\sigma}_{v} n^{-1/3}$ is adopted in local linear smoothing according to Zhou and Liang ^[6], where $\hat{\sigma}_{v}$ is sample root variance of the variable $V_{i}$. We use the standard normal density for $\phi(t)$ in the modal regression and the bandwidth $h$ is selected through minimizing (4.1). We use the cubic B-spline basis function in our simulation. The sample size $n$ is 100,250 and 500, and run simulation is 200.

To evaluate the performance of the proposed estimation procedure in the parametric component, we compared the modal regression (MR) estimators to the profile least squares (PLS) estimators proposed by Zhou and Liang ^[6]. Three types of estimators include the naive estimators in which $\mathit{\boldsymbol{\xi}}_{i}$ is simply replaced by $\mathit{\boldsymbol{\eta}}_{i}$ (denoted by MR-N and PLS-N), our proposed calibrated estimators (denoted by MR-C and PLS-C), and the benchmark estimators which assumes that $\mathit{\boldsymbol{\xi}}_{i}$ is known before (denoted by MR-B and PLS-B). Tables 1–6 show the bias and standard deviation (denoted by SD in the brackets) of the estimators of parameters $\mathit{\boldsymbol{\beta}}$ and $\mathit{\boldsymbol{\theta}}$.

According to Tables 1–6, we draw the following conclusions: (1) The estimators of parameters are obtained under six estimation approaches. Their bias and SD become smaller with the increasing of the sample size in all situations. This indicates that all estimation procedures are reasonable. (2) The bias and SD of the PLS estimator are smaller than MR estimator while the model error yielding to normal distribution. However, when the model error is the $t$ or MN distribution, the proposed MR estimator gives smaller bias and SD than the PLS estimator. It illustrates that the PLS estimator is sensitive to the non normal distributed model error. (3) Compared the estimators among MR-N, MR-C and MR-B, we find that our proposed MR-C estimator is superior to MR-N estimator and it performs just a little bit worse than MR-B. And the PLS-C estimator has the similar performance to the MR-C estimator, which demonstrates that it is necessary to calibrate measurement errors in the calculation and we should not ignore the effect of the measurement errors. (4) Given the sample size, with the increase of the measurement error, the SD of PLS-C and MR-C become smaller. This reveals that the calibration of the error-prone covaritates is effective.

Furthermore, to evaluate the estimation results on the nonparametric functions, we present the fitting curve of $\alpha_1(u)$ and $\alpha_2(u)$. Figures 1 and 2 show the plot of estimated curve with the dotted, the dashed and the solid lines which denote the PLS-C estimator, MR-C estimator and the real function, respectively. We only give the results of the measurement error $\sigma_{e}^{2} = 0.4^2$ and $\sigma_{e}^{2} = 0.6^2$ under $n = 250$ with model error $\varepsilon_{i}\sim t(3)$. From Figures 1 and 2, we can see that two estimation procedures of PLS-C and MR-C estimator are all closed to the real curve, and the MR-C is much closer to the real curve. This indicates that our proposed estimation procedure is more robust when the data follows non normal distributions.

In addition, we apply the root of mean square error(RMSE) to illustrate the effectiveness of the nonparametric estimators, and the form of the RMSE is

where $M = 200$, $u_{i}$ is the equidistant points on $(0, 1)$. The boxplots of the RMSE for the nonparametric functions are given in Figures 3 and 4 by PLS-C and MR-C estimators when the measurement error is $\sigma_{e}^{2} = 0.4^2$ and $\sigma_{e}^{2} = 0.6^2$, respectively, and the sample size $n = 250$ and $n = 500$ under the model error $\varepsilon_{i}\sim t(3)$. From Figures 3 and 4, the MR-C estimator performs smaller RMSE than PLS-C estimator. It further confirms that our proposed estimation procedure is robust when the model error is non normal. Meanwhile, the RMSEs of two procedures are decreasing as the sample size increases.

6.   A real example

We adopt a diabetes data to our proposed estimation procedure in this section. The data set contains 10 baseline variables, including age, sex, body mass index, average blood pressure, six blood serum measurements obtained from 442 observations, and the response of interest, a quantitative measurement of disease progression 1 year after baseline. A detailed description of each variable can be seen by Efron et al. ^[21]. Similar to Sun et al. ^[11], the following varying-coefficient partially linear measurement error model is given as

Low-density level (ldl) and high-density level(hdl) are the covariates $\mathit{\boldsymbol{\xi}}$ which are measured with error. The exactly observed covariates $\textbf{W}$ are triglycerides(tc), total cholesterol(tch) and tension glaucoma level(ltg). The covariates $\textbf{Z}$ are the glucose concentration(glu) and average blood pressure(bp) and the variable $U$ is age. Further, the variable body mass index(bmi) is considered as the auxiliary instrumental variable $V$.

We presented two methods with PLS-C and IMR-C to estimate parameter and nonparametric functions in model (6.1). The parametric estimating results are given in Table 7. From Table 7, the calibrated methods based on auxiliary instrumental variables are available in practice and the parametric estimators with PLS-C and MR-C have similar effects on response variable. In addition, the fitting curves of the coefficient functions are given in Figure 5, where the solid line and the dashed line describe the PLS-C method and the MR-C method, respectively. From Figure 5, we can see that the fitting trend of the nonparametric functions with PLS-C is close to MR-C. In a word, our proposed estimation procedures are effective in both parametric and nonparametric aspects.

7.   Conclusions

We have presented a robust estimation procedure for varying-coefficient partially linear model when some covariates in parametric part are error-prone. We assume that some linear covariates in the model can not be observed directly, and do not specify the model structure of measurement error. Firstly, we calibrate the covariates with the auxiliary instrumental variable by the local linear smoothing technique. Then with the B-spline basis function and modal regression, the components of parameters and coefficient functions are simultaneously estimated. The proposed estimation procedure not only can attenuate the effect of the measurement errors, but also is robust against outliers or heavy-tail distributed data. Under some mild conditions, the consistency and the asymptotic normality of the resulting estimators are established. We also give the strategy to choose optimal bandwidth and specific algorithm in practice. Some simulation studies and a real data analysis are carried out to demonstrate the proposed estimation methods are satisfactory.

Acknowledgments

This work is supported by the Natural Science Basic Research Plan in Shaanxi Province of China (2022JM-027)

Conflict of interest

There are no conflicts of interest as declared by the authors.

Appendix: Proofs of theorems

Lemma 1. Under the Assumptions of the theoretical properties in Section 3, we have

uniformly on $v\in V$, where $V$ is the compact support of $v$ and ${e}_{ik}$ is the $k$th entry of $e_i$. $\mu_{L_2} = \int v^2 L(v) {\mathit{\text{d}}}v$ and $f_V (v)$ is the density function of $V$.

Proof. This result is directly followed by Zhou and Liang ^[6]. □

Proof of Theorem 1. Let $\delta = n^{-r/(2r+1)}$. Define ${\bf{\Theta}} = {\bf{\Theta}}_0+\delta {\bf v}_1$ and $\mathit{\boldsymbol{\gamma}} = \mathit{\boldsymbol{\gamma}}_0+\delta {\bf v}_2$, and ${\bf v} = ({\bf v}_1^T, {\bf v}_2^T)^T$, where ${\bf{\Theta}}_0 = (\mathit{\boldsymbol{\beta}}_0^T, \mathit{\boldsymbol{\theta}}_0^T)^T$ and $\mathit{\boldsymbol{\gamma}}_0$ are the true values of ${\bf{\Theta}}$ and $\mathit{\boldsymbol{\gamma}}$, ${\bf v}_1$ and ${\bf v}_2$ are vector with $p$ and $qL$ dimension, respectively.

We show that, for any given $\varepsilon > 0$, there exists a large enough constant $C$ such that

where $\hat L({\bf{\Theta}}, \mathit{\boldsymbol{\gamma}})$ is defined in (2.7).

Let $R_k(u) = \alpha_k(u)-\mathit{\boldsymbol{B}}(u)^T\mathit{\boldsymbol{\gamma}}_{0k}$, ${R\mathit{\boldsymbol{}}}(u) = (R_1(u), R_2(u), \ldots, R_q(u))^T$. Adopting Taylor expansion to (2.7), we have

where $\hat{{\bf X}}_i = (\hat{\mathit{\boldsymbol{\xi}}}_i^T, \mathit{\boldsymbol{W}}_i^T)^T$, $\zeta_{i}$ lies in $\varepsilon_i+\mathit{\boldsymbol{Z}}_i^T \mathit{\boldsymbol{R}}(U_i)+(\mathit{\boldsymbol{\xi}}_i-\hat {\mathit{\boldsymbol{\xi}}}_i)^T\mathit{\boldsymbol{\beta}}_0$ and $\varepsilon_i+\mathit{\boldsymbol{Z}}_i^T \mathit{\boldsymbol{R}}(U_i)+(\mathit{\boldsymbol{\xi}}_i-\hat {\mathit{\boldsymbol{\xi}}}_i)^T\mathit{\boldsymbol{\beta}}_0-\delta(\hat{{\bf X}}_i^T{\bf v}_1+{\bf \Pi}_i^T{\bf v}_2)$.

Using Taylor expanding $J_1$ around $\varepsilon_i$, we have

where $\varepsilon^*_{i}$ is between $\varepsilon_{i}$ and $\varepsilon_{i}+\mathit{\boldsymbol{Z}}_i^T \mathit{\boldsymbol{R}}(U_i)+(\mathit{\boldsymbol{\xi}}_i-\hat {\mathit{\boldsymbol{\xi}}}_i)^T\mathit{\boldsymbol{\beta}}_0$.

By Corollary 6.21 in Schumaker ^[19], we have $||\mathit{\boldsymbol{R}}(U_i)|| = O_p(K^{-r}) = O_p(\delta)$. In addition, $\hat{\mathit{\boldsymbol{\xi}}}_i$ is the local linear estimators of $\mathit{\boldsymbol{\xi}}_i$, invoking condition A5 and A9 and the simple calculation yields

Similarly, we also have

Therefore, by choosing a sufficient large $C$, $J_2$ dominates $J_1$ uniformly $||{\bf v}|| = C$. Similarly, we can prove that ${J}_{3} = O_{p}(n\delta^{3}||{\bf v}||^{3}).$ Since $\delta\rightarrow 0$, it follows that $\delta ||{\bf v}||\rightarrow 0$ with $||{\bf v}|| = C$, which lead to $J_3 = o_p(J_2)$. Therefor, $J_3$ is also dominated by $J_2$ in $||{\bf v}|| = C$.

Hence, Eq (7.1) holds and there exists local maximizers $\hat{{\bf{\Theta}}}$ and $\hat{\mathit{\boldsymbol{\gamma}}}$ such that $||\hat{{\bf{\Theta}}}-{\bf{\Theta}}_0|| = O_p(\delta)$ and $||\hat{\mathit{\boldsymbol{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0|| = O_p(\delta)$. This completes the proof of part (ⅰ) in Theorem 1.

where $\textbf{H} = \int _0^1 \mathit{\boldsymbol{B}}(u) \mathit{\boldsymbol{B}}(u)^T du$.

Then invoking $||\textbf{H}|| = O(1)$ and $||\hat{\mathit{\boldsymbol{\gamma}}}-\mathit{\boldsymbol{\gamma}}_0|| = O_p(\delta)$, we have $(\hat{\mathit{\boldsymbol{\gamma}}}_k-{\mathit{\boldsymbol{\gamma}}}_{0k})^T \textbf{H} (\hat{\mathit{\boldsymbol{\gamma}}}_k-{\mathit{\boldsymbol{\gamma}}}_{0k}) = O_p(\delta^2).$ Moreover, we can get that $\int _0^1 R_k(u)^2 du = O_p(\delta^2)$. Then $||\hat \alpha_{k}(u)-\alpha_{0k}(u)|| = O_p(\delta)$ is obtained. This completes the Theorem 1. □

Proof of Theorem 2. Since ${(\hat{\mathit{\boldsymbol{\beta}}}, \hat{\mathit{\boldsymbol{\theta}}}, \hat{\mathit{\boldsymbol{\gamma}}})}$ is the solution of maximizing the objective function (2.7), we have

and

Note that

Thus, from (7.9) and Taylor expression, we have

Thus, from Lemma A.1 and (7.10), we can obtain that (7.8) has the following asymptotic expression,

By the condition A2 and A6, and some calculations based on (7.11), it follows that

where

Then, by the law of large numbers, we have $\Gamma_n\stackrel{p}\longrightarrow\Gamma$ and $\Phi_n\stackrel{p}\longrightarrow\Phi$, where $\Gamma = E[F(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h)\Pi \Pi^T]$ and $\Phi = E[F(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h)\Pi \mathit{\boldsymbol{X}}^T]$.

Substituting (7.12) into (7.7), and similar to (7.11), we obtain that

Note that

and

Let $\tilde{{\bf X}}_i = {{\bf X}}_i-\Phi_n^T\Gamma_n^{-1}{\bf \Pi}_i$, (7.13) can be expressed by

Since $\frac{1}{n}\sum_{i = 1}^n \phi_h''(\varepsilon_i)\tilde{{\bf X}}_i {\bf \Pi}_i^T = \frac{1}{n}\sum_{i = 1}^n \phi_h''(\varepsilon_i){{\bf X}}_i {\bf \Pi}_i^T- \Phi_n^T \Gamma_n^{-1}\frac{1}{n}\sum_{i = 1}^n \phi_h''(\varepsilon_i){\bf \Pi}_i{\bf \Pi}_i^T = 0$, then $I_3 = 0.$ Together with $||\mathit{\boldsymbol{R}}(U_i)|| = O(K^{-r}) = o(1)$, a simple calculation yields $I_2 = o_p(1).$

For $I_4$,

As $nb^{4}\rightarrow 0$, we have $I_{41} = o_p(1)$. Then through the similar argument of Zhang et al. ^[9], we have

Then, (7.14) can be written as

By the Slutsky's theorem and (7.15), it follows that

where ${\mathit{\boldsymbol{\Omega}}} = \mathrm{E}(G(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h){\bf \tilde{X}}{\bf \tilde{X}}^T)+\mathit{\boldsymbol{\beta}}_0^T\Sigma_e\mathit{\boldsymbol{\beta}}_0 \mathrm{E}(F(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Z}}, U, h)^2 E({\bf \tilde{X}}{\bf \tilde{X}}^T))$. This completes the proof of Theorem 2.

□