Multiple robust estimation of parameters in varying-coefficient partially linear model with response missing at random

1.   Introduction

The model considered in this paper is a classical semi-parametric model, varying-coefficient partially linear model, and it has the following form

where $Y$ is the response variable, $X$, $Z$ and $T$ are $q-$dimensional, $p-$dimensional and one-dimensional covariates, respectively. $\beta = (\beta_1, \cdots, \beta_p)^T$ is a $p-$dimensional unknown parameter vector, $\theta(\cdot) = (\theta_1(\cdot), \cdots, \theta_q(\cdot))^T$ is a $q-$dimensional unknown non-parametric function vector, $\varepsilon$ is the random error and satisfies $E(\varepsilon|X, Z, T) = 0$. Model (1.1) has been well studied by many statisticians, see the literatures, for example, Fan and Huang ^[1], You and Zhou ^[2], Huang and Zhang ^[3], Zhao ^[4], Feng, Zhang and Lu ^[5] among others.

In practical applications, missing data problems are frequently encountered in almost all research areas, such as psychological sciences, medical studies, industrial and agricultural production. The complete-case (CC) method will lose the estimation efficiency due to the disregard of the information from the missing values, and may result in biased results if the data is not missing completely at random. For details, see Little and Rubin ^[6]. The inverse probability weighted (IPW) method is another frequently used method dated back to Horvitz and Thompson ^[7] that can be applied to the case of missing covariates. This method is to take the inverse of the selection probability as the weight to the fully observed data, and under missing at random (MAR) assumption this method is unbiased. It has attracted much attention in statistical analysis with missing data, but still doesn't make full use of the incomplete data. The imputation method is a popularly method to deal with missing responses in many studies which was introduced by Yates ^[8]. The concept of imputation is to fill in each missing data with a suitable value, and then use the observed value and the imputed value for statistical inference by the standard method. This method can improve the efficiency of the resulted estimators, see the literatures, for example, Cheng ^[9], Wang and Rao ^[10,11], Wang, Linton and Hardle ^[12], and so forth. In order to further improve the efficiency of estimation, Robins, Rotnitzky and Zhao ^[13] propose an augmented inverse probability weighted (AIPW) method. This method has the double robustness, that is, if the selection probability and the conditional expectation function are both correctly specified, the resulted estimator will reach the semi-parametric effective bound, and if either of the two assumed models is correctly specified, the estimator is consistent, see the details in Robins and Rotnitzky ^[14] and Scharfstein, Rotnitzky and Robins ^[15]. In subsequent ten years, the doubly robust estimation has been well studied, see for example, Kang and Schafer ^[16], Qin, Shao and Zhang ^[17], Cao, Tsiatis and Davidian ^[18], Han ^[19], and Rotnitzky et al. ^[20].

However, double robustness does not provide sufficient protection for estimation consistency, since it allows only one model for the selection probability and one for the conditional expectation function. It is often risky to assume that one of these two models is correctly specified with an unknown data generating process. Noticed this, Han and Wang ^[21] propose multiple robust estimator for the population mean when the response variable is subject to ignorable missingness. They suggest multiple models for both the selection probability function and the outcome regression model, and the resulted estimator is consistent if any of the multiple models is correctly specified, and attains the semi-parametric efficiency bound when one selection probability and one outcome regression model are correctly specified, without requiring knowledge of which models are correct. For the details please resort to Han and Wang ^[21]. Subsequently, Han ^[22] studies the multiple robust estimator for the linear regression model. He discusses the numerical implementation of the proposed method through a modified Newton-Raphson algorithm, derives the asymptotic distribution of the resulted estimator and provides some ways to improve the estimation efficiency. Later, Sun, Wang and Han ^[23] propose multiple robust kernel estimating equations (MRKEEs) for nonparametric regression, demonstrate its multiple robustness, and show that the resulted estimator achieves the optimal efficiency within the class of augmented inverse propensity weighted (AIPW) kernel estimators when including correctly specified models for both the missingness mechanism and the outcome regression. Please refer to Sun, Wang and Han ^[23] for more discussion. In addition, the multiple robust estimation with nonignorably missing data has been studied recently, and here we just list some literatures, see for example, Han ^[24] and Li, Yang and Han ^[25].

To the best of our knowledge, the multiple robust estimation for the parameters of the varying-coefficient partially linear model with response missing at random has not been studied. So in this paper, applying the idea of Han ^[22] and Sun, Wang and Han ^[23], we consider the multiple robust estimation method for the parameters of the varying-coefficient partially linear model with missing response, and the proposed method is demonstrated superior over the existing competitors via simulation studies.

This paper is organized as follows. The proposed estimation technique and its multiple robustness are presented in Section 2. Numerical simulation studies are conducted in Section 3 in order to examine the performance of the proposed method. The technical proofs are also provided in Section 4. Conclusions are summarized in Section 5.

2.   The proposed estimator

Suppose the available incomplete data $\left\{(R_{i}, Y_{i}, X_{i}, Z_{i}, T_{i}), i = 1, 2, \ldots, n\right\}$ is a random sample from model (1.1), that is

where $R_{i}$ is an indicator variable, when $Y_{i}$ can be observed, then $R_{i} = 1$, and when $Y_{i}$ is missing with $R_{i} = 0$. The covariate $X_{i}, Z_{i}$ and $T_{i}$ are all observed. Following Han ^[22] and Sun, Wang and Han ^[23], we also suppose the auxiliary variables $S_{i}$ relate to $(R_{i}, Y_{i}, X_{i}, Z_{i}, T_{i})$ is available. Just as Han ^[22] points out that the auxiliary variables do not enter the regression model and are not of direct statistical interest, but they can reduce the impact of missing data on estimation and improve the estimation efficiency. Let $V_{i} = (X_{i}^{T}, Z_{i}^{T})^{T}$ denote the covariates. The missing mechanism we assume in this paper is MAR mechanism that commonly used in practice. Specifically, given the covariates $V_{i}$, $T_{i}$ and the available auxiliary variables $S_{i}$, the missing of $Y_{i}$ is independent of $Y_{i}$, that is,

Here we assume that $\pi(\cdot)$ is only related to $V$ and $S$.

We first carry out the estimator of the varying coefficient functions $\theta(\cdot)$. For any $t$ in a small neighborhood of $t_{0}$, using the local linear fitting for $\theta_{j}(t), j = 1, 2, \ldots, q$, we have

Suppose the parameter $\beta$ is known, and then minimizing the following objective function

about $(a_{j}, b_{j}), j = 1, 2, \ldots, q$, we can obtain the estimator of $\theta(t)$ at $t_{0}$, where $K_{h}(\cdot) = h^{-1}k(\cdot/h)$, $k(\cdot)$ is a kernel function, and $h$ is the bandwidth. Let

and

then the estimator of the coefficient functions $\theta(t)$ at $t_{0}$ is given by

Substituting (2.3) into (2.1), we obtain

where $\tilde{Y}_{i} = Y_{i}-X_{i}^{T}\hat{g}(T_{i})$, $\tilde{Z}_{i} = Z_{i}-\hat{\mu}^{T}(T_{i})X_{i}$ with $\hat{g}(t) = \sum_{k = 1}^{n}S_{k}(t)Y_{k}$ and $\hat{\mu}(t) = \sum_{k = 1}^{n}S_{k}(t)Z_{k}^{T}$.

For model (2.4), using the complete data, the CC estimator of $\beta$ can be obtained by solving the following estimation equation

where

From Little and Rubin ^[6] we know that the CC estimator maybe biased unless the missing mechanism is missing completely at random. So following the works of Robins, Rotnitzky and Zhao ^[13], the doubly robust estimator $\hat{\beta}_{AIPW}$ of $\beta$ can be defined by

where $\hat{\pi}(V_{i}, S_{i})$ is some estimated value of $\pi(V_{i}, S_{i})$, $\eta_{i}(\beta) = E[\hat{\xi}_{i}(\beta)|V_{i}, T_{i}, S_{i}]$. $\hat{\beta}_{AIPW}$ has been improved in terms of consistency, but in practice it is still a great risk to assume that one of the two assumed models is correctly specified. So inspired by Han ^[22] and Sun, Wang and Han ^[23], next we shall give the multiple robust estimation for $\beta$.

Suppose there are $J$ and $K$ models used to estimate $\pi(V, S)$ and $E(Y|V, T, S)$. Let $\mathscr{P} = \{\pi^j(\alpha^j): j = 1, \cdots, J\}$ and $\mathscr{F} = \{a^k(\gamma^k): k = 1, \cdots, K\}$ denote the set of these two models respectively, where $\alpha^j$ and $\gamma^k$ are the corresponding parameters.

Let $\hat{\alpha}^j$, $\hat{\gamma}^k$ be the estimator of $\alpha^j$, $\gamma^k$ respectively. Usually, $\hat{\alpha}^j$ can be obtained by maximizing the binomial likelihood

According to the property of MAR assumption, it can be seen that $Y$ and $R$ are conditionally independent with respect to $(V, T, S)$, that is, $E(Y|V, T, S) = E(Y|R = 1, V, T, S)$. Therefore, using the complete observation data to fit the model $a^k(\gamma^k)$, we can obtain $\hat{\gamma}^k$. Let $\hat{\beta}^k$ be the solution of

Obviously, $\hat{\beta}^k$ is an estimated value of $\beta$.

Next, let $m = \sum^{n}_{i = 1}R_i$ represents the number of the observable response variables. Without loss of generality, $R_1 = \cdots = R_m = 1$, $R_{m+1} = \cdots = R_n = 0$. Let $\omega(V, S) = \frac{1}{\pi(V, S)}$, similar to Han ^[22], the following formulas hold

where $j = 1, \cdots, J$, $k = 1, \cdots, K$, $U^k(\beta, \gamma^k) = \{Z-\mu^{T}(T)X\}\{a^k(\gamma^k)-X^T\theta(T)-Z^T\beta\}$. Therefore, the weights $\omega_i, i = 1, \cdots, m$ can be defined by

where

Based on the empirical likelihood method, under the above constraints, the Lagrange multiplier method is used to solve the maximum value problem of $\prod^m_{i = 1}\omega_i$, and we use the solution as the weight $\omega_i(i = 1, \cdots, m)$ to estimate the parameter $\beta$. For ease of presentation, let $\hat{\alpha}^T = \{(\hat{\alpha}^1)^T, \cdots, (\hat{\alpha}^J)^T\}$, $\hat{\beta}^T = \{(\hat{\beta}^1)^T, \cdots, (\hat{\beta}^K)^T\}$, $\hat{\gamma}^T = \{(\hat{\gamma}^1)^T, \cdots, (\hat{\gamma}^K)^T\}$, and $\hat{g}_i(\hat{\alpha}, \hat{\beta}, \hat{\gamma})^T = [\pi_i^1(\hat{\alpha}^1)-\nu^1(\hat{\alpha}^1), \cdots, \pi_i^J(\hat{\alpha}^J)-\nu^J(\hat{\alpha}^J),$ $\{\hat{U}_i^1(\hat{\beta}^1, \hat{\gamma}^1)-\eta^1(\hat{\beta}^1, \hat{\gamma}^1)\}^T, \cdots, \{\hat{U}_i^K(\hat{\beta}^K, \hat{\gamma}^K)-\eta^K(\hat{\beta}^K, \hat{\gamma}^K)\}^T]$. Based on the empirical likelihood theory, we have

where $\hat{\rho}^T = (\hat{\rho}_1, \cdots, \hat{\rho}_{J+pK})$ is the $(J+pK)$-dimension Lagrange multiplier, and is the solution of

Due to the non-negativity of the weight $\hat{\omega}_i$, $\hat{\rho}$ satisfies

So we can solve the equation

to obtain the multiple robust estimator of the parameter $\beta$, denoted by $\hat{\beta}_{MR}$.

In calculation of the weight $\hat{\omega}_{i}$, the Lagrange multiplier $\hat{\rho}$ is essential. The calculation algorithm we used is similar to Han ^[22], for the details please refer to Han ^[22], here we omit.

The multiple robustness of $\hat{\beta}_{MR}$ is given by the following theorem.

Theorem 2.1. Suppose that the conditions C1–C5 in Section 4 hold, and if $\mathscr{P}$ contains a model that correctly specifies $\pi(V, S)$, or $\mathscr{F}$ contains a correctly specified model for $E(Y|V, T, S)$, then $\sum_{i = 1}^{m}\hat{\omega}_{i}\hat{\xi}_{i}(\beta) \stackrel{P}\longrightarrow 0$ with $n\rightarrow \infty$.

3.   Simulation study

In this section, we conduct some numerical simulations to evaluate the feasibility of the above method and the finite sample performance of the proposed estimator $\hat{\beta}_{MR}$. Several indices of multiple robust estimates, inverse probability weighted estimates, and augmented inverse probability weighted estimates are compared and analyzed under different sample sizes.

We consider five mutually independent covariates, namely: $X\sim N(0, 1), T\sim U(0, 1), Z_{1}\sim N(1, 5), Z_{2}\sim B(0.5, 1), Z_{3}\sim N(0, 1)$. The response variable is generated by the model $Y = X^T\theta(T)+Z^T\beta+\varepsilon$, where $\theta(t) = \sin(\pi t)$ and $\beta = (1, 1, 2)^{T}$. In addition, We consider three auxiliary variable, namely $S^{(1)} = 1+Z^{(1)}-Z^{(2)}+\varepsilon_1$, $S^{(2)} = I\{S^{(1)}+0.4\varepsilon_2 > 2.8\}$, $S^{(3)} = \exp[\{S^{(1)}/9\}^2]+\varepsilon_3$, where $I(\cdot)$ is an indicator function. $(\varepsilon, \varepsilon_1, \varepsilon_2, \varepsilon_3)^T\sim N(0, \Sigma)$. The diagonal elements of the matrix $\Sigma$ are $1, 0.5, 1, 2$, the elements at positions $(1, 2)$ and $(2, 1)$ are $0.5$, and the remaining elements are all $0$. The probability of selection is $\rm{logit}\{\pi(V, S)\} = 3.5-5S^{(2)}$, under which there are approximately $34\%$ of the subjects with missing $Y$. The models for correctly estimating $\pi(V, S)$ and $E(Y|V, T, S)$ are $\rm{logit}\{\pi^1(\alpha^1)\} = \alpha_1^1+\alpha_2^1S^{(2)}$ and $a^1(\gamma^1) = X^{T}\theta(T)+\gamma_1^1Z_{1}+\gamma_2^1Z_{2}+\gamma_3^1Z_{3}+\gamma_4^1S^{(3)}$ respectively. In addition, we also use two incorrect models in the simulation process, namely $\rm{logit}\{\pi^2(\alpha^2)\} = \alpha_1^2+\alpha_2^2Z_{1}+\alpha_3^2Z_{2}+\alpha_4^2Z_{3}$, $a^2(\gamma^2) = X^{T}(-4T^{2}+4T)+\gamma_1^2Z_{1}+\gamma_2^2Z_{2}+\gamma_3^2Z_{3}+\gamma_4^2S^{(3)}$. For simplicity, we use the Rule of Thumb method to obtain the optimal bandwidth when estimating the nonparametric functions, that is, $h = 1.06*\{min(qr, sig)\}*n^{-1/5}$, where $sig$ is the standard deviation of covariate $T$, $qr = (Q_{3}-Q_{1})/1.34$, $Q_{1}$ and $Q_{3}$ are the first and third quartile, respectively. In simulation, we generate random samples with $n = 200$ and $n = 500$ respectively, and repeat the process $500$ times to calculate the average biases, mean squared errors (MSEs), the root of mean squared errors (RMSEs) and median absolute error (MAEs).

In order to verify the superiority of the multiple robust estimation method, we give the calculated indices of the parameter $\beta$ under different estimation methods, which are the inverse probability weighted estimates $\hat{\beta}_{IPW}$, and the augmented inverse probability weighted estimates $\hat{\beta}_{AIPW}$ and multiple robust estimates $\hat{\beta}_{MR}$. To distinguish all the estimators constructed based on different methods and models, each estimator is assigned a name with the form "Method-$0000$", where each digit of the four-digit number, from left to right, indicates whether $\pi^1(\alpha^1), \pi^2(\alpha^2), a^1(\gamma^1), a^2(\gamma^2)$ is used in the construction ($1$ means yes, $0$ means no), respectively. The simulation results are reported in Table 1 and Table 2 with the sample size $n = 200$ and $n = 500$.

It can be seen from the two tables that regardless of the estimation method, the larger the sample size, the better the estimation effect. And when the models for estimating the selection probability and the conditional expectation are all specified correctly, the estimated results obtained by the multiple robust estimation method, the inverse probability weighted estimation method and the augmented inverse probability weighted estimation method are not much different, but the effect of multiple robust estimation is better in terms of MSE. When all the models for estimating the selection probability and the conditional expectation are specified incorrectly, the $AIPW-0101$ has unsatisfactory effects, the resulted estimators have larger deviations, but our proposed $MRE-0101$, despite using two incorrect models, can generate better estimators. The interesting observation that $\hat{\beta}_{MR}$ seems to still provide a reasonable (at least not too bad) estimate of $\beta$ even if there is no model correctly specified is similar to Han ^[22]. In a word, it is obvious that our proposed multiple robust estimation method is better than the two competitors.

4.   Proofs

Before we give the proof of Theorem 2.1, some notations and interpretations are presented firstly.

Let $\Phi(t) = E[RXZ^{T}|T = t]$, $\Psi(t) = E[RXX^{T}|T = t]$, then

Substituting (4.1) into (2.1), we obtain

where $\check{Y}_{i} = Y_{i}-X_{i}^{T}g(T_{i})$, $\check{Z}_{i} = Z_{i}-\mu^{T}(T_{i})X_{i}$, with $g(T_{i}) = \{\Psi(T_{i})\}^{-1}E[R_{i}X_{i}Y_{i}|T_{i}]$, $\mu(T_{i}) = \{\Psi(T_{i})\}^{-1}\Phi(T_{i})$. From model (4.2), using the complete data, the CC estimator of $\beta$ can be obtained by solving the following estimation equation

where $\xi_{i}(\beta) = \check{Z}_{i}(\check{Y}_{i}-\check{Z}_{i}^{T}\beta) = (Z_{i}-\mu^{T}(T_{i})X_{i})[Y_{i}-X_{i}^{T}g(T_{i})-(Z_{i}-\mu^{T}(T_{i})X_{i})^{T}\beta]$, and $E[\xi_{i}(\beta)] = 0.$

Suppose $C$ to be a positive constant which can represent different values, and assume the following conditions C1–C5 hold.

C1 The bandwidth $h$ satisfies $h = Cn^{-1/5}$, that is $h\rightarrow0$ and $nh\rightarrow \infty$ as $n\rightarrow \infty$, where $C > 0$ is a given positive constant.

C2 The kernel function $K(\cdot)$ is a symmetric probability kernel function, and $\int t^{2}K(t)dt \neq 0$, $\int t^{4}K(t)dt < \infty$.

C3 For each $t\in(0, 1)$, $f(t), \Phi(t), \Psi(t)$ and $\theta(t)$ are twice continuous differentiable at point $t$, where $f(t)$ is the density function of the variable T.

C4 $\sup \limits_{0\leq t\leq 1}E[\varepsilon_{i}^{4}|T_{i} = t] < \infty$, $\sup \limits_{0\leq t\leq 1}E[X_{ir}^{4}|T_{i} = t] < \infty$, and they are continuous about $t$, where $X_{ir}$ is the $r$-th component of $X_{i}$, $i = 1, \cdots, n$, $r = 1, \cdots, q$.

C5 For a given $t$, $\Psi(t)$ is a positive definite matrix.

Next, a Lemma is needed in proof of Theorem 2.1, and the proof can be found in Zhao ^[4].

Lemma 4.1. Suppose conditions C1–C5 hold, then we have

where $C_{n} = h^{2}+(\frac{\log(1/h)}{nh})^{1/2}$.

Proof of Theorem 2.1: Assuming that $\mathscr{P}$ contains a model that correctly specifies $\pi(V, S)$, without loss of generality, let $\pi^1(\alpha^1)$ be the model, $\alpha_0^1$ represents the truth value of $\alpha^1$, that is $\pi^1(\alpha_0^1) = \pi(V, S)$. Next, we combine the theory of empirical likelihood to prove that $\hat{\beta}_{MR}$ is a consistent estimator of $\beta$.

Referring to the method in Han ^[22] to establish the relationship between the weight $\hat{\omega}_i$ and the empirical likelihood on the biased sample. Let $p_i$ represent the conditional empirical probability on the biased sample $(Y_i, X_i, Z_i, T_i, S_i), R_i = 1, i = 1, \cdots, m$, based on $(2.8), (2.9)$ and $\omega(V, S) = \frac{1}{\pi^1(\alpha_0^1)}$, a more reasonable value of $p_i$ can be given by the following constrained optimization problem:

Using the Lagrange multiplier method again, we get

where $\hat{\lambda}^T = (\hat{\lambda}_1, \cdots, \hat{\lambda}_{J+pK})$ is the $(J+pK)$-dimensional Lagrange multiplier, and satisfies

Due to the non-negativity of $\hat{p}_i$, $\hat{\lambda}$ satisfies $1+\hat{\lambda}^T\hat{g}_i(\hat{\alpha}, \hat{\beta}, \hat{\gamma})/\pi_i^1(\hat{\alpha}^1) > 0, i = 1, \cdots, m.$ Since

then the solution of $(2.11)$, $\hat{\rho}$, can be written as $\hat{\rho}_1 = (\hat{\lambda}_1+1)/\nu^1(\hat{\alpha}^1)$ and $\hat{\rho}_l = \hat{\lambda}_l/\nu^1(\hat{\alpha}^1), l = 2, \cdots, J+pK.$ Therefore

Just like White ^[24], let $\alpha_\ast^j$, $\beta_\ast^k$ and $\gamma_\ast^k$ are the minimum points of the corresponding Kullback-Leibler distance respectively, then we have $\hat{\alpha}^j\stackrel{P}\longrightarrow\alpha_\ast^j, \hat{\beta}^k\stackrel{P}\longrightarrow\beta_\ast^k, \hat{\gamma}^k\stackrel{P}\longrightarrow\gamma_\ast^k$, and $n^{1/2}(\hat{\alpha}^j-\alpha_\ast^j)$, $n^{1/2}(\hat{\beta}^k-\beta_\ast^k)$ and $n^{1/2}(\hat{\gamma}^k-\gamma_\ast^k)$ are bounded by probability. At the same time, $\nu^j(\hat{\alpha}^j)\stackrel{P}\longrightarrow\nu_\ast^j$, $\eta^k(\hat{\beta}^k, \hat{\lambda}^k)\stackrel{P}\longrightarrow\mu_\ast^k$, where $\nu_\ast^j = E[\pi^j(\alpha_\ast^j)]$, $\mu_\ast^k = E[U^k(\beta_\ast^k, \gamma_\ast^k)]$. Generally speaking, when the model $\pi^j(\alpha^j)$ for $\pi(V, S)$ is correctly specified, we have $\pi^j(\alpha_\ast^j) = \pi(V, S)$, and when the model $a^k(\gamma^k)$ for $E(Y|V, T, S)$ is correctly specified, we have $a^k(\gamma_\ast^k) = E(Y|V, T, S)$. Let $\alpha_\ast^T = \{(\alpha_\ast^1)^T, \cdots, (\alpha_\ast^J)^T\}$, $\beta_\ast^T = \{(\beta_\ast^1)^T, \cdots, (\beta_\ast^K)^T\}$, $\gamma_\ast^T = \{(\gamma_\ast^1)^T, \cdots, (\gamma_\ast^K)^T\}$, and suppose $\hat{\rho}\stackrel{P}\longrightarrow\rho_\ast$.

Based on the empirical likelihood theory, it can be known that $\hat{\lambda}\stackrel{P}\longrightarrow0$. According to the appendix in Han ^[22], $\hat{\lambda} = O_p(n^{-1/2})$ holds. Since the model $\pi^1(\alpha^1)$ is correct, then we have $\frac{m}{n}\stackrel{P}\longrightarrow\nu_\ast^1$, and

Refer to Zhao ^[4], since

and $E[X_{i}\varepsilon_{i}] = 0$, $E[(Z_{i}-\mu^{T}(T_{i})X_{i})X_{i})^{T}] = 0$, we have

Combine conditions C1, C4, C5 and Lemma 4.1, we have

That is, $\frac{1}{n}\sum\limits_{i = 1}^{n}\frac{R_i}{\pi_i^1(\alpha_\ast^1)}\hat{\xi}_i(\beta)\stackrel{P}\longrightarrow\frac{1}{n}\sum\limits_{i = 1}^{n}\frac{R_i}{\pi_i^1(\alpha_\ast^1)}\xi_i(\beta)$. Then we have

Therefore, when $n\rightarrow \infty$, $\beta$ is the solution of the formula $(2.13)$, which shows that $\hat{\beta}_{MR}$ is a consistent estimator of $\beta$.

Next, suppose that $\mathscr{F}$ contains a model that correctly specifies $E(Y|V, T, S)$. Without loss of generality, let $a^1(\gamma^1)$ be the true model and $\gamma_0^1$ be the true value of $\gamma^1$, that is $a^1(\gamma_0^1) = E(Y|V, T, S)$, and $\gamma_\ast^1 = \gamma_0^1$. A previous constraint is actually

and $\hat{\beta}^1\stackrel{P}\longrightarrow\beta_\ast^1 = \beta$, so we get $\frac{1}{n}\sum\limits_{i = 1}^{n}\hat{U}_i^1(\hat{\beta}^{1}, \hat{\gamma}^{1})\stackrel{P}\longrightarrow0.$

Let $g(\alpha_\ast, \beta_\ast, \gamma_\ast)^T = [\pi^1(\alpha_\ast^1)-\nu_\ast^1, \cdots, \pi^J(\alpha_\ast^J)-\nu_\ast^J, \{U^1(\beta_\ast^1, \gamma_\ast^1)-\eta_\ast^1\}^T, \cdots, \{U^K(\beta_\ast^K, \gamma_\ast^K)-\eta_\ast^K\}^T],$ due to $\frac{1}{n}\sum\limits_{i = 1}^{n}\hat{U}_i^1(\hat{\beta}^1, \hat{\gamma}^1)\stackrel{P}\longrightarrow\frac{1}{n}\sum\limits_{i = 1}^{n}U_i^1(\beta, \gamma_0^1), \|\frac{1}{n}\sum\limits_{i = 1}^{n}[\hat{\xi}_{i}(\beta)-\xi_i(\beta)]\|\stackrel{P}\longrightarrow0,$ and $E[U^1(\beta, \gamma_0^1)] = 0$, then we have

This shows that $\hat{\beta}_{MR}$ is a consistent estimator of $\beta$.

So the proof of Theorem 2.1 is completed.

5.   Conclusions

In this article, we have proposed the multiple robust estimators for parameters in varying-coefficient partially linear model with missing response at random, and the multiple robustness of our proposals has been shown theoretically under some regular conditions. Our simulation studies fully demonstrate the superiority of our multiple robust estimation method through Table 1 and Table 2. Finally, we point out some problems for the future researches. First, we only discuss the multiple robust estimation process of parameters, and the fitting of nonparametric function curves can be expanded in the future studies. Next, based on the model in this article, if the missing mechanism is nonignorable missing, how to obtain the robust estimation of parameters is also worth studying.

Acknowledgments

The research is supported by NSF projects (ZR2021MA077 and ZR2019MA016) of Shandong Province of China.

Conflict of interest

The authors declare that they have no conflicts of interest to this work.