纵向数据是指对同一组个体在不同时间点上重复观测得到的数据，在经济学、生物学、医学等领域有着广泛的应用.半参数模型既有参数分量，又含有非参数分量，能够刻画纵向数据内部的相关性.Zeger和Diggle^[1]用迭代算法讨论半参数模型的估计问题，He^[2]采用回归样条估计非参数部分，柴根象等^[3]利用二阶段估计方法讨论该模型的估计问题，赵江^[4]研究了半参数回归函数的收敛性，田萍^[5-6]等讨论了纵向数据下半参数模型的相合性及渐近性质.本文建立纵向半参数负二项模型，对此模型参数估计，并得出估计的性质，最后并利用R软件随机模拟.

1 模型介绍

考虑以下纵向数据半参数模型

$ {y_{ij}} = x_{ij}^T{\beta _p} + V\left( {{t_{ij}}} \right) + {\varepsilon _{ij}},i = 1,2, \cdots ,m,j = 1,2, \cdots ,{n_i} $

(1)

其中β是p维未知参数，X_ij是协变量，X_ijβ_j是参数固定效应部分，t_ij为随机变量且独立同分布，ε_ij为模型误差，E(ε_ij)=0，Var(ε_ij)=σ² < ∞，y_ij是第i个个体在时间t_ij的响应值，N=mn_i.

2 参数估计

对变量x_ij=(x_ij1, x_ij2, …, x_ijp)给定m×n_i个观测样本，同一般线性模型一样，利用观测数据对模型中的参数β_p及函数V(t_ij)做出估计.用π(x_ij, t_ij)表示当x=x_ij, t=t_ij时事件发生的概率，则Y_ij是服从参数n和π(x_ij, t_ij)的负二项分布.V(t_ij)是关于时间t_ij的函数.

设

$ V\left( {{t_{ij}}} \right) = {\beta _{p + 1}}{T_1}\left( {{t_{ij}}} \right) + {\beta _{p + 2}}{T_2}\left( {{t_{ij}}} \right) + \cdots + {\beta _{p + q}}{T_q}\left( {{t_{ij}}} \right) = T_{ij}^T{\beta _q} $

(2)

其中T_ij=(T₁(t_ij), T₂(t_ij), …, T_q(t_ij))^T是基函数向量.W_ij=(x_ij^T, T_ij^T)^Tβ=(β_p^T, β_q^T)^T,

则相应得半参数负二项模型为：

$ {y_{ij}} = x_{ij}^T{\beta _p} + T_{ij}^T{\beta _q} + {\varepsilon _{ij}} = W_{ij}^T\beta + {\varepsilon _{ij}},i = ,1,2, \cdots ,m,j = 1,2, \cdots ,{n_i} $

则

(Y₁₁=y₁₁, Y₁₂=y₁₂, …)的对数似然函数为

$ \begin{array}{l} L\left( {\beta ,{y_{11}},{y_{12}}, \cdots ,{y_{m{n_i}}}} \right)\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {C_{n + {y_{ij}} - 1}^{n - 1}\pi {{\left( {{x_{ij}},{t_{ij}}} \right)}^n}\left( {1 - \pi {{\left( {{x_{ij}},{t_{ij}}} \right)}^{{y_{ij}}}}} \right)} } \\ \ln L\left[ {P\left( {{Y_{ij}} = {y_{ij}}} \right)} \right] = \ln \left( {\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {C_{n + {y_{ij}} - 1}^{n - 1}} } } \right) + n\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {\ln \pi \left( {{x_{ij}},{t_{ij}}} \right)} } + {y_{ij}}\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {\ln \left( {1 - \pi \left( {{x_{ij}},{t_{ij}}} \right)} \right)} } \\ = \ln \left( {\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {C_{n + {y_{ij}} - 1}^{n - 1}} } } \right) + n\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {\ln \left( {\frac{{\exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}{{1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}} \right)} } + \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {\ln \left( {\frac{1}{{1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}} \right)} } {y_{ij}}\\ = \ln \left( {\prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {C_{n + {y_{ij}} - 1}^{n - 1}} } } \right) + n\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {\ln \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right) - \left( {n + 1} \right)} } + \prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {\ln \left( {1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)} \right)} } + \prod\limits_{i = 1}^m {\prod\limits_{j = 1}^{{n_i}} {\ln {y_{ij}}} } \end{array} $

(3)

故对数似然函数对β进行一阶求导得出如下结果：

$ \begin{array}{l} \frac{{\partial \ln L\left( {\beta ,{y_{11}},{y_{12}}, \cdots ,{y_{1n}},{y_{21}},{y_{22}}, \cdots ,{y_{mni}}} \right)}}{{{\partial _\beta }}}\\ = n\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {W_{ij}^T - \left( {n + 1} \right)} } \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {\frac{{\exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)W_{ij}^T}}{{1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}} } \\ = n\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {W_{ij}^T - \left( {n + 1} \right)} } \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {W_{ij}^T\left( {1 - \frac{1}{{1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}} \right)} } \\ = \left( {n + 1} \right)\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {\frac{{W_{ij}^T}}{{1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)}}} } - \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {W_{ij}^T} } \end{array} $

(4)

关于β求解方程，得到β的极大似然估计 ${\hat \beta }$ ，但是方程(4) 是非线性方程，需要用迭代方法求解.

3 Fisher信息阵

下面在(3) 式基础上对似然函数进行参数的二阶偏导数，并给出Fisher信息阵.

$ \begin{array}{l} \frac{{{\partial ^2}\ln L\left( {\beta ,{y_{11}}, \cdots ,{y_{n1}},{y_{21}}, \cdots ,{y_{m{n_i}}}} \right)}}{{\partial \beta \times \partial {\beta ^T}}}\\ = - \left( {n + 1} \right)\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {{{\left[ {1 + \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right)} \right]}^{ - 2}} \cdot \exp \left( {W_{ij}^T\beta + {\varepsilon _{ij}}} \right){W_{ij}}W_{ij}^T} } \\ = - \left( {n + 1} \right)\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {{W_{ij}}W_{ij}^T\pi \left( {{x_{ij}},{t_{ij}}} \right)\left( {1 - \pi \left( {{x_{ij}},{t_{ij}}} \right)} \right)} } \end{array} $

(5)

由Fisher信息阵的性质可知I_ij(θ)= ${E_\theta }\left( {\frac{{\partial {\rm{In}}{f_\theta }\left( x \right)}}{{\partial {\theta _i}}}\frac{{\partial {\rm{In}}{f_\theta }\left( x \right)}}{{\partial {\theta _j}}}} \right)$，则参数β的Fisher信息阵可以表示为

$ \begin{array}{l} I\left( \beta \right) = E{\left( { - \frac{{{\partial ^2}\ln L\left( {\beta ,{y_{11}}, \cdots ,{y_{n1}},{y_{21}}, \cdots ,{y_{m{n_i}}}} \right)}}{{\partial \beta \times \partial {\beta ^T}}}} \right)_{\left( {p + q} \right) \times \left( {p + q} \right)}}\\ = {\left[ {\left( {n + 1} \right)\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^{{n_i}} {{W_{ij}}W_{ij}^T\pi \left( {{x_{ij}},{t_{ij}}} \right)\left( {1 - \pi \left( {{x_{ij}},{t_{ij}}} \right)} \right)} } } \right]_{\left( {p + q} \right) \times \left( {p + q} \right)}} \end{array} $

4 似然函数的Newton-Raphson迭代方法

对数似然函数ln(β, y₁₁, y₁₂, …, y_1n1, …, y_mni)为目标函数，当β=β^(t)时，

$ {u^{\left( t \right)}} = \frac{{\partial \ln L\left( {\beta ,{y_{11}},{y_{12}}, \cdots ,{y_{1n}},{y_{21}},{y_{22}}, \cdots ,{y_{mni}}} \right)}}{{{\partial _\beta }}}, $

$ {v^{\left( t \right)}} = \frac{{{\partial ^2}\ln L\left( {\beta ,{y_{11}}, \cdots ,{y_{n1}},{y_{21}}, \cdots ,{y_{m{n_i}}}} \right)}}{{\partial \beta \times \partial {\beta ^T}}}, $

$ {U^{\left( t \right)}} = {\left( {u_0^{\left( t \right)},u_1^{\left( t \right)}, \cdots ,u_{p + q - 1}^{\left( t \right)}} \right)^T},{V^{\left( t \right)}} = {\left( {{v^{\left( t \right)}}} \right)_{\left( {p + q} \right) \times \left( {p + q} \right)}} $

则可以得到负二项模型参数β的最大似然估计的Newton-Raphson迭代公式是

$ {\beta ^{\left( {t + 1} \right)}} = {\beta ^{\left( t \right)}} - {\left( {{V^{\left( t \right)}}} \right)^{ - 1}}{U^{\left( t \right)}}. $

5 MLE的渐近性质

参数的极大似然估计(MLE)具有大样本性质，即有相合性和渐近正态性.根据文献[7]我们令

θ=(β，σ²)^T，设θ的极大似然估计${\hat \theta }$是由最大化目标函数

$ {G_N}\left( {y,\theta } \right) = \frac{1}{N}\sum\limits_{i = 1}^N {g\left( {{y_i},\theta } \right)} $

得到，则θ的估计方程是${{G'}_N}\left( {y, \theta } \right) = \frac{1}{N}\sum\limits_{i = 1}^N {g'} \left( {{y_i}, \theta } \right) = 0, g'\left( {y, \theta } \right) = \frac{{\partial g\left( {{y_i}, \theta } \right)}}{{\partial \theta }}$.

假设θ₀是${\overline G _N}\left( \theta \right)$ =E[G_N(y, θ)]的最大值，从而E[G_N(y, θ)]=0将${{G'}_N}\left( {y, {\theta _N}} \right)$在θ₀处展开得${{G'}_N}\left( {y, \theta } \right) \approx {{G'}_N}\left( {y, {\theta _0}} \right) + {{G''}_N}\left( {y, \tilde \theta } \right)\left( {{\theta _N} - {\theta _0}} \right)$.其中$\tilde \theta = \alpha '{\theta _N} + \left( {1 - \alpha '} \right){\theta _0}, 0 < \alpha ' < 1$.由于${{G'}_N}\left( {y, \theta } \right) = 0$，故

$ \sqrt N \left( {{{\hat \theta }_N} - {\theta _0}} \right) = {\left[ { - {{G''}_N}\left( {y,\tilde \theta } \right)} \right]^{ - 1}}\left( {\sqrt N {{G'}_N}\left( {y,{\theta _0}} \right)} \right). $

定理1  在一定的正则条件下，根据文献[8]有

1) $\hat \theta \to {\theta _0}\;\;{\rm{a}}{\rm{.s}}{\rm{.}}$

2) $\sqrt N \left( {\hat \theta - {\theta _0}} \right)\mathop L\limits_ \to N\left( {0, J_{_N}^{ - 1}\left( {{\theta _0}} \right)} \right)$.

其中${J_{_N}^{ - 1}\left( {{\theta _0}} \right)}$ =M_N^-1Q_NM_N^-1, M_N= $ - {{\bar G''}_N}\left( {{\theta _0}} \right),{{\bar G''}_N}\left( {{\theta _0}} \right) = E\left[ {{{G''}_N}\left( {y,\theta } \right)} \right],{Q_N} = {\mathop{\rm var}} \left[ {\sqrt N {{G'}_N}\left( {y,{\theta _0}} \right)} \right]$.