贝叶斯半参数回归解析：多分类时空数据的建模新思路

Bayesian Semiparametric Regression Analysis of Multicategorical Time-Space Data

贝叶斯半参数回归解析：多分类时空数据的建模新思路

https://epub.ub.uni-muenchen.de/1592/1/paper_202.pdf

摘要
我们提出了一种基于马尔可夫随机场先验的统一半参数贝叶斯方法，用于分析多类别响应变量与时间、空间及其他协变量之间的依赖关系。该一般模型通过引入空间效应以及以灵活的半参数形式表示的度量协变量的非线性效应，扩展了用于分类时间序列和纵向数据的动态模型（或状态空间模型）。趋势成分、季节性成分、不同类型协变量以及空间效应都在同一通用框架下进行处理，只需为它们分配具有不同形式和平滑程度的适当先验分布。推理过程完全基于贝叶斯方法，并使用MCMC（马尔可夫链蒙特卡罗）技术进行后验分析。我们提供了两种方法：第一种方法基于对观测似然函数的直接评估；第二种方法基于潜在半参数效用模型，特别适用于Probit模型。这些方法通过失业数据和森林损害调查的应用案例进行了说明。
关键词 ：分类时空数据、森林损害、马尔可夫随机场、MCMC、半参数贝叶斯推理、失业。

1 引言

多类别纵向数据由观测值 (Yit; xit)，i = 1, ..., n，t = 1, ..., T 组成，表示对一组 n 个个体随时间的观察结果，其中响应变量 Y 被观测为有序或无序的类别 r ∈ {1, ..., k}。协变量可以是时间恒定的，也可以是随时间变化的。当时间点数 T 相对于 n 较小时，广义估计方程（GEE）方法是数据分析中的常用选择。而对于中等或较大的 T，动态模型或状态空间模型则是一个有用的替代方案，参见例如 Fahrmeir 和 Tutz（1994, 2000, 第8章）。在特殊情况下（n=1），即分类时间序列，动态广义线性模型已成为进行近似或完整贝叶斯推理的一种成熟工具。

在本文中，我们考虑的是多类别的时空数据，其中每个个体的空间位置或站点 s 在一个空间网格 {1, ..., s, ..., S} 上给出，作为额外的信息。我们还区分了度量型协变量 xt = (xt1, ..., xtp)'，其效应将采用非参数方式建模和估计；以及另一个协变量向量 wt，其效应将以通常的线性形式进行参数化建模。由此，n 个个体或单位的多类别时空数据包括如下观测值：

一个典型的例子是德国劳动局1980年至1995年间的月度登记数据 ，其中 Yit 表示个体 i 在月份 t 的就业状况（例如：失业、兼职、全职），而 si 是个体 i 所在的德国行政区。森林健康调查数据是另一个例子：在年度调查中，树木 i 在年份 t 的健康状态 Yit 通过其叶面积损失程度来表示，该状态被划分为有序类别（从严重受损到无损），si 则表示该树木在一个格网地图上的位置。在这两个例子中，协变量可以是分类变量或连续变量，也可能随时间变化。

一般来说，这种类型的时空数据无法通过现有的非参数或传统参数方法进行充分分析。我们提出了一种统一的半参数贝叶斯框架，用于联合建模和分析时间、空间以及不同类型协变量对分类响应的影响。趋势成分、季节性成分、空间效应、具有非线性效应的度量型协变量以及具有固定效应的普通协变量，都在同一个通用框架下进行处理，只需为它们分配具有不同形式和平滑程度的适当先验分布。这一广泛的模型类别包含了以往研究中所考虑的用于分类时间序列的状态空间模型作为特例。推理过程完全基于贝叶斯方法，并使用了最新的MCMC（马尔可夫链蒙特卡罗）技术。我们提出了两种MCMC推理方法：

第一种方法适用于在给定协变量和未知参数的情况下，数据的似然函数容易计算的情形，如累积模型或多项式Logistic模型。此时可以通过 Fahrmeir 和 Lang（1999）为单变量响应开发的Metropolis-Hastings算法的扩展来生成马尔可夫链样本。

第二种方法基于潜在变量，其中可观测的分类响应通过阈值机制或效用机制生成。对于潜在高斯变量，这将导致多类别Probit模型。参见 Albert 和 Chib（1993）关于线性预测器的简化情形，以及 Yau、Kohn 和 Wong（2000）关于使用基函数进行非参数回归的情形。在MCMC推理中，将高斯潜在变量视为未知的附加“参数”，并在Gibbs抽样方案中与其它参数一起生成。高效地从高维高斯马尔可夫随机场中抽样的方法被整合进来，作为主要构建模块。

第2节描述了我们针对在时间和空间上观测到的分类响应变量所提出的贝叶斯半参数回归模型，这些响应变量依赖于未知函数和参数。第3节介绍了MCMC算法。第4节将这些方法应用于基于（非）就业状态的分类时空数据对再就业机会的研究，以及森林健康调查数据的分析。

2 多类别时空数据的半参数贝叶斯模型

分类响应模型可以从潜在变量的角度来构建。这种方法不仅有助于模型的建立，也有助于贝叶斯推理，即将潜在变量视为额外的未知“参数”。

对于具有无序类别 1, ..., k 的名义响应 Y 的情况，设 Ur 是与第 r 个类别相关联的潜在变量或效用值。假设 Ur 满足：

最后，我们注意到，在公式 (2.6) 中，某个特定协变量的效应可能只对响应变量的某些类别 k 存在。此外，我们可能会观察到一些只与响应变量的某一特定类别相关的协变量，即所谓的类别特定协变量。这导致了类别特定协变量和全局协变量之间的区别。尽管对公式 (2.6) 进行修改并估计这类模型是直接的（并且已经实现），但在此我们不深入讨论，因为在这篇论文的应用中，包含类别特定协变量并不是必要的。更多细节可以在 Fahrmeir 和 Tutz (1994, 2000, 第3章) 及其中的参考文献中找到。

其中，Ns 表示相邻区域的数量，记号 u ∈ ∂s 表示区域 u 是区域 s 的邻居。因此，fstr(s) 的（条件）均值是相邻区域的 fstr(u) 值的平均。同样地，方差 τ_str² 控制着平滑程度。这一先验将在我们关于失业持续时间的第一个应用中使用。在某些其他应用中，如我们关于森林损害数据的第二个例子中，采用一种更一般的先验设定似乎更为合适。在此应用中，我们假设若两个站点 si 和 sj 之间的距离在某个特定距离 d 之内，则它们是相邻的。此外，我们现在假设 fstr(s) 的条件均值是相邻站点 fstr(u) 值的一个加权平均，而不是像式（2.11）中那样的未加权平均。权重被设定为与相邻站点到站点 s 的距离成反比。以 wsu 表示权重，一般的空间先验可以定义如下：

与自回归先验类似，(2.12) 可以写成形式 (2.10)，其中罚矩阵 K的元素由以下给出：

3 通过 MCMC 进行后验分析

以下中，f表示所有函数评估值的向量，包括时间 t的趋势项和季节项，以及结构化和非结构化的空间效应；τ表示所有方差的向量；β=β对于名义模型而言，而对于有序模型则 β=(β,ξ)。

对于名义逻辑回归模型（nominal logit model）或累积逻辑回归模型（cumulative logit model），给定参数的情况下，观测值 Yit对数据似然函数 p(Y∣f,β)的贡献可以很容易地计算出来。

然后，贝叶斯推理可以基于如下后验分布进行：

通过 Metropolis-Hastings (MH) 步骤，从非标准化的完整条件分布 p(fj[u;v]∣⋅)中生成用于 fj[u;v]的马尔可夫链样本，其中提议（proposal）使用基于条件先验的分布，如 Knorr-Held（1999）所提出的方法。

从 p(β∣⋅)以及非结构化空间效应 funstr（或其他非结构效应）中抽样，可以通过 Gamerman（1997）提出的加权最小二乘提议方法或其轻微改进版本实现。

方差参数的更新通过 Gibbs 步骤完成，直接从逆伽马分布中抽样。

对于单变量响应变量的情形，Fahrmeir 和 Lang（1999）描述了这些更新方案的详细内容。

由此产生的用于生成后验样本的混合 MCMC 方案是通过从以下完整条件分布中进行抽样来定义的：

抽样方案 1：

(i) 使用带有加权最小二乘提议的 MH 步骤对“固定效应”参数 β进行抽样。对于累积模型（2.5），阈值必须满足顺序限制（order restrictions）。

4 应用

我们考虑两个实际应用。在第一个应用中，关于失业持续时间，我们分析来自德国联邦劳工局（“Bundesanstalt für Arbeit”）的失业数据。这是一个庞大的数据集，包含约280,000条观测记录，展示了我们的方法即使在非常大的数据集上也具有良好的实用性。在第二个应用中，我们分析了1983年至1997年间在巴伐利亚北部的Rothenbuch森林区域收集的森林健康纵向数据。所有计算均使用 BayesX 完成，这是一个用于贝叶斯推理的软件包，由我们系开发。该程序可供公众使用，网址为 http://www.stat.uni-muenchen.de/~lang/。有关 BayesX 功能的详细描述，请参见 Lang 和 Brezger（2000）。

4.1 再就业机会

在我们的第一个应用中，我们分析了德国联邦劳工局1980年至1995年间的月度失业数据。我们的分析仅限于前西德地区（不包括柏林）以及女性个体的数据。对于每位个体，数据提供了其在第 t个月的就业状况、居住地区以及若干个人特征的信息。由于我们关注的是再就业机会的分析，并区分全职与兼职工作，我们将三分类响应变量 Yit定义为事件指示变量：

在 Fahrmeir 和 Lang（2000）中可以找到使用多项逻辑回归模型以及抽样方案 1 来从完整条件分布中抽取样本的类似预测变量的分析。

图2展示了获得全职工作（左列）和兼职工作（右列）的持续时间效应、日历时间趋势项和季节项的估计结果。持续时间效应呈现出其他研究中也观察到的典型模式：在2-3个月后达到峰值，随后逐渐下降。全职和兼职工作的日历时间趋势显示出相似的整体模式：1982年之前下降，之后缓慢上升至1990年（即德国统一一年后），随后再次明显下降，并伴有中期回升。这与德国在此期间劳动力市场的经济趋势相符。估计的季节效应在整个时期基本稳定，尽管其幅度有所变化。为了更深入了解，图f)和g)展示了1992年的估计效应片段，显示春季和秋季有典型高峰，12月则为整体最低点。

年龄效应如图1所示。女性约30岁时存在局部最小值，这可能是“家庭”效应所致。50岁以上人群的失业概率急剧下降尤其显著。60岁之后的上升可能由边界效应导致，因此不作解释。请注意，寻找全职工作的女性其年龄效应比寻找兼职工作的女性更为显著。

结构化区域效应如图3和图4所示。图3展示了估计的后验均值，图4则显示了“概率图”，其中颜色等级对应于“显著负向”（黑色）、“不显著”（灰色，即估计值的置信区间包含零）以及“显著正向”（白色）。要解释结构化效应，还需考虑非结构化效应。因此表4总结了不同地区的非结构化空间效应的估计后验均值。

我们观察到，获得全职工作的结构化效应比获得兼职工作的更强。更重要的是，兼职工作的非结构化效应明显超过其结构化效应，这与全职工作的估计效应形成对比。尽管图3b)中兼职工作的结构化效应后验均值显示出一些空间变异，但图4清楚地表明，在后验概率意义上并不存在“显著”的变异。相反，获得全职工作的结构化效应显示出“显著”的空间变异，南部地区的再就业机会优于中部和北部地区。

图3a)和图4b)中的两个深色区域标记了上世纪八九十年代以结构性经济问题著称的地区。

获得全职工作的固定效应估计见表2，获得兼职工作的见表3。表2和表3的结果确认了一些通过传统方法已知的事实：相比接受过职业培训或未接受过职业培训的女性，拥有大学学位的德国女性更容易再就业。这两个效应在兼职工作中更为显著。以往失业期的数量被用作劳动力市场经验的代理变量：以往失业次数越多，失业持续时间越短的概率越高。失业补助的影响显著为负，而失业救济金的影响则显著为正，这似乎与关于失业救济金具有负面副作用的普遍看法相矛盾。然而，这可能是因为“失业救济金”这一变量也充当了一个替代变量，代表那些过去曾工作过并定期缴纳保险费的人群。

4.2 森林健康

在这项关于树木健康状况的纵向研究中，我们分析了日历时间、树木年龄、冠层密度以及样地位置对山毛榉（beech）落叶程度（defoliation degree）的影响。数据来源于1983年至1997年间在巴伐利亚北部Rothenbuch森林区域进行的年度森林损害调查。该区域内共有80个观测点，分布范围约东-西方向15公里，南-北方向10公里，见图5。

在预测项中没有引入非结构化空间效应，原因如下两点：

首先，观察观测点的地图（图5）可以发现，有些样地只有一个邻居，这使得区分结构化效应和非结构化效应变得非常困难甚至不可能。其次，对于这80个样地中的每一个，只有15次对同一棵树的观测数据，且响应类别的变化很小。事实上，仅有少数几个样地观察到了全部三类响应结果。因此，在模型中引入非结构化效应会导致结构化与非结构化效应之间严重的不可识别问题，这一点可以通过检查参数的抽样路径来观察到。

我们最初采用了标准参数化的有序概率模型进行分析。然而，在抽样方案2的第(ii)步中，阈值 ξ1和 ξ2的后验样本混合性并不理想，见图6所示。根据 Chen 和 Dey（2000）的建议，我们对模型进行了重新参数化。首先，在式（4.1）中引入常数项 β0，从而可以将 ξ1=0。其次，由于潜在高斯模型预测项中的参数只能在乘法因子下被识别，我们假设误差项 εit服从 N(0,σ2)分布，其中方差 σ2未知。这使得我们可以设定 ξ2=1。

参数 β0与固定效应参数 β1和 β2同时进行抽样。对于 σ2，我们指定其先验为逆伽马分布，从而其完整条件分布也为逆伽马分布，便于生成后验样本。

关于估计结果的解释需注意以下几点：根据我们对模型（2.3）至（2.5）的定义，预测项（4.1）中（或其中各项效应）的较高（较低）值对应于树木更健康（较差）的状态。

表1给出了 β0以及年龄效应的估计结果。正如我们所预期的那样，较年轻的树木比年老的树木状态更健康。图7展示了日历时间趋势和冠层密度效应的后验均值估计结果。可以看出，树木在1986年前后的糟糕年份之后有所恢复，但自1994年后健康状况再次下降到较低水平。冠层密度超过30%后效应呈现出明显的单调递增趋势，表明在冠层密度较高的林分中，山毛榉能够更好地抵御不良环境因素的影响。

图8以后验概率的形式展示了估计的（结构化）空间效应，其中黑色斑块表示具有严格负向可信区间的区域，即该区域中较多树木处于较差状态。这些黑色标注的样地大多位于森林区域中海拔较高的地方，相比其他区域，这些地区的营养物质含量和土壤质量较差。

5 结论

这些应用表明，所发展的贝叶斯方法是在现实复杂分类回归模型中进行推理的有用且灵活的工具。通过修改或推广观测模型、预测项和光滑先验，可以实现多种扩展。例如，基于潜在效用的probit模型可以通过考虑潜在的多变量半参数高斯模型，推广到相关分类响应或混合连续-分类响应的情形。通过引入协变量之间的非参数交互效应（如Clayton（1996）和Knorr-Held（2000）所建议的），可以使预测项更加灵活。

将高斯先验替换为重尾分布，则可以用于处理不平滑的回归函数。

原文链接：https://epub.ub.uni-muenchen.de/1592/1/paper_202.pdf