数学物理学报, 2020, 40(2): 460-474 doi:

论文

缺失数据下部分非线性变系数EV模型的统计推断

马奕佳,, 薛留根, 芦飞

Statistical Inference in Partially Nonlinear Varying-Coefficient Errors-in-Variables Models with Missing Responses

Ma Yijia,, Xue Liugen, Lu Fei

通讯作者: 马奕佳, E-mail: redbamboo12138@126.com

收稿日期: 2018-05-18  

基金资助: 国家自然科学基金.  11971001
北京市自然科学基金.  1182002

Received: 2018-05-18  

Fund supported: the NSFC.  11971001
the Beijing Natural Science Foundation.  1182002

摘要

该文研究了响应变量缺失下半参数部分非线性变系数EV模型的统计推断问题,利用逆概率加权局部纠偏profile最小二乘法构造了模型中非参数分量和参数分量的估计,证明了估计量的渐近正态性.通过数值模拟和实际数据分析,验证了所提出的估计方法是有效的.

关键词: 部分非线性变系数模型 ; 缺失数据 ; 局部纠偏 ; 测量误差 ; 渐近正态性

Abstract

This paper considers about the estimation of varying-coefficient partial nonlinear errors-in-variables models with missing responses. Firstly, we develop inverse probability weighted approaches and local bias-corrected restricted profile least squares estimators. Asymptotic normality of estimators is established. Moreover, both simulation results and a real data show that local bias-corrected restricted profile least squares estimated approach are better than the performance ignoring the measurement error.

Keywords: Varying-coefficient partial nonlinear model ; Missing data ; Local bias-corrected restricted profile least-squares approach ; Measurement error ; Asymptotic normality

PDF (601KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

马奕佳, 薛留根, 芦飞. 缺失数据下部分非线性变系数EV模型的统计推断. 数学物理学报[J], 2020, 40(2): 460-474 doi:

Ma Yijia, Xue Liugen, Lu Fei. Statistical Inference in Partially Nonlinear Varying-Coefficient Errors-in-Variables Models with Missing Responses. Acta Mathematica Scientia[J], 2020, 40(2): 460-474 doi:

1 引言

考虑半参数部分非线性变系数模型

$ \begin{equation} \left. Y = X^{T} a(U) + g(Z, \beta) + \varepsilon, \right. \end{equation} $

其中$ Y $是响应变量, $ \beta = (\beta_1, \cdots , \beta_r)^{T} $是未知的$ r $维参数向量,变系数部分$ a(\cdot) = (a_1(\cdot), \cdots, $$ a_q(\cdot))^{T} $$ q\times1 $维未知非参数函数向量, $ X, Z, U $均为协变量, $ \varepsilon $是均值为0的模型误差.为免于维数灾祸,假定$ U $为一维变量. $ g({Z}, {\beta}) $是已知的非线性函数向量, $ Z $$ \beta $是未知的有限维参数,它们的维数不一定相同,为方便起见,证明过程中假设两者维数相同, $ \varepsilon $是均值为0且方差为$ \sigma^2 $的随机误差.在实际应用中, $ Y $$ Z $一般不是线性关系,所以将模型拓展到半参数部分非线性变系数模型.当非线性函数部分为0时,是经典的变系数模型,详见Fan和Zhang[1].当$ q = 1 $, $ X = 1 $时,模型成为部分非线性模型. Li和Nie[2]在分析生态学中真实数据集时提出此模型.当没有协变量$ X $时,模型降为著名的非线性回归模型.许多例子和应用可以在Bates和Watts[3]中有所发现.

就作者所知,目前将响应变量缺失,协变量带测量误差这一复杂数据情况应用到部分非线性变系数模型的结果还较少.然而在实际问题中,经常会遇到数据缺失现象,缺失数据的处理一直受到统计学家的重视.另外,在多数情况下很可能得不到协变量的准确观测值,只能得到含有误差的观测值.也就是说变量带测量误差的模型在响应变量缺失下的研究是很有必要的,带测量误差的模型统称为EV模型.针对缺失数据, Zhou等[4]研究了一类存在缺失数据时的借补估计方程方法,这些方法可以用到部分线性及非线性模型中来; Xue和Li[5]构造了响应变量有缺失数据的填补估计,证明了估计量的渐近性质,并在响应变量缺失下构造了线性模型中参数的经验似然置信域.针对测量误差数据, Carroll和Ruppert[6]详细讨论了非线性EV模型;冯三营等[7]研究了未知系数函数的经验似然问题,利用纠偏方法得到了似然比统计量,并给出了相关的性质、证明及模拟.针对缺失数据和测量误差数据同时存在的情况, Liang等[8]讨论了协变量含测量误差和响应变量存在缺失时的部分线性模型的统计推断问题; Wei[9]针对部分线性变系数模型研究了响应变量缺失和协变量带有测量误差情况下的估计问题.基于以上原因,本文研究响应变量缺失且协变量带有测量误差情况下的半参数部分非线性变系数模型的统计推断问题,采用局部线性回归技术(详见文献[10])来估计非参数函数分量$ a(u) $.然后,对其进行逆概率加权和纠偏的处理,最终得到了$ \beta $的profile最小二乘估计及其良好的渐近性质.另外,在选择概率未知的情况下,使用logistic模型估计选择概率并结合局部纠偏逆概率加权最小二乘法运用到部分非线性变系数EV模型中,实现该方法的推广.最后,通过数值模拟和实例分析验证了该模型强大的适应性.

2 估计方法与主要结果

假设$ \{(Y_i, X_i, Z_i, U_i), i = 1, 2, \cdots, n\} $为来自总体模型的独立同分布样本,则部分非线性变系数模型具有如下形式:

$ \begin{equation} \left\{ \begin{array}{ll} Y_i & = X_i^{T} a(U_i) + g(Z_i, \beta) + \varepsilon_i, \\ W_i& = X_i+e_i, \\ \end{array} \right. \end{equation} $

其中$ \{\varepsilon_i, i = 1, \cdots, n\} $是独立同分布的随机误差,满足$ E(\varepsilon_i|X_i, U_i, Z_i) = 0 $,方差为$ \sigma^2<\infty $. $ W_i $$ X_i $的替代变量,测量误差$ \{e_i, i = 1, \cdots, n\} $独立同分布,与$ (X_i^{T}, Z_i^{T}, U_i, \varepsilon_i)^{T} $相互独立且均值为0.假定测量误差的协方差阵$ \Sigma_e $已知.另外还需要讨论测量误差对缺失机制类型的影响.若$ X_i $可以被精确观测,则响应变量$ Y $是随机缺失的.但由于$ X_i $带有测量误差,其精确值不能得到.所以假设缺失机制为$ \rm MAR $来处理测量误差带来的缺失机制的变化, $ \delta = 1 $表示$ Y $值可观测, $ \delta = 0 $表示$ Y $缺失,此机制表明选择概率只依赖于观测到的数据,即

$ \begin{equation} {\Bbb P}(\delta_i = 1|W_i, U_i, Z_i, Y_i) = {\Bbb P}(\delta_i = 1|W_i, U_i, Z_i)\equiv\pi(W_i, U_i, Z_i) = :\pi(V_i), \end{equation} $

其中$ V_i = (W_i^{T}, U_i, Z_i^{T})^{T} $.假设$ \beta $已知,则模型的第一个式子可以写为如下形式:

$ \begin{equation} Y_i-g(Z_i, \beta) = X_i^{T} a(U_i)+\varepsilon_i, \; \; \; \; i = 1, \cdots, n. \end{equation} $

模型(2.3)可以看成是传统的变系数模型.所以可以利用局部线性方法估计变系数函数$ a(\cdot) $.$ U $的邻域内任意一点$ u $,对$ a_j(U) $进行泰勒展开,用一个线性函数局部逼近$ a_j(u) $,即: $ a_j(U)\approx a_j(u)+a'_j(u)(U-u), $转化为极小化下面的加权最小二乘问题来估计系数函数:

$ \begin{equation} \sum^n\limits_{i = 1}\Big\{Y_i-g(Z_i, \beta)-\sum^q\limits_{j = 1}[a_j(u)+a'_j(u)(U_i-u)]X_{ij}\Big\}^2 \frac{\delta_i}{\pi(V_i)}K_h(U_i-u), \end{equation} $

其中$ K_h(\cdot) = K(\cdot/h)/h $, $ K(\cdot) $是核函数, $ h $是带宽且为收敛于0的常数.引入记号:

极小化$ (2.4) $式,不难得到$ \psi(u) $的估计是

$ \begin{equation} {\bf \hat{\psi}}(u;\beta) = \{(D^{X}_{u})^{T}\omega_\delta(D^{X}_{u})\}^{-1} (D^{X}_{u})^{T}\omega_\delta[Y-g(\beta)]. \end{equation} $

由于$ X_i $不可观测,若直接用$ W_i $代替$ X_i $,则得到的估计$ (2.5) $不是相合估计.为了消除测量误差对估计造成的影响,采用文献[11]的方法,对$ (2.5) $式进行局部纠偏,得到

$ \begin{equation} {\bf \hat{\psi}}(u;\beta) = \{(D^{W}_{u})^{T}{\omega}_\delta(D^{W}_{u})-{\Omega}(u)\}^{-1} (D^{W}_{u})^{T}{\omega}_\delta[Y-g(\beta)]. \end{equation} $

于是当给定$ \beta $时,非参数函数的估计为

$ \begin{equation} \hat{a}(u;\beta) = (I_q, 0_{q\times q})\{(D^{W}_{u})^{T}{\omega}_\delta(D^{W}_{u})-{\Omega}(u)\}^{-1} (D^{W}_{u})^{T}{\omega}_\delta[Y-g(\beta)], \end{equation} $

其中$ \Omega(u) = \sum^n\limits_{i = 1} \Sigma_e \otimes\left(\begin{array}{cc} 1 & \frac{U_i-u}{h}\\\frac{U_i-u}{h}&(\frac{U_i-u}{h})^2\end{array}\right)K_h(U_i-u)\left[\frac{\delta_i}{\pi(V_i)}\right] $, $ \otimes $表示Kronecker乘积.极小化

$ \begin{equation} R(\beta) = \sum^n\limits_{i = 1}\{ [Y_i-{W_i^{T}}\hat{{a}}(U_i, \beta) -g({Z}_i, \beta)]^2 -\hat{{a}}(U_i, \beta)^{T}\Sigma_e \hat{{a}}(U_i, \beta)\}\frac{\delta_i}{\pi(V_i)}K_h(U_i-u), \end{equation} $

可得到$ \beta $的逆概率加权局部纠偏profile最小二乘估计,估计方程和算法步骤如下:

$ \begin{equation} \sum^n\limits_{i = 1}\frac{\delta_i}{\pi(V_i)}J_i(\beta)(Y_i-f({V_i}, \beta))- \sum^n\limits_{i = 1}\frac{\delta_i}{\pi(V_i)}\frac{\partial\hat{{a}}(U_i, \beta)} {\partial\beta}\Sigma_e\hat{{a}}(U_i, \beta) = 0, \end{equation} $

其中$ J_i(\beta) = \partial f(V_i, \beta)/\partial\beta, f({V_i}, \beta) = {W_i^{T}}\hat{{a}}(U_i, \beta)+g({Z}_i, \beta). $迭代公式为$ \beta^{(k)} = \beta^{(k-1)}+(J^{T}\Delta J+B)^{-1}(J^{T}\Delta\varepsilon+A), k $为迭代次数. $ J = (J_1, \cdots, J_n)^{T}, A = \sum\limits_{i = 1}^n\frac{\partial\hat{{\bf a}}^{T}(U_i, {\beta})}{\partial\beta}\Delta\otimes\Sigma_e\hat{{a}}(U_i, {\beta}), $$ B = \sum\limits_{i = 1}^ng'(Z_i, \beta)Q_h^{T}\Delta\otimes\Sigma_e Q_hg'(Z_i, \beta), Q_h = (I_q, 0_{q\times q})\{(D^{W}_{u})^{T}{\omega}_\delta(D^{W}_{u})-{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}{\omega}_\delta, $$ \varepsilon = (Y_1-f({V}_1, \beta), \cdots, Y_n-f({V}_n, \beta))^{T}. $

步骤1   给定初始估计$ {\beta}^{(0)} $,可得到非参数函数的初始估计为

步骤2  由$ \hat{a}^{(0)}(u) $,代入(2.8)式并极小化$ R(\beta) $可得$ \beta $的估计$ {\beta}^{(k)}. $

步骤3  令$ {\beta}^{(0)} = {\beta}^{(k)}, \hat{a}^{(0)}(u) = \hat{a}^{(k)}(u;\beta) $,重复计算步骤1和步骤2直到收敛,得到最终估计为$ \hat{\beta}, \hat{a}(u;\hat{\beta}) $,记$ \hat{a}(u) = \hat{a}(u;\hat{\beta}), \hat{\psi}(u) = \hat{\psi}(u;\hat{\beta}) $.

考虑真实参数$ \beta $和函数$ a(u) $,引入记号$ A^{\otimes2} = AA^{T} $,其中$ A $是向量或矩阵, $ {g}'(Z, \beta)_{r\times 1} = \partial g(Z, \beta)/\partial \beta $, $ {\bf g}''(Z, \beta)_{r\times r} = \partial^{2}{g}(Z, \beta)/\partial \beta\partial\beta^{T} $, $ \mu_j = \int u^jK(u){\rm d}u, \; \nu_j = \int u^jK(u){\rm d}u $, $ \Gamma(U) = E[XX^{T}|U], \; \Phi(U) = E[X{g}'({Z}, \beta)^{T}|U] $, $ C_n = \{\frac{\log(1/h)}{nh}\}^{1/2}+h^{2}, \Sigma = E[{\bf g}'(Z, \beta) $$ {\bf g}'(Z, \beta)^{T}]-E[\Phi(U)^{T}\Gamma(U)^{-1}\Phi(U)]. $为了证明估计量的大样本性质,首先给出一些正则条件.

C1:随机变量$ U $具有有界支撑$ {\Omega} $,它的密度函数$ f(u) $是二次连续可微的,并在其支撑集$ {\Omega} $下有界且远离0.

C2: $ \Gamma(u) $对任何$ u\in{\Omega} $是非奇异的, $ \Gamma(u) $$ \Phi(u) $均具有连续的二阶导数且$ \Gamma(u)^{-1} $是Lipschitz连续的.

C3:存在正数$ s>2.5 $使得$ E(||{W}||^{2s})<\infty $, $ E(||{X}||^{2s})<\infty $, $ E(||{g}'(Z, \beta)||^{2s})<\infty $, $ E({\bf \varepsilon}^{2s})<\infty $,且存在某个$ 0<\delta<2-s^{-1} $使得$ n^{2\delta-1}h\rightarrow\infty $.

C4:非参数函数$ a(u) $具有连续的二阶导数.

C5:核权函数$ K(\cdot) $是具有紧支撑的对称密度函数,满足Lipschitz连续,函数$ u^3K(u) $, $ u^3K'(u) $有界且$ \int u^4K(u){\rm d}u<\infty $.

C6:当$ n\rightarrow\infty $时, $ nh^5\rightarrow0, nh^2/{{\rm log}^{2}(1/h)}\rightarrow\infty $.

C7:对任何$ z $, $ g(z, \beta) $连续且具有连续二阶导数, $ \beta\in{ß} $, $ {ß} $是一个紧集.

C8: $ E[{g}'(Z, \beta)^{\otimes2}], E[E({g}'(Z, \beta)|U)^{\otimes2}], \; E[E( $Vech$ \{{\bf g}''(Z, \beta)\}|U)^{\otimes2}] $$ \beta $邻域内有界.

C9:矩阵$ E[{g}'(Z, \beta){g}'(Z, \beta)^{T}]-E[\Phi(U)^{T}\Gamma(U)^{-1}\Phi(U)] $是正定矩阵.

C10:选择概率$ \pi(V, \theta) $关于$ \theta $二阶可导且存在常数$ C_0 $,使$ \rm {inf}\{\pi(V, \theta)\}\geq{ \text sl C_0}>0 $.

注2.1  条件C1, C3–C4, C7–C8是常用于部分线性变系数模型和非线性函数中的条件.条件$ \rm C2, C9 $保证参数估计$ \hat{\beta} $渐近方差的存在性,条件$ \rm C5–C6 $是限制核函数和带宽的常规条件,条件$ \rm C10 $为逆概率加权法的使用提供保证.

定理2.1  在条件$ \rm C1–C10 $且缺失概率已知情况下,若$ \beta $是参数真实值,则

其中$ \Lambda = E[{\pi(V)^{-1}}((g'(Z, \beta)-\Phi(U)^{T}\Gamma(U)^{-1}X)(\varepsilon-e^{T}a(U)))^{\otimes2}]+\sigma^{2}E[{\pi(V)^{-1}}\Phi(U)^{T} \Gamma(U)^{-1}\\ \Sigma_e\Gamma(U)^{-1}\Phi(U)]+E[{\pi(V)^{-1}}(\Phi(U)^{T}\Gamma(U)^{-1}(ee^{T}-\Sigma_e)a(U))^{\otimes2}]; $$ \stackrel{D}\longrightarrow $表示依分布收敛.

定理2.2  在条件$ \rm C1–C10 $且缺失概率已知情况下,如果$ a(u) $是系数函数真值,则

其中$ \psi(u) = (a(u)^{T}, ha'(u)^{T})^{T}, \; \eta = \Sigma_e-ee^{T}-Xe^{T}, \Sigma^{*} = \Gamma^{-1}(u)[\sigma^{2}(E({\pi(V)^{-1}}XX^{T}|U = u)+E({\pi(V)^{-1}}ee^{T}|U = u))+E({\pi(V)^{-1}}\eta a(u)a(u)^{T}\eta^{T}|U = u)]\Gamma^{-1}(u). $进而有

注2.2  定理$ \rm 2.1 $表明利用逆概率加权局部纠偏$ \rm profile $最小二乘法得到的参数估计具有渐近正态性.定理$ \rm 2.2 $表明系数函数的渐近性质.值得注意的是,当非线性函数$ g(Z, \beta) $取为$ Z^{T}\beta $时且选择概率$ \pi(V) = I $时, $ \hat{\beta} $与冯三营[7]提出的估计量具有相同的渐近方差; $ \hat{a}(u) $与冯三营[7]以及You, Zhou和Chen[11]提出的估计量具有相同的渐近方差.

下面给出选择概率未知时估计的大样本性质.为减少“维数灾祸”现象的影响,本文考虑用$ \rm logistic $模型作为缺失机制,即$ \pi(V) = \frac{\exp(\theta^{T}(1, V))}{1+\exp(\theta^{T}(1, V))} $,其中$ V = (W^{T}, U, Z^{T})^{T} $, $ \theta = (\theta_0, \theta_1^{T}, \theta_2, \theta_3^{T})^{T} $是未知参数向量.文献[15]给出了$ \theta $的极大似然估计$ \hat{\theta} $,从而得到选择概率的估计$ \pi(V, \hat{\theta}) $,以及$ {\hat{\omega}}_\delta = \omega(u)\hat{\Delta}, \; \hat{\Omega}(u) = \Omega(u, \hat{\theta}), \; \hat{S}_h = ({W}^{T}, 0)\{(D^{W}_{u})^{T}{\hat{\omega}}_\delta(D^{W}_{u})-\hat{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}{\hat{\omega}}_\delta $, $ {\hat{Q}}_h = (I_q, 0_{q\times q})\{(D^{W}_{u})^{T}{\hat{\omega}}_\delta(D^{W}_{u})-{\hat{\Omega}}(u)\}^{-1}(D^{W}_{u})^{T}{\hat{\omega}}_\delta $.进而可以得到$ a(u) $的估计$ \tilde{a}(u) $, $ \beta $的估计$ \tilde{\beta} $,其中$ \tilde{a}(u) = \tilde{a}(u;\beta) = (I_q, 0_{q\times q})\{(D^{W}_{u})^{T}\hat{{\omega}}_\delta(D^{W}_{u})-\hat{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}\hat{{\omega}}_\delta[Y-{\bf g}(\beta)], $然后极小化相应目标函数可得$ \tilde{\beta} $.

定理2.3  设条件$ \rm C1–C10 $成立,在用$ \rm logistic $模型估计缺失概率的情况下,若$ \beta $是参数真值,则$ \sqrt{n}(\tilde{\beta}-\beta)\stackrel{D}\longrightarrow N( {\bf 0}, \Sigma^{-1}\Lambda\Sigma^{-1}). $

定理2.4  设条件$ \rm C1–C10 $成立,在用$ \rm logistic $模型估计缺失概率的情况下,若$ a(u) $是系数函数真值,则$ \sqrt{nh}[\tilde{a}(u)-a(u)-2^{-1}h^2\mu_2 a''(u)]\stackrel{D}\longrightarrow N({\bf 0}, \nu_0f(u)^{-1}\Sigma^{*}). $

3 模拟研究

我们通过数值模拟来验证上一节提出的估计方法,基于模型$ 2.1 $来产生数据,其中$ \beta = (\beta_1, \beta_2)^{T} = (1, 1.5)^{T} $为真值, $ a_{1}(U) = \sin(2\pi U) $, $ a_{2}(U) = 3.5[\exp(-(4U-1)^2)+\exp(-(4U-3)^2)]-1.5 $, $ g(Z_1, Z_2;\beta_1, \beta_2) = \exp(\beta_1Z_1+\beta_2Z_2) $, $ X_1\sim N(1, 1) $, $ X_2\sim U(0, 3) $, $ \varepsilon\sim N(0, 0.25) $, $ Z_1\sim N(1, 1), Z_2\sim N(0, 1) $,测量误差$ e $服从均值为0,方差为$ \sigma_e^{2} $的正态分布,这里取两种测量误差,即$ \sigma_e^2 = 0.01, 0.25 $,重复模拟500次,分别取样本容量100, 200, 400.平均选择概率选为0.5,即

模拟中选取高斯核函数$ K(u) = \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}} $,采用$ \rm GCV $方法选取带宽.由表 1可以得到结论:缺失概率固定时,随着样本量的增加,两种参数估计得到的均值都越来越接近于真实值,它们的SD和MSE均在减小.另一方面,测量误差大小影响估计精度, $ \sigma_e^{2} $越大估计效果越差.

表 1   基于完全情况数据下参数$\beta$的估计结果

$\sigma_e^{2} = 0.1^{2}$$\sigma_e^{2} = 0.5^{2}$
$n$$\beta$BiasSDMSEBiasSDMSE
$100$$\hat{\beta}_1$-0.00160.01141.334e-04-0.01700.10450.0112
$\hat{\beta}_2$ 0.00080.01462.148e-04 0.01410.10530.0113
$200$$\hat{\beta}_1$ 0.00010.00502.546e-05-0.00050.00512.697e-05
$\hat{\beta}_2$-0.00040.00664.418e-05-0.00080.00816.685e-05
$400$$\hat{\beta}_1$-0.00010.00224.916e-06-0.00040.00266.698e-06
$\hat{\beta}_2$-0.00010.00277.188e-06-0.00010.00256.395e-06

新窗口打开| 下载CSV


表 2中将利用逆概率加权局部纠偏法得到的估计记为$ \hat{\beta}^{HT} $,忽略测量误差仅用逆概率加权法得到的估计记为$ \hat{\beta}^{NE} $.从表中结果可以进一步看出,在其他条件一样的情况下,局部纠偏得到的估计结果优于直接忽略测量误差得到$ \beta $的估计结果.两表联系起来观察可以发现,综合局部偏差纠正和逆概率加权方法要优于只针对完全数据进行局部纠偏或者利用逆概率加权法却忽略测量误差这两种方法得到的估计结果.说明了本文提出方法的有效性.通过图 1可以更加直观地看出本文提出的逆概率加权局部纠偏profile最小二乘法良好的拟合结果.因为采用的是估计的均值来拟合的曲线,所以图中结果显得更为光滑.

表 2   基于逆概率加权法考虑测量误差和忽略测量误差情况下得到的参数估计的表现

$\hat{\beta}^{HT}$$\hat{\beta}^{NE}$
$\beta$$\sigma_e^{2}$$n$BiasSDMSEBiasSDMSE
$\beta_1$$0.1^{2}$$100$-0.00130.00918.558e-050.01050.01423.112e-04
$200$-0.00070.00452.117e-050.00680.00577.865e-05
$400$0.00030.00162.614e-06-0.00080.00278.244e-06
$0.5^{2}$$100$0.00190.00939.010e-050.01140.01473.461e-04
$ 200$-0.00080.00522.843e-050.00620.00658.049e-05
$400$0.00020.00183.443e-060.00480.00514.916e-05
$\beta_2$$0.1^{2}$$100$0.00010.01161.345e-040.00220.01552.451e-04
$200$0.00030.00654.240e-05-0.00090.00715.114e-05
$400$-0.00040.00193.881e-060.00030.00391.557e-05
$0.5^{2}$$100$-0.00160.01381.930e-040.00260.02747.575e-04
$200$-0.00070.00654.382e-050.00140.01121.271e-04
$400$-0.00010.00267.035e-060.00040.00482.353e-05

新窗口打开| 下载CSV


图 1

图 1   $n = 400$时的模拟结果.在每个图中,实线为真实函数曲线,点线为局部纠偏的非参数函数估计曲线,虚线为忽略测量误差的非参数函数估计曲线


4 实例分析

这里选取一个心脏病数据来验证提出的方法在实际中的可用性.在此次实际数据分析中选取最大心率作为响应变量($ Y $),受检者年龄($ \rm age $)、性别($ \rm sex $)、胆固醇($ \rm mg/dl $)、静止血压($ \rm resting blood pressure $)作为与其相关的协变量,其中$ X_1 $代表血压, $ X_2 $代表胆固醇, $ U $代表年龄, $ Z $代表性别.由Qian[12],假设$ Z, \beta $有如下关系: $ g(Z, \beta) = \frac{\beta_1Z}{\beta_2+Z} $,为更好的审视估计结果,我们认为真实数据是不带测量误差的,利用普通profile最小二乘得到的估计曲线作为基准曲线,然后假定协变量$ X = (X_1, X_2) $带有测量误差$ e_1, e_2 $,再局部纠偏得到系数函数估计曲线,根据查阅相关资料,假定$ e_1\sim N(0, 3^2), e_2\sim N(0, 1^2) $,将带有测量误差下利用局部偏差纠正方法和仅基于profile最小二乘的得到的估计曲线与基准曲线进行对比.核函数和带宽选取同数值模拟中方法.由于实际数据量较少,平均缺失概率仅选取为0.02,采用部分非线性变系数模型$ Y_i = \sum\limits^{2}_{i = 1}a_i(U)X_i+\frac{\beta_1Z_i}{\beta_2+Z_i}+\varepsilon_i $来拟合给定的数据,参数估计值$ (\hat{\beta}_1, \hat{\beta}_2) = (-0.589, 0.991). $图 2为非参数系数函数$ a_1(u), a_2(u) $的估计曲线, $ a_1(u) $为血压协变量的系数, $ a_2(u) $为胆固醇协变量的系数.一方面通过计算可决系数$ R^2 $的值,未纠偏时为$ 0.456 $,纠偏后为$ 0.692 $,说明拟合效果有所提升.另一方面,由图同样可以看出局部偏差纠正的profile最小二乘估计接近基于真实数据得到的基准估计.特别是受检者的年龄较小时,在有测量误差的情况下,仅依赖普通的profile最小二乘估计得到的估计曲线与局部纠偏方法得到的估计结果差距明显,说明本文所提出的局部偏差纠正方法是有效的.在青年人群和老年人群中,胆固醇和血压与最大心率之间的关系比较密切,而中年人群则与此两者关系并不明显,且年龄越小或年龄越大两者对于心脏健康的影响越为明显.由图中可以看到纵坐标的值还是比较小的,所以未来可以进一步研究非参数系数函数的检验问题,从而更加精准地判断其显著性.

图 2

图 2   基于局部纠偏法的实际数据非参数函数的估计曲线


5 定理的证明

引理5.1  设$ (X_1, Y_1), (X_2, Y_2), \cdots, (X_n, Y_n) $是独立同分布的随机向量,其中$ f(x, y) $表示$ (X, Y) $联合概率密度函数, $ K(\cdot) $是满足Lipschitz条件且具有有界支撑的有界正函数,给定$ \varepsilon<1-s^{-1} $,使得$ n^{2\varepsilon-1}-h\rightarrow\infty, $$ E|Y|^{s}<+\infty, \sup\limits_x\int |y|^{s}f(x, y){\rm d}y<\infty $,则有

$ \begin{equation} \sup\limits_{x}\left|n^{-1}\sum^n\limits_{i = 1}\{K_h(X_i-x)Y_i-E[K_h(X_i-x)Y_i]\}\right| = O_p\left(\left\{\frac{\rm {log} (1/h)}{nh}\right\}^{1/2}\right). \end{equation} $

  证明详见文献[13].

引理5.2  设$ T_1, \cdots, T_n $为独立同分布的随机变量,若$ s>1 $, $ E|T_i|^s $有界,则$ \max\limits_{1\leq i\leq n}|T_i| = o(n^{1/s}), $ a.s..

  证明详见文献[14].

引理5.3  在条件$ \sup\limits_{i\geq1}\|\tau_i\|<\infty $$ \lambda_{\min}\rightarrow\infty $下,其中$ \tau_i = (1, X_i^{T}, U_i, Z_i^{T}) $, $ A = E(\tau_1\tau_1^{T}\pi_1 (1-\pi_1)), $$ \lambda_{\min} $表示$ \sum\limits^n_{i = 1}\tau_i\tau_i^{T} $最小特征根,缺失机制中参数$ \theta $的拟似然估计$ \hat{\theta} = (\hat{\theta}_0, \hat{\theta}_1, \hat{\theta}_2, \hat{\theta}_3)^{T} $

$ \begin{equation} \sqrt{n}(\hat{\theta}-\theta) = A^{-1}n^{-\frac{1}{2}}\sum^n\limits_{i = 1}\tau_i(\delta_i-\pi_i)+o_p(1). \end{equation} $

  证明详见文献[15].

引理5.4  设$ {\bf g}'(\beta) = (g'(Z_1, \beta), g'(Z_2, \beta), \cdots, g'(Z_n, \beta))^{T}) $,则在条件C1–C9下,有

$ \begin{equation} \{(D^{W}_{u})^{T}{\omega}_\delta D^{W}_{u}-{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}{\omega}_\delta {\bf g}' (\beta) = \Gamma^{-1}(u)\Phi(u)\otimes(1, 0)^{T}+O_p(C_n), \end{equation} $

$ \begin{equation} \{(D^{W}_{u})^{T}{\omega}_\delta D^{W}_{u}-{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}{\omega}_\delta M^X = a(u)\otimes(1, 0)^{T}+O_p(C_n), \end{equation} $

$ \begin{equation} \{(D^{W}_{u})^{T}{\omega}_\delta D^{W}_{u}-{\Omega}(u)\}^{-1}(D^{W}_{u})^{T}{\omega}_\delta\varepsilon = (1, 0)^{T}O_p(C_n). \end{equation} $

  由于三式的证明类似,所以这里仅给出(5.3)式的证明.简单计算可得

利用引理5.1,条件期望的平滑性并结合文献[11,定理2]的证明过程,易证(5.3)式成立.

引理5.5  在条件$ \rm C1–C9 $下,有

$ \begin{equation} n^{-1}{\bf g}'(\beta)^{T}[(I_n-S_h)^{T}\Delta(I_n-S_h)-Q_h^{T}\Delta\otimes\Sigma_eQ_h]{\bf g}'(\beta) = \Sigma+o_p(1), \end{equation} $

$ \begin{equation} n^{-\frac{1}{2}}{\bf g}'(\beta)^{T}[(I_n-S_h)^{T}\Delta(I_n-S_h)-Q_h^{T}\Delta\otimes\Sigma_e Q_h](Y-g(\beta)) = \sqrt{n}\xi_n+o_p(1), \end{equation} $

其中

  结合引理(5.1)–(5.3)并利用大数定律,类似Li和Mei[16]可证(5.6)式成立.下证(5.7)式.利用引理5.1,引理5.3和引理5.4简单计算可得

引理5.5证毕.

引理5.6  设C1, C5, C6, C10条件成立,则当$ n\rightarrow\infty $时,有

$ \begin{equation} \max\limits_{1\leq i\leq n} \Big|\frac{\delta_i}{\pi(V_i, \hat{\theta})}- \frac{\delta_i}{\pi(V_i, \theta)}\Big| = o_p(n^{-\frac{1}{2}+\frac{1}{2s}}). \end{equation} $

  根据引理$ \rm 5.3 $,将$ \frac{1}{\pi(V_i, \hat{\theta})} $在点$ \theta $处一阶Taylor展开得

结合条件$ \rm C10 $和引理$ \rm 5.2 $

引理5.6证毕.

引理5.7  在条件$ \rm C1–C10 $成立的情况下,有

$ \begin{equation} n^{-1}(\hat{\Omega}(u)-\Omega(u)) = o_p(n^{-\frac{1}{2}+\frac{1}{2s}}) \end{equation} $

  结合引理$ \rm 5.1 $,引理$ \rm 5.6 $和条件$ \rm C3 $简单运算可得,此处省略.

引理5.8  在条件$ \rm C1–C10 $成立的情况下,有

$ \begin{equation} n^{-1}\{[(D^{W}_{u})^{T}\hat{\omega}_\delta D^{W}_{u}-\hat{\Omega}(u)]-[(D^{W}_{u})^{T}\omega_\delta D^{W}_{u}-\Omega(u)]\} = o_p(n^{-\frac{1}{2}+\frac{1}{2s}}). \end{equation} $

  结合与引理$ \rm 5.1 $,引理$ \rm 5.4 $,引理$ \rm 5.6 $类似的证明方法和条件$ \rm C1–C6 $,并利用引理$ \rm 5.7 $的结果得证.

引理5.9  在条件$ \rm C1–C10 $成立的情况下,有

$ \begin{equation} \hat{S}_h-S_h = o_p(n^{-\frac{1}{2}+\frac{1}{2s}}). \end{equation} $

  结合条件C3, C5,引理$ \rm 5.6 $并利用同引理$ \rm 5.7 $和引理$ \rm 5.4 $相似的方法可以证明本引理成立.

引理5.10  设$ {\bf g}'(\beta) = (g'(Z_1, \beta), \cdots , g'(Z_n, \beta))^{T}) $,在$ \rm C1–C10 $条件成立下有

$ \begin{equation} n^{-1}{\bf g}'(\beta)^{T}[(I_n-\hat{S}_h)^{T}\hat{\Delta}(I_n-\hat{S}_h)-\hat{Q}_h^{T}\hat{\Delta}\otimes\Sigma_e\hat{Q}_h]{\bf g}'(\beta) = \Sigma+o_p(1), \end{equation} $

$ \begin{equation} n^{-\frac{1}{2}}{\bf g}'(\beta)^{T}[(I_n-\hat{S}_h)^{T}\hat{\Delta}(I_n-\hat{S}_h)-{\hat{Q}_h}^{T}\hat{\Delta}\otimes\Sigma_e \hat{Q}_h](Y-g(\beta)) = \sqrt{n}\xi_n+o_p(1). \end{equation} $

  将(5.12)式拆解成2项,简单计算可得

下面对$ R_1 $$ R_2 $部分分别进行展开计算.首先考虑

下证$ R_{1i} (i\neq 5) = o_p(1) $.$ R_{11}, R_{12} $为例.记$ A_{r\times n} = [(I_n-S_h){\bf g}'(\beta)]^{T} $,由于

因此,根据引理5.2、5.6,条件C3、C7、C8有

同理设$ B_{n\times r} = (S_h-\hat{S}_h){\bf g}'(\beta) $,根据引理5.2、5.9和条件C3、C7、C8容易知

类似计算$ R_{1i}\ (i = 3, 4, 6, 7, 8) = o_p(1) $.其次考虑

下证$ R_{2i}\ (i\neq5) = o_p(1) $.$ R_{21} $为例,由引理5.2和引理5.6和条件C3、C7、C8以及矩阵范数的相容性可以得到

类似计算$ R_{2i}\ (i = 2, 3, 4, 6, 7, 8) = o_p(1) $.由引理$ 5.5 $$ R_{15}-R_{25} = \Sigma+o_p(1) $,从而可得

类似方法将(5.13)式拆解为4项,即

接着将$ T_{11} $$ T_{12} $$ T_{21} $$ T_{22} $分别展开计算.首先考虑

下面证明$ T_{11i}\ (i\neq5) = o_p(1) $.$ T_{111} $为例,记$ G_{n\times 1} = (I_n-S_h)M^{X} $,那么有

取上面矩阵中的任意一个元素来分析,由条件C2–C8知$ E(A^{2}_{ij}G^{2}_{j}) $有界设为$ C $.再结合引理5.6,有

类似可以计算$ T_{11i}(i = 2, 3, 4, 6, 7, 8) = o_p(1) $.其次考虑

下面证明$ T_{12i}(i\neq5) = o_p(1) $,以$ T_{121} $为例,结合条件C2–C8和引理5.2,引理5.6可得

类似可以计算$ T_{12i}\ (i = 2, 3, 4, 6, 7, 8) = o_p(1) $.按照类似的思路计算$ T_{21}, T_{22} $可得$ T_{2ji}\ (j = 1, 2; $$ i = 2, 3, 4, 6, 7, 8) = o_p(1) $.又由引理$ 5.5 $$ T_{115}-T_{125}+T_{215}-T_{225} = \sqrt{n}\xi_n+o_p(1) $,可得

引理5.10得证.

定理2.1的证明  结合$ R(\beta) = (Y-g(\beta))^{T}[(I_n-S_h)^{T}\Delta(I_n-S_h)-Q^{T}_h\Delta\otimes\Sigma_eQ_h](Y-g(\beta)), $考虑对$ R(\beta+t) $$ \beta $泰勒展开,有

$ \begin{equation} R(\beta+t) = R(\beta)+R'(\beta)^{T}t+\frac{1}{2}t^{T}R''(\beta^{*})t, \end{equation} $

其中$ \beta^{*}\in(\beta, \beta+t) $$ R'(\beta) $$ R(\beta) $$ \beta $处的一阶导,利用引理5.5中的(5.7)式的计算过程,简单计算可得

$ G(\beta) = t^{T}{\bf g}''(\beta)t $.因为

再结合引理5.5和条件C7, C8,由于$ M^{X}, \varepsilon $不含变量$ \beta $,可以看做是有限常数,从而只考虑$ d $这部分即可.由$ \beta_{0}\in(\beta, \beta^{*}) $,利用中值定理可得

由于$ \|t\| $可以任意小,所以$ t^{T}R''(\beta^{*})t $的大小和正负由$ nt^{T}\Sigma\{1+o_p(1)\}t $决定,所以有上式成立.又因为$ \Sigma $是正定矩阵,所以$ t^{T}R''(\beta^{*})t>0 $.由前面的计算可知$ R'(\beta)^{T}t = O_p(\sqrt{n}\|t\|) $,即$ R'(\beta)^{T}t $至少与$ \sqrt{n}\|t\| $同阶,而$ t^{T}R''(\beta^{*})t\triangleq2nt^{T}\Sigma\{1+o_p(1)\}t+O_p(||t||^3) $至少与$ (\sqrt{n}\|t\|)^{2} $同阶,进而$ R(\beta+t)-R(\beta) = R'(\beta)^{T}t+\frac{1}{2}t^{T}R''(\beta^{*})t $,其大小和正负仍然由$ \frac{1}{2}t^{T}R''(\beta^{*})t $决定.故$ R(\beta+t)>R(\beta) $$ \|t\| = a $的邻域里几乎处处成立,换句话说$ R(\beta) $的局部极小值存在,所以$ R'(\hat{\beta}) = 0 $.在具备以上结果之后,就可以进一步证明$ \hat{\beta} $的渐近正态性了.首先利用中值定理有

$ \begin{equation} R'(\hat{\beta}) = R'(\beta)+R''(\beta^{*})^{T}(\hat{\beta}-\beta) = 0, \end{equation} $

其中$ \beta^{*} $满足$ |\beta^{*}-\beta|\leq |\hat{\beta}-\beta| $.又因为$ n^{-1}R'(\beta) = -2(\xi_n+O_p(C_n)O_p(n^{-\frac{1}{2}})). $联系(5.15)式,简单运算可以得到$ \Sigma\{1+o_p(1)\}(\hat{\beta}-\beta) = \sqrt{n}\xi_n+O_p(C_n) = \sqrt{n}\xi_n+o_p(1). $因为$ \xi_n $是独立随机变量和的形式,由中心极限定理知其渐近期望为0,渐近方差为

根据Slutsky定理有$ \sqrt{n}(\hat\beta-\beta)\stackrel{D}\longrightarrow\; N( {\bf 0}, \Sigma^{-1}\Lambda\Sigma^{-1}). $

定理2.2的证明  由$ Y-g(\hat{\beta}) = M^{X}+\varepsilon-g(\hat{\beta})+g(\beta) $,简单计算得

首先考虑$ I_3 $,根据定理2.1得到的$ \hat{\beta} $$ \sqrt{n} $相合性结果并结合泰勒展开方法,有$ I_3 = O_p(n^{-\frac{1}{2}}). $然后考虑$ I_1 $,仍然利用泰勒展开,可得

将上式代入$ I_1 $,利用引理5.4的相关结果可以推出

最后考虑$ I_{2} $,首先讨论$ (D^{W}_{u})^{T}\omega_\delta\varepsilon $.经计算可得其渐近方差为

其次计算

$ \eta = \Sigma_e-ee^{T}-Xe^{T} $,用类似于$ I_{1}, I_{2} $方法并结合中心极限定理有

然后计算$ P_{22} $,通过换元积分结合引理5.1和条件C3–C6,可得$ \frac{\sqrt{nh}}{n}P_{22} = o_p(1). $将上面结果代回到$ \hat{\psi}(u) $的表达式,有

利用Slutsky整理可得

其中$ \psi(u) = (a(u)^{T}, h a'(u)^{T})^{T}, \eta = \Sigma_e-ee^{T}-Xe^{T}, \Sigma^{*} = \Gamma^{-1}(u)[\sigma^{2}(E({\pi(V)}^{-1}XX^{T}|U = u)+E({\pi(V)}^{-1}ee^{T}|U = u))+E({\pi(V)}^{-1}\eta a(u)a(u)^{T}\eta^{T}|U = u)]\Gamma^{-1}(u). $进而有

定理2.2得证.

定理2.3和定理2.4的证明  应用定理2.1和2.2和引理5.3–5.10可以将缺失概率未知时参数$ \beta $和非参数函数$ a(U) $的估计写为

然后根据Slutsky定理,便可证明定理2.3和定理2.4.此处不再做详细证明过程.

参考文献

Fan J Q , Zhang W Y .

Statistical methods with varying coefficient models

Statistics and Its Interface, 2008, 1 (1): 179- 195

DOI:10.4310/SII.2008.v1.n1.a15      [本文引用: 1]

Li R Z , Nie L .

Efficient statistical inference procedures for partially nonlinear models and their applications

Biometrics, 2008, 64 (3): 904- 911

DOI:10.1111/j.1541-0420.2007.00937.x      [本文引用: 1]

Bates D M, Watts D G. Nonlinear Regression Analysis and Its Applications. New York: Wiley, 1988

[本文引用: 1]

Zhou Y , Wan A T K , Wang X .

Estimating equations inference with missing data

Journal of the American Statistical Association, 2008, 103, 1187- 1199

DOI:10.1198/016214508000000535      [本文引用: 1]

Xue L G .

Empirical likelihood for linear models with missing responses

Journal of Multivariate Analysis, 2009, 100 (7): 1353- 1366

DOI:10.1016/j.jmva.2008.12.009      [本文引用: 1]

Carroll R J , Ruppert D , Stefanski L A .

Measurement error in nonlinear models

Metrika, 1997, 45 (3): 182- 183

URL     [本文引用: 1]

冯三营, 牛惠芳.

部分线性变系数EV模型估计的渐近正态性

河南科技大学学报, 2011, 32 (2): 83- 87

DOI:10.3969/j.issn.1672-6871.2011.02.021      [本文引用: 3]

Feng S Y , Niu H F .

Asymptotic normality of estimators for partially linear varying coefficient errors-in-variables Models

Journal of Henan University of Science and Technology, 2011, 32 (2): 83- 87

DOI:10.3969/j.issn.1672-6871.2011.02.021      [本文引用: 3]

Liang H , Wang S , Carroll R J .

Partially linear models with missing response variables and error-prone covariates

Biometrika, 2007, 94 (1): 185- 198

URL     [本文引用: 1]

Wei C H .

Estimation in varying-coefficient errors-in-variables models with missing response variables

Communications in Statistics-Simulation and Computation, 2011, 40 (3): 383- 393

DOI:10.1080/03610918.2010.542846      [本文引用: 1]

Fan J Q, Gijbels I. Local Polynomial Modeling and Its Applications. London: Chapman and Hall, 1996

[本文引用: 1]

You J , Zhou Y , Chen G .

Corrected local polynomial estimation in varying-coefficient models with measurement errors

The Canadian Journal of Statistics, 2006, 34 (3): 391- 410

DOI:10.1002/cjs.5550340303      [本文引用: 3]

Qian Y Y .

Statistical inference for a varying-coefficient partially nonlinear model with measurement errors

Statistical Methodology, 2016, 32, 122- 130

DOI:10.1016/j.stamet.2016.05.004      [本文引用: 1]

Mark Y P , Silverman B W .

Weak and strong uniform consistency of kernel regression estimates

Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 1982, 61 (3): 405- 415

DOI:10.1007/BF00539840      [本文引用: 1]

Shi J , Lau T S .

Empirical likelihood for partially linear models

Journal of Multivariate Analysis, 2000, 72 (1): 132- 148

DOI:10.1006/jmva.1999.1866      [本文引用: 1]

李志强, 薛留根.

协变量随机缺失下的广义半参数模型

北京工业大学学报, 2007, 7, 761- 765

URL     [本文引用: 2]

Li Z Q , Xue L G .

The generalized semiparatric models with missing covariates

Journal of Beijing University of Technology, 2007, 7, 761- 765

URL     [本文引用: 2]

Li T Z , Mei C L .

Estimation and inference for varying coefficient partially nonlinear models

Journal of the Statistical Planning and Inference, 2013, 143, 2023- 2037

DOI:10.1016/j.jspi.2013.05.011      [本文引用: 1]

Fan J Q , Huang T .

Profile likelihood inferences on semiparametric varying-coefficient partially linear models

Bernoulli, 2005, 11 (6): 1031- 1057

DOI:10.3150/bj/1137421639     

Zhou Y , Liang H .

Statistical inference for semiparametric varying-coefficient partially linear models with error-prone linear covariates

The Annals of Statistics, 2009, 37 (1): 427- 458

DOI:10.1214/07-AOS561     

/