Robust Estimation and Statistical Inference for Longitudinal Multi-Kink Regression Model Based on Exponential Square Loss

考虑折点个数K及折点位置τ = (τ₁, …, τ_K)^T均未知的纵向多折点回归模型:

其中: n为个体数, m_i为第i个个体的重复观测次数, (a)₊ =a ·I (a ≥0), I为示性函数, Y_ij为响应变量观测值, Z_ij为p维协变量观测值, X_ij为有界门限变量, e_ij为随机误差, 记e_i =(e_i1, …, e_imi)^T. 由模型(1) 可知, 变量X_ij的回归系数在折点τ =(τ₁, …, τ _K)^T处会发生变化. 记b =(b₁, …, b_K)^T, ϑ = (a₀, a₁, b ^T, β^T)^T, 以及参数向量θ = (ϑ^T, τ^T)^T, 其中$\boldsymbol{\vartheta} \in G \subset \mathbb{R}^{2+K+p}$, τ ∈ T ⊂ D^K. G, T均为紧集, D为X_ij的支撑集. 下面, 首先对假设折点个数K已知后的参数估计算法进行介绍, 随后给出用于确定折点个数K的模型选择方法.

在模型(1) 中, 对于k =1, …, K, b_k (X_ij -τ_k)₊关于τ_k不可导, 使用文献[3] 所提的基于局部线性平滑的快速迭代法, 将b_k (X_ij -τ_k)₊在与折点真实位置相近的τ_k⁽⁰⁾附近进行一阶泰勒展开

进而有如下的近似回归模型

其中$ \mu_k=b_k\left(\tau_k-\tau_k^{(0)}\right).$记$\boldsymbol{\mu}=\left(\mu_1, \cdots, \mu_K\right)^{\mathrm{T}}, $ $ \boldsymbol{\eta}=\left(\boldsymbol{\vartheta}^{\mathrm{T}}, \boldsymbol{\mu}^{\mathrm{T}}\right)^{\mathrm{T}}$为未知系数向量. 为便于表达, 记$\tilde{\boldsymbol\chi}_{i j}=(1, \left.X_{i j}, \left(X_{i j}-\tau_1^{(0)}\right)_{+}, \cdots, \left(X_{i j}-\tau_k^{(0)}\right)_{+}, Z_{i j}^{\mathrm{T}}, -I\left(X_{i j}>\tau_1^{(0)}\right), \cdots, -I\left(X_{i j}>\tau_K^{(0)}\right)\right)^{\mathrm{T}}, \tilde{\boldsymbol\chi}_i=\left(\tilde{\boldsymbol\chi}_{i 1}, \cdots, \right.\left.\tilde{\boldsymbol\chi}_{i m_i}\right)^{\mathrm{T}}, \tilde{\boldsymbol{e}}_i=\left(\tilde{e}_{i 1}^{(0)}, \cdots, \tilde{e}_{i m_i}^{(0)}\right)^{\mathrm{T}}.$基于文献[12], 本文基于指数平方损失函数对模型(2) 的未知参数进行估计, 即关于η最小化如下目标函数

其中$\varphi_\gamma(x)=1-\exp \left(-\frac{x^2}{\gamma}\right)$为指数平方损失函数, γ称为调谐参数. 不难发现, 当γ取值较小时, 可以为大的绝对值残差赋予较小的权重, 从而降低异常值对参数估计的影响. 当γ → ∞时该损失函数与最小二乘的平方损失函数近似. 因此在指数平方损失函数中, 可以通过选择合适的γ实现模型参数的自适应稳健估计. 考虑纵向数据个体重复观测之间的相关性, 结合(4) 式, 为提高参数估计效率, 本文提出求解如下广义估计方程

以获得未知参数向量η的估计. 其中$\boldsymbol{Y}_i=\left(Y_{i 1}, \cdots, Y_{i m_i}\right)^{\mathrm{T}}, \psi_\gamma(x)=\frac{2 x}{\gamma} \exp \left(-\frac{x^2}{\gamma}\right) \cdot \boldsymbol{V}_i$ 用于刻画个体残差之间的相关性，理论上最优的矩阵$\boldsymbol{V}_i=\operatorname{Cov}\left(\psi_\gamma\left(\boldsymbol{Y}_i-\tilde{\boldsymbol\chi}_i \boldsymbol{\eta}\right)\right)$然而由于个体内部相关性无法观测，V _i的具体形式亦难以确定. 文献[7]提出QIF法，该方法使用基矩阵的线性形式替代V_i^-1，但所能刻画的仍是等相关和AR(1)等一些具有特殊结构的矩阵. 本文使用文献[8]提出的修正Cholesky分解方法，该方法不局限于特殊相关结构，具有更广泛的适用性. 具体来说，同文献[15]类似，将V _i分解为如下的矩阵乘积形式

其中Φ _i是主对角线元素全为1的下三角矩阵，第(j，k)个元素是如下自回归方程系数ϕ_ijk, _γ的相反数

特别地，当j=1时，$\varphi_\gamma\left(Y_{i j}-\tilde{\boldsymbol\chi}_{i j} \boldsymbol{\eta}\right)=\varepsilon_{i j}. \boldsymbol{D}_i=\operatorname{diag}\left(d_{i 1}^2, \cdots, d_{i m_i}^2\right)$, 其中$d_{i j}^2=\operatorname{Var}\left(\varepsilon_{i j}\right)$为新息方差. 对ϕ_ijk及d_ij²建立广义线性方程

其中：$\varrho^{}$与r_ij为q维向量，ρ与ω_ijk为p维向量，r_ij及ω _ijk通常与观测时间相关. 参考文献[15]，在本文的模拟与实证中均设定$\boldsymbol{r}_{i j}=\left(1, t_{i j}, \cdots, t_{i j}^{q-1}\right)^{\mathrm{T}}, \boldsymbol{\omega}_{i j k}=\left(1, t_{i j}-t_{i k}, \cdots, \left(t_{i j}-t_{i k}\right)^{p-1}\right)^{\mathrm{T}}, $其中t_ij为时间变量，记录第i名受试者第j次观测的时间. 记$\boldsymbol{\xi}=\left(\boldsymbol{\rho}^{\mathrm{T}}, \boldsymbol{\varrho}^{\mathrm{T}}, \boldsymbol{\eta}^{\mathrm{T}}\right)^{\mathrm{T}}$关于ξ求解如下广义估计方程组

其中

$\boldsymbol{T}_i^{\mathrm{T}}=-\frac{\partial \boldsymbol{\varepsilon}_i^{\mathrm{T}}}{\partial \boldsymbol{\rho}}, \boldsymbol{\varepsilon}_i=\left(\varepsilon_{i 1}, \cdots, \varepsilon_{i m_i}\right)^{\mathrm{T}}, \boldsymbol{d}_i^2=\left(d_{i 1}^2, \cdots, d_{i m_i}^2\right)^{\mathrm{T}}, \boldsymbol{R}_i=\left(r_{i 1}, \cdots, r_{i m_i}\right)^{\mathrm{T}}, \boldsymbol{W}_i=\operatorname{Cov}\left(\boldsymbol{\varepsilon}_i^2\right)$参考文献[15]，在本文的模拟与实证中均设定$\boldsymbol{W}_i=2 \operatorname{diag}\left(d_{i 1}^4, \cdots, d_{i m_i}^4\right) .$

本文使用Newton-Raphson迭代算法求解方程组(8)以保证估计的精度. 下面给出在指定调谐参数γ下，基于指数平方损失的纵向多折点回归模型的参数估计算法：

步骤1  给定初始折点位置τ ⁽⁰⁾，计算模型(3)的普通最小二乘估计作为η的初始值η ⁽⁰⁾.

步骤2  基于τ ^(s)及η ^(s)，应用修正的Cholesky分解法估计协方差阵V _i，使用Newton-Raphson算法迭代求解方程组(8)获得η ^(s+1)，具体步骤如下：

步骤2.1  设定V_i^(s+1，0)为m_i×m_i的单位阵，ρ ^(s+1，0)= 0，$\varrho^{(s+1, 0)}=\bf{0}$给定折点位置τ ^(s)，以模型(3)的普通最小二乘估计作为η ^(s+1，0).

步骤2.2  基于η ^(s+1，r)，$e^{(s+1, r)}$以及ρ ^(s+1，r)，求解η ^(s+1，r+1)，具体步骤如下：

步骤2.2.1  基于η ^(s+1，r)及$e^{(s+1, r)}$通过下式迭代至收敛得到ρ ^(s+1，r+1)：$\boldsymbol{\rho}^{\left(s+1, r+1, \iota_1+1\right)}=\boldsymbol{\rho}^{\left(s+1, r+1, \iota_1\right)}+\left.\left\{\left(\sum_{i=1}^n \boldsymbol{T}_i^{\mathrm{T}} \boldsymbol{D}_i^{-1} \boldsymbol{T}_i\right)^{-1} \times \sum_{i=1}^n \boldsymbol{T}_i^{\mathrm{T}} \boldsymbol{D}_i^{-1} \boldsymbol{\varepsilon}_i\right\}\right|_{\eta=\eta(s+1, r), \boldsymbol{\rho}=\boldsymbol{\rho}(s+1, r+1, l1), \varrho=e^{(s+1, r)}}$

步骤2.2.2  基于η ^(s+1，r)及ρ ^(s+1，r+1)，通过下式迭代至收敛得到$e^{(x+1, r+1)}$：$\varrho^{\left(s+1, r+1, l_2+1\right)}=\varrho^{\left(s+1, r+1, l_2\right)}+\left.\left\{\left(\sum_{i=1}^n \boldsymbol{R}_i^{\mathrm{T}} \boldsymbol{D}_i \boldsymbol{W}_i^{-1} \boldsymbol{D}_i \boldsymbol{R}_i\right)^{-1} \times \sum_{i=1}^n \boldsymbol{R}_i^{\mathrm{T}} \boldsymbol{D}_i \boldsymbol{W}_i^{-1}\left(\boldsymbol{\varepsilon}_i^2-d_i^2\right)\right\}\right|_{\eta=\eta(s+1, r), \boldsymbol{\rho}=\rho(s+1, r+1), \varrho=\varrho^{\left(s+1, r+1, l^2\right)}}$

步骤2.2.3  根据(6)式及(7)式，可得V_i，r+1. 通过下式迭代至收敛得到η ^(s+1，r+1)：$\boldsymbol{\eta}^{\left(s+1, r+1, l_3+1\right)}=\boldsymbol{\eta}^{\left(s+1, r+1, l_3\right)}+$ $\left\{\left(\sum_{i=1}^n \tilde{\boldsymbol\chi}_i^{\mathrm{T}} \boldsymbol{V}_{i, r+1}^{-1} \boldsymbol{\Lambda}_i(\boldsymbol{\eta}) \tilde{\boldsymbol\chi}_i\right)^{-1} \times \sum_{i=1}^n \tilde{\boldsymbol\chi}_i^{\mathrm{T}} \boldsymbol{V}_{i, r+1}^{-1} \psi_\gamma\left(\boldsymbol{Y}_i-\tilde{\boldsymbol\chi}_i \boldsymbol{\eta}\right)\right\}\boldsymbol{\eta}=\boldsymbol{\eta}^{(s+1, r+1, l 3)}, \boldsymbol{\rho}=\boldsymbol{\rho}^{(s+1, r+1)}, \boldsymbol{\varrho}=\boldsymbol{\varrho}^{(s+1, r+1)}$

其中$\boldsymbol{\Lambda}_i(\boldsymbol{\eta})=\operatorname{diag}\left\{\Lambda_{i 1}(\boldsymbol{\eta}), \cdots, \Lambda_{i m_i}(\boldsymbol{\eta})\right\}, \Lambda_{i j}(\boldsymbol{\eta})=\psi_\gamma^{\prime}\left(Y_{i j}-\tilde{\boldsymbol\chi}_{i j} \boldsymbol{\eta}\right), \psi_\gamma^{\prime}(x)=\exp \left(-\frac{x^2}{\gamma}\right)\left(\frac{2}{\gamma}-\frac{4 x^2}{\gamma^2}\right)$

步骤2.3  重复步骤2.2直至参数收敛，可得η ^(s+1).

步骤3  更新折点位置τ ^(s+1). 通过

步骤4  重复步骤2-步骤3直至参数收敛.

注1  当V _i为单位阵时，方程(5)和最小化目标函数(4)的解等价. 因此，在上述步骤中省略步骤2.2.1及步骤2.2.2，设置步骤2.2.3中V _i，r+1为单位阵，即可获得独立工作矩阵情况的参数估计.

在指数平方损失函数$\varphi_\gamma(x)=1-\exp \left(-\frac{x^2}{\gamma}\right)$中，如何选择调谐参数γ实现回归参数的最佳稳健估计是一个重要问题. 参考文献[15]，使用网格搜索法，选择使得回归参数估计值的渐近协方差阵的行列式最小的γ作为最优调谐参数γ_opt，具体实现方法见第2节中注 2.

1.2节给出了当折点个数给定的时候模型参数的估计算法. 然而实际问题中，折点个数真值K₀通常是未知的. 参考文献[4]，提出如下基于指数平方损失的贝叶斯准则以确定折点个数：

其中$\begin{aligned} & \wedge \\ & \boldsymbol{\eta}_k \end{aligned}$为指定折点个数为k时获得的参数估计，P_k为此时模型中未知参数个数，S_n的表达式见(4)式，C_n > 0是与样本量有关的常数. 给定折点个数的最大值K^*，依次设定折点个数k=0，1，…，K^*，选择使得(10)式最小的k为实际折点个数的估计值，记作$\hat{K} .$

记参数$\boldsymbol{\xi}=\left(\boldsymbol{\rho}^{\mathrm{T}}, \boldsymbol{\varrho}^{\mathrm{T}}, \boldsymbol{\vartheta}^{\mathrm{T}}, \boldsymbol{\mu}^{\mathrm{T}}\right)^{\mathrm{T}}$的真值为$\boldsymbol{\xi}_0=\left(\boldsymbol{\rho}_0^{\mathrm{T}}, \boldsymbol{\varrho}_0^{\mathrm{T}}, \boldsymbol{\vartheta}_0^{\mathrm{T}}, \boldsymbol{\mu}_0^{\mathrm{T}}\right)^{\mathrm{T}}, $参数的估计值记作$\hat{\boldsymbol{\xi}}=\left(\hat{\boldsymbol{\rho}}{ }^{\mathrm{T}}, \hat{\varrho}^{\mathrm{T}}, \hat{\boldsymbol{\vartheta}}^{\mathrm{T}}, \right.\left.\hat{\boldsymbol{\mu}}^{\mathrm{T}}\right)^{\mathrm{T}}$定义$\tilde{\boldsymbol\chi}_i=\left(\tilde{\boldsymbol\chi}_{1, i}, \tilde{\boldsymbol\chi}_{2, i}\right), \tilde{\boldsymbol\chi}_{1, i}=\left(\tilde{\boldsymbol\chi}_{1, i 1}, \cdots, \tilde{\boldsymbol\chi}_{1, i m_i}\right)^{\mathrm{T}}, \tilde{\boldsymbol\chi}_{1, i j}=\left(1, X_{i j}, \left(X_{i j}-\tau_1^{(0)}\right)_{+}, \cdots, \left(X_{i j}-\right.\right.$ $\left.\left.\tau_K^{(0)}\right)_{+}, \boldsymbol{Z}_{i j}^{\mathrm{T}}\right), \tilde{\boldsymbol\chi}_{2, i}=\left(\tilde{\boldsymbol\chi}_{2, i 1}, \cdots, \tilde{\boldsymbol\chi}_{2, i m_i}\right)^{\mathrm{T}}, \tilde{\boldsymbol\chi}_{2, i j}=\left(-I\left(X_{i j}>\tau_1^{(0)}\right), \cdots, -I\left(X_{i j}>\tau_K^{(0)}\right)\right)$因此，

则$\boldsymbol{U}(\boldsymbol{\xi})=\left(\boldsymbol{U}_1^{\mathrm{T}}(\boldsymbol{\xi}), \boldsymbol{U}_2^{\mathrm{T}}(\boldsymbol{\xi}), \boldsymbol{U}_{1, 3}^{\mathrm{T}}(\boldsymbol{\xi}), \boldsymbol{U}_{2, 3}^{\mathrm{T}}(\boldsymbol{\xi})\right)^{\mathrm{T}} .$记$\boldsymbol{\Gamma}_n=\left(\boldsymbol{\Gamma}_n^{k l}\right)_{k, l=1, 2, 3, 4}$为$\frac{\boldsymbol{U}\left(\boldsymbol{\xi}_0\right)}{\sqrt{n}}$的协方差阵. 当$n \rightarrow \infty$时$\boldsymbol{\Gamma}_n^{k l} \stackrel{p}{\longrightarrow} \boldsymbol{\Gamma}^{k l}$记$ \boldsymbol{\Gamma}=\left(\boldsymbol{\Gamma}^{k l}\right)_{k, l=1, 2, 3, 4}.$根据1.2节的讨论，协方差阵V _i与参数ξ相关，为方便表示，本节将V _i重记为V _i(ξ). 定义如下矩阵

给出如下条件：

(ⅰ) E(ψ_γ(e_ij))=0，E(ψ′ _γ(e_ij))>0，且对于任意γ>0，E(ψ_γ(e_ij)²)有界.

(ⅱ) 折点的真值K₀及协变量Z_ij的维数p是固定的，个体的观测次数m_i一致有界，且个体数n趋于无穷.

(ⅲ) X_ij在定义域上有连续的密度函数. 当$n \rightarrow \infty, \frac{1}{n} \sum_{i j}\left\|\boldsymbol{Z}_{i j}\right\|^2=O$(1). 另外，w_ijk，r_ij及矩阵W_i均有界.

(ⅳ) 记Σ _i，0=Cov(ψ_γ(e_i))，矩阵Σ _i，0及V_i^-1(ξ₀) Σ _i，0 V_i^-1(ξ₀)的最大特征值均有界.

(ⅴ) 对于k=1，…，K，折点参数的初始值τ_k⁽⁰⁾满足$\tau_{k, 0}-\tau_k^{(0)}=O_p\left(n^{-\frac{1}{2}}\right) .$

(ⅵ) 参数空间$\mathit{\Theta}\in \mathbb{R}^{2+p+2 K}$是紧集，参数真值η₀是参数空间Θ的内点，S_n(η)在η ₀处唯一取得全局最小值.

(ⅶ) 当$n \rightarrow \infty$时，$\frac{C_n \log n}{\sqrt{n}} \rightarrow 0$

定理1  假设真实折点个数K₀已知，且条件(ⅰ)-(ⅴ)成立，则当$n \rightarrow \infty$

其中$\boldsymbol{A}=\operatorname{diag}\left(\boldsymbol{A}^{11}, \boldsymbol{A}^{22}, \boldsymbol{A}^{33}, \boldsymbol{A}^{44}\right), \boldsymbol{A}^{k k}=\lim \limits_{n \rightarrow \infty} \boldsymbol{A}_n^{k k}, k=1, 2, 3, 4 .$

注2  由定理1，回归系数估计$\hat{\boldsymbol{\eta}}=\left(\hat{\boldsymbol{\vartheta}}^{\mathrm{T}}, \hat{\boldsymbol{\mu}}^{\mathrm{T}}\right)^{\mathrm{T}}$的渐近协方差阵可通过下式对其进行估计，即

其中

选择使参数估计值的渐近协方差阵行列式，即$\operatorname{det}(\overbrace{\operatorname{Cov}(\hat{\boldsymbol{\eta}})})$最小的γ为最优调谐参数γ_opt. 在本文后续的模拟与实证中，采用网格搜索法寻找γ_opt，在2~50之间以步长3进行搜索.

定理2  假设真实折点个数K₀已知，且条件(ⅰ)-(ⅴ)成立，则当$n \rightarrow \infty$时

其中$\boldsymbol{B}=\operatorname{diag}\left(\frac{1}{b_{1, 0}}, \cdots, \frac{1}{b_{K, 0}}\right).$

定理3  假设条件(ⅰ)-(ⅶ)成立，则对于$\stackrel{\wedge}{K}=\arg\min _{k=0, \cdots, K^*} \mathrm{BIC}(k)$当$n \rightarrow \infty$时

对于多折点回归模型的折点存在性检验，文献[4]基于无折点的原假设H₀：b_k=0，k=1，…，K提出CUSUM(cumulative summation)型统计量. 本文沿用文献[4]所提方法，提出基于指数平方损失的纵向折点回归模型的折点存在性检验方法，具体步骤如下：

步骤1  计算H₀成立时系数参数$\zeta=\left(a_0, a_1, \boldsymbol{\beta}^{\mathrm{T}}\right)^{\mathrm{T}}$的估计$\hat{\zeta}, $，即关于ζ最小化如下目标函数

其中$\boldsymbol{M}_{i j}=\left(1, X_{i j}, \boldsymbol{Z}_{i j}^{\mathrm{T}}\right)^{\mathrm{T}}, \varphi_\gamma(x)=1-\exp \left(-\frac{x^2}{\gamma}\right)$为指数平方损失函数.

步骤2  基于原始样本数据计算$T_n(\gamma)=\sup _{\tau \in T}\left|R_n(\tau, \gamma, \hat{\zeta})\right|$其中

$\psi_\gamma(x)=\frac{2 x}{\gamma} \exp \left(-\frac{x^2}{\gamma}\right), (a)_{-}=a \cdot I(a <0), $ T是折点参数的取值范围.

步骤3  生成服从标准正态分布的随机样本{v₁，…，v_n}，计算$T_n^*(\gamma)=\sup _{\tau \in T}\left|R_n^*(\tau, \gamma, \hat{\zeta})\right|$其中

步骤4  重复步骤3L次，计算

作为统计量T_n (γ)的p值. 一般情况，若$\begin{aligned} & \wedge \\ & p_n \end{aligned}$小于显著性水平0.05，则拒绝原假设.

注3  在步骤2和步骤3中$R_n^*(\tau, \gamma, \hat{\zeta})$与$R_n(\tau, \gamma, \hat{\zeta})$渐近等价，其证明过程与文献[6]类似.

根据如下模型生成数据

其中$Z_{i j} \stackrel{\mathrm{i} \mathrm{idd}}{\sim} N(5, 1) .$对于门限变量X_ij，时间变量t_ij以及随机误差e_ij，考虑以下3种情形：

情形1  考虑时间t_ij为门限变量，残差服从正态分布的情况. 设置t_ij=X_ij~U(1，10)，m_i=10，i=1，…，n. e_i=(e_i1，…，e_imi)^T服从多元正态分布$N\left(\bf{0}, \boldsymbol{\Sigma}_i\right), \boldsymbol{\Sigma}_i=\nabla_i^{-1} \boldsymbol{B}_i\left(\nabla_i^{\mathrm{T}}\right)^{-1}, $其中B _i对角元素为$\exp \left(-0.5+0.5 u_{i j}\right), u_{i j} \sim N(0, 2) . \nabla_i$是单位下三角矩阵，j行k列元素$\delta_{j k}^{(i)}=-\left(0.2+0.8\left(t_{i j}-t_{i k}\right)\right), k <j .$

情形2  考虑门限变量不为时间，且数据存在异常值的情况. 设置X_ij~U(-5，5)，t_ij={1，2，…，10}，随后每个观测以20%的概率缺失以生成不平衡数据. e_i服从自由度为3，协方差阵为情形1中Σ_i的多元t分布，随后随机选取10%的e_ij服从标准柯西分布.

情形3  为公平起见，以等相关结构生成残差来说明修正的Cholesky分解法估计参数的有效性. 设置e _i服从自由度为3，协方差阵为$\boldsymbol{\Sigma}_i=\boldsymbol{A}_i^{\frac{1}{2}} \boldsymbol{C}_i \boldsymbol{A}_i^{\frac{1}{2}}$的多元正态分布，其中$\boldsymbol{A}_i=\operatorname{diag}\left(\sigma_{i 1}^2, \cdots, \sigma_{i m_i}^2\right), \sigma_{i j}^2=2 t_{i j}^2, \boldsymbol{C}_i$是相关系数为0.85的等相关结构. 随后随机选取10%的e_ij服从标准柯西分布. 其他设置与情形2相同.

每种情形均考虑3种折线效应，对于情形1，设定：①K=1，τ =7，ϑ =(-2，1，-3，1)^T；②K=2，τ =(5，7)^T，ϑ =(2，1，-4，5，1)^T；③K=3，τ =(2，5，7)^T，ϑ =(0，1，-3，2，-6，1)^T. 对于情形2、3，系数ϑ设置与情形1相同，折点个数和位置设置为：①K=1，τ =-1；②K=2，τ =(-1，1)^T；③K=3，τ =(-4，-1，4)^T. 设定样本量n=400，重复模拟次数为100次. 使用1.2节所提算法求解最小化目标函数(4)及方程组(8)，分别记为“ESL.IND”和“ESL.CHO”，即独立工作矩阵结构及基于Cholesky分解的模型参数的估计方法.

表 1给出3种情形下所有模拟中选择的最优参数γ_opt的均值. 可以看到，无异常值的情形下，γ_opt均值较大；情形2，3的γ_opt较小，以降低异常值产生的影响. 这说明在指数平方损失函数中，可以通过选择合适的调谐参数实现回归系数的自适应稳健估计.

表 2给出了真实折点个数为2时，不同C_n在3种情形下的折点个数的正确选择率. 可以看到，所有情形下，本文所提出的折点个数选择方法均具有较高的正确选择率，并且“ESL.CHO”能够实现更高的选择正确率. 而且样本量不变时，较大的C_n能略微提高正确选择折点数的概率.

为说明所提方法优势，与以下估计量进行比较：

1) 方程(5)的工作协方差矩阵V_i为等相关及AR(1)的特定结构，利用文献[7]提出的QIF方法对(5)式进行求解，同时分别考虑指数平方损失及经典的平方损失两种损失函数，记所得估计量分别为“ESL.EXCH”“OLS.EXCH”“ESL.AR1”“OLS.AR1”；

2) 文献[4]提出的基于分位数回归的纵向多折点模型的估计量，指定分位水平为0.5，记为“MKQR”，使用R程序包MultiKink实现；

3) 文献[3]提出的多折点模型的最小二乘估计量，记为“SEG”，使用R程序包segmented实现.

表 3展示了情形2在K=2时的参数估计结果，汇总了100次模拟中估计量的平均偏差、标准差及均方误差. 其余情形的参数估计的结果已上传至Github(https://github.com/Tangming-hub). 可以发现：

1) 当残差向量服从正态分布时(情形1)，几种方法所得到的估计量的估计效果相近. 然而当数据中存在异常值时(情形2和情形3)，基于指数平方损失和分位数回归的估计量“ESL.IND”“ESL.CHO”“ESL.EXCH”“ESL.AR1”和“MKQR”均能提供回归参数的有效估计，其中相对于分位数回归方法，基于指数平方损失函数的估计量表现更佳，而基于经典的平方损失函数的估计量的估计效果均较差；

2) 仅考虑指数平方损失函数的估计方法时，相对于独立工作矩阵的估计量“ESL.IND”，融合组内相关性的估计量“ESL.CHO”“ESL.EXCH”“ESL.AR1”均具有更优良的估计效果，且本文所提出的“ESL.CHO”在各指标上表现最佳. 这说明有效考虑纵向数据个体重复观测有利于提高回归参数的估计效率. 综上，本文所提出的估计方法，可以为纵向多折点回归模型的参数提供更为稳健的、有效的估计.

此外，不同情形下K=2时，3种估计方法“ESL.CHO”“ESL.IND”“MKQR”的标准差、标准误、95%的Wald型置信区间的平均长度及经验覆盖率亦上传至Github. 由结果可见，3种估计量的标准差与经验标准误差接近，并且所构造的置信区间的经验覆盖率均在置信水平95%左右. 并且，相对于“MKQR”方法，本文所提方法的置信区间的长度更短.

为研究本文第3节所提出的折线效应存在性检验统计量T_n的有限样本性质，考虑折点个数为2时检验统计量的功效. 对于情形1，设置-b₁=b₂=0，0.1，0.2，0.3；对于情形2和情形3，设置-b₁=b₂=0，0.15，0.3，0.45，其中当系数b₁和b₂都等于0时，不存在折线效应. 设置L=300，显著性水平α=0.05. 表 4展示所得检验统计量均值(Mean of T_n，简记为Mean-T_n)及经验p值(Power). 根据统计量的定义，当原假设成立时，统计量T_n应接近于0，模拟结果与理论一致. 注意到，当折线效应不存在时，经验P值接近名义上的显著性水平0.05，而随着折线效应增强，检验功效增加，并趋近于1，这说明本文提出的检验统计量能够有效识别折线效应.

本文亦使用所提方法“ESL.IND”“ESL.CHO”分析文献[4]的纵向黄体酮数据，参数的估计值、平均绝对误差、置信区间，以及拟合曲线(结果见Github). 实证结果显示，相较“SEG”“MKQR”方法，所提估计具有明显竞争优势.

为了证明定理1，我们给出如下引理.

引理1  若真实折点个数K₀已知，且条件(ⅰ)-(ⅴ)成立，则当$n \rightarrow \infty$时$\boldsymbol{\eta} \stackrel{\mathrm{p}}{\longrightarrow} \boldsymbol{\eta}_0 .$

引理1的证明  根据文献[16]结论，证$\hat{\eta}$为$\boldsymbol{\eta}_0$的一致估计量，只需证明存在常数C>0，对任意ε，

(16) 式等价于

下面证明$\hat{\vartheta}$的渐近一致性，$\hat{\rho}, \hat{\varrho}, \hat{\mu}$的性质证明可类似进行. 使用中值定理，经计算，有

其中$\tilde{e}_i^*$位于$\tilde{e}_i$及$e_{i, 0}$之间，$\boldsymbol{E}_i=\operatorname{diag}\left(\psi^{\prime}{ }_\gamma\left(\tilde{e}_i^*\right)\right)$对于$\hat{\boldsymbol{\vartheta}}-\boldsymbol{\vartheta}_0=C O_p\left(n^{-\frac{1}{2}}\right)$条件(ⅰ)成立时，显然有$\mathrm{E}\left(I_1\right)=0, \mathrm{E}\left(I_2\right)=0.$且在条件(ⅰ)-(ⅳ)成立时，有

以及

另一方面，由于$\hat{\boldsymbol{\vartheta}}-\boldsymbol{\vartheta}_0=C O_p\left(n^{-\frac{1}{2}}\right), $因此，$I_1=C O_p(1), I_2=C^2 O_p\left(n^{-\frac{1}{2}}\right)$注意到，$\tilde{e}_{i j}^{(0)}-e_{i j, 0}=\tilde{\boldsymbol\chi}_{1, i j}\left(\boldsymbol{\vartheta}_0-\hat{\boldsymbol{\vartheta}}\right)+\sum_{k=1}^K b_{k, 0}\left\{\left(X_{i j}-\tau_{k, 0}\right)_{+}-\left(X_{i j}-\tau_k^{(0)}\right)_{+}+\left(\tau_{k, 0}-\tau_k^{(0)}\right) I\left(X_{i j}>\tau_k^{(0)}\right)\right\}, $记$b_{k, 0}\left\{\left(X_{i j}-\tau_{k, 0}\right)_{+}-\right.\left.\left(X_{i j}-\tau_k^{(0)}\right)_{+}+\left(\tau_{k, 0}-\tau_k^{(0)}\right) I\left(X_{i j}>\tau_k^{(0)}\right)\right\}$为$\Delta_{i j, k} \text {, }$因此，

其中$\Delta_{i, k}=\left(\sum_{k=1}^K \Delta_{i 1, k}, \cdots, \sum_{k=1}^K \Delta_{i m i, k}\right)^{\mathrm{T}} .$当条件(ⅰ)-(ⅳ)成立时，有$I_{3.1}=-C^2 O_p(1).$经计算，$\left|\Delta_{i j}\right| <\left|b_{k, 0}\right|\left|\tau_{k, 0}-\tau_k^{(0)}\right|$因此，当条件(ⅴ)成立时，$\Delta_{i j, k}=O_p\left(n^{-\frac{1}{2}}\right)$由条件(ⅰ)-(ⅴ)可得I_3，2=CO_p(1). 因此，I₃=-C²O_p(1). 类似可得，$I_4=C^3 O_p\left(n^{-\frac{1}{2}}\right)$综合上述结论，可得

因此，存在常数C>0使得$\left(\hat{\boldsymbol{\vartheta}}-\boldsymbol{\vartheta}_0\right)^{\mathrm{T}} \boldsymbol{U}_{1, 3}(\hat{\boldsymbol{\xi}}) <0.$故而引理1可证.

定理1的证明  下面证明$\hat{\vartheta}$渐近服从正态分布，$\hat{\boldsymbol{\mu}}, \hat{\boldsymbol{\rho}}, \hat{\varrho}$的证明可类似进行. 由中值定理，根据估计方程(8)所得的估计量$\hat{\vartheta}$满足

其中ξ ^*位于$\hat{\xi}$与ξ ₀之间，因此，据引理1，$n \rightarrow \infty$时$\xi^* \stackrel{\mathrm{p}}{\longrightarrow} \xi_{0 .}$由连续映射定理可知，

记$\boldsymbol{U}_{1, 3}^{(i)}\left(\boldsymbol{\xi}_0\right)=\left.\tilde{\boldsymbol\chi}_{1, i}^{\mathrm{T}} \boldsymbol{V}_i^{-1} \psi_\gamma\left(\boldsymbol{Y}_i-\tilde{\boldsymbol\chi}_{1, i} \boldsymbol{\eta}\right)\right|_{\xi=\xi_0}, $条件(ⅰ)成立时，其期望$\mathrm{E} \boldsymbol{U}_{1, 3}^{(i)}\left(\xi_0\right)=\bf{0}, $协方差阵

条件(ⅱ)-(ⅲ)的成立保证了存在常数κ₀使得$\operatorname{Cov}\left(\boldsymbol{U}_{1, 3}^{(i)}\left(\xi_0\right)\right) \leqslant \kappa_0 \boldsymbol{I}_{m_i \times m_i}$因此$\sum_{i=1}^{\infty} \frac{\operatorname{Cov}\left(\boldsymbol{U}_{1, 3}^{(i)}\left(\boldsymbol{\xi}_0\right)\right)}{i^2} <\infty$应用李雅普洛夫中心极限定理，有

结合(17) 式及(18) 式,

故而定理1可证.

定理2的证明  记折点位置真值为$\tau_0=\left(\tau_{1, 0}, \cdots, \tau_{K, 0}\right)^{\mathrm{T}} .$由(9)式可得等式

根据定理1，$\begin{gathered} \wedge \\ \mu_k \end{gathered}$与$\hat{b}_k$分别为$ \mu_{k, 0}, b_{k, 0}$的$\sqrt{n}$相合估计，因此

根据Slutsky定理可得$\sqrt{n}\left(\hat{\tau}_k-\tau_{k, 0}\right)$与$\frac{1}{b_{k, 0}} \sqrt{n}\left(\hat{\mu}_k-\mu_{k, 0}\right)$具有相同的渐近分布. 定理2可证.

定理3的证明  为证$\stackrel{\wedge}{K}$为折点个数真值K₀的相合估计，我们只需证明当$n \rightarrow \infty$时

经计算，

根据条件(ⅵ)，有$S_n\left(\hat{\boldsymbol{\eta}}_K\right)-S_n\left(\boldsymbol{\eta}_0\right) \geqslant 0 .$对$S_n\left(\boldsymbol{\eta}_0\right)-S_n\left(\hat{\boldsymbol{\eta}}_{K_0}\right)$使用泰勒展开，

由定理1知$\hat{\boldsymbol{\eta}}_{K 0}$为η₀的$\sqrt{n}$相合估计，因此，$S_n\left(\boldsymbol{\eta}_0\right)-S_n\left(\hat{\boldsymbol{\eta}}_{K_0}\right)=o\left(n^{\frac{1}{2}}\right)$计算可得

考虑K > K₀和K < K₀两种情况. K>K₀时，由于K-K₀>0，$n \rightarrow \infty$时显然有BIC(K)-BIC(K₀) >0依概率趋于1；对于K < K₀，条件(ⅶ)成立时，当$n \rightarrow \infty$时

定理3得证.

本文基于指数平方损失提出纵向多折点回归模型参数估计和统计推断方法. 为处理折点回归模型，本文首先基于局部线性平滑方法将折点回归模型转化为普通的纵向线性模型. 然后为融合纵向数据中重复观测间的相关性，本文利用修正的Cholesky分解方法对重复观测间的协方差阵进行建模，以提高回归模型参数的估计效率. 本文讨论了参数估计的大样本性质，并同时讨论了指数平方损失函数中调谐参数的选择、折点个数的确定方法和折线效应的检验问题等. 数值模拟和实证分析结果显示本文所提方法可以为纵向多折点回归模型的参数提供更为稳健的、有效的估计.