基于分组数据的时空平滑、插值和分布预测-7-26

2024_04_27_9364f8eb1b0ca10a3b62g

然而，由于最高收入阶层是开放式的，粗略估计的构建相当随意，其不确定性的量化可能是不可能的。

大多数市镇的第三收入阶层（200 万至 300 万日元）和第四收入阶层（300 万至 400 万日元）居民所占比例最大。远离市区的市町村，如日本西南部和北部的市町村，第一收入阶层（100 万日元以下）和第二收入阶层（100 万至 200 万日元）的居民所占比例较大。其他几轮 HLS 的图表见补充材料。

家庭收入调查每五年进行一次，至少一年后才会公布结果。为了在时间和空间上灵活、快速地制定细粒度的政策，最好能定期（如每年）提供和更新所有城市的收入状况信息。

总之，尽管 HLS 收集了日本全国家庭收入的一些详细数据，但分析人员仍将面临以下挑战：

这些数据将以分组数据的形式出现，无法轻易从中获得收入和贫困衡量标准；
数据将包括许多未被抽样调查的城市；
尽管人们希望了解最新的收入状况，例如每年一次，但调查将每五年进行一次。

下一节提出的时空混合模型可以克服这些挑战。

700 万至 1000 万日元

1,000 万至 1,500 万日元

图 1：2018 年 HLS 数据中九个收入等级的比例。

3 时空混合物建模

3.1 模式

假设我们感兴趣的是

地区在

期间的收入分布情况，用

表示，密度函数为

，

和

。我们还假设无法观测到单个家庭的收入，只能获得分组数据。具体而言，设

为

第三期的

收入等级，其中

和

通常分别设为 0 和

。因此，每个区间的住户数都是观测值，用

表示。与 HLS 数据一样（见第 5 节），我们假定所有的

的区域不采样。在不失一般性的前提下，对第一个

地区进行抽样，其余

地区不进行抽样，

。此外，分组数据中第

个地区和第

个时期的住户数用

表示。

为了灵活地模拟基本的时空收入分布，我们考虑使用对数正态分布的混合物，如下所示：

其中，

是与收入数据相关的某些外部信息（如消费水平或税收水平）的

维向量，

和

分别是

th 分量的系数和标度参数，

表示带有

和

参数的对数正态分布的密度函数。

(1) 中成分分布的协变量

至关重要。关于非采样城市的空间插值，

，必须对所有城市进行观测。此外，在时间预测方面，还假定

与

的最新观测数据比 HLS 更频繁获得，例如每年一次。我们的 HLS 数据分析基于每个纳税人的应税收入（见第 5.1 节）。

众所周知，仅采用对数正态分布对收入数据的拟合效果很差。然而，对数正态分布的混合分布却能很好地拟合数据（Lubrano 和 Ndoye，2016 年）。此外，(1) 只是针对未观察到的连续家庭收入的假设。对于分组收入数据，

th 地区在

th 期间的似然贡献采用多叉形式：

为

和

，其中

是

第

期面积的分布函数，其分布函数为对数正态分布，

。

混合比例

模拟为

其中，

表示总体混合比例，

表示特定区域的空间效应，

表示时间效应。我们将所有

地区的行标准化邻接矩阵表示为

。利用

，我们独立假设

的同步自回归（SAR）模型为

。即

取样参数

和相关参数

。取样区域空间效应的边际分布

为

。我们假设时间效应，

，遵循随机漫步过程，

，为

，其中

为初始位置参数。

3.2 以前的分配

在这里，我们指定了参数的先验分布。由于这些先验分布在数据增强后是有条件共轭的，因此便于吉布斯抽样算法的开发。关于分量分布的参数，

和

，我们假设

和

为

。由于这些分量分布在空间和时间上是共通的，可以稳定估计，因此默认的超参数选择为

和

，这样它们就是弥散的。关于

，假定正态先验为

，默认值为

。通过使用行标准化邻接矩阵，空间相关参数的先验值由

给出。对于合成孔径雷达模型的精度参数

，我们假定

为默认值

，这样就大大涵盖了可能的值。关于随机漫步过程的初始位置和方差，我们假设

和

为

。默认超参数选择为

, 和

。

3.3 后验推理

后验推断基于马尔科夫链蒙特卡罗（MCMC）方法。我们开发的吉布斯采样算法无需调整，因为它不需要任何 Metropolis-Hastings (MH) 步骤。具体细节见补充材料，我们在此简要介绍一下构建吉布斯采样器的策略。首先，为了便于对

中的变量进行采样，我们采用了 Polson 等人（2013 年）的 Pólya-gamma 数据扩增法，以便从正态分布中对这些变量进行采样。其次，我们引入了一个潜变量，它表示

第三收入阶层中属于

第三混合物成分的家庭数量。第三，引入潜在的个体家庭收入，以方便对

和

进行采样。根据

第 1 个收入类别和

第 1 个混合物成分中的数据和家庭数量，从截断正态分布中生成单个家庭收入的对数。最后，从

和

的全条件分布中取样，根据条件共轭，它们分别是正态分布和逆伽马分布。

根据 MCMC 的输出结果，可以得到相关数量的后验推断。在这种情况下，我们对

和

的平均收入、

、收入中位数和基尼指数感兴趣。在混合模型下，收入中位数和基尼指数无法通过分析获得。关于收入中位数，

，对

进行数值求解。对于基尼指数（由

定义），按照 Lubrano 和 Ndoye（2016 年）的方法对积分进行数值计算。

考虑到参数和潜变量的 MCMC 抽样，

从 SAR 模型的条件分布中抽样：

。关于时间预测，

取自

。

在实践中，成分的数量

是未知的。然而，它却会影响估算性能或对估算结果的解释。根据 Celeux（1999 年）、Fürwirth-Schnatter 和 Pyne（2010 年）以及 Malsiner-Walli 等人（2016 年）的研究，

-均值聚类应用于

的 MCMC 输出。假设我们有

MCMC 抽样

，表示为

，为

和

。然后，

, s 作为

-means聚类的输入，为

和

的每个

分配标签

。计算

-means聚类的聚类标签与

的排列相匹配的 MCMC 抽样分数。在实践中，匹配分数的计算公式为

，其中

表示指标函数。例如，当

聚类分隔良好，且

等同于

的集合时，

。如果给定

的 MCMC 抽样的匹配分数小于 1，则表明混合物过度拟合；因此，首选较小的

。

4 模拟研究

首先使用模拟数据演示了所提出的方法。二维空间上的

区域单位是通过从

中抽取每个坐标生成的。一个区域的邻域定义为坐标距离在 0.2 欧几里得范围内的区域，平均邻域数为 5.77。在每个区域，

由

生成，共 21 期。各组分别由

、

和

与

代表所有

。这些分组与后三轮 HLS 中的分组相对应（表 1）。在本次模拟研究中，分组数据由混合模型（1）生成，

。分量分布的参数设置如下：

。每个协变量（

）均由均匀分布生成。

使用相同的成分分布，我们考虑了以下两种混合比例设置。在第一种情况下，混合比例由

生成，

。在第二种设置中，空间和时间效应由以下确定性序列决定：

。此外，二维空间

被划分为

个方形块。同一区块内的区域对混合比例的影响相同：

。在

在这种情况下，区域

是第

个区块的子区域。因此，

区域对应的

区块的混合比例与

对应的

成比例。每个区块的子区域平均为 8 个。

首先，我们对第 3 节所述的时空混合物模型进行拟合，拟合时使用不同的成分数：

。建议的吉布斯采样器运行 20,000 次迭代，初始消耗期为 10,000 次迭代。图 2 显示了两种设置下建议的混合物模型的匹配分数。在两种设置下，只要

，匹配分数都是 1。然而，图中显示，在

, 和 6 , 的过度指定模型中，匹配分数变得小于 1，表明混合模型的过度拟合。从图中可以看出，从

的较小值开始监测匹配分数，并在扭结处选择适当的

是合理的方法。

图 2：MCMC 抽样中

-均值聚类的输出与

的排列相匹配的百分比。

接下来，我们将拟议的时空混合模型与四个替代模型进行比较。第一个替代模型是基于（1）和（3）的混合模型。混合比例用代表空间和时间异质性的双向独立随机效应建模：

分别表示空间和时间的异质性。由于该模型具有空间和时间效应，因此可用于空间插值和时间预测。第二个替代模型是对数正态分布的混合模型，其混合比例仅包括

和

。这是一个纯粹的空间模型，对每个

进行独立拟合，并用于空间插值。

只是。该模型可视为 Kawakubo 和 Kobayashi（2023 年）提出的分组数据小区域模型的空间混合扩展。第三种选择是分组数据的简单对数正态模型。该模型使用最大似然法对每个地区和时期进行独立估计。第四种选择是 SAE 模型。更具体地说，我们考虑了 Torabi 和 Rao（2014 年）提出的线性混合模型，这也是其中一位审稿人的建议。该模型如下

其中，

是该地区的分地区数，

，这样，

是下文所述的粗平均收入，

是地区效应，

是分地区效应，

是抽样误差。平均收入是根据

计算的，其中

是收入等级的中点。最高收入阶层

的区间是开放的，

设为 20。该模型在 "设置 2 "中拟合，而没有

的模型在 "设置 1 "中拟合，因为真正的模型没有分地区结构。小区域模型在每个时期独立估算。

地区和

期间的数据用于模型估计。其余

地区和最后一个时期的数据则用于空间插值和时间预测。

首先根据样本内估计、空间插值和时间预测中平均收入的均方根误差（RMSE）和

可信区间或预测区间的覆盖率对模型进行比较。例如，在样本内估计的情况下，均方根误差由

定义，覆盖率由

定义，其中

是真相，

是抽样地区的后验均值或最大似然估计值，

和

是贝叶斯模型下

可信区间的端点。关于空间插值和时间预测情况，

对应于

和

的非采样区域的后验预测均值；它也对应于

的所有区域的后验预测均值。同样，

和

与

预测区间的端点相对应。

表 2 列出了五种模型下平均收入的均方根误差和覆盖率。总体而言，拟议模型的均方根误差优于其他模型。覆盖率在目标值 0.95 左右。至于双向模型，虽然样本内估计的均方根误差在两种情况下都是第二小的，但在提出的模型之后，它在空间插值和时间预测方面产生了较高的均方根误差。此外，覆盖率也低于 0.95 的目标值。空间模型在空间插值方面的均方根误差在两种情况下都是第二小的，紧随拟议模型，覆盖率约为 0.95。这些结果表明了时空效应的有效性，而包含这两种效应的拟议模型则表现出更优越的性能。由于缺乏跨时空的借用信息，通过最大似然法分别估计每个地区和时期的分组数据模型会导致样本内估计的均方根误差很高。采用的 SAE 方法不是专门为分组数据或时空环境设计的。因此，它在样本内估计和空间插值中产生的均方根误差最大。

最后，根据 RMSE 和住户数的覆盖率，比较了建议模型、双向模型和空间模型。表 3 列出了两种设置的结果。具有空间和时间效应的拟议模型表现良好，除了设置 2 中的样本内估计和目标 0.95 附近的覆盖率外，均方根误差最小。在所有情况下，双向模型产生的均方根误差最大。在空间插值方面，空间模型产生的均方根误差小于双向模型。不过，借用时间信息的拟议模型产生的均方根误差小于空间模型。

表 2：使用提议的时空方法（Proposed）、双向独立时空方法（Two-way）、独立空间方法（Spatial）、SAE 和独立分组最大似然方法（GML）进行样本内估计、非样本区域空间插值和时间预测的平均收入均方根误差和覆盖率。

			Proposed	Two-way	Spatial	SAE	GML
Setting 1	RMSE	In-sample	0.035	0.085	0.185	0.345	0.217
		Spatial	1.094	1.400	1.188	1.420	-
		Temporal	0.493	0.597	-	-	-
	Coverage	In-sample	0.952	0.553	0.992	0.929	-
		Spatial	0.976	0.810	0.953	0.888	-
		Temporal	0.990	0.910	-	-	-
Setting 2	RMSE	In-sample	0.048	0.049	0.131	0.417	0.254
		Spatial	0.375	1.032	0.471	0.762	-
		Temporal	0.254	0.875	-	-	-
	Coverage	In-sample	0.920	0.908	0.947	0.958	-
		Spatial	0.985	0.862	0.978	0.874	-
		Temporal	0.975	0.875	-	-	-

表 3：在拟议时空模型、双向时空模型和单独空间模型下，样本内估算、非样本区域空间插值和时间预测的住户数的均方根误差和覆盖率

			Proposed	Two-way	Spatial
Setting 1	RMSE	In-sample	5.027	5.264	5.229
		Spatial	20.446	25.184	21.552
		Temporal	6.721	8.059	-
	Coverage	In-sample	0.968	0.959	0.970
		Spatial	0.963	0.906	0.951
		Temporal	0.973	0.964	-
Setting 2	RMSE	In-sample	5.065	5.074	4.681
		Spatial	7.798	18.382	9.130
		Temporal	5.104	11.721	-
	Coverage	In-sample	0.968	0.968	0.981
		Spatial	0.977	0.938	0.980
		Temporal	0.979	0.962	—

5 HLS 数据分析

5.1 设置

在分析中，我们使用了家庭收入调查的一般家庭年度总收入数据，包括工资收入和转移性收入。我们使用 800 个城市的数据进行模型估计，其中的观测值在所有轮次中都可用。其余的观测数据被保留下来，以比较模型的性能（第 5.2 节）。由于数据跨度长达 20 年，因此收入等级的终点（表 1）以 2018 年的消费物价指数为标准进行了调整。此后，可以对收入分布进行时间比较。

关于协变量，我们采用纳税人人均应纳税收入（单位：百万日元），用税率

表示，并设置

。该协变量取自税法。

日本总务省的调查数据

。日本所有市镇每年都有这些数据，可以根据 HLS 合理预测收入。与收入等级一样，应纳税收入也根据消费价格指数进行了调整。每个纳税人应纳税收入的空间分布见补充材料。

值得注意的是，虽然这些辅助信息每年都可在市一级获得，但却不能单独用于推断收入分布。不过，为了准确地进行空间插值和时间预测，我们建议的模型必须将这一重要信息作为协变量纳入其中。

5.2 模型比较

建议的吉布斯采样器运行 30,000 次迭代，初始消耗期为 10,000 次迭代。图 3 显示了

-means聚类的输出与

的排列相匹配的 MCMC 抽样的分数，混合成分的数量在 2 到 8 之间。在拟议模型下，

的匹配分数为 1。至于

，匹配分数低于 1，表明过度拟合。图 3 显示了双向独立效应模型的结果；我们观察到的模式与建议模型类似。

使用空间插值法预测了被排除在外的城市收入组中的家庭数量。建议的模型与 Sugasawa 等人（2020）提出的分组数据双向和参数空间收入模型进行了比较。在 Sugasawa 等人（2020 年）的模型中，参数收入分布的参数通过空间随机效应的适当变换随空间而变化。我们将 Singh-Maddala 和 Dagum 分布视为参数模型，因为它们常用于收入建模（Kleiber 和 Kotz，2003 年）。我们根据

预测区间的覆盖率，以及

所定义的户数的相对绝对误差（RAE）的 0.5、0.9 和 0.99 量级，对四个模型的性能进行了比较，其中

表示保留市的后验预测均值。表 4 显示，拟议模型和双向模型

图 3：MCMC 抽样的匹配分数，其

- 均值聚类的输出与 HLS 数据

的排列相匹配。

在这两个模型中，Singh-Maddala 模型，特别是 Dagum 模型，产生了约 0.95 的理想覆盖率。Singh-Maddala 模型，特别是 Dagum 模型，导致覆盖率不足。虽然拟议模型和达古姆模型产生的 RAE 中位数相当，但达古姆模型产生的 RAE 的 0.9 和 0.99 百分位数大于拟议模型产生的 RAE 百分位数。Singh-Maddala 模型产生的 RAE 量值最大。总体而言，下表显示了拟议模型的最佳性能。

根据图 3 和表 4 中的观察结果，在接下来的分析中，我们采用

作为拟议的时空混合模型。

的参数后验分布见补充材料。

表 4：拟议时空模型（拟议）、双向时空模型（双向）以及 Dagum 和 Sing-Maddala 模型下的保留市镇的覆盖率以及 RAE 的 0.5、0.9 和 0.99 百分位数

		RAE 的定量
Model	Coverage	0.5 th	0.9 th	0.99 th
Proposed	0.955	0.136	0.411	1.105
Two-way	0.950	0.161	0.469	1.259
Dagum	0.553	0.130	0.445	1.383
Singh-Maddala	0.869	0.196	0.875	3.654

5.3 日本收入地图

我们提供了日本在人口普查年份以及未进行人口普查的 2020 年的完整收入地图。在调查年份，对未抽样的市镇进行了空间插值。在 2020 年，所有市町村都采用了时间插值法。

预测。2020 年的时间效应由

生成，因为

的增量相当于五年的增量。

图 4 显示了日本平均收入的完整地图。颜色较深的城市表示收入较高的地区。随着时间的推移，地图的整体色调有变浅的趋势。例如，本岛中部（本州）的市町村比西部岛屿（四国和九州）的市町村颜色要深得多。然而，在 2018 年，高收入市镇大多只集中在大都市地区。2020 年，预测的平均收入呈现出与 2018 年相似的模式。我们还观察到收入中位数的类似模式；结果见补充材料。

图 5 显示了基尼系数的后验均值和后验预测均值。深色表示基尼系数相对较高的城市。20 年来，全国各地基尼系数的差异可能有所缩小。例如，1998 年，日本北部特别是西南部的市镇用深色表示，而本岛中部的市镇则用浅色表示。从图中可以看出，2003 年日本全国的基尼系数都很高。然而，在 2018 年，虽然日本西部和北部地区的市镇仍然表现出相对较高的基尼指数，但地图的整体色调似乎更加均匀。

图 6 展示了 1998 至 2020 年间日本全国平均收入、中位数收入和基尼系数分布的变化。总体平均收入和中位收入呈现出明显的下降趋势。图中还显示，基尼系数的变化随着时间的推移而减小。

图 4：已完成的平均收入地图

图 5：已完成的基尼系数地图

图 6：平均收入、收入中位数和基尼系数的分布随时间的变化。垂直线表示最近一次 HLS 于 2018 年进行

5.4 按地区划分的收入分布

在此，我们重点关注任意选取的八个城市的地区收入分布情况。所选城市的位置如图 7 所示。抽样城市和非抽样城市分别用蓝色和红色表示。

图 8 显示了调查年份和 2020 年未调查年份的收入分布后验均值。上下两行分别对应抽样城市和非抽样城市。通过所提出的时空混合模型，可以得到不同的收入分布形状。例如，作为东京最富裕地区之一的港区，与其他市町村相比，表现出较长的右尾部。与此相反，本州最北部的二户、萨摩和中土佐，以及本州最西部和四国南部分别表现出较短的右尾。这些分布的主体主要集中在年收入 1000 万日元以下的地区。

有趣的是，忍野町的收入分布显示出较长的右尾部，这与港町相似。忍野是一个小镇，没有在 HLS 中抽样。然而，它也是一家大型电气机械制造商的总部和工厂所在地，因此收入水平较高。结果表明，所提出的模型可以通过协变量和时空混合比例等辅助信息推断出该地区的收入分布情况。在图中列出的所有市镇中，除忍野市外，收入分布的右尾

在过去的 20 年中，低收入地区的人口密度略有下降，分布密度有所增加。

图 9 显示了选定地区的平均收入、收入中位数和基尼系数随时间变化的后验均值和后验预测均值。从图中可以看出，除忍野市外，所有市镇的平均收入和收入中位数都随着时间的推移而下降。随着时间的推移，港区的基尼系数约为 0.41，而忍野区的基尼系数则从 0.38 逐步上升到约 0.4。其余地区的基尼系数变化在过去 20 年中有所减小，与全国的结果类似。

图 7：任意选择的城市位置（蓝色：取样城市，红色：未取样城市）

图 8：所选抽样城市（上排）和非抽样城市（下排）收入分布的后验均值和后验预测均值随时间的变化情况

图 9：选定地区的平均收入、收入中位数和基尼系数随时间的变化。垂直虚线表示最近一次人口动态统计调查是在 2018 年进行的

5.5 事先敏感性检查

最后，我们评估了规模参数先验规范的敏感性。图 10 显示了在默认的反伽马先验和替代先验下平均收入、收入中位数和基尼指数的后验均值散点图。在其他先验设定中，标准差

和

遵循均值为 0.2 的指数分布，这有利于较小的标准差，先验为

标准差超过 0.5 的概率约为 0.08。从图中可以看出，这些相关的后验量对于先验规范是稳健的。补充材料提供了这些先验设置下标准偏差参数后验分布的比较。

图 10：默认逆伽马先验（IG）和替代指数先验（Exp）下的平均收入、收入中位数和基尼指数的后验均值散点图。

6 讨论

在本研究中，我们提出了一种基于分组数据分析收入分布的新型时空混合模型。通过协变量信息和时空异质性效应，所提出的模型可以推断和预测任何时间点上采样地区和非采样地区的收入分布。如果有比分组收入数据更高频率的适当协变量，就可以通过预测对最新的收入状况进行深入了解。通过将建议的方法应用于 HLS 数据，我们得到了平均收入、收入中位数和基尼指数的完整地图。我们发现，在过去 20 年中，总体平均收入和中位收入均有所下降，同期基尼系数的变化也有所减小。我们还证明，采用适当的协变量对揭示各地区的收入状况至关重要。对于非抽样地区来说尤其如此，因为协变量信息是唯一可用来推断收入状况的信息。

建议的方法可以揭示收入分布和相关的贫困指标。

只要协变量信息可用，就可以对任何地区和任何时期的收入水平进行估算。虽然本文没有展示，但也可以估算出每个地区收入低于任何给定收入水平临界值的人口比例。通过探索所提议方法的结果，政策制定者可以确定需要进行政治干预的领域，以通过有效的政策实施和资源分配来缓解全国范围内的贫困和不平等现象。这种政治干预措施可以根据收入分布的最新信息来实施，而拟议模型的空间插值和时间预测则阐明了这一点。

本研究采用全 MCMC 方法进行后验推断。建议的吉布斯采样器的优点是无需对算法进行调整。由于包含许多潜变量可能会导致效率低下，我们可能会认为需要另一种计算策略。事实上，我们也曾尝试实现一种

算法及其自适应版本对

和

进行采样，而不产生单个家庭的收入。换句话说，(2) 直接用于贡献似然函数。然而，MH 算法未能探索参数空间，大多停留在初始值附近。其他一些更快的后验计算方法，如集成嵌套拉普拉斯近似（Rue 等人，2009 年）和变异贝叶斯近似（Blei 等人，2017 年），都是很有吸引力的替代方法。关于 HLS 数据的应用，基于苹果 M1 Ultra 芯片上未经优化的 R 代码，在没有任何并行化的情况下，计算时间（例如

）约为 40 小时。由于我们更倾向于确保后验推断的质量，而且计算时间并不是我们应用中的主要问题，因此我们认为目前的计算方法是可以接受的。不过，开发计算效率更高的算法也是一个必须探索的相关研究方向。

鸣谢

本研究得到了日本学术振兴会（#22K01421、#21K01421、#21H00699、#20H00080 和 #18H03628）和日本经济研究中心的支持。

税务调查数据见 http://www5.cao.go.jp/keizai-shimon/kaigi/special/ future/keizai-jinkou_data/csv/file11.csv