基于纵向数部分线性模型的分位数模型平均

Quantile Model Averaging Based on Longitudinal Partial Linear Model

摘要: 本文针对纵向数据部分线性回归模型中的参数估计与非参数估计问题，基于分位数回归估计方法提出了一种稳健的模型平均估计量。为了提高估计效率，采用工作相关矩阵分解和估计方程平滑法处理纵向数据的组内相关性，并通过局部线性估计方法处理模型的非参数部分，给出了模型参数与非参数估计的Newton-Raphson迭代算法。数值模拟表明，新的估计方法具有良好的估计性能。将新估计方法应用到空气质量数据的预测分析中，证明了该方法在实际应用中也具有可行性。

Abstract: Based on the problem of parameter estimation and non-parametric estimation in the partial linear regression model of longitudinal data, this paper proposes a robust model average estimator based on the quantile regression estimation method. In order to improve the estimation efficiency, the working correlation matrix decomposition and estimation equation smoothing method are used to deal with the intra-group correlation of longitudinal data, and the non-parametric part of the model is processed by the local linear estimation method. The Newton-Raphson iterative algorithm for model parameter and non-parametric estimation is given. Numerical simulation shows that the new estimation method has good estimation performance. The new estimation method is applied to the prediction and analysis of air quality data, which proves that the method is also feasible in practical applications.

1. 引言

假设有n个观测样本，记$T_{ij}$ 为第i个样本第j次观测的时间变量。$X_{ij}$ 与$Y_{ij}$ 分别为第i个样本第j次观测的协变量矩阵和响应变量，$m_i$ 为第i个样本重复观测的次数，则有纵向数据样本集$\left\{\left(X_{ij},T_{ij},Y_{ij}\right),i=1,\cdots ,n;j=1,\cdots ,m_i\right\}$ 。一类重要的纵向数据部分线性回归模型为下列形式：

$Y_{ij}=X_{ij}^{\text{T}}\beta +g\left(T_{ij}\right)+{\epsilon }_{ij}$ (1)

其中$\beta$ 为p维分位数回归系数向量，$g(\cdot )$ 是定义在有界闭区间$\left[0,1\right]$ 上的未知光滑函数，${\epsilon }_{ij}$ 为模型的连续误差项。

Zeger和Diggle [1]率先对模型(1)进行了研究，并将其应用于HIV研究中CD4细胞数的时间趋势的估计问题中。Lin和Carroll [2]提出了profile-kernel估计方法，估计兴趣参数。Hu等[3]将Zeger和Diggle提出的后移算法与Lin和Carroll提出的profile-kernel估计方法进行了比较。Fan和Li [4]使用局部多项式回归方法估计非参数函数，提出了两种简单有效的参数分量的估计方法：差分估计和profile最小二乘估计。田萍[5]详细总结了前人关于模型(1)在随机点列情况下所做过的工作，并讨论了其在固定设计点列下的情况。Tian和Xue [6]在部分线性回归模型框架内对纵向数据进行经验似然推断。王明辉和尹居良[7]研究了模型(1)的分位数估计问题，并证明了参数估计量的渐近正态性。刘会明[8]在模型(1)框架下，对局部多项式进行改进并提出了全局拟似然方法，给出了回归系数和非参数函数的拟似然估计。

在纵向数据分析中，分位数回归也是一种被广泛使用的估计方法。为了降低分位数回归推断的效率损失，Jung [9]首次提出了一种中位数回归的准似然方法，该方法将纵向数据重复测量之间的相关性纳入其中。Fu和Wang [10]结合重复测量之间的相关性，引入了诱导平滑方法来获得参数估计及其方差的方法。Leng和Zhang [11]通过组合多组无偏估计方程，提出了一种在考虑重复测量之间的相关性的同时，产生更有效估计的分位数回归模型。Fu和Wang、Leng和Zhang都将诱导平滑方法扩展到了分位数回归，从而获得了允许应用Newton-Raphson迭代的平滑目标函数，后者还给出了参数的估计及其协方差矩阵的Sandwich估计。Lu和Fan [12]在他们的基础上，提出了一个结合了纵向数据重复测量之间的相关性结构的纵向数据分位数回归模型。模型平均方法能够有效平衡估计的方差和偏差。在分位数回归框架下，Lu和Su [13]将Jackknife模型平均(JMA)方法应用到线性分位数回归，证明了该估计量的渐近最优性。对于部分线性模型的模型平均预测分析。Zhang和Wang [14]提出了具有独立误差的部分线性模型最佳模型平均。Fang [15]提出了一种用于二分反应的半参数模型平均预测方法。胡国治和曾婕[16]构造了基于部分线性分位数回归模型的模型平均估计量并探究了其大样本性质。这些半参数模型平均的例子都是在独立假设下进行的。在纵向数据中，Hu等[17]研究了纵向数据变系数部分线性模型下模型平均估计量的渐近分布。Li等[18]考虑过具有相关误差的半参数模型的模型平均，提供了一种最佳模型平均方法来改进纵向数据部分线性模型的预测。

受上述文献启发，针对模型(1)的参数与非参数估计，本文的创新点在于将分位数回归方法和模型平均估计拓展到其中的，采用局部线性估计对模型(1)的非参数部分拟合，对模型(1)的参数估计部分，提出一种稳健的分位数回归估计量，给出了估计的Newton-Raphson迭代算法及协方差矩阵的Sandwich估计提高了估计效率，利用JMA估计模型平均的最优权重。并通过数值模拟与实例分析，证明该方法优于纵向数据部分线性回归模型分位数估计、纵向数据线性回归模型分位数估计、纵向数据线性回归模型、分位数模型平均估计。

2. 分位数平均估计模型模型

2.1. 构建候选子模型

定义纵向数据样本集${\left\{\left(X_i,T_i,Y_i\right)\right\}}_{i=1}^n$ ，其中$Y_i={\left(Y_{i1},\cdots ,Y_{i{m}_i}\right)}^{\text{T}}$ 是第i个个体$m_i\times 1$ 维响应变量，$X_i={\left(X_{i1},\cdots ,X_{i{m}_i}\right)}^{\text{T}}$ 是第i个个体$m_i\times p$ 维协变量矩阵，$X_{ij}={\left(X_{ij1},\cdots ,X_{ijp}\right)}^{\text{T}}$ 是第i个个体第j次观测的p维协变量，$j=1,\cdots ,m_i$ ，$T_i={\left(T_{i1},\cdots ,T_{i{m}_i}\right)}^{\text{T}}$ 是第i个个体$m_i\times 1$ 维时间变量，为不失一般性，假定$T_{ij}\in \left[0,1\right]$ 。本文主要考虑下列纵向数据部分线性分位数回归模型：

${\mu }_{\tau ij}=Q_{\tau }\left(Y_{ij}|X_{ij}\right)=X_{ij}^{\text{T}}{\beta }_{\tau }+g_{\tau }\left(T_{ij}\right)$ (2)

其中${\beta }_{\tau }$ 为p维分位数回归系数向量，$g_{\tau }(\cdot )$ 是定义在有界闭区间$\left[0,1\right]$ 上的未知光滑函数，$\tau \in \left(0,1\right)$ 为感兴趣的分位数水平。记$g_{\tau }\left(T_i\right)={\left(g_{\tau }\left(T_{i1}\right),\cdots ,g_{\tau }\left(T_{i{m}_i}\right)\right)}^{\text{T}}$ ，${\mu }_{\tau i}={\left({\mu }_{\tau i1},\cdots ,{\mu }_{\tau i{m}_i}\right)}^{\text{T}}$ ，则有下列纵向数据部分线性分位数回归模型：

${\mu }_{\tau i}=Q_{\tau }\left(Y_i|X_i\right)=X_i{\beta }_{\tau }+g_{\tau }\left(T_i\right)$ (3)

为探究上述模型的模型平均估计量，我们采用嵌套的方式构建候选子模型，记$X_{ij}^{\left(s\right)}={\left(X_{ij1},\cdots ,X_{ijs}\right)}^{\text{T}}$ 是第s个候选子模型第i个个体第j次观测的协变量，$s=1,\cdots ,p$ ，此时第s个候选子模型$M_s$ 为：

${\mu }_{\tau i}^{\left(s\right)}=Q_{\tau }\left(Y_i|X_i^{\left(s\right)}\right)=X_i^{\left(s\right)}{\beta }_{\tau }^{\left(s\right)}+g_{\tau }^{\left(s\right)}\left(T_i\right)$ (4)

其中$X_i^{\left(s\right)}={\left(X_{i1}^{\left(s\right)},\cdots ,X_{i{m}_i}^{\left(s\right)}\right)}^{\text{T}}$ 是$m_i\times s$ 维协变量矩阵，${\beta }_{\tau }^{\left(s\right)}$ 为s维分位数回归系数向量，$g_{\tau }^{\left(s\right)}(\cdot )$ 是定义在有界闭区间$\left[0,1\right]$ 上的未知光滑函数。

2.2. 子模型的参数估计

为获得第s个候选子模型$M_s$ 中未知回归系数向量${\beta }_{\tau }^{\left(s\right)}$ 与未知光滑函数$g_{\tau }^{\left(s\right)}(\cdot )$ 的估计值，我们采用两阶段估计法。第一阶段，不考虑纵向数据组内相关性，获得${\beta }_{\tau }^{\left(s\right)}$ 与$g_{\tau }^{\left(s\right)}(\cdot )$ 的初始估计值${\tilde{\beta }}_{\tau }^{\left(s\right)}$ 与${\tilde{g}}_{\tau }^{\left(s\right)}(\cdot )$ 。

对于候选子模型中未知光滑函数$g_{\tau }^{\left(s\right)}(\cdot )$ 的估计，我们采用局部线性估计的非参数估计方法，在$T_{ij}$ 的一个小区域内对$g_{\tau }^{\left(s\right)}\left(T_{ij}\right)$ 进行一阶泰勒展开，即

$g_{\tau }^{\left(s\right)}\left(T_{ij}\right)\approx g_{\tau }^{\left(s\right)}\left(t\right)+{\dot{g}}_{\tau }^{\left(s\right)}\left(t\right)\left(T_{ij}-t\right)=Z_{ij}^{\text{T}}{\theta }_{\tau }^{\left(s\right)}$

其中${\dot{g}}_{\tau }^{\left(s\right)}(\cdot )$ 是$g_{\tau }^{\left(s\right)}(\cdot )$ 的一阶导数，$Z_{ij}={\left(1,T_{ij}-t\right)}^{\text{T}}$ 是局部线性估计中的设计矩阵，${\theta }_{\tau }^{\left(s\right)}={\left(g_{\tau }^{\left(s\right)}\left(t\right),{\dot{g}}_{\tau }^{\left(s\right)}\left(t\right)\right)}^{\text{T}}$是相应的参数向量。记$D_{ij}^{\left(s\right)}={\left(X_{ij}^{\left(s\right)\text{T}},Z_{ij}^{\text{T}}\right)}^{\text{T}}$ ，${\eta }_{\tau }^{\left(s\right)}={\left({\beta }_{\tau }^{\left(s\right)\text{T}},{\theta }_{\tau }^{\left(s\right)\text{T}}\right)}^{\text{T}}$ ，可以通过最小化下列目标函数得到${\eta }_{\tau }^{\left(s\right)}$ 的估计值

${\tilde{\eta }}_{\tau }^{\left(s\right)}=\underset{{\eta }_{\tau }^{\left(s\right)}}{\text{argmin}}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\rho }_{\tau }}}\left(Y_{ij}-D_{ij}^{\left(s\right)\text{T}}{\eta }_{\tau }^{\left(s\right)}\right)K_{h_1}\left(T_{ij}-t\right)$

其中，${\rho }_{\tau }\left(u\right)=u\left(\tau -{\text{I}}_{\left(u\le 0\right)}\right)$ 为给定$\tau$ 分位数下模型的损失函数，$K_{h_1}(\cdot )=K\left(\cdot /h_1\right)$ ，$K(\cdot )$ 是核函数，$h_1$ 为给定的窗宽。从而可以得到$g_{\tau }^{\left(s\right)}\left(t\right)$ 的初始估计值${\tilde{g}}_{\tau }^{\left(s\right)}\left(t\right)=\left(0_{1\times s},1,0\right){\tilde{\eta }}_{\tau }^{\left(s\right)}$ 。

尽管在上述估计中我们得到了${\beta }_{\tau }^{\left(s\right)}$ 的估计值，但由于该估计值是局部最优估计，其收敛速率未达到$\sqrt{n}$ 。因此，我们利用$g_{\tau }^{\left(s\right)}\left(T_{ij}\right)$ 的初始估计值${\tilde{g}}_{\tau }^{\left(s\right)}\left(T_{ij}\right)$ ，得到${\beta }_{\tau }^{\left(s\right)}$ 的全局最优初始估计

${\tilde{\beta }}_{\tau }^{\left(s\right)}=\underset{{\beta }_{\tau }^{\left(s\right)}}{\text{argmin}}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\rho }_{\tau }}}\left[Y_{ij}-{\tilde{g}}_{\tau }^{\left(s\right)}\left(T_{ij}\right)-X_{ij}^{\left(s\right)\text{T}}{\beta }_{\tau }^{\left(s\right)}\right]$

第二阶段，在考虑纵向数据组内相关性的基础上，利用${\beta }_{\tau }^{\left(s\right)}$ 与$g_{\tau }^{\left(s\right)}(\cdot )$ 在第一阶段获得的初始估计值${\tilde{\beta }}_{\tau }^{\left(s\right)}$ 与${\dot{g}}_{\tau }^{\left(s\right)}(\cdot )$ ，得到其最终估计值。令${\epsilon }_i^{\left(s\right)}={\left({\epsilon }_{i1}^{\left(s\right)},\cdots ,{\epsilon }_{i{m}_i}^{\left(s\right)}\right)}^{\text{T}}$ ，其中${\epsilon }_{ij}^{\left(s\right)}=Y_{ij}-g_{\tau }^{\left(s\right)}\left(T_{ij}\right)-X_{ij}^{\left(s\right)\text{T}}{\beta }_{\tau }^{\left(s\right)}$ 是候选子模型的连续误差项，具有未知条件密度函数$f_{ij}^{\left(s\right)}(\cdot )$ 。将${\phi }_{\tau }\left({\epsilon }_i^{\left(s\right)}\right)$ 的协方差矩阵表示为

$V_i^{\left(s\right)}=\text{Cov}\left({\phi }_{\tau }\left({\epsilon }_i^{\left(s\right)}\right)\right)=\text{Cov}\left(\begin{array}{c}\tau -{\text{I}}_{\left({\epsilon }_{i1}^{\left(s\right)}\le 0\right)}\\ ⋮\\ \tau -{\text{I}}_{\left({\epsilon }_{i{m}_i}^{\left(s\right)}\le 0\right)}\end{array}\right)$

且令${\Gamma }_i^{\left(s\right)}=\text{diag}\left[f_{i1}^{\left(s\right)}\left(0\right),\cdots ,f_{i{m}_i}^{\left(s\right)}\left(0\right)\right]$ 。Jung [9]提出了下列估计方程，用于计算考虑响应变量组内相关性后${\beta }_{\tau }^{\left(s\right)}$ 的估计值

$U_0\left({\beta }_{\tau }^{\left(s\right)}\right)={\displaystyle \sum _{i=1}^n{X}_i^{\left(s\right)\text{T}}}{\Gamma }_i^{\left(s\right)}{\left(V_i^{\left(s\right)}\right)}^{-1}{\phi }_{\tau }\left(Y_i-g_{\tau }^{\left(s\right)}\left(T_i\right)-X_i^{\left(s\right)}{\beta }_{\tau }^{\left(s\right)}\right)=0$ (5)

其中${\phi }_{\tau }\left(u\right)={\rho }_{{\tau }^{\prime }}\left(u\right)=\tau -1_{\left(u\le 0\right)}$ 。在求解估计方程(5)时${\Gamma }_i^{\left(s\right)}$ 的元素$f_{ij}^{\left(s\right)}\left(0\right)$ 的估计值${\widehat{f}}_{ij}^{\left(s\right)}\left(0\right)$ 可以参考Hendricks和Knenker [19]得到

${\widehat{f}}_{ij}^{\left(s\right)}\left(0\right)=2h_n{\left[X_{ij}^{\left(s\right)\text{T}}\left({\widehat{\beta }}_{\tau +h_n}^{\left(s\right)}-{\widehat{\beta }}_{\tau -h_n}^{\left(s\right)}\right)\right]}^{-1}$

当然，在某些情况下，当$f_{ij}^{\left(s\right)}\left(0\right)$ 的值难以估计时，${\Gamma }_i^{\left(s\right)}$ 可简单视为一个单位矩阵，此时${\beta }_{\tau }^{\left(s\right)}$ 的估计效率略有损失，本文就将${\Gamma }_i^{\left(s\right)}$ 视为一个单位矩阵。

由于在进行参数估计时，协方差矩阵$V_i^{\left(s\right)}$ 的估计非常复杂，很难正确指定，于是Lu和Fan [12]提出了以下估计方程

$U_1\left({\beta }_{\tau }^{\left(s\right)}\right)={\displaystyle \sum _{i=1}^n{X}_i^{\left(s\right)\text{T}}}{\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\phi }_{\tau }\left(Y_i-g_{\tau }^{\left(s\right)}\left(T_i\right)-X_i^{\left(s\right)}{\beta }_{\tau }^{\left(s\right)}\right)=0$ (6)

其中${\Sigma }_i^{\left(s\right)}\left(\rho \right)$ 是${\phi }_{\tau }\left({\epsilon }_i^{\left(s\right)}\right)$ 的工作协方差矩阵，${\Sigma }_i^{\left(s\right)}\left(\rho \right)=A_i^{\left(s\right)1/2}{C}_i^{\left(s\right)}\left(\rho \right)A_i^{\left(s\right)1/2}$ ，$A_i^{\left(s\right)}=\text{diag}\left[{\sigma }_{i11}^{\left(s\right)},\cdots ,{\sigma }_{i{m}_i{m}_i}^{\left(s\right)}\right]$ 是一个$m_i\times m_i$ 维对角矩阵，${\sigma }_{ijj}^{\left(s\right)}=\text{Var}\left({\phi }_{\tau }\left({\epsilon }_{ij}^{\left(s\right)}\right)\right)$ ，$C_i^{\left(s\right)}\left(\rho \right)$ 是${\phi }_{\tau }\left({\epsilon }_i^{\left(s\right)}\right)$ 的工作相关矩阵，$\rho$ 是相关指数参数。假定工作协方差矩阵${\Sigma }_i^{\left(s\right)}\left(\rho \right)$ 具有一般的平稳自相关结构，则工作相关矩阵$C_i^{\left(s\right)}\left(\rho \right)$ 具有如下形式

$C_i^{\left(s\right)}\left(\rho \right)=\left(\begin{array}{ccccc}1& {\rho }_1& {\rho }_2& \cdots & {\rho }_{m_i-1}\\ {\rho }_1& 1& {\rho }_2& \cdots & {\rho }_{m_i-2}\\ ⋮& ⋮& ⋮& & ⋮\\ {\rho }_{m_i-1}& {\rho }_{m_i-2}& {\rho }_{m_i-3}& \cdots & 1\end{array}\right)$

对于所有$i=1,\cdots ,n$ ，${\rho }_l$ 的估计值

${\widehat{\rho }}_l=\frac{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i-l}{\tilde{y}}_{ij}^{\left(s\right)}{\tilde{y}}_{i,j+l}^{\left(s\right)}/n\left(m_i-l\right)}}}{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{m_i}{\tilde{y}}_{ij}^{\left(s\right)2}/n{m}_i}}}$

其中$l=1,\cdots ,m_i-1$ ，Hendricks和Knenker [19]定义

${\tilde{y}}_{ij}^{\left(s\right)}=\frac{{\phi }_{\tau }\left(Y_{ij}-g_{\tau }^{\left(s\right)}\left(T_{ij}\right)-X_{ij}^{\left(s\right)\text{T}}{\beta }_{\tau }^{\left(s\right)}\right)}{\sqrt{{\delta }_{ijj}^{\left(s\right)}}}$

为了估计${\sigma }_{ijj}^{\left(s\right)}=\text{Var}\left({\phi }_{\tau }\left({\epsilon }_{ij}^{\left(s\right)}\right)\right)$ ，我们令${\tilde{\epsilon }}_i^{\left(s\right)}={\left({\tilde{\epsilon }}_{i1}^{\left(s\right)},\cdots ,{\tilde{\epsilon }}_{i{m}_i}^{\left(s\right)}\right)}^{\text{T}}$，${\tilde{\epsilon }}_{ij}^{\left(s\right)}=Y_{ij}-{\tilde{g}}_{\tau }^{\left(s\right)}\left(T_{ij}\right)-X_{ij}^{\left(s\right)\text{T}}{\beta }_{\tau }^{\left(s\right)}$。此时应用${\phi }_{\tau }\left({\tilde{\epsilon }}_i^{\left(s\right)}\right)={\left({\phi }_{\tau }\left({\tilde{\epsilon }}_{i1}^{\left(s\right)}\right),\cdots ,{\phi }_{\tau }\left({\tilde{\epsilon }}_{i{m}_i}^{\left(s\right)}\right)\right)}^{\text{T}}$，${\phi }_{\tau }\left({\tilde{\epsilon }}_{ij}^{\left(s\right)}\right)=\tau -1_{\left({\tilde{\epsilon }}_{ij}^{\left(s\right)}\le 0\right)}$ ，可得到

${\widehat{\delta }}_{ijj}^{\left(s\right)}=P\left({\tilde{\epsilon }}_{ij}^{\left(s\right)}\le 0\right)\left(1-P\left({\tilde{\epsilon }}_{ij}^{\left(s\right)}\le 0\right)\right)$

求解估计方程(6)的主要困难在于不可微的非凸和非连续目标函数$U_1\left({\beta }_{\tau }^{\left(s\right)}\right)$ ，为了克服这些问题，Fu和Wang [10]将诱导平滑方法扩展到对估计方程(7)的求解中。令${\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)=E_{H^{\left(s\right)}}\left[U_1\left({\beta }_{\tau }^{\left(s\right)}\right)+{\Omega }_1^{\left(s\right)1/2}{H}^{\left(s\right)}\right]$，其中$H^{\left(s\right)}~N\left(0,I_s\right)$ ，$I_s$ 为$s\times s$ 维单位矩阵，${\Omega }_1^{\left(s\right)}$ 为参数估计器更新后的协方差矩阵的估计。经过一系列代数计算，可得到下列平滑估计方程

${\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)={\displaystyle \sum _{i=1}^n{X}_i^{\left(s\right)\text{T}}}{\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\tilde{\phi }}_{\tau }\left(Y_i-{\tilde{g}}_{\tau }^{\left(s\right)}\left(T_i\right)-X_i^{\left(s\right)}{\beta }_{\tau }^{\left(s\right)}\right)=0$ (8)

其中，${\tilde{\phi }}_{\tau }\left({\tilde{\epsilon }}_i^{\left(s\right)}\right)={\left(\tau -1+\Phi \left(\frac{{\tilde{\epsilon }}_{i1}^{\left(s\right)}}{{\tilde{r}}_{i1}}\right),\cdots ,\tau -1+\Phi \left(\frac{{\tilde{\epsilon }}_{i{m}_i}^{\left(s\right)}}{{\tilde{r}}_{i{m}_i}}\right)\right)}^{\text{T}}$，${\tilde{r}}_{ij}=\sqrt{X_{ij}^{\left(s\right)\text{T}}{\Omega }_1^{\left(s\right)}{X}_{ij}^{\left(s\right)}}$。此时${\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)$ 对${\beta }_{\tau }^{\left(s\right)}$ 的微分可以很容易的计算出来，并且可以用$\frac{\partial {\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}$ 近似地代替$\frac{\partial U_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}$ 为

$\frac{\partial {\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}=-{\displaystyle \sum _{i=1}^n{X}_i^{\left(s\right)\text{T}}}{\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\tilde{\Lambda }}_i^{\left(s\right)}{X}_i^{\left(s\right)}$

其中${\tilde{\Lambda }}_i^{\left(s\right)}$ 为$m_i\times m_i$ 维对角矩阵，第j个对角元为$\varphi \left({\tilde{\epsilon }}_{ij}^{\left(s\right)}/{\tilde{r}}_{ij}\right)/{\tilde{r}}_{ij}$ 。

此时，${\beta }_{\tau }^{\left(s\right)}$ 及其协方差矩阵${\Omega }_1^{\left(s\right)}$ 的平滑估计量可以从下列Newton-Raphson迭代算法中获得：

步骤1：给定${\beta }_{\tau }^{\left(s\right)}$ 和协方差矩阵${\Omega }_1^{\left(s\right)}$ 的初始值，分别为${\widehat{\beta }}_{\tau }^{\left(s\right)}\left(0\right)={\tilde{\beta }}_{\tau }^{\left(s\right)}$ ，${\widehat{\Omega }}_1^{\left(s\right)}\left(0\right)=\frac{1}{n}{I}_s$。

步骤2：使用第r次迭代出的${\widehat{\beta }}_{\tau }^{\left(s\right)}\left(r\right)$ 和${\widehat{\Omega }}_1^{\left(s\right)}\left(r\right)$，通过下列公式更新${\widehat{\beta }}_{\tau }^{\left(s\right)}\left(r+1\right)$ 和${\widehat{\Omega }}_1^{\left(s\right)}\left(r+1\right)$

${\widehat{\beta }}_{\tau }^{\left(s\right)}\left(r+1\right)={\widehat{\beta }}_{\tau }^{\left(s\right)}\left(r\right)+{\left[\frac{\partial {\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}\times {\left[{\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)\right]}_r$

${\widehat{\Omega }}_1^{\left(s\right)}\left(r+1\right)={\left[-\frac{\partial {\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}\times {\left[\text{Cov}\left({\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)\right)\right]}_r\times {\left[-\frac{\partial {\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)}{\partial {\beta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}$

其中，$\text{Cov}\left({\tilde{U}}_1\left({\beta }_{\tau }^{\left(s\right)}\right)\right)={\displaystyle \sum _{i=1}^n{X}_i^{\left(s\right)\text{T}}}{\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\tilde{\phi }}_{\tau }\left({\tilde{\epsilon }}_i^{\left(s\right)}\right){\tilde{\phi }}_{\tau }^{\text{T}}\left({\tilde{\epsilon }}_i^{\left(s\right)}\right){\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{X}_i^{\left(s\right)}$，${\left[\right]}_r$ 表示方括号内的表达式在${\beta }_{\tau }^{\left(s\right)}={\widehat{\beta }}_{\tau }^{\left(s\right)}\left(r\right)$ 的计算结果。

步骤3：重复步骤2直至收敛，并记最终的收敛值为${\widehat{\beta }}_{\tau }^{\left(s\right)}$ ，${\widehat{\Omega }}_1^{\left(s\right)}$。

上述估计方法提供了${\beta }_{\tau }^{\left(s\right)}$ 及其协方差矩阵${\Omega }_1^{\left(s\right)}$ 的有效估计${\widehat{\beta }}_{\tau }^{\left(s\right)}$ 和${\widehat{\Omega }}_1^{\left(s\right)}$。接下来，利用${\widehat{\beta }}_{\tau }^{\left(s\right)}$ 估计候选子

模型中的$g_{\tau }^{\left(s\right)}(\cdot )$ 。令$\check{\epsilon }{}_i^{\left(s\right)}=\left(\check{\epsilon }{}_{i1}^{\left(s\right)},\cdots ,\check{\epsilon }{}_{i{m}_i}^{\left(s\right)}\right)$，$\check{\epsilon }{}_{ij}^{\left(s\right)}=Y_{ij}-X_{ij}^{\left(s\right)\text{T}}{\widehat{\beta }}_{\tau }^{\left(s\right)}-g_{\tau }^{\left(s\right)}\left(T_{ij}\right)$，利用局部线性估计得到的一阶

泰勒展开，参考Lv等[20]的估计方程(4)，提出下列平滑估计方程用于估计${\theta }_{\tau }^{\left(s\right)}$

${\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)={\displaystyle \sum _{i=1}^n{Z}_i^{\text{T}}}{K}_i\left(t;h_2\right){\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\tilde{\phi }}_{\tau }\left(Y_i-X_i^{\left(s\right)}{\widehat{\beta }}_{\tau }^{\left(s\right)}-Z_i{\theta }_{\tau }^{\left(s\right)}\right)=0$ (9)

其中${\tilde{\phi }}_{\tau }\left(\check{\epsilon }{}_i^{\left(s\right)}\right)={\left(\tau -1+\Phi \left(\frac{\check{\epsilon }{}_{i1}^{\left(s\right)}}{\check{r}{}_{i1}}\right),\cdots ,\tau -1+\Phi \left(\frac{\check{\epsilon }{}_{i{m}_i}^{\left(s\right)}}{\check{r}{}_{i{m}_i}}\right)\right)}^{\text{T}}$，$\check{r}{}_{ij}=\sqrt{Z_{ij}^{\text{T}}{\Omega }_2^{\left(s\right)}{Z}_{ij}}$，$Z_i={\left(Z_{i1},\cdots ,Z_{i{m}_i}\right)}^{\text{T}}$ ，$K_i\left(t;h_2\right)=\text{diag}\left(K_{h_2}\left(T_{i1}-t\right),\cdots ,K_{h_2}\left(T_{i{m}_i}-t\right)\right)$ 。同理，我们可以得到${\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)$ 关于${\theta }_{\tau }^{\left(s\right)}$ 的一阶微分为

$\frac{\partial {\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)}{\partial {\theta }_{\tau }^{\left(s\right)}}=-{\displaystyle \sum _{i=1}^n{Z}_i^{\text{T}}}{K}_i\left(t;h_2\right){\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}\check{\Lambda }{}_i^{\left(s\right)}{Z}_i$

其中$\check{\Lambda }{}_i^{\left(s\right)}$为$m_i\times m_i$ 维对角矩阵，第j个对角元为$\varphi \left(\check{\epsilon }{}_{ij}^{\left(s\right)}/\check{r}{}_{ij}\right)/\check{r}{}_{ij}$。

此时，${\theta }_{\tau }^{\left(s\right)}$ 及其工作协方差矩阵${\Omega }_2^{\left(s\right)}$ 的光滑估计量可以从下列Newton-Raphson迭代算法中获得：

步骤1：给定${\theta }_{\tau }^{\left(s\right)}$ 和协方差矩阵${\Omega }_2^{\left(s\right)}$ 的初始值，分别为${\widehat{\theta }}_{\tau }^{\left(s\right)}\left(0\right)={\left(\left(0_{1\times s},1,0\right){\tilde{\eta }}_{\tau }^{\left(s\right)},\left(0_{1\times s},0,1\right){\tilde{\eta }}_{\tau }^{\left(s\right)}\right)}^{\text{T}}$，

${\widehat{\Omega }}_2^{\left(s\right)}\left(0\right)=\frac{1}{n}{I}_2$。

步骤2：使用第r次迭代出的${\widehat{\theta }}_{\tau }^{\left(s\right)}\left(r\right)$ 和${\widehat{\Omega }}_2^{\left(s\right)}\left(r\right)$，通过下列公式更新${\widehat{\theta }}_{\tau }^{\left(s\right)}\left(r+1\right)$ 和${\widehat{\Omega }}_2^{\left(s\right)}\left(r+1\right)$

${\widehat{\theta }}_{\tau }^{\left(s\right)}\left(r+1\right)={\widehat{\theta }}_{\tau }^{\left(s\right)}\left(r\right)+{\left[-\frac{\partial {\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)}{\partial {\theta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}\times {\left[{\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)\right]}_r$

${\widehat{\Omega }}_2^{\left(s\right)}\left(r+1\right)={\left[-\frac{\partial {\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)}{\partial {\theta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}\times {\left[\text{Cov}\left({\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)\right)\right]}_r\times {\left[-\frac{\partial {\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)}{\partial {\theta }_{\tau }^{\left(s\right)}}\right]}_r^{-1}$

其中，$\text{Cov}\left({\tilde{U}}_2\left({\theta }_{\tau }^{\left(s\right)}\right)\right)={\displaystyle \sum _{i=1}^n{Z}_i^{\text{T}}}{K}_i\left(t;h_2\right){\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{\tilde{\phi }}_{\tau }\left(\check{\epsilon }{}_i^{\left(s\right)}\right){\tilde{\phi }}_{\tau }^{\text{T}}\left(\check{\epsilon }{}_i^{\left(s\right)}\right){\left[{\Sigma }_i^{\left(s\right)}\left(\rho \right)\right]}^{-1}{K}_i\left(t;h_2\right)Z_i$。

步骤3：重复步骤2直至收敛，并记最终的收敛值为${\widehat{\theta }}_{\tau }^{\left(s\right)}$ ，${\widehat{\Omega }}_2^{\left(s\right)}$。

最后，我们得到${\widehat{g}}_{\tau }^{\left(s\right)}\left(t\right)=\left(1,0\right){\widehat{\theta }}_{\tau }^{\left(s\right)}$ ，那么第s个候选子模型$M_s$ 下${\widehat{\mu }}_{\tau i}^{\left(s\right)}$ 的估计值

${\widehat{\mu }}_{\tau i}^{\left(s\right)}=X_i^{\left(s\right)}{\widehat{\beta }}_{\tau }^{\left(s\right)}+{\widehat{g}}_{\tau }^{\left(s\right)}\left(T_i\right)$ (10)

其中${\widehat{g}}_{\tau }^{\left(s\right)}\left(T_i\right)={\left({\widehat{g}}_{\tau }^{\left(s\right)}\left(T_{i1}\right),\cdots ,{\widehat{g}}_{\tau }^{\left(s\right)}\left(T_{i{m}_i}\right)\right)}^{\text{T}}$。

2.3. Jackknife权重选择

设$\omega ={\left({\omega }_1,\cdots ,{\omega }_p\right)}^{\text{T}}$ 是集合$\mathcal{W}=\left\{\omega \in {\left[0,1\right]}^p:0\le {\omega }_s\le 1,{\displaystyle \sum _{s=1}^p{\omega }_s}=1\right\}$ 中的权重向量，则${\mu }_{\tau i}$ 的模型平均估计量

${\widehat{\mu }}_{\tau i}\left(\omega \right)={\displaystyle \sum _{s=1}^p{\omega }_s}{\widehat{\mu }}_{\tau i}^{\left(s\right)}$ (11)

设${\widehat{\beta }}_{\tau i}^{\left(s\right)}$ ，${\widehat{g}}_{\tau i}^{\left(s\right)}(\cdot )$ 表示第s个子模型$M_s$ 中删除第i个个体的观测数据$\left\{X_i,T_i,Y_i\right\}$ 后${\mu }_{\tau i}^{\left(s\right)}$ 的Jackknife估计值，则将留一交叉验证的准则函数定义为

$C{V}_n\left(\omega \right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\rho }_{\tau }}\left(Y_i-{\displaystyle \sum _{s=1}^p{\omega }_s}\left[X_i^{\left(s\right)}{\widehat{\beta }}_{\tau i}^{\left(s\right)}+{\widehat{g}}_{\tau i}^{\left(s\right)}\left(T_i\right)\right]\right)$ (12)

此时，通过选择$\omega \in \mathcal{W}$ 以最小化上述准则函数，获得权重向量$\omega$ 的Jackknife选择$\widehat{\omega }=\left({\widehat{\omega }}_1,\cdots ,{\widehat{\omega }}_p\right)$ ，即

$\widehat{\omega }=\underset{\omega \in \mathcal{W}}{\text{argmin}}C{V}_n\left(\omega \right)$ (13)

最后，通过$\widehat{\omega }$ 可以得到${\mu }_{\tau i}$ 的Jackknife模型平均(JMA)估计量

${\widehat{\mu }}_{\tau i}\left(\widehat{\omega }\right)={\displaystyle \sum _{s=1}^p{\widehat{\omega }}_s}{\widehat{\mu }}_{\tau i}^{\left(s\right)}$ (14)

3. 数值模拟

3.1. 数据生成

本节利用R软件，根据模型(15)生成三种不同类型的数据，对比Lu和Fan [12]提出的纵向数据线性回归模型的分位数回归估计(LQR)、基于Lu和Fan [12]提出的纵向数据线性回归模型分位数回归估计的Jackknife模型平均估计(MLQR)、基于本文子模型估计方法的纵向数据部分线性回归模型分位数回归估计(BQR)、本文提出的纵向数据部分线性回归模型分位数回归模型平均估计(MBQR)这四种估计方法在样本外预测误差方面的表现情况。

将非参数核函数固定为Epanechnikov核，即$K\left(u\right)=0.75{\left(1-u^2\right)}_{+}$ ，最优窗宽h₁和h₂的选择采用五折交叉验证标准。以h₁为例，设$T-T^v$ 和$T^v$ 分别为$v=1,\cdots ,5$ 的交叉验证训练集和测试集，其中T是完整的数据集。对于每个h₁和v，我们使用训练集$T-T^v$ 找到${\eta }_{\tau }^{\left(s\right)}$ 的估计${\tilde{\eta }}_{\tau }^{\left(s\right)}$ ，然后，利用下列五折交叉验证标准：

$\text{CV}\left(h_1\right)={\displaystyle \sum _{v=1}^5\sum _{\left(Y_{ij},D_{ij}^{\left(s\right)}\right)\in T^v}}{\rho }_{\tau }\left(Y_{ij}-D_{ij}^{\left(s\right)\text{T}}{\tilde{\eta }}_{\tau }^{\left(s\right)}\right)$

通过使用网格搜索方法，我们可以找到最优窗宽$h_{1opt}={\text{min}}_{h_1}\text{CV}({\text{h}}_{\text{1}})$ 。

对模型(16)式纵向数据部分线性分位数回归模型中参数向量$\beta$ 与非参数平滑函数$g\left(t\right)$ 的设定如下：选择$n=30,50,100$ 代表小样本、中等样本和大样本量，数据模拟次数$M=100$ 。$p=5$ ，$m_i$ 服从二项分布${\rm B}\left(10,0.8\right)$ ，参数向量$\beta ={\left(-2,-1,1,3,0\right)}^{\text{T}}$ ，非参数平滑函数$g\left(t\right)$ 分别取$\sin \left(2\pi t\right)$ ，$\cos \left(5\pi t\right)$ ，${\text{e}}^{-2t}$ ，$5t^2+1$ ，

解释变量$X_{ij}$ 服从多元正态分布$N\left(0,\Sigma \right)$ ，其中，$\Sigma ={\left({\Sigma }_{ij}\right)}_{p\times p}$ ，${\Sigma }_{i,j}={\rho }^{\left|i-j\right|}$ ，$1\le i,j\le p$ ，$\rho =0.5$ ，时间变

量$T_{ij}$ 服从$\left[0,1\right]$ 上的均匀分布。在此设定下，根据模型(1.1)中误差项的分布，生成下述三种模拟数据：

模拟1：误差项${\epsilon }_i~N\left({\mu }_i,{\Phi }_i\right)$ 。其中，${\mu }_i$ 是一个$m_i$ 维向量，每个元素值都等于负标准正态分布的$\tau$ 分位数；${\Phi }_i$ 为一阶自回归AR(1)矩阵，相关系数${\rho }_{ij}=0.85$ 。

模拟2：误差项${\epsilon }_i~t\left(v_i,{\mu }_i,{\Phi }_i\right)$ 。其中，自由度$v_i=3$ ；${\mu }_i$ 是一个$m_i$ 维向量，每个元素值都等于负标准t分布的$\tau$ 分位数；${\Phi }_i$ 与模拟1相同。

模拟3：在模拟1生成机制的基础上，随机将样本数据中百分之五的响应变量用异常值进行替换，异常值为原值的0.1倍。

选取分位数$\tau =0.5,0.75$ ，使用R软件中的Mass、quantreg、sampling等程序包可以实现上述模拟数据的生成。根据四种估计方法在三种不同数据模型下的样本外预测误差来展示估计方法的表现情况，样本外预测误差的计算公式如下：

$\text{FPE}\left({\widehat{\mu }}_{\tau i}\right)=\frac{1}{M}{\displaystyle \sum _{r=1}^M\frac{1}{n}}{\displaystyle \sum _{i=1}^n\frac{1}{m_i}}{\displaystyle \sum _{j=1}^{m_i}{\rho }_{\tau }}\left(Y_{ij}-{\widehat{\mu }}_{\tau ij}\right)$

其中${\widehat{\mu }}_{\tau ij}$ 代表${\mu }_{\tau ij}=Q_{\tau }\left(Y_{ij}|X_{ij}\right)$ 的估计量，数据集$\left\{\left(X_{ij},Y_{ij},T_{ij}\right),i=1,\cdots ,n;j=1,\cdots ,m_i\right\}$ 代表预测集。

3.2. 结果分析

利用R软件自行编写算法程序，将上述模拟数据带入程序中，得到模拟结果如下表1~3所示：

表1四种估计方法在模拟1下针对不同的非参数平滑函数$g\left(t\right)$ 和样本容量n的预测误差。

Table 1. Out-of-sample forecast error under simulation 1

表1. 模拟1下样本外预测误差

$g\left(t\right)$	n	FPE ($\tau =0.5$ )				FPE ($\tau =0.75$ )
		MBQR	BQR	MLQR	LQR	MBQR	BQR	MLQR	LQR
$\sin \left(2\pi t\right)$	30	0.5224^*	0.5339	0.5624	0.5609	0.4107^*	0.4116	0.4296	0.4309
	50	0.6117^*	0.6156	0.6357	0.6367	0.4381^*	0.4396	0.4596	0.4609
	100	0.5503^*	0.5606	0.5767	0.5745	0.5025^*	0.5048	0.5320	0.5286
$\cos \left(5\pi t\right)$	30	0.5538^*	0.5740	0.6143	0.6517	0.5194^*	0.5206	0.5413	0.5478
	50	0.6325^*	0.6381	0.6345	0.6389	0.4558^*	0.4563	0.4723	0.4740
	100	0.5791^*	0.5802	0.5873	0.5886	0.5071^*	0.5092	0.5404	0.5399
${\text{e}}^{-2t}$	30	0.6139^*	0.6287	1.7130	1.7247	0.4798^*	0.4854	2.2677	2.2773
	50	0.5781^*	0.5830	1.5226	1.5236	0.4601^*	0.4706	1.9695	1.9773
	100	0.6393^*	0.6495	1.6368	1.6377	0.4876^*	0.6011	2.0626	2.0652
$5t^2+1$	30	0.5392^*	0.5531	1.4624	1.4663	0.4738^*	0.5070	1.5009	1.5012
	50	0.7352^*	0.8391	1.3237	1.3261	0.4095^*	0.4202	1.5545	1.5597
	100	0.6546^*	0.6649	1.3695	1.3665	0.4712^*	0.4795	1.6653	1.6717

注：^*表示四种估计方法中预测误差的最小值。

表2四种估计方法在模拟2下针对不同的非参数平滑函数$g\left(t\right)$ 和样本容量n的预测误差。

Table 2. Out-of-sample forecast error under simulation 2

表2. 模拟2下样本外预测误差

$g\left(t\right)$	n	FPE ($\tau =0.5$ )				FPE ($\tau =0.75$ )
		MBQR	BQR	MLQR	LQR	MBQR	BQR	MLQR	LQR
$\sin \left(2\pi t\right)$	30	0.7480^*	0.7554	0.7710	0.7727	0.5669^*	0.5676	0.5933	0.6043
	50	0.9442^*	0.9593	0.9768	0.9799	0.5751^*	0.5771	0.5874	0.5885
	100	0.7355^*	0.7378	0.7679	0.7698	0.6406^*	0.6444	0.6583	0.6586
$\cos \left(5\pi t\right)$	30	0.6530^*	0.6551	0.6670	0.6623	0.5927^*	0.6087	0.6168	0.6309
	50	0.8775^*	0.8783	0.8953	0.8940	0.5192^*	0.5196	0.5222	0.5206
	100	0.8327^*	0.8405	0.8426	0.8428	0.6976^*	0.6999	0.7478	0.7469
${\text{e}}^{-2t}$	30	0.9116^*	0.9137	1.6703	1.6777	0.6090^*	0.6145	1.9911	1.9982
	50	0.8357^*	0.8557	1.7444	1.7472	0.6564^*	0.6752	1.8838	1.8796
	100	0.8457^*	0.8698	1.7498	1.7514	0.6737^*	0.6823	2.1490	2.1453
$5t^2+1$	30	0.9263^*	0.9531	1.6178	1.6295	0.6872^*	0.7040	1.3731	1.3829
	50	0.7935^*	0.8615	1.4238	1.4211	0.5417^*	0.5717	1.7032	1.6960
	100	0.8614^*	0.8684	1.5538	1.5530	0.5779^*	0.5869	1.6428	1.6380

注：^*表示四种估计方法中预测误差的最小值。

表3四种估计方法在模拟3下针对不同的非参数平滑函数$g\left(t\right)$ 和样本容量n的预测误差。

Table 3. Out-of-sample forecast error under simulation 3

表3. 模拟3下样本外预测误差

$g\left(t\right)$	n	FPE ($\tau =0.5$ )				FPE ($\tau =0.75$ )
		MBQR	BQR	MLQR	LQR	MBQR	BQR	MLQR	LQR
$\sin \left(2\pi t\right)$	30	0.6307^*	0.6359	0.6670	0.6673	0.5980^*	0.6055	0.6090	0.6112
	50	0.6474^*	0.6691	0.6782	0.6784	0.4589^*	0.4643	0.4832	0.4842
	100	0.6332^*	0.6547	0.6618	0.6622	0.4975^*	0.4976	0.4985	0.4993
$\cos \left(5\pi t\right)$	30	0.6501^*	0.6666	0.6722	0.6809	0.5288^*	0.5355	0.5573	0.5742
	50	0.6249^*	0.6260	0.6297	0.6294	0.6299^*	0.6424	0.6788	0.6775
	100	0.6813^*	0.7117	0.7213	0.7205	0.5386^*	0.5391	0.5527	0.5526
${\text{e}}^{-2t}$	30	0.6173^*	0.6222	1.5898	1.5989	0.5783^*	0.5868	1.8352	1.8387
	50	0.7603^*	0.7629	1.7070	1.7090	0.4736^*	0.5139	1.8773	1.8813
	100	0.7722^*	0.7906	1.6655	1.6642	0.5410^*	0.5606	2.0703	2.0668
$5t^2+1$	30	0.6870^*	0.7053	1.3539	1.3505	0.6128^*	0.6255	1.8885	1.8715
	50	0.7890^*	0.8474	1.4400	1.4387	0.4714^*	0.4809	1.5493	1.5523
	100	0.6687^*	0.6705	1.4034	1.4013	0.4693^*	0.4763	1.5419	1.5421

注：^*表示四种估计方法中预测误差的最小值。

分析上述模拟结果，以表1中$\tau =0.5$ ，$n=100$ 时四种估计方法的样本外预测误差为例，可以发现在四种非参数光滑函数下，BQR的样本外预测误差值优于MLQR、LQR的样本外预测误差值，说明部分线性回归模型分位数回归估计方法的预测精度优于线性回归模型下的分位数回归估计与分位数回归Jackknife模型评价估计方法。而MBQR的样本外预测误差值优于BQR的样本外预测误差值，说明了提出的纵向数据部分线性回归模型分位数回归模型平均估计方法在一定条件下是可以达到提高样本预测精度的目的。

分析表1~3的全部表格结果可以发现，三角函数、指数函数和幂函数三种不同类型的非参数函数在随机误差分布分别为正态分布、t分布和数据中存在异常值的情况下的研究结果表明，本文提出的纵向数据部分线性回归模型分位数回归模型平均估计方法(MBQR)在数值模拟中的样本外预测误差都表现良好。

以模拟1中$\tau =0.5$ ，$n=100$ 为例，绘制出了在四种非参数光滑函数下MBQR估计方法得出的子模型权重分布情况如图1所示：

Figure 1. Weight distribution of sub-models under four non-parametric smooth functions

图1. 四种非参数光滑函数下子模型的权重分布情况图

分析图1可以发现在四种非参数光滑函数下高权重主要分布在第4和第5个子模型中，其余子模型权重都较低，说明这两个子模型能够涵盖大部分模型信息，根据我们在数据生成过程中将参数向量设置为$\beta ={\left(-2,-1,1,3,0\right)}^{\text{T}}$ 可以发现，MBQR下得到的权重分布情况是合理且符合实际的。总体而言，本文提出的估计方法在有限样本情况下具有良好的效果，这说明本文估计方法是可行的。

4. 实例分析

随着我国经济的迅猛发展，工业化、城市化进程的加速，带来了大量的污染气体排放和资源消耗，空气质量问题成为城市居民日常生活中不可忽视的困扰。我们通过收集和分析空气质量监测数据，可以更准确地评估不同地区、不同时间段的空气质量状况，了解污染物的排放来源和影响范围。同时，利用空气质量指数(AQI)等综合性评价指标，预测分析空气质量问题的变化趋势，为相关政策制定和环境治理提供科学依据。因此，对空气质量数据的分析是解决空气污染问题的重要手段。

本节将本文提出的方法应用到来自“中国环境检测总站”的空气质量面板数据，数据选取了北京、上海等50个城市在2020年1月到2020年12月的空气质量月度数据，包含600组观测值以及七个变量。表4列出了空气质量数据集中的变量名称及其含义：

表4空气质量数据的变量信息。

Table 4. Variable information of instance data

表4. 实例数据的变量信息

变量	名称	含义
$Y$	AQI指数	空气质量指数
$X_1$	PM2.5	空气中细颗粒物质量浓度
$X_2$	PM10	空气中可吸入颗粒物质量浓度
$X_3$	SO2	空气中二氧化硫质量浓度
$X_4$	NO2	空气中二氧化氮质量浓度
$X_5$	CO	空气中一氧化碳质量浓度
$X_6$	O3	空气中臭氧质量浓度

通常一个地区的AQI空气质量指数是由当地一段时间内空气中PM_2.5质量浓度、PM₁₀质量浓度、SO₂质量浓度、NO₂质量浓度、CO质量浓度、O₃质量浓度等六项参数通过计算而得出的空气污染程度及空气质量状况的表述。在本文的实证分析中，我们将AQI空气质量指数作为模型(1)中的响应变量Y，将其余六项参数作为模型(1.1)中的协变量矩阵X或T。为了确定模型(1)中参数部分的协变量矩阵X与非参数部分的协变量T，分别作出AQI空气质量指数与其余六项参数的散点图，如图2所示。

分析图2可以发现，AQI空气质量指数与PM_2.5质量浓度、PM₁₀质量浓度、NO₂质量浓度、CO质量浓度、O₃质量浓度有明显的线性关系，与SO₂质量浓度没有明显线性关系。因此，将PM_2.5质量浓度、PM₁₀质量浓度、NO₂质量浓度、CO质量浓度、O₃质量浓度作为参数部分的协变量X，SO₂质量浓度作为非参数部分的协变量T。将50个城市每月的空气质量数据作为一个观测点，为了探究数据点之间的相关性，我们从50个城市中随机选取了10个城市绘制出来每个城市之间AQI空气质量指数的相关系数图如图3所示。

Figure 2. Scatter plot of AQI air quality data and other six parameters

图2. AQI空气质量数据与其余六项参数数据的散点图

Figure 3. Correlation coefficient diagram of AQI air quality index between cities

图3. 城市与城市之间AQI空气质量指数的相关系数图

分析图3可以发现大多数城市之间的AQI空气质量指数相关性并不高，极少一部分城市因为地域相近所以AQI空气质量指数之间存在相关性。由此，我们针对50个城市的空气质量数据假定50个城市单个城市月与月之间的数据具有相关性，城市与城市之间的数据相互独立，适用于本文提出的纵向数据部分线性回归模型的估计方法。

为了体现本文方法与其它常用方法在实际应用中的差异，将本文提出的纵向数据部分线性回归模型分位数回归模型平均估计(MBQR)与Lu和Fan [14]提出的纵向数据线性回归模型的分位数回归估计(LQR)、基于Lu和Fan [14]提出的纵向数据线性回归模型分位数回归估计的Jackknife模型平均估计(MLQR)、基于本文提出的子模型估计方法的纵向数据部分线性回归模型分位数回归估计(BQR)三种估计方法同时应用到空气质量数据的估计中。在空气质量数据线性回归模型的建模中，我们将AQI空气质量指数作为线性回归模型中的响应变量Y，其余六项参数作为线性回归模型中的协变量矩阵X。将50个城市的空气质量数据以7:3的比例选择城市分别作为训练集和测试集，进行100次随机选择，得到四种估计方法分别在分位数$\tau =0.5,0.75$ 下的样本外预测误差如表5所示。

表5空气质量数据下四种不同估计方法的预测误差。

Table 5. Prediction error of example data

表5. 实例数据的预测误差

估计方法	MBQR	BQR	MLQR	LQR
FPE ($\tau =0.5$ )	2.7020^*	2.7178	2.7434	2.7698
FPE ($\tau =0.75$ )	2.2418^*	2.2518	2.2823	2.3166

注：^*表示四种估计方法中预测误差的最小值。

分析表5中的样本外预测误差可以看出，本文估计方法MBQR相较于其余三种估计方法在不同分位数情况下的预测误差都最小，说明本文提出的估计方法在实际数据的预测分析中也能表现良好。我们还绘制出了上述空气质量数据在$\tau =0.5,0.75$ 时，MBQR估计方法得出的子模型权重分布情况如图4所示：

Figure 4. Weight distribution of sub-models under MBQR estimation method

图4. MBQR估计方法下子模型的权重分布情况图

分析图4可以发现，高权重主要分布在第5个子模型中，其余子模型权重都较低，说明第5子模型能够涵盖大部分模型信息，根据我们对空气质量数据七个指标的实际情况下分析AQI空气质量指数是由当地一段时间内空气中PM_2.5质量浓度、PM₁₀质量浓度、SO₂质量浓度、NO₂质量浓度、CO质量浓度、O₃质量浓度等六项参数综合计算而得出的关于空气污染程度及空气质量状况的表述，MBQR下得到的权重分布情况是合理且符合实际的。总体而言，本文提出的MBQR估计方法为实际数据中的纵向数据分析提供了一种更有效的新途径。

5. 结语

本文的创新点在于将模型平均方法与分位数回归估计方法相结合，提出了针对纵向数据部分线性回归模型的新估计方法，用工作相关矩阵分解和估计方程平滑法来处理纵向数据的组内相关性，用局部线性估计来处理部分线性回归模型中的非参数部分，给出了模型参数估计的Newton-Raphson迭代算法。并通过数值模拟和实例分析验证了该估计方法的优良性。