纵向零膨胀计数数据的多参数联合回归模型

Multi-parameter joint regression model for longitudinal zero-inflated count data

摘要

纵向零膨胀计数数据是对一系列试验个体在多个时间点进行测量得到的数据，兼具纵向数据、计数数据及零膨胀数据的特点，数据中包含的零值数量远超泊松、负二项分布等经典离散分布能够随机产生的零值个数。在对纵向零膨胀计数数据进行回归分析时，一般认为数据由零膨胀数据与随机采样数据构成，均值、零膨胀率与离散度是数据的主要特征。目前多数零膨胀数据回归模型都仅考虑协变量对均值和零膨胀率的影响而忽略离散度，或将其设定为固定值，因而无法适用于不同观测时间点的离散度随时间或其他协变量动态变化情形。对于服从零膨胀负二项分布的数据，本文设计了包含全部3个参数与协变量的广义回归模型，用EM（Expectation Maximization）算法得到了回归系数的极大似然估计。理论分析表明，相较双参数回归模型估计，本文的估计具有相合性和渐近正态性，能够同时精准估计各回归系数。最后，模拟分析结果显示，相较于未考虑时变离散度的回归模型，本文模型更加准确有效。

Abstract

Longitudinal zero-inflated count data are gathered from a series of measurements of experimental individuals at multiple time points and possess the very characteristics of longitudinal data， count data and zero-inflated data simultaneously.The number of zeros contained in count data far exceeds that randomly generated by classical discrete distributions such as Poisson and negative binomial distributions.In the regression analysis of longitudinal zero-inflated count data， the data are generally assumed to consist of zero-inflated components and random sampling components， and the mean， zero-inflation rate and dispersion are three main characteristic parameters of the data.Nowadays， most regression models for zero-inflated data only consider the influence of covariates on the mean and zero-inflation rate and the dispersion is unfortunately igored or set to a fixed value.As a result， these models cannot be applied to the situations where dispersion at different observation time points dynamically changes with time or other covariates.For the count data following the zero-inflated negative binomial distribution， a generalized regression model is develpoed to describe the relationship between all three data characteristics and the covariates.The maximum likelihood estimate of the regression coefficients is obtained by using the EM （Expectation-Maximization） algorithm to simultaneously the three regression coefficients.Theoretical analysis shows that， in comparison with the two-parameter regression model estimation， the proposed likelihood estimator is consistent and asymptotically normal， enabling accurate simultaneous estimation of the three regression coefficients.Finally， simulation results demonstrate that the proposed model is more accurate and effective than the regression models that do not consider time-varying dispersion.

Graphical abstract

关键词

纵向零膨胀计数数据 / 联合建模 / 零膨胀负二项回归 / EM算法

Key words

longitudinal zero-inflated count data / joint modeling / zero-inflated negative binomial regression / EM algorithm

零膨胀计数数据常见于农业、医疗、公共卫生及金融等领域，数据中零值的个数远超经典离散分布（如泊松、负二项分布等）所能随机产生的零值个数。一般认为，零膨胀计数数据由结构零和随机零组成。例如，在某地区居民因哮喘每年就诊次数数据中，许多居民的就诊次数为零，则其中无哮喘病居民的就诊次数必然为零（结构零），同时虽有哮喘但病情控制良好或因其他原因没有就诊的居民的就诊次数也是零（随机零）。

零膨胀计数数据的分布通常呈现双峰特征，一个峰在零值附近，另一个峰在非零值区域。对零膨胀计数数据进行建模主要使用零膨胀泊松（Zero-inflated Poisson， ZIP）分布和零膨胀负二项（Zero-inflated Negative Binomial， ZINB）分布。Mullahy^［1］首先提出了具有混合比例的ZIP分布。Lambert^［2］假设零膨胀数据以概率p服从退化分布，以概率\mathrm{1}-p服从泊松分布，对ZIP分布中的零膨胀率进行协变量回归建模得到了ZIP回归模型。Yau等^［3］的研究发现，ZINB分布对零膨胀数据的拟合效果更好，并考虑用ZINB分布进行建模，同时对零膨胀率进行协变量回归，提出了ZINB回归模型。目前，ZINB回归模型已在医学^［4-5］、生态学^［6］、社会学^［7-8］等领域得到广泛应用。

随着研究的不断深入，ZINB回归模型和ZIP回归模型也被不断拓展。Weng等^［9］研究了一种基于最大似然回归树的ZINB回归模型，用于预测航运事故死亡率，同时检验影响航运事故中人员生命损失的相关因素。Cho等^［10］提出了一种灵活的双变量ZINB回归模型，用于研究存在两个相关计数随机变量的联合分布，以满足生物医学研究的复杂分析需求。胡亚南等^［11］基于ZIP分布研究了零膨胀计数数据的联合建模及变量选择问题。王芝皓等^［12］使用函数型变系数零膨胀模型来研究零膨胀计数数据。这些研究主要关注截面零膨胀数据，数据来自于不同个体的单次观察。

近年来，在生物和医学等领域，纵向零膨胀计数数据逐渐得到重视^［13-14］。纵向计数数据是在多个观测时刻对一系列试验个体进行重复测量所获得的计数数据。对于每个个体，纵向计数数据记录其在不同时刻的测量结果。若随着时间变化对试验个体进行观测所得到的计数数据为零膨胀计数数据，则称其为纵向零膨胀计数数据。相较截面数据，纵向零膨胀计数数据引入了时间维度，数据特征（零膨胀率、离散度和均值等）可能随时间变化，且存在个体差异。Samani等^［15］提出了一种零膨胀幂级数回归模型。Zhu等^［16］利用ZIP分布和ZINB分布假设对均值、零膨胀率建立联合模型，并在其中引入了随机效应。Zhu等^［17］研究了一种新的可用于处理带缺失数据的纵向零膨胀计数数据的半参数方法。Wang等^［18］提出了一个两阶段混合效应模型。Zhang等^［19］利用Copula函数对纵向零膨胀数据进行了均值-相关性回归分析。

现有的纵向零膨胀计数数据回归分析多关注均值、零膨胀率随时间和协变量的变化，并假设离散度是固定的，很少考虑离散度随时间与某些协变量变化的情形。实际上，在许多纵向零膨胀计数数据问题中，离散度会随着时间与某些协变量变化。例如，随哮喘发病率上升或医疗资源投入增加，某地区居民哮喘的年治疗次数的离散度也会发生变化。因此，为了更好地适应数据的复杂性，扩大纵向零膨胀计数数据回归分析方法的应用场景，本文假定纵向零膨胀计数数据的均值、零膨胀率、离散度均受到时间与其他协变量的影响，并假设数据服从ZINB分布。本文对数据的均值、离散度和零膨胀率等参数进行广义线性回归建模，并设计了EM算法同时对各参数的广义线性回归系数进行极大似然估计。最后，本文通过模拟实验分析了模型的性能。

1 联合回归模型

零膨胀计数数据的零值超过标准计数分布（例如泊松分布、负二项（Negative Binomial， NB）分布）能够解释的范围，本文用退化零分布与NB分布的混合分布（即ZINB分布）来描述零膨胀计数数据Y，

P(Y=y)=\left\{\begin{array}{ll}\phi +(\mathrm{1}-\phi )f(\mathrm{0};\mu ,\alpha ),& y=\mathrm{0},\\ (\mathrm{1}-\phi )f(y;\mu ,\alpha ),& y\mathrm{0}\end{array}\right.（1）

其中，\phi为零膨胀率，即数据以概率\phi来自退化零分布（数据为结构零），以概率\mathrm{1}-\phi来自NB分布； \mu \mathrm{0}为均值； \alpha \mathrm{0}为离散度； f(·)为NB概率密度函数^［2］，

\begin{array}{l}f(y;\mu ,\alpha )=\frac{\mathrm{\Gamma }(y+\frac{\mathrm{1}}{\alpha })}{\mathrm{\Gamma }(y+\mathrm{1})\mathrm{\Gamma }(\frac{\mathrm{1}}{\alpha })}{(\frac{\mathrm{1}}{\mathrm{1}+\alpha \mu })}^{\frac{\mathrm{1}}{\alpha }},\\ y=\mathrm{0,1},\cdots \end{array}（2）

这里的\mathrm{\Gamma }(·)表示伽马函数。当\alpha \to \mathrm{0}时，ZINB分布退化为ZIP分布；当\phi =\mathrm{0}时，ZINB分布退化为NB分布。

1.1 极大似然估计模型

记在时刻集{\mathit{t}}_i=(t_{i\mathrm{1}},t_{i\mathrm{2}},\cdots ,t_{i{n}_i}{)}^{\mathrm{T}}对第i个个体的观测为{\mathit{y}}_i=(y_{i\mathrm{1}},y_{i\mathrm{2}},\cdots ,y_{i{n}_i}{)}^{\mathrm{T}}，其中y_{ij}为在时刻t_{ij}对第i个个体的观测，i=\mathrm{1},\cdots ,m,j=\mathrm{1},\cdots ,n_i。设y_{ij}由服从ZINB分布的Y_{ij}产生，

P(Y_{ij}=y_{ij})=\left\{\begin{array}{cc}{\phi }_{ij}+(\mathrm{1}-{\phi }_{ij})f(\mathrm{0};{\mu }_{ij},{\alpha }_{ij})),& y_{ij}=\mathrm{0},\\ (\mathrm{1}-{\phi }_{ij})f(y_{ij};{\mu }_{ij},{\alpha }_{ij}),& y_{ij}\mathrm{0}\end{array}\right.（3）

其中，{\mu }_{ij}为均值，{\alpha }_{ij}表示均值和方差的离散程度，{\phi }_{ij}为数据的零膨胀率。记在时刻集{\mathit{t}}_{\mathrm{1}},{\mathit{t}}_{\mathrm{2}},\cdots ,{\mathit{t}}_m分别观测m个个体得到的总数据集为\mathit{Y}=\{y_{ij}\}。

纵向零膨胀计数数据引入了时间维度，数据特征（例如零膨胀率、离散度和均值等）随时间与某些协变量变化，且个体间存在差异。本文为均值{\mu }_{ij}，离散度{\alpha }_{ij}和零膨胀概率参数{\phi }_{ij}建立如下的广义线性回归模型：

\left\{\begin{array}{l}\mathrm{l}\mathrm{o}\mathrm{g}{\mu }_{ij}={\mathit{x}}_{ij}^{\mathrm{T}}\mathit{\beta },\\ \mathrm{l}\mathrm{o}\mathrm{g}{\alpha }_{ij}={\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\gamma },\\ \mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}{\phi }_{ij}=\mathrm{l}\mathrm{o}\mathrm{g}\frac{{\phi }_{ij}}{\mathrm{1}-{\phi }_{ij}}={\mathit{w}}_{ij}^{\mathrm{T}}\mathit{\tau }\end{array}\right.（4）

其中，参数\mathit{\beta }=({\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\cdots ,{\beta }_b{)}^{\mathrm{T}}、\mathit{\gamma }=({\gamma }_{\mathrm{1}},{\gamma }_{\mathrm{2}},\cdots ,{\gamma }_q{)}^{\mathrm{T}}和\mathit{\tau }=({\tau }_{\mathrm{1}},{\tau }_{\mathrm{2}},\cdots ,{\tau }_p{)}^{\mathrm{T}}为回归系数，{\mathit{x}}_{ij}、{\mathit{z}}_{ij}和{\mathit{w}}_{ij}分别为对应的b维、q维和p维协变量，\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(·)为Logit函数。联合回归模型为零膨胀纵向计数数据建立了基于3参数联合建模的零膨胀负二项回归模型，后文将其简称为ZINB-3模型。

为了用极大似然估计法对回归模型（4）中的回归系数\mathit{\beta }\mathrm{、}\mathit{\gamma }\mathrm{、}\mathit{\tau }进行估计，本文引入指示变量b_{ij}、b_{ij}=\mathrm{0}时表示y_{ij}\mathrm{0}；b_{ij}=\mathrm{1}时表示y_{ij}=\mathrm{0}。在ZINB分布假设下，观测数据的对数似然函数记为：

\begin{array}{l}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau };\mathit{Y})={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}{b}_{ij}\mathrm{l}\mathrm{o}\mathrm{g}[{\phi }_{ij}+\\ (\mathrm{1}-{\phi }_{ij})f(\mathrm{0};{\mu }_{ij},{\alpha }_{ij})]+{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-b_{ij})\mathrm{l}\mathrm{o}\mathrm{g}[(\mathrm{1}-{\phi }_{ij})f(y_{ij};{\mu }_{ij},{\alpha }_{ij})]\end{array}（5）

其中f(y_{ij};{\mu }_{ij},{\alpha }_{ij})由式（2）给出。由（4）式可得：

\left\{\begin{array}{l}{\mu }_{ij}={\mathrm{e}}^{{\mathit{x}}_{ij}^{\mathrm{T}}\mathit{\beta }},\\ {\alpha }_{ij}={\mathrm{e}}^{{\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\gamma }},\\ {\phi }_{ij}={\mathrm{e}}^{{{\mathit{w}}_{ij}}^{\mathrm{T}}\mathit{\tau }}/(\mathrm{1}+{\mathrm{e}}^{{{\mathit{w}}_{ij}}^{\mathrm{T}}\mathit{\tau }})\end{array}\right.（6）

即{\mu }_{ij},{\alpha }_{ij},{\phi }_{ij}是待定回归系数\mathit{\beta },\mathit{\gamma },\mathit{\tau }的函数。这样，回归系数\mathit{\beta },\mathit{\gamma },\mathit{\tau }估计的极大似然优化问题为：

\underset{\mathit{\beta },\mathit{\gamma },\mathit{\tau }}{\mathrm{m}\mathrm{a}\mathrm{x}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau };\mathit{Y})（7）

1.2 EM算法

对数似然函数（5）较为复杂，一般来说极大似然优化问题（7）也难以求解。为此，本文引入潜变量{\mathit{d}}_i=(d_{i\mathrm{1}},\cdots ,d_{i{n}_i})，其中d_{ij}=\mathrm{1}表示y_{ij}由退化零分布产生，d_{ij}=\mathrm{0}表示y_{ij}由NB分布产生。潜变量d_{ij}可看成来自一次试验成功概率为{\phi }_{ij}的两点分布的随机变量。所有观测数据\mathit{Y}和随机潜变量\mathit{D}=({\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\cdots ,{\mathit{d}}_m^{\mathrm{T}}{)}^{\mathrm{T}}构成完全数据。相应的对数似然函数（5）可以改写为：

\begin{array}{l}{l}_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau };\mathit{Y},\mathit{D})=\\ l_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau };\mathit{Y},\mathit{D})+l_{\mathrm{c}}(\mathit{\tau };\mathit{Y},\mathit{D})\end{array}（8）

其中，

l_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma };\mathit{Y},\mathit{D})={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-d_{ij})\mathrm{l}\mathrm{o}\mathrm{g}f(y_{ij};{\mu }_{ij},{\alpha }_{ij})（9）

\begin{array}{l}{l}_{\mathrm{c}}(\mathit{\tau };\mathit{Y},\mathit{D})={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}{d}_{ij}\mathrm{l}\mathrm{o}\mathrm{g}[\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{w}}_{ij}^{\mathrm{T}}\mathit{\tau })]-\\ {\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}\mathrm{l}\mathrm{o}\mathrm{g}[\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{w}}_{ij}^{\mathrm{T}}\mathit{\tau })]\end{array}（10）

本文使用EM算法^［20］求解如下完全数据的对数极大似然问题：

\underset{\mathit{\beta },\mathit{\gamma },\mathit{\tau }}{\mathrm{m}\mathrm{a}\mathrm{x}}{l}_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau };\mathit{Y},\mathit{D})（11）

通过E步和M步迭代交替求解，可以得到参数的估计\widehat{\mathit{\theta }}=({\widehat{\mathit{\beta }}}^{\mathrm{T}},{\widehat{\mathit{\gamma }}}^{\mathrm{T}},{\widehat{\mathit{\tau }}}^{\mathrm{T}}{)}^{\mathrm{T}}，其中 E步和M步的具体流程如下。

E步： 已知观测数据与M步计算得到的参数估计值\tilde{\mathit{\theta }}=({\tilde{\mathit{\beta }}}^{\mathrm{T}},{\tilde{\mathit{\gamma }}}^{\mathrm{T}},{\tilde{\mathit{\tau }}}^{\mathrm{T}}{)}^{\mathrm{T}}（初始迭代时给定参数的初始值），计算潜变量d_{ij}的条件期望，

\begin{array}{l}{\widehat{d}}_{ij}=E(d_{ij}|\tilde{\mathit{\beta }},\tilde{\mathit{\gamma }},\tilde{\mathit{\tau }};\mathit{Y})=\\ P(d_{ij}=\mathrm{1}|\tilde{\mathit{\beta }},\tilde{\mathit{\gamma }},\tilde{\mathit{\tau }};\mathit{Y})=\\ \frac{P(d_{ij}=\mathrm{1},Y_{ij}=y_{ij}|\tilde{\mathit{\theta }})}{P(Y_{ij}=y_{ij}|\tilde{\mathit{\theta }})}=\\ \frac{P(Y_{ij}=y_{ij}|d_{ij}=\mathrm{1},\tilde{\mathit{\theta }})P(d_{ij}=\mathrm{1}|\tilde{\mathit{\theta }})}{\sum _{k=\mathrm{0,1}}P(Y_{ij}=y_{ij}|d_{ij}=k,\tilde{\mathit{\theta }})P(d_{ij}=k|\tilde{\mathit{\theta }})}=\\ \frac{{\tilde{\phi }}_{ij}\mathrm{1}(y_{ij}=\mathrm{0})}{{\tilde{\phi }}_{ij}\mathrm{1}(y_{ij}=\mathrm{0})+(\mathrm{1}-{\tilde{\phi }}_{ij})P(Y_{ij}=y_{ij})}=\\ \frac{\mathrm{1}(y_{ij}=\mathrm{0})}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(-{\mathit{w}}_{ij}^{\mathrm{T}}\tilde{\mathit{\tau }})f(\mathrm{0};{\tilde{\mu }}_{ij},{\tilde{\alpha }}_{ij})}\end{array}（12）

其中{\tilde{\phi }}_{ij}\mathrm{、}{\tilde{\mu }}_{ij}\mathrm{、}{\tilde{\alpha }}_{ij}为M步的零膨胀概率、均值及离散度估计，\mathrm{1}(·)为示性函数。将似然函数（8）中的潜变量d_{ij}替换为它的条件期望{\widehat{d}}_{ij}，有

Q_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau })=Q_{\mathrm{1}}(\mathit{\beta },\mathit{\gamma })+Q_{\mathrm{2}}(\mathit{\tau })（13）

其中，

Q_{\mathrm{1}}(\mathit{\beta },\mathit{\gamma })={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-{\widehat{d}}_{ij})\mathrm{l}\mathrm{o}\mathrm{g}f(y_{ij};{\mu }_{ij},{\alpha }_{ij})（14）

\begin{array}{l}{Q}_{\mathrm{2}}(\mathit{\tau })={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}{\widehat{d}}_{ij}\mathrm{l}\mathrm{o}\mathrm{g}[\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{w}}_{ij}^{\mathrm{T}}\mathit{\tau })]-\\ {\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}\mathrm{l}\mathrm{o}\mathrm{g}[\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{w}}_{ij}^{\mathrm{T}}\mathit{\tau })]\end{array}（15）

M步： 本文利用牛顿迭代法求解优化问题：

\tilde{\mathit{\theta }}=\underset{\mathit{\beta },\mathit{\gamma },\mathit{\tau }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}{Q}_{\mathrm{c}}(\mathit{\beta },\mathit{\gamma },\mathit{\tau })（16）

其中，\tilde{\mathit{\theta }}=({\tilde{\mathit{\beta }}}^{\mathrm{T}},{\tilde{\mathit{\gamma }}}^{\mathrm{T}},{\tilde{\mathit{\tau }}}^{\mathrm{T}}{)}^{\mathrm{T}}为M步给出的参数估计。牛顿迭代法求解的具体流程如下。

1） 计算Q_{\mathrm{1}}(\mathit{\beta },\mathit{\gamma })\mathrm{、}{Q}_{\mathrm{2}}(\mathit{\tau })的1阶偏导S_{\mathit{\beta }}\mathrm{、}{S}_{\mathit{\gamma }}\mathrm{、}{S}_{\mathit{\tau }}和2阶偏导（Hessian矩阵）H_{\mathit{\beta }}\mathrm{、}{H}_{\mathit{\gamma }}\mathrm{、}{H}_{\mathit{\tau }}。Q_{\mathrm{1}}(\mathit{\beta },\mathit{\gamma })关于\mathit{\beta }的1阶和2阶偏导数为：

\begin{array}{l}{S}_{\mathit{\beta }}={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-{\widehat{d}}_{ij})\left(y_{ij}{\mathit{x}}_{ij}-\frac{{\mu }_{ij}{\mathit{x}}_{ij}}{\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}}-\frac{y_{ij}{\mu }_{ij}{\alpha }_{ij}{\mathit{x}}_{ij}}{\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}}\right),\\ H_{\mathit{\beta }}={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-{\widehat{d}}_{ij})\left[\frac{-{\mu }_{ij}{\mathit{x}}_{ij}{\mathit{x}}_{ij}^{\mathrm{T}}}{(\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}{)}^{\mathrm{2}}}-\frac{y_{ij}{\mu }_{ij}{\alpha }_{ij}{\mathit{x}}_{ij}{\mathit{x}}_{ij}^{\mathrm{T}}}{(\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}{)}^{\mathrm{2}}}\right],\end{array}

Q_{\mathrm{1}}(\mathit{\beta },\mathit{\gamma })关于\mathit{\gamma }的1阶和2阶偏导数为：

\begin{array}{l}{S}_{\mathit{\gamma }}={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-{\widehat{d}}_{ij})\cdot \\ [-{\mathrm{\Gamma }}^{\text{'}}(y_{ij}+{\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{1}}{\mathit{z}}_{ij}+{\mathrm{\Gamma }}^{\text{'}}({\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{1}}{\mathit{z}}_{ij}+\\ {\alpha }_{ij}^{-\mathrm{1}}{\mathit{z}}_{ij}\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+{\mu }_{ij}{\alpha }_{ij})-\frac{{\mu }_{ij}{\mathit{z}}_{ij}}{\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}}+\\ y_{ij}{\mathit{z}}_{ij}-\frac{y_{ij}{\mu }_{ij}{\alpha }_{ij}{\mathit{z}}_{ij}}{\mathrm{1}+{\mu }_{ij}{\alpha }_{ij}}],\end{array}

\begin{array}{l}{H}_{\mathit{\gamma }}={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}(\mathrm{1}-{\widehat{d}}_{ij})\cdot \\ [{\mathrm{\Gamma }}^{″}(y_{ij}+{\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{2}}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}+{\mathrm{\Gamma }}^{\text{'}}(y_{ij}+\\ {\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{1}}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}-\\ {\mathrm{\Gamma }}^{″}({\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{2}}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}-{\mathrm{\Gamma }}^{\text{'}}({\alpha }_{ij}^{-\mathrm{1}}){\alpha }_{ij}^{-\mathrm{1}}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}-\\ {\alpha }_{ij}^{-\mathrm{1}}\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+{\alpha }_{ij}{\mu }_{ij}){\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}+\frac{{\mu }_{ij}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}}{\mathrm{1}+{\alpha }_{ij}{\mu }_{ij}}+\\ \frac{{\mu }_{ij}{\alpha }_{ij}{\mu }_{ij}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}-y_{ij}{\alpha }_{ij}{\mu }_{ij}{\mathit{z}}_{ij}{\mathit{z}}_{ij}^{\mathrm{T}}}{(\mathrm{1}+{\alpha }_{ij}{\mu }_{ij}{)}^{\mathrm{2}}}],\end{array}

Q_{\mathrm{2}}(\mathit{\tau })关于\mathit{\tau }的1阶和2阶偏导数为：

\begin{array}{l}{S}_{\tau }={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}({\widehat{d}}_{ij}{\mathit{w}}_{ij}-{\phi }_{ij}{\mathit{w}}_{ij}),\\ H_{\mathit{\tau }}={\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^{n_i}}\frac{{\mathit{w}}_{ij}{\mathit{w}}_{ij}^{\mathrm{T}}{e}^{-{\mathit{w}}_{ij}^{\mathrm{T}}\tau }}{[\mathrm{1}+e^{-{\mathit{w}}_{ij}^{\mathrm{T}}\tau }{]}^{\mathrm{2}}},\end{array}

其中，式（6）中的{\mu }_{ij}\mathrm{、}{\alpha }_{ij}\mathrm{、}{\phi }_{ij}是待定回归系数\mathit{\beta }\mathrm{、}\mathit{\gamma }\mathrm{、}\mathit{\tau }的函数，{\mathrm{\Gamma }}^{\text{'}}(·)和{\mathrm{\Gamma }}^{″}(·)分别表示伽马函数的1阶和2阶导数。

2） 使用牛顿迭代法更新参数，

\left\{\begin{array}{c}{\mathit{\beta }}^{(t+\mathrm{1})}={\mathit{\beta }}^{(t)}-{\mathit{H}}_{{\mathit{\beta }}^{(t)}}^{-\mathrm{1}}{S}_{{\mathit{\beta }}^{(t)}},\\ {\mathit{\gamma }}^{(t+\mathrm{1})}={\mathit{\gamma }}^{(t)}-{\mathit{H}}_{{\mathit{\gamma }}^{(t)}}^{-\mathrm{1}}{S}_{{\mathit{\gamma }}^{(t)}},\\ {\mathit{\tau }}^{(t+\mathrm{1})}={\mathit{\tau }}^{(t)}-{\mathit{H}}_{{\mathit{\tau }}^{(t)}}^{-\mathrm{1}}{S}_{{\mathit{\tau }}^{(t)}}\end{array}\right.（17）

其中，{(·)}^{(t)}为第t次迭代的参数估计值，{\mathit{H}}_{(·)}^{-\mathrm{1}}为Hessian矩阵的逆。

进行参数估计时，若所选取的参数初值与真实值偏差过大，很可能导致迭代次数增加，加大计算时间成本。本文用参数\mathit{\beta }\mathrm{、}\mathit{\tau }的传统估计作为初值，这些估计可利用R软件的pscl包^［21］中的zeroinfl函数得到。对于\mathit{\gamma }的初值，本文首先基于纵向数据均值与方差的关系^［22］反推离散参数{\alpha }_{ij}，即：

{\alpha }_{ij}=[D(y_{ij})-E(y_{ij})]/[{\mu }_{ij}^{\mathrm{2}}(\mathrm{1}-{\phi }_{ij})]-{\phi }_{ij}，

然后用模型\mathrm{l}\mathrm{o}\mathrm{g}{\alpha }_{ij}={\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\gamma }的回归系数的最小二乘估计作为\mathit{\gamma }的初始值。

2 回归系数估计的统计性质

通过EM算法求解极大似然问题（11）得到参数\mathit{\theta }=({\mathit{\beta }}^{\mathrm{T}},{\mathit{\gamma }}^{\mathrm{T}},{\mathit{\tau }}^{\mathrm{T}})的似然估计后，记该参数似然估计为\widehat{\mathit{\theta }}=({\widehat{\mathit{\beta }}}^{\mathrm{T}},{\widehat{\mathit{\gamma }}}^{\mathrm{T}},{\widehat{\mathit{\tau }}}^{\mathrm{T}}{)}^{\mathrm{T}}。本节基于以下的假设来讨论\widehat{\mathit{\theta }}的渐近性质：

（H1） \mathit{\theta }所在的参数空间\mathrm{\Theta }是R^{b+q+p}中的紧子集；

（H2） \mathit{\theta }的真实值{\mathit{\theta }}_{\mathrm{0}}在参数空间\mathrm{\Theta }的内部；

（H3） 观测y_{ij}独立同分布且服从ZINB分布，对数似然函数l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})满足：

E(\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})}{\partial \mathit{\theta }})=\mathrm{0}（18）

Fisher信息矩阵I(\mathit{\theta })=E(-\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}})正定有界且满足：

\frac{\mathrm{1}}{\sqrt[]{m}}\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})}{\partial \mathit{\theta }}{\left|\right.}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\stackrel{d}{\to }N(\mathrm{0},\overline{I}({\mathit{\theta }}_{\mathrm{0}}))（19）

其中\overline{I}({\mathit{\theta }}_{\mathrm{0}})=\frac{\mathrm{1}}{m}I(\mathit{\theta })。

（H4） 假设\mathrm{𝒩}为{\mathit{\theta }}_{\mathrm{0}}的邻域，存在函数M(\mathit{Y})，使得对于\mathit{\theta }\in \mathrm{𝒩}有

\left|\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})}{\partial {\theta }_j\partial {\theta }_c\partial {\theta }_k}\right|\le M(\mathit{Y})（20）

其中， {\theta }_j\mathrm{、}{\theta }_c\mathrm{、}{\theta }_k为参数\mathit{\theta }的满足E_{{\mathit{\theta }}_{\mathrm{0}}}(M(\mathit{Y}))\mathrm{\infty }的任意分量。

定理2.1 假设（H1）~（H4）成立。则极大似然估计\widehat{\mathit{\theta }}是真实参数{\mathit{\theta }}_{\mathrm{0}}的\sqrt[]{m}-相合估计，

\Vert \widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}{\Vert }_{\mathrm{2}}=O_p(m^{-\frac{\mathrm{1}}{\mathrm{2}}})（21）

且当m\to \mathrm{\infty }时服从渐近正态分布，

\sqrt[]{m}(\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}})\stackrel{d}{\to }N(\mathrm{0},\overline{I}({\mathit{\theta }}_{\mathrm{0}}{)}^{-\mathrm{1}})（22）

证明 首先证明相合性。为记号方便，将l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta };\mathit{Y})简写为l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })。相合性证明需要证明对任意给定的\epsilon \mathrm{0}，存在常数M_{\epsilon }\mathrm{0}，使得：

\begin{array}{l}\underset{m\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}P\{\underset{\Vert \mathit{v}{\Vert }_{\mathrm{2}}=M_{\epsilon }}{\mathrm{s}\mathrm{u}\mathrm{p}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+m^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v})\\ l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})\}\ge \mathrm{1}-\epsilon \end{array}（23）

从而当m\to \mathrm{\infty }时在球体\{{\mathit{\theta }}_{\mathrm{0}}+m^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v}:\Vert \mathit{v}{\Vert }_{\mathrm{2}}\le M_{\epsilon }\}中存在局部极大似然估计\widehat{\mathit{\theta }}，使得\Vert \widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}{\Vert }_{\mathrm{2}}=O_p(m^{-\frac{\mathrm{1}}{\mathrm{2}}})大概率成立。

为证明式（23），对似然函数进行泰勒展开，得：

\begin{array}{l}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+m^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v})-l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})=\\ \nabla l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}{)}^{\mathrm{T}}{m}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v}-\\ \left(-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}^{\mathrm{T}}\left({\nabla }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})/m\right)\mathit{v}\right)+{\mathit{v}}^{\mathrm{T}}{o}_p(\mathrm{1})\mathit{v},\end{array}

其中，

\begin{array}{l}\nabla l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})=\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}},\\ {\nabla }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})=\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\end{array}。

根据给定的正则化条件（H3），当m\to \mathrm{\infty }时，以下结果成立：

\left\{\begin{array}{l}\nabla l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})=O_p(\sqrt[]{m}),\\ -{\nabla }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})/m=\overline{I}({\mathit{\theta }}_{\mathrm{0}})+o_p(\mathrm{1}),\\ \overline{I}({\mathit{\theta }}_{\mathrm{0}})=\mathrm{E}(-{\nabla }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})/m)=O_p(\mathrm{1})\end{array}\right.（24）

从而当m\to \mathrm{\infty }时，l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+m^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v})-l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})的展开项主要由-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}^{\mathrm{T}}\{\overline{I}({\mathit{\theta }}_{\mathrm{0}})+o_p(\mathrm{1})\}\mathit{v}控制。因此，对任意给定的\epsilon \mathrm{0}，存在一个充分大的常数M_{\epsilon }\mathrm{0}，使得：

P\{\underset{\Vert \mathit{v}{\Vert }_{\mathrm{2}}=M_{\epsilon }}{\mathrm{s}\mathrm{u}\mathrm{p}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+m^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathit{v})-l_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}})\mathrm{0}\}\ge \mathrm{1}-\epsilon \mathrm{。}

综上，\widehat{\mathit{\theta }}是{\mathit{\theta }}_{\mathrm{0}}的\sqrt[]{m}相合估计。

接下来证明渐近正态性。已知\widehat{\mathit{\theta }}是{\mathit{\theta }}_{\mathrm{0}}的\sqrt[]{m}-相合估计，则\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}=\mathrm{0}\mathrm{。} 由泰勒展开，得：

\begin{array}{l}\mathrm{0}=\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{\left|\right.}_{\mathit{\theta }=\widehat{\mathit{\theta }}}=\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{\left|\right.}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}-\\ [-\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}{\left|\right.}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\times (\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}})]+R({\mathit{\theta }}^{\mathrm{*}})\end{array}（25）

其中{\mathit{\theta }}^{\mathrm{*}}在{\mathit{\theta }}_{\mathrm{0}}和\widehat{\mathit{\theta }}之间，

R({\mathit{\theta }}^{\mathrm{*}})=\frac{\mathrm{1}}{\mathrm{2}}{\nabla }^{\mathrm{3}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}({\mathit{\theta }}^{\mathrm{*}})(\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}{)}^{\mathrm{T}}(\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}})

为展开式的余项。由于\Vert \widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}{\Vert }_{\mathrm{2}}=O_p(m^{-\frac{\mathrm{1}}{\mathrm{2}}})，结合给定的正则化条件（H4）可得R({\mathit{\theta }}^{\mathrm{*}})=O_p(\mathrm{1})。由式（25）有

\begin{array}{l}\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}={\left[-\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\right]}^{-\mathrm{1}}\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}+\\ {\left[-\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\right]}^{-\mathrm{1}}R({\mathit{\theta }}^{\mathrm{*}})\mathrm{。}\end{array}

根据弱大数定律，可以证明：

-\frac{\mathrm{1}}{m}\frac{{\partial }^{\mathrm{2}}{l}_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}\stackrel{p}{\to }\overline{I}({\mathit{\theta }}_{\mathrm{0}})\mathrm{。}

因此，当m\to \mathrm{\infty }时，结合Slutsky定理，有

\begin{array}{l}\sqrt[]{m}(\widehat{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}})=\overline{I}({\mathit{\theta }}_{\mathrm{0}}{)}^{-\mathrm{1}}\frac{\mathrm{1}}{\sqrt[]{m}}\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}+\\ m^{-\frac{\mathrm{1}}{\mathrm{2}}}\overline{I}({\mathit{\theta }}_{\mathrm{0}}{)}^{-\mathrm{1}}R({\mathit{\theta }}^{\mathrm{*}})=\\ \overline{I}({\mathit{\theta }}_{\mathrm{0}}{)}^{-\mathrm{1}}\frac{\mathrm{1}}{\sqrt[]{m}}\frac{\partial l_{\mathrm{o}\mathrm{b}\mathrm{s}}(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}}+\\ O_p(m^{-\frac{\mathrm{1}}{\mathrm{2}}})\stackrel{d}{\to }N(\mathrm{0},\overline{I}({\mathit{\theta }}_{\mathrm{0}}{)}^{-\mathrm{1}})\mathrm{。}\end{array}

证毕。

3 模拟实验

为验证方法性能，本文随机生成200个服从ZINB分布的独立样本，每个样本的观察时刻为{\mathit{t}}_i=(t_{ij}{)}_{j=\mathrm{1}}^{\mathrm{10}}={(\mathrm{0.5,1},\mathrm{1.5},\cdots ,\mathrm{5})}^{\mathrm{T}}。考虑协变量：

{\mathit{x}}_{ij}=(\mathrm{1},t_{ij},t_{ij}^{\mathrm{2}},\cdots ,t_{ij}^{p-\mathrm{1}}{)}^{\mathrm{T}},

{\mathit{z}}_{ij}=(\mathrm{1},t_{ij},t_{ij}^{\mathrm{2}},\cdots ,t_{ij}^{b-\mathrm{1}}{)}^{\mathrm{T}},

{\mathit{w}}_{ij}=(\mathrm{1},t_{ij},h_i{)}^{\mathrm{T}}，

其中，h_i为0-1变量，h_i=\mathrm{0}表示第i个个体接受治疗，h_i=\mathrm{1}表示第i个个体不接受治疗。设p=\mathrm{3}，b=\mathrm{4}。考虑估计方法在两组不同零膨胀率和离散度场景中的回归系数估计问题：

场景1 零膨胀率为20%，协变量回归系数：

\begin{array}{l}\mathit{\beta }={(\mathrm{3,0.5,0.05,0.001})}^{\mathrm{T}},\\ \mathit{\gamma }={(-\mathrm{0.5},-\mathrm{0.3},-\mathrm{0.05})}^{\mathrm{T}},\\ \mathit{\tau }={(-\mathrm{0.2},-\mathrm{0.5,0.1})}^{\mathrm{T}};\end{array}

场景2 零膨胀率为70%，协变量回归系数：

\begin{array}{l}\mathit{\beta }={(\mathrm{3,0.5,0.05,0.001})}^{\mathrm{T}},\\ \mathit{\gamma }={(-\mathrm{0.5},-\mathrm{0.3},-\mathrm{0.05})}^{\mathrm{T}},\\ \mathit{\tau }={(-\mathrm{0.8,0.7,0.1})}^{\mathrm{T}}\mathrm{。}\end{array}

为了评估ZINB-3模型在不同场景中的表现，本文将其与ZIP模型^［2］与NB模型^［23］进行比较。ZIP模型假设数据服从ZIP分布，对零膨胀率进行广义回归建模，并给出回归系数的极大似然估计。NB模型假设数据服从NB分布，对均值和离散度进行广义回归建模，并给出回归系数的极大似然估计。

表1给出了3种模型的参数估计的平均绝对误差（Mean Absolute Error， MAE）值和相应的标准差，所有结果均由500次重复实验得到，参数\mathit{\beta }的MAE计算公式为：

S_{\mathrm{M}\mathrm{A}\mathrm{E},\mathit{\beta }}=\frac{\mathrm{1}}{b}{\displaystyle \sum _{k=\mathrm{1}}^b}|{\widehat{\beta }}_k-{\beta }_k|（26）

其中，{\beta }_k与{\widehat{\beta }}_k分别是向量\mathit{\beta }与参数估计\widehat{\mathit{\beta }}的第k个分量，参数\mathit{\gamma }和\mathit{\tau }的MAE计算公式类似可得。

由表1可知，相较于ZIP和NB模型，ZINB-3模型的估计误差更小，精度更高，说明均值、离散度和零膨胀率的联合建模能够同时提升模型对各参数的估计精度。此外，相较于ZIP和NB模型，ZINB-3模型的均值系数估计的标准差较小，说明均值、离散度和零膨胀率联合建模能够提升均值模型估计的有效性。

基于ZINB-3模型的估计结果可以得到均值参数和离散度的变化趋势。为了直观展示均值和离散度的变化，本文以零膨胀率为20%情况为例，分别在图1a、b中展示了均值和离散度的变化趋势。可以看到，数据的均值呈上升趋势，而离散度则持续下降，表明随着时间推移，尽管响应变量的均值不断增加，但数据的分散程度却在逐渐减小。此外，均值和离散度的置信带都比较窄，表明本文的模型是有效的。

4 结论

为了解释纵向零膨胀计数数据的均值、离散度和零膨胀率对时间与协变量的依赖关系，本文对这3个参数进行了广义线性回归联合建模，提出了基于EM算法的回归系数极大似然估计方法。模拟实验结果表明，相较于两参数回归模型，本文的模型具有更好的估计效果，能够更好地分析纵向零膨胀计数数据。

参考文献

[1] Mullahy J.Specification and testing of some modified count data models ［J］.J Econom， 1986， 33（3）： 341-365.

[2] Lambert D.Zero-inflated Poisson regression， with an application to defects in manufacturing ［J］.Technometrics， 1992， 34（1）： 1-14.

[3] Yau K K W， Wang K， Lee A H.Zero-inflated negative binomial mixed regression modeling of over-dispersed count data with extra zeros ［J］.Biom J， 2003， 45（4）： 437-452.

[4] Fontanive V N， Broch B， Seoane M， et al.Oral impact of dental prosthesis use in different patterns of tooth loss ［J］.J Prosthet Dent， 2025， 134（5）： 1781-1786.

[5] Cho H， Liu C， Preisser J S， et al.A bivariate zero-inflated negative binomial model and its applications to biomedical settings ［J］.Stat Methods Med Res， 2023， 32（7）： 1300-1317.

[6] Ozor J， Ronde P， Tung S， et al.The “Middleground” as a catalyst for the dynamics of innovation in an ecosystem？ The case of Eurasanté in Hauts-de-France ［J］.Technol Forecasting Social Change， 2024， 209： 123731.

[7] Caniglia M.Economic insecurity in young adulthood among children maltreated during early adolescence ［J］.Child Abuse Neglect， 2024， 157： 107047.

[8] Kang K I， Kang K， Kim C.Risk factors influencing cyberbullying perpetration among middle school students in Korea： Analysis using the zero-inflated negative binomial regression model ［J］.Int J Environ Res Public Health， 2021， 18（5）： 2224.

[9] Weng J， Yang D， Qian T， et al.Combining zero-inflated negative binomial regression with MLRT techniques： An approach to evaluating ship accident casualties ［J］.Ocean Eng， 2018， 166： 135-144.

[10] Cho H， Liu C， Preisser J S， et al.A bivariate zero-inflated negative binomial model and its applications to biomedical settings ［J］.Stat Methods Med Res， 2023， 32（7）： 1300-1317.

[11] Hu Y N， Tian M Z.Joint modeling and variable selection from zero-inflated count data ［J］.Statistical Research， 2019， 36（1）： 104-114.

[12] 胡亚南， 田茂再.零膨胀计数数据的联合建模及变量选择［J］.统计研究， 2019， 36（1）： 104-114.

[13] Wang Z H， Liu Y X， Tian M Z， et al.Functional partially varying coefficient models for zero-inflated count data ［J］.Statistical Research， 2021， 38（7）： 127-139.

[14] 王芝皓， 刘艳霞， 田茂再， 等.零膨胀计数数据函数型部分变系数模型［J］.统计研究， 2021， 38（7）： 127-139.

[15] Alfò M， Maruotti A.Two-part regression models for longitudinal zero-inflated count data ［J］.Can J Stat， 2010， 38（2）： 197-216.

[16] Atkins D C， Baldwin S A， Zheng C， et al.A tutorial on count regression and zero-altered count models for longitudinal substance use data ［J］.Psychol Addict Behav， 2013， 27（1）： 166-177.

[17] Samani E B， Amirian Y， Ganjali M.Likelihood estimation for longitudinal zero-inflated power series regression models ［J］.J Appl Stat， 2012， 39（9）： 1965-1974.

[18] Zhu H， Luo S， DeSantis S M.Zero-inflated count models for longitudinal measurements with heterogeneous random effects ［J］.Stat Methods Med Res， 2017， 26（4）： 1774-1786.

[19] Zhu H， DeSantis S M， Luo S.Joint modeling of longitudinal zero-inflated count and time-to-event data： A Bayesian perspective ［J］.Stat Methods Med Res， 2018， 27（4）： 1258-1270.

[20] Wang J， Reyes-Gibby C， Shete S.An approach to analyze longitudinal zero-inflated microbiome count data using two-stage mixed effects models ［J］.Stat Biosci， 2021， 13： 267-290.

[21] Zhang W， Qian J W F， Chen Y.A joint mean-correlation modeling approach for longitudinal zero-inflated count data ［J］.Braz J Probab Stat， 2020， 34（1）： 35-50.

[22] McLachlan G J， Krishnan T.The EM algorithm and extensions ［M］.Hoboken： John Wiley Sons， 2008.

[23] Zeileis A， Kleiber C， Jackman S.Regression models for count data in R ［J］.J Stat Softw， 2008， 27（1）： 1-25.

[24] Cameron A C， Pravin K T.Regression analysis of count data ［M］.Cambridge： Cambridge University Press， 2013.

[25] Hilbe J M.Negative binomial regression ［M］.Cambridge： Cambridge University Press， 2011.

基金资助

重庆市教委科学技术研究项目(KJQN202301156)

重庆理工大学研究生教育高质量发展项目(gzlcx20245278)

四川省科技计划重点研发项目(2025YFHZ0025)