南京信息工程大学数学与统计学院，江苏 南京

收稿日期：2021年11月23日；录用日期：2021年12月17日；发布日期：2021年12月24日

摘要

受实验观测时长或观测条件导致的数据双侧截断在生存分析等领域十分常见，故而对该类数据在不同适用模型下的讨论具有研究意义。考虑双侧截断数据下的部分线性模型，针对观测到数据的逆概率加权，利用核估计对模型线性部分参数给出估计，并给出估计量的偏差与方差的渐进表达式，与常见的模型参数估计方法相比估计值偏差与标准差有一定程度的提升且更准确。通过数值模拟和AIDS病毒感染年龄与潜伏期的相关数据，验证了该方法的实用性。

关键词

双侧截断，部分线性模型，加权核估计

Parametric Estimation for Partial Linear Models with Double Truncated Data

Zhezi Shen, Peng Lai^*

School of Mathematics and Statistics, Nanjing University of Information Science & Technology, Nanjing Jiangsu

Received: Nov. 23^rd, 2021; accepted: Dec. 17^th, 2021; published: Dec. 24^th, 2021

ABSTRACT

Double truncated data due to the limit of observing time and other objective reasons of the experiments are quite common in the field of survival analysis. So the discussion of this type of data under different models is of great significance. Considering the partial linear model under the two-sided truncated data, based on the inverse probability weighting of the observed data, kernel estimation is used to estimate the linear parameter. Asymptotic expressions for the bias and variance of the estimator are given. Compared with the common method, the bias and variance are improved. The real data analysis of AIDS virus infection age and incubation proves the practicability of the proposed method.

Keywords:Double Truncation, Partial Linear Model, Weighted Kernel Estimation

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

双侧截断是一种常出现在生存分析和工业检测上的现象，双侧截断数据指的是数据的两侧存在由于随机出现的截断点导致响应变量无法被观测，从而导致样本点缺失的数据集。即当且仅当响应变量Y落在截断区间 $\left(L,R\right)$ 内时，该样本才可以被观测记录。由于响应变量缺失导致的对应协变量的缺失并无规律可循，而忽视这些有一定占比的缺失数据可能导致估计回归模型时参数有偏，且拟合效果会随着截断概率的增加而下降，故而如何在数据存在双侧截断的条件下更好地估计回归模型是一个有实际意义的课题。双侧截断数据的研究最早在天体观测和生存分析领域展开，从参数法与非参数法出发，得出了一些有效的估计方法。Efron和Petrosian (1999) [1] 针对双侧截断数据的生存分布函数提出了非参数最大似然估计(NPMLE)，Shen (2010) [2] 进而研究了NPMLE的渐近特性，并引入了截断分布函数的非参数估计量。Moreira和de Uña-Álvarez (2010) [3] 证明了bootstrap重采样可用于近似NPMLE的抽样分布。当假设的参数模型正确时，参数法的估计效果要显著强于非参数方法。故而近年来对双侧截断条件下的半参数模型的讨论居多，Shen (2013) [4] 针对半参数变换模型利用EM算法给出回归参数估计。进一步的，Shen (2017) [5]，Mandel et al. (2018) [6] 提出逆概率加权估计法拟合Cox回归模型。Rennert和Xie (2018) [7] 提出了双截断数据的Cox回归模型，并使用加权估计方程方法来估计回归参数。Ying et al. (2019) [8] 提出扩展的Mann-Whitney型损失函数，在响应变量双侧截断时估计线性回归模型中的回归参数。目前已有的研究多集中在生存分析领域常用模型在双侧截断数据情况下的扩展，但对于一般的常见模型的讨论较少。

部分线性模型是最初由Engle (1986) [9] 提出的半参数模型，综合了线性回归的可解释性和非参数回归的稳健性，有着极广的应用范围。对部分线性回归的估计方法大致有两类：一是回归样条，即使用样条函数逼近模型中的非参数部分，再利用最小二乘法同时估计线性部分参数及非参数部分的样条函数；二是Backfitting估计，在给定线性部分参数初值条件下估计非参数部分，再通过极小化估计标准误求得线性部分参数。然而对于双侧截断数据而言，已有的估计方法都会因为响应变量的截断导致估计结果有偏。由于双侧截断数据常出现在适用于部分线性模型常见的生存分析，工业观测等领域，故而引入适用该类数据的部分线性模型回归方法是有实际价值的。

2. 基于双侧截断的部分线性回归模型的参数估计与拟合

本文针对响应变量双侧截断的部分线性回归模型提出一种参数估计方法。

在不考虑数据截断情况下假设全部样本为： $\tilde{Z}=\left\{{\tilde{X}}_i,{\tilde{T}}_i,{\tilde{Y}}_i,L_i,R_i,i=1,\cdots ,N\right\}$，针对部分线性模型：

$\tilde{Y}={\tilde{X}}^{\text{T}}\beta +g\left(\tilde{T}\right)+\tilde{\epsilon }$. (1)

对(1)式两侧求期望后直接计算可得：

$\tilde{Y}-E\left(\tilde{Y}|T\right)={\beta }^{\text{T}}\left[\tilde{X}-E\left(\tilde{X}|T\right)\right]+\tilde{\epsilon }$. (2)

记 $E\left(\tilde{Y}|T\right)=m_{\tilde{Y}}\left(T\right)$，$E\left(\tilde{X}|T\right)=m_{\tilde{X}}\left(T\right)$，并且分别将其估计值记为 ${\hat{m}}_{\tilde{Y}}\left(T\right)$，${\hat{m}}_{\tilde{X}}\left(T\right)$。

由于响应变量Y存在双侧截断，当且仅当满足 $L_i<{\tilde{Y}}_i<R_i$ 时协变量及响应变量可被观测，记被观测到的样本为 $Z=\left\{X_i,T_i,Y_i,U_i,V_i,i=1,\cdots ,n\right\}$，其中 $U,V$ 为观测到的样本的双侧截断点，满足 $U_i\le Y_i\le V_i$。截断后的数据满足：

$Y={\beta }^{\text{T}}X+g\left(T\right)+\epsilon$. (3)

考虑截断后 $E\left(Y|T\right)$ 和 $E\left(X|T\right)$ 的非参数估计。记 $E\left(Y|T\right)=m_Y\left(T\right)$，$E\left(X|T\right)=m_X\left(T\right)$。

记 $F_1\left(\cdot |t\right)$ 为给定 $T=t$ 时 $\tilde{Y}$ 的条件分布函数， $F_2\left(\cdot |t\right)$ 为给定 $T=t$ 时 $\tilde{X}$ 的条件分布函数，则有 $m_{\tilde{Y}}\left(t\right)={\displaystyle {\int }_{-\infty }^{+\infty }s\text{d}{F}_1\left(s|t\right)}$，$m_{\tilde{X}}\left(t\right)={\displaystyle {\int }_{-\infty }^{+\infty }s\text{d}{F}_2\left(s|t\right)}$。给定 $T=t$ 时 $\tilde{Y}$ 落在未被截断区间的概率记为 $\alpha \left(t\right)$，则有 $\alpha \left(t\right)=P\left(L\le \tilde{Y}\le R|T=t\right)={\displaystyle {\int }_{-\infty }^{+\infty }H\left(s\right)\text{d}{F}_1\left(s|t\right)}$，代表无截断情况下响应变量可以被观测到的条件概率，其中响应变量可被观测到的概率记为 $H(\cdot )$，有 $H\left(y\right)=P\left(L\le y\le R\right)$。参照Moreira C.等(2014) [10] 推得

截断后响应变量的条件分布函数 $F\left(y|t\right)=\alpha {\left(t\right)}^{-1}{\displaystyle {\int }_{-\infty }^{y}H\left(s\right)\text{d}F\left(s|t\right)}$。再记 ${\alpha }^{*}\left(t\right)=\alpha {\left(t\right)}^{-1}$，则有 ${\alpha }^{*}\left(t\right)={\displaystyle {\int }_{-\infty }^{+\infty }H{\left(s\right)}^{-1}\text{d}F\left(s|t\right)}=E\left(H{\left(Y\right)}^{-1}|T=t\right)$。

记 $m_Y^{\ast }\left(t\right)=E\left(YH{\left(Y\right)}^{-1}|T=t\right)$，故而 $m_Y\left(t\right)$ 可以表示为： $m_Y\left(t\right)=\frac{m_Y^{\ast }\left(t\right)}{{\alpha }_Y^{\ast }\left(t\right)}$。由于 $m^{*}\left(t\right)$ 与 ${\alpha }^{*}\left(t\right)$ 均为观测变量的条件均值，故二者可以通过NW核估计给出估计值。

在双侧截断条件下二者的核估计值分别为：

${\hat{m}}_X^{\ast }\left(T\right)=\frac{{\displaystyle {\sum }_{i=1}^{n}K\left(\frac{T_i-T}{h}\right)X_i{H}_n\left(Y_i^{-1}\right)}}{{\displaystyle {\sum }_{i=1}^{\tilde{n}}K\left(\frac{T_i-T}{h}\right)}}$，${\hat{m}}_Y^{\ast }\left(T\right)=\frac{{\displaystyle {\sum }_{i=1}^{n}K\left(\frac{T_i-T}{h}\right)Y_i{H}_n\left(Y_i^{-1}\right)}}{{\displaystyle {\sum }_{i=1}^{\tilde{n}}K\left(\frac{T_i-T}{h}\right)}}$ (4)

${\hat{\alpha }}_X^{\ast }\left(T\right)=\frac{{\displaystyle {\sum }_{i=1}^{n}K\left(\frac{T_i-T}{h}\right)H_n\left(Y_i^{-1}\right)}}{{\displaystyle {\sum }_{i=1}^{\tilde{n}}K\left(\frac{T_i-T}{h}\right)}}$，${\hat{\alpha }}_Y^{\ast }\left(T\right)=\frac{{\displaystyle {\sum }_{i=1}^{n}K\left(\frac{T_i-T}{h}\right)H_n\left(Y_i^{-1}\right)}}{{\displaystyle {\sum }_{i=1}^{\tilde{n}}K\left(\frac{T_i-T}{h}\right)}}$ (5)

其中 $K(\cdot )$ 为核函数，h为窗宽， $H_n$ 为响应变量被观测到的概率H的非参数估计。

对于响应变量被观测到概率H的非参数估计 $H_n$，Shen (2010) [2] 提出 $H_n\left(y\right)={\displaystyle \int S_n\left(\text{d}L,\text{d}R\right)}$，其中 $S_n$ 为双侧截断点的联合分布函数的NPMLE，并给出了 $S_n$ 的迭代估计步骤：

Step 1. 设初始估计为 ${\varphi }_i^{\left(0\right)}=\stackrel{n}{\underset{m=1}{{{\displaystyle \sum }}^{\text{}}}}{\phi }_m^{\left(0\right)}I\left(L_i\le Y_m\le R_i\right)$，其中 ${\phi }_m^{\left(0\right)}=\frac{1}{n},1\le m\le n$ ；

Step 2. 重新对S给出估计：令 ${\psi }_j^{\left(1\right)}={\left[\stackrel{n}{\underset{i=1}{{{\displaystyle \sum }}^{\text{}}}}\frac{1}{{\varphi }_i^{\left(0\right)}}\right]}^{-1}\frac{1}{{\varphi }_j^{\left(0\right)}},1\le j\le n$ 并记其累加和为 ${\psi }_i^{\left(1\right)}=\stackrel{n}{\underset{m=1}{{{\displaystyle \sum }}^{\text{}}}}{\psi }_m^{\left(1\right)}I\left(L_m\le Y_i\le R_m\right)$ ；

Step 3. 更新对F的估计值：令 ${\phi }_j^{\left(1\right)}={\left[\stackrel{n}{\underset{i=1}{{{\displaystyle \sum }}^{\text{}}}}\frac{1}{{\psi }_i^{\left(1\right)}}\right]}^{-1}\frac{1}{{\psi }_j^{\left(1\right)}},1\le j\le n$ 并记其累加和为 ${\varphi }_i^{\left(1\right)}=\stackrel{n}{\underset{m=1}{{{\displaystyle \sum }}^{\text{}}}}{\phi }_m^{\left(1\right)}I\left(L_i\le Y_m\le R_i\right)$ ；

Step 4. 规定收敛条件为 $\underset{1\le j\le n}{\max }\left|{\phi }_j^{\left(k-1\right)}-{\phi }_j^{\left(k\right)}\right|\le {10}^{-6}$，重复以上两个迭代步骤直至达到收敛条件。

设在第k次迭代结束够满足收敛条件，则 $S_n$ 的估计值为：

$S_n\left(L,R\right)=\stackrel{n}{\underset{j=1}{{{\displaystyle \sum }}^{\text{}}}}{\psi }_j^{\left(k\right)}I\left(L_j\le L,R_j\le R\right)$. (6)

由此可知 $H_n\left(y\right)={\displaystyle \int S_n\left(\text{d}L,\text{d}R\right)}=\stackrel{n}{\underset{j=1}{{{\displaystyle \sum }}^{\text{}}}}{\psi }_j^{\left(k\right)}I\left(L_j\le y\le R_j\right)$。

至此可计算出 ${\hat{m}}_Y\left(T\right)$，${\hat{m}}_X\left(T\right)$ 分别为：

${\hat{m}}_X\left(T\right)=\frac{{\hat{m}}_X^{\ast }\left(T\right)}{{\hat{\alpha }}_X^{\ast }\left(T\right)}$，${\hat{m}}_Y\left(T\right)=\frac{{\hat{m}}_Y^{\ast }\left(T\right)}{{\hat{\alpha }}_Y^{\ast }\left(T\right)}$ (7)

进一步地，记 $\bar{Y}=Y-{\hat{m}}_Y\left(T\right)$，$\bar{X}=X-{\hat{m}}_X\left(T\right)$，则有：

$\bar{Y}={\beta }^{\text{T}}\bar{X}+\epsilon$. (8)

由于双侧截断的缺失机制，直接对回归系数给出估计存在较大偏差，故通过逆概率加权对其纠偏，利用最小二乘法对线性参数 $\beta$ 作出估计。

记 $Y_i^{*}=\frac{{\bar{Y}}_i}{H_n\left(Y_i\right)}$，可求得参数 $\beta$ 的估计量为：

$\hat{\beta }={\left({\displaystyle {\sum }_{i=1}^n{\bar{X}}_i{\bar{X}}_i^{\text{T}}}\right)}^{-1}{\displaystyle {\sum }_{i=1}^n{\bar{X}}_i{Y}_i^{\ast }}$. (9)

得到线性回归参数的估计量 $\hat{\beta }$ 之后，代入(3)式可以得到非线性部分的拟合值。

为简化表述，将下述 ${\hat{m}}_Y\left(t\right)$ 与 ${\hat{m}}_X\left(t\right)$ 统称为 $\hat{m}\left(t\right)$，给出 $\hat{m}\left(t\right)$ 的偏差和方差的渐近表达式。定义协变量 的密度函数为 $f\left(t\right)$，$\alpha =P\left(R\le Y^{*}\le L\right)$，$\alpha >0$ 为无条件下的观测概率，同时给定 $\mu \left(K\right)={\displaystyle \int t^{2}K\left(t\right)\text{d}t}$，$R\left(K\right)={\displaystyle \int K{\left(t\right)}^2\text{d}t}$。为简化表述令 $a_n~b_n$ 表示实数列 $a_n,b_n$ 在 $n\to \infty$ 时满足 $\frac{a_n}{b_n}\to 1$。引入：

${\sigma }^2\left(x\right)=E\left[{\left(\tilde{Y}-m\left(\tilde{T}\right)\right)}^{2}H{\left(\tilde{Y}\right)}^{-1}|T=t\right]$. (10)

给定如下假设：

(C1). $m\left(t\right)$ 的二阶导 $m^{″}$ 连续且在t的领域内有界；

(C2). 核函数 $K(\cdot )$ 为一关于零对称的密度函数且存在其紧支集；

(C3). $n\to \infty$ 的同时 $h\to 0$，且 $nh\to 0$ ；

(C4). 条件期望 $E\left[{\left(\tilde{Y}-m\left(\tilde{T}\right)\right)}^{d}H{\left(\tilde{Y}\right)}^{-1}|T=t\right],d=0,1,2$ 有界。

基于假设(C1)-(C4)，类似Moreira C.等(2014) [10] 证明过程可以给出 $\hat{m}\left(t\right)$ 的渐进偏差和方差：

$Bias\left(\hat{m}\left(t\right)\right)~\frac{1}{2}\mu \left(K\right)h^2\left[m^{″}\left(t\right)f\left(t\right)+2m^{\prime }\left(t\right)f^{\prime }\left(t\right)\right]/f\left(t\right)\equiv h^{2}B\left(t\right)$， (11)

$Var\left(\hat{m}\left(t\right)\right)~\frac{1}{nh}R\left(K\right)\alpha {\sigma }^2\left(t\right)/f\left(t\right)\equiv \frac{1}{nh}V\left(t\right)$. (12)

3. 数值模拟

本节模拟研究文中所提估计方法在有限样本下的实际表现，比较在响应变量双侧截断的条件下本文提出的估计方法与原始的部分线性模型回归方法的回归效果。对模拟给定的两个模型分别取截断概率为0.2，0.4，可被观测到的样本量分别取 $n=100,200,400$。Bootstrap重抽样数B = 1000。记忽略截断条件的参数估计量为 ${\hat{\beta }}_{PLM}$，本文提出的非参数估计法得到的估计量为 ${\hat{\beta }}_{NP}$。对于响应变量双侧截断点的选取分别令其服从均匀分布，正态分布进行模拟。具体参数设置如下：

模型一：考虑协变量均为一维变量，通过模型 $\tilde{Y}=3\tilde{X}+\sin \left(2\pi \tilde{T}\right)+\epsilon$ 产生观测数据，其中 $\tilde{X}~U\left(-1,1\right)$，$\tilde{T}~U\left(0,1\right)$，$\epsilon ~N\left(0,0.01\right)$，随机生成双侧截断点： $L~U\left(-4,-1\right)$，$R~U\left(1,4\right)$。

模型二：考虑多元协变量，通过模型 $\tilde{Y}={\beta }^{\text{T}}\tilde{X}+\cos \left(2\pi \tilde{T}\right)+\epsilon$ 产生观测数据，其中 $\tilde{X}=\left({\tilde{X}}_1,{\tilde{X}}_2,{\tilde{X}}_3\right)$，${\tilde{X}}_i~U\left(-1,1\right)$，$\left(i=1,2,3\right)$，$\tilde{T}~U\left(0,1\right)$，$\beta ={\left(1,2,-1\right)}^{\text{T}}$，$\epsilon ~N\left(0,0.01\right)$，随机生成双侧截断点： $L~N\left(-2,1\right)$，$R~N\left(2,1\right)$。

由表1与表2可知，在相同样本量以及相同截断概率的条件下，本文提出的估计方法相较于传统的估计方法有更小的偏差和较小的估计标准误差，且95%置信区间在真实值附近更加精确。由模型一的参数部分估计结果可以看出，在双侧截断数据条件下，已有的部分线性模型估计方法给出的线性部分参数估计结果并没有随着样本数的增加而逼近真实值，可见数据截断对于该方法的准确性造成一定影响。反观本文提出的非参数方法给出的估计值均为样本数越大，估计越准确。图1与图2分别展示了两个模型的非参数部分回归曲线与原始数据曲线(实线)的相近程度。横轴上的线段分布表示了变量的截断情况，线段越稀疏代表该处被截断从而无法观测到的数据越多。可以看出随着截断后样本数量n的增加以及样本被观测概率的增长，传统回归方法(虚线)对原始数据拟合效果提升。而本文提出的对非参数部分的加权核回归(点虚线)几乎与原始数据(实线)重合，效果明显优于传统方法。

表1. 模型一的参数估计

图1. 不同 $\left(n,p\right)$ 取值下两种方法对模型一非参数部分的拟合

表2. 模型二的参数估计

4. 实例分析

现将本文提出的方法应用于Kalbfeisch (1989)报告的AIDS孵化数据共295条。数据信息包含患者年龄(AGE)，感染AIDS病毒的时长(INF)以及潜伏期(ICUB)。研究目的为判断患者年龄是否会对其病毒潜伏期构成影响。该数据为双侧截断数据，左侧截断是由于1982年之前HIV疾病未知，右侧截断则是实验终止导致。响应变量ICUB的截断区间为 $\left(L,R\right)=\left(54-\left(\text{INF}-0.5\right),138-\left(\text{INF}-0.5\right)\right)$。

图2. 不同 $\left(n,p\right)$ 取值下两种方法对模型二非参数部分的拟合

根据患者年龄将数据划分为三类分别求其年龄与潜伏期的关系，分别划分为0~18为儿童，19~59为

成人，60及以上为老年人，考虑模型： $Y={\beta }^{\text{T}}X+g\left(T\right)+\epsilon$。类似Liu (2019)构造变量： $T=\frac{\text{INF}-E\left(\text{INF}\right)}{\sqrt{Var\left(\text{INF}\right)}}$。利用本文提出的双侧截断数据模型参数估计方法，给出不同年龄段的病毒潜伏期与年龄的线性回归参数分别为： ${\beta }_1=3.5$，${\beta }_2=3.11$，${\beta }_3=1.49$。可得结论：不同年龄段患者的潜伏期时长受其年龄影响，且随患者感染年龄的增长该影响逐渐减小。

5. 结论

本文基于双侧截断数据，使用逆概率加权核估计对部分线性模型给出线性部分参数估计，并给出估计量的偏差与方差的渐进表达式。从数值模拟部分的结果可以看出，该方法给出的估计值具有更小的偏差与更强的准确性，且后续对非线性部分的拟合效果有显著的提升，实例分析进一步证明了该方法的可行性。故而本文提出的针对双侧截断数据下部分线性模型的估计方法在生存分析等领域的适用数据具有实际价值。

基金项目

国家自然科学基金面上项目：超高维复杂数据统计降维研究(11771215)，2018.1~2021.12。

文章引用

沈哲子,来 鹏. 基于双侧截断的部分线性回归模型的参数估计
Parametric Estimation for Partial Linear Models with Double Truncated Data[J]. 应用数学进展, 2021, 10(12): 4365-4372. https://doi.org/10.12677/AAM.2021.1012464

参考文献