中国县域人口少子化的预测与解释：机器学习的视角
Mengtaoa,Lanlanb
a 云南民族大学经济与管理学院 教授 b 云南民族大学经济与管理学院 副教授 通讯作者dunzi21@qq.com
撰写于2024-12-28
Prediction and Explanation of Low Fertility in China’s County Populations: A Machine Learning Perspective
摘要
2024 年中央一号文件强调农村生育支持，凸显农村少子化问题。中国农村生育模式变化显著，少子化引发劳动力、人口结构、家庭养老及教育资源等多方面问题，全球多地亦面临此困境，中国生育率已低于稳定水平。学界虽对少子化成因有研究，但中国农村少子化研究存在解释缺深度、预测缺基础等问题。本文利用 2000 - 2020 年人口普查分县数据库及相关数据，采用机器学习方法研究农村少子化。经特征工程与筛选确定 11 个特征建模，选 CatBoost 算法并调优。研究发现部分人口学与地理区位特征重要，特征间非线性与交互效应显著、个体异质性强。 边际贡献体现在：解释上，揭示人口转变路径依赖，中国农村少子化进程快，经济因素未入选可能因特征相关性及家庭生育决策自适应；预测上，考虑非线性与相关性，SHAP 方法的局部解释为分县预测提供基础，优于传统建模。但研究存在变量受限、数据未涵盖未来变化、模型研究不全面及未涉城市人口变迁等局限。
Abstract In 2024, the No. 1 Central Document emphasized support for rural fertility, highlighting the issue of rural low fertility. The fertility pattern in rural China has changed significantly. Low fertility has led to problems in multiple aspects such as labor supply, population structure, family elderly care, and educational resources. Many regions around the world are also facing this dilemma, and China’s fertility rate has fallen below the level required for population stability. Although the academic community has conducted research on the causes of low fertility, studies on rural low fertility in China have problems such as lack of in-depth explanation and prediction foundation. This paper uses the county-level database of the 2000 - 2020 population census and related data, and adopts the machine learning method to study rural low fertility. After feature engineering and screening, 11 features are determined for modeling, and the CatBoost algorithm is selected and optimized. The study finds that some demographic and geographical location features are important, and there are significant non-linear and interaction effects among features, as well as strong individual heterogeneity. The marginal contributions are as follows: in terms of explanation, it reveals the path dependence of population transformation. The process of rural low fertility in China is rapid, and the non-selection of economic factors may be due to feature correlation and the adaptive adjustment of family fertility decisions; in terms of prediction, considering non-linearity and correlation, the local explanation of the SHAP method provides a basis for county-level prediction, which is superior to traditional modeling. However, this study has limitations such as limited variables, data not covering future changes, incomplete model research, and not involving urban population changes.
关键词
农村少子化；机器学习；人口预测；影响因素；模型解释
Keywords rural low fertility; machine learning; population prediction; influencing factors; model interpretation
1 引言
2024 年，中央一号文件首次提出加强农村生育支持政策。此政策的出台，标志着我国对农村人口问题的重视程度提升至新的高度，农村少子化问题已然凸显。
中国农村生育率虽始终高于全国水平，但在 21 世纪，其生育模式发生急剧变化。具体表现为生育率下降与生育推迟并存，生育高峰后移且更为分散；乡村与镇的生育模式逐渐趋同；农村年轻女性的多孩生育比例上升，区域差异显著。
(吴, 2024; 周 et al., 2024; 徐 and 张, 2021)。
农村地区的少子化现象在当前形势下，极有可能带来以下几个方面的问题。
劳动力供给呈现不足态势：少子化现象预示着未来农村劳动力人口数量的减少。伴随时间的逐步推移，此情况将对农村经济发展产生重大影响。尤其在农业生产等领域，极有可能面临劳动力短缺之问题，进而对粮食安全等重要方面造成影响。甚至可能致使部分土地无人耕种或者耕种不及时，从而影响农业的可持续发展。
人口结构失衡态势加剧：少子化与老龄化相互作用，致使农村人口老龄化程度愈发加深，使得老年人口占比进一步提升，从而增加了社会养老负担。例如，2020 年，我国农村老年人口数量已达 1.21 亿，相当于平均每四个农村人口中约有一位 60 岁及以上老人。在此情形下，农村养老保障面临着更大的压力，也对农村社会养老服务体系的建设与完善提出了更高的要求。
家庭养老功能呈削弱态势：子女数量有所减少，家庭规模渐趋变小，传统的家庭养老模式遭受冲击。农村老人的养老保障面临严峻挑战，亟需更多的社会养老服务予以补充。例如，部分农村老人的子女在城市工作生活，难以在其身边进行照顾，老人在生活照料、精神慰藉等方面存在问题，进而影响老人的生活质量与幸福感(刘静丽. et al., 2024; 马, 2024)。
教育资源利用不足。随着农村儿童数量的减少，部分农村学校出现了生源不足的情况。这致使教育资源被闲置浪费，教育规模效益下滑，进而影响了农村教育质量的提升以及教育事业的可持续发展(何 and 黄, 2024; 刘升., 2024; 吕, 2024)。
20世纪以来，少子化问题陆续在欧洲、美洲、东亚等地出现。日本自1992年起就面临严重的少子化问题。2020年，其生育率为1.34，65岁及以上老年人口占总人口的比例已超过28%，人口持续减少，一些地方甚至出现城镇荒废的情况。欧洲的匈牙利、保加利亚等国也曾有过较长时间的人口负增长历史。
一份最新的权威研究成果表明：在 1950 年至 2021 年这长达七十年的时间区间内，全球总生育率（TFR）降幅超过一半。1950 年之际，TFR 为 4.84（即每一名女性平均生育约五个孩子），而到了 2021 年，TFR 降至 2.23。据预测，未来生育率将在全球范围内持续下降，至 2050 年全球总生育率预计将达 1.83，到 2100 年则将达到 1.59。(Bhattacharjee et al., 2024)。
我国目前已步入老龄化与少子化并存的社会状态。少儿人口比重呈持续下降态势，2022 年出生人口仅为 956 万，此数据创下自 1949 年以来的最低纪录。生育率已然跌破 1.1，不仅远低于维持人口总体规模稳定所需的 2.1 生育率水平，亦低于国际社会通行的 1.5 警戒线(我国2022年出生人口956万，三孩及以上占比为15%, n.d.)。
学界在少子化成因方面已然取得众多成果，且呈现出典型的跨学科研究特质。Roser (2024) 对迄今为止的文献进行了系统总结。其指出少子化的原因主要包含三个方面。其一为经济因素。当前，生活成本持续上升，特别是在住房、教育、医疗等领域，支出显著增加。这使得众多年轻人面临较大的经济压力，难以承担养育孩子所需的费用。其二是社会文化因素。在现代社会中，个人主义观念日益盛行，人们更加注重自我价值的实现以及个人生活品质的提升，生育意愿相对降低。同时，女性社会地位的提高与教育程度的提升，促使其更加追求职业发展和个人独立，晚婚晚育乃至不婚不育的现象愈发普遍。此外，传统家庭观念逐渐淡化，如大家庭的解体、亲属关系的疏远等，致使年轻人在养育孩子方面缺乏充足的家庭支持。其三为政策因素。部分国家的政策尚不完善，未能有效缓解年轻人的生育压力。例如，育儿津贴、教育补贴等福利政策不足，无法满足家庭养育孩子的实际需求。另外，税收政策、劳动保障政策等方面也可能存在不利于生育的因素，如税收负担过重、产假和陪产假制度不够完善等。
中国学者指出，农村生育率下降的原因在很大程度上与此相同。除经济因素之外，中国还面临城市化进程加速以及劳动力迁移等问题。与此同时，农村地区在观念转变方面也面临着一定的挑战(刘升., 2024; 刘静丽. et al., 2024; 李婷 et al., 2019; 王增文 et al., 2023; 解 and 林, 2023; 闫 and 霍, 2024)，育儿成本被视作最为主要的原因之一。(梁建章，黄文政，何亚福, 2024)，张孝栋 et al. (2021) 的综述为国内研究进展提供了更为详尽的全景展示。最新的一项权威研究表明，生育意愿下降以及延迟生育乃是致使生育率下降的主要因素，在医学层面体现为育龄人口生育能力的降低(Wang et al., 2024)。
文献检索结果表明：在中国农村地区，少子化问题的研究仍存在一定的欠缺。其一，解释力度欠佳。当前，虽对农村地区少子化现状有所描述，且该问题已成为普遍共识，亦进行了部分原因推测，然而缺乏更多跨学科实证研究的证据支撑。已有文献显示，影响农村少子化的因素众多，除人口学、经济学、政策因素之外，地理区位、环境污染、地缘政治风险、贸易条件等亦具有一定影响。故而，需要具备跨学科视角，纳入尽可能多的解释变量。其二，预测基础不足。原因分析固然重要，但对农村少子化未来状况的预测则更为关键，唯有如此，方能为政策制定提供有效依据。并且，这种预测不应局限于人口学领域，也不能仅仅是平均趋势预测，而应细化至具体的县级区域（或地级区域），以便进行具有针对性的分类指导。
鉴于此，本文以 2020 年第七次人口普查县域数据为基础，结合县域其他经济社会指标、地理交通状况、土地利用情况、灯光污染等可获取之指标，运用机器学习方法对中国农村地区少子化的影响因素及其重要性进行分析，并对今后一段时期农村地区的县域少子化情况予以解释与预测。
本文之所以选用机器学习而非传统的线性统计建模方法，缘由如下：其一，人口少子化的影响因素繁多，绝非仅几个因素便可进行阐释。在建模层面，几乎难以获取精确的理论指引；其二，非线性特征显著。姑且不论单个影响因素的幂次设定，仅影响因素之间的相互关系便极为复杂。例如，女性的教育与经济条件变化之间存有明显的交互作用，类似的相互作用不胜枚举，传统建模对此无能为力；其三，影响因素之间无疑存在相关性。即便共线性尚可容忍，但对参数的解释却捉襟见肘。因为传统模型的参数含义是在假定“其他条件不变”（Ceteris Paribus）的前提下做出的解释。鉴于影响因素存在相关性，那么“其他条件不变”在现实中乃是根本不存在的虚假前提。(James et al., 2013)。
本文以下论述将分为四个部分予以展开：第2节为研究使用的数据与方法介绍；第3节为机器学习回归过程的建模，涵盖特征工程以及特征选择、算法挑选与比较、超参数调优三个环节进行阐释；第4节为机器学习的模型解释；第5节为全文的总结与评述。
2 文献/数据与方法
2.1 文献回顾
关于机器学习解释与人口出生率的文献数量尚为有限，然而已有相关的有益探索。M. Li and Xu (2022) 依托中国综合社会调查（CGSS）的微观数据，采用机器学习及逻辑回归方法，对影响中国居民二胎生育意愿的因素进行了系统分析。主要研究发现如下：年龄乃是影响二胎生育意愿的最为关键因素；家庭收入对生育意愿有所影响，但在所有变量中仅位居第六；兄弟姐妹数量越多，生育意愿则越强；地区对生育意愿存在影响；健康状况良好的人群生育意愿更高。同时，性别差异普遍存在。Zhuang et al. (2021) 采用深度学习方法，结合 Landsat 图像对中国 1985 年至 2010 年多时间尺度人口分布展开研究。通过获取中国 2010 年城镇级人口普查数据以及 WorldPop 1km 分辨率的网格化人口数据，对城镇人口予以重新分配。经研究发现，高人口密度区域的面积与人口比例显著增加；中等人口密度区域的面积和人口比例均有一定程度减少；低密度区域面积略有扩大，然而人口比例降低。此现象或许与中国的扶贫及山区移民政策存在关联。Farahani (2024) 着重探讨了人工智能于社会科学问题中的应用，尤其在预测人口变化方面的案例研究。选取了 OECD 国家的 18 个经济指标作为输入变量，以反映经济发展的各个层面对于人口变化的影响。在对神经网络与遗传算法进行比较后发现：与传统的统计和计量模型相比，AI 技术能够提升人口预测的准确性，更优地适应复杂的数据结构及问题。Tzitiridou-Chatzopoulou et al. (2024) 主要就生育率与地缘政治风险之间的关系展开了深入探讨，并运用多种机器学习方法进行分析。经研究发现，非线性机器学习回归模型，诸如随机森林、额外树以及 XGBoost 等，在对苏格兰未来出生人数进行预测时展现出较为良好的能力，然而线性回归的预测效果则相对较差。总体而言，bagging 树模型在多步预测中表现出较强的竞争力，LGBM 算法亦取得了较为理想的结果，同时该研究着重强调了超参数调优的重要性。
其他若干研究，其中包括运用小样本针对伊朗女性生育倾向所进行的研究。 (Moulaei et al., 2024)；以国家作为研究区域，运用联合国数据库收集 1950 年至 2015 年期间 205 个国家的七个人口统计特征数据，就其在预测国家人口增长率方面的性能展开比较研究。(Otoom, 2021)；基于人口统计学的机器学习方法，对 G7 国家以及韩国的国家养老基金资产变化展开预测的研究(Song et al., 2023)等，强调了多种机器学习算法的综合运用进行比较与挑选，同时尽可能多地纳入各个领域的解释变量，以实现对出生率的预测。相较而言，S. Li et al. (2018) 仅采用了最小二乘支持向量机（LS-SVM）与莱斯利矩阵模型相结合的方式对中国人口结构进行预测。该方法所需假定过多，且仅将性别比作为解释特征，其他经济社会变量均未涵盖。
已有研究在本文的应用方面给予了强有力的指导，主要涵盖三个方面：其一，需对多个机器学习算法的表现进行比较，尤其要重视对新型算法的比较，例如袋装法及梯度速降算法；其二，应包含尽可能多的特征（解释变量）；其三，需进行细致的超参数调优以确定模型。然而，已有研究主要存在解释方面的不足，具体表现为仅体现了全局（Global）解释的重要性，却极少涉及局部（Local）解释，且对于交互项和异质性的说明有所缺位。而这些正是机器学习被诟病为“黑箱”的主要原因所在。
2.2 数据介绍
本文的基础数据源自 2000 至 2020 年中国人口普查分县数据库，涵盖第五次人口普查（2000 年）、第六次人口普查（2010 年）以及第七次人口普查（2020 年）这三次人口普查（以下简称为“五普、六普、七普”）的分县全部指标。本文聚焦于农村地区的少子化问题，故而删除了城市数据（在本文中，将“区”界定为城市）。本文的样本地区包括县、县级市、旗等（以下均简称为“县”）。
基础数据库与其他公开渠道实现对接，可搜集到县域经济、社会、地理、人口、交通、土地利用、灯光亮度以及污染指标（pm2.5）等信息。有效样本保留五普/六普/七普没有行政区划变更的县，保证三次普查的可比性。以确保三次普查具备可比性。最终保留的有效样本县数量为 1857 个。
在考量人口、经济、交通等时变变量于时序上的依存问题时，尤其需注意本文的被解释变量为第七次全国人口普查数据0-14 岁人口比例，其极大程度地依赖于前期的人口学特征。因此，被解释变量以及部分人口结构变量，包括第五次、第六次全国人口普查的变量（即滞后两期和滞后一期的人口学变量）。0-14 岁人口比例乃是衡量人口结构是否合理的重要指标之一。维持适度的该年龄段人口占比，有助于优化人口结构，避免诸如人口老龄化过快、少子化等问题的出现。中国的人口普查每十年开展一次，采用（参考文献：bili）作为地区少子化水平的测量指标，相较于出生率更为适宜，类似的文献亦采用此种方式测量地区少子化水平(丁 et al., 2023; 许 et al., 2024; 马, 2024)。
其他时变的社会经济、交通以及污染等变量均仅选取 2020 年的数值。被解释变量与人口学变量滞后项的引入，实际上充分考量了各县的非观察异质性，那些影响过去且无法观测到的因素同样也对现在产生影响。最终的数据库呈现为以第七次全国人口普查数据为核心的截面数据形式，其中包含 67 个特征（在机器学习中，解释变量集合的术语称为“特征”，即features；被解释变量称为“目标”，即target；全部的特征与目标清单请参见6.3）。观测值的缺失比例为1.1%，比例很小（参见13），在后续的特征工程过程中将进行多重查补处理。
在保留的样本县中，国家级贫困县达 770 个1，贫困县与非贫困县作为重要的分组标志，代表着经济社会、地理区位等方面的广泛异质性。依据此分组对三次普查0-14 岁人口比例，结果呈现在图 1 ，具有三个基本特点：其一，0-14 岁人口比例第五次普查结果最高，第六次普查下降最为明显，第七次普查与第六次普查相比变化不大，并未呈现逐次下降的趋势；其二，五普/六普/七普都是贫困县比例显著高于非贫困县，二者的差距在三次普查并未拉大或者缩小（即所谓的 Difference in Difference）；其三，从分布情况来看，三次普查逐渐由相对集中的分布转向相对分散的分布。就非贫困县而言，第五次/第六次普查的分布集中在中位数附近略高的位置，到了第七次普查则出现向低分位数集中的情况（低于 25 分位附近）。贫困县则从第五次普查在 75 分位数附近的集中转变为第七次普查的非常分散形态，均匀分布在四分位距区间上。
分组比较的检验结果表明，无论是组间差异（即三次普查之间的两两比较，经 Bonferroni 校正后的结果），还是组内比较（同一次普查中贫困县与非贫困县之间的比较），均在统计上具有高度显著性（p 值<0.001）。在图 1 绘制了参考线，用以标识七普全国总体0-14 岁人口比例 （17.95%），由此可以看出，非贫困县的情况与全国水平几乎完全一致，而贫困县的0-14 岁人口比例则明显高于全国水平。换言之，贫困县的（ref:bili）高于全国平均水平。
Figure 1: 五普/六普/七普0-14 岁人口比例 贫困县与非贫困县分组比较
2.3 机器学习方法说明
2.3.1 数据预处理与特征工程
机器学习需历经特定的数据整理过程，该过程被称为“特征工程” (Johnson, n.d.-b)。唯有将数据整理妥善，方可进入机器学习的模型构建阶段。[本文中数据的特征工程依照 mlr3 流程 Bischl et al. (2023) 的管道算子予以处理。特征工程主要涵盖三个环节：
针对缺失值查补事宜，因子变量（即离散变量，当前数据中并无缺失的因子变量）的缺失值查补采用 oor 方法进行。数值（连续）变量则采用直方图查补与随机森林相结合的处理手段，具体为采用 10 次多重查补后取均值的办法。同时，已对查补后的结果进行验证，确认其与原始数据无显著差异。
因子编码可分为以下三种情形进行处理：情形一，对于因子取值超过 10 个以上的特征，采用特征编码方案。例如，本文中的省份/地级市因子均为超过 10 个取值的特征编码。情形二，当因子具有 2 至 10 个取值时，采用独热编码（One-hot encode），即将每个取值展开成为哑变量，但不删除对照组。例如，具有 4 个取值的东中西和东北地区分组可展开成为 4 个独热编码。情形三，二分变量直接作为哑变量。
进行标准化处理。为确保特征具有可比性，对所有连续变量特征实施标准化处理，以使其均值达到 0，标准差变为 1。
在实际操作过程中，上述环节通过采用管道算符进行链接，以构建图算子运行(Binder et al., 2021)。
经过特征工程的处理，获取了符合机器学习建模规范的数据库。随后转入正式建模阶段。在机器学习中，根据目标的不同属性，建模任务分为分类与回归两种。本文所涉及的特征（参考来源：bili）为连续变量，故采用的建模任务类型为回归任务（标识为 regr）。
2.3.2 机器学习回归任务建模技术路线
依照惯例，将 70%的样本设定为训练集，30%的样本用作测试集。基本的建模流程与技术路线在图 2中予以呈现。
Figure 2: 机器学习回归任务建模的技术路线
3 建模过程
图 2 展示了正式的机器学习技术路线。首先排除第五次、第六次人口普查样本，保留第七次人口普查样本，从而获得 1857 个有效观测值。其中，1269 个县被用作训练集，588 个县作为测试集，且训练集与测试集不存在任何交集。基于训练集进行建模，基于测试集进行模型验证与评价，主要测量指标为 RMSE。在确立任务并划分训练集/测试集之后，分两步进行操作。
第一步为特征筛选。调用 mlr3 的集成自动挑选流程，采用三次嵌套的五折交叉验证方式，运用递归特征消除（RFE）结合嵌入方法实施特征选择。其原理在于持续带入并挑选特征，将那些对均方根误差（RMSE）的下降几无帮助的特征予以删除，而保留下来的特征即为建模纳入的最终特征(Johnson, n.d.-a; 黄, 2024年7月)。最终结果实现了最小的均方根误差（RMSE）。从总计 67 个特征中，最终确定了 11 个特征用于建模。
第二步为算法选择。基准算法采用普通线性回归，常规算法选取了 kknn（k 近邻）、svm（支持向量机）、nnet（神经网络）以及 ranger（随机森林）这四种算法；同时还有三种表现较为优异的袋装梯度速降算法，即 LightGBM、XGBoost 以及最新的 CatBoost 这三种算法。基于十折交叉验证，RMSE分布比较结果在3呈现。其中表现最为出色的为 CatBoost 算法，其 RMSE 的中位数小于 2，而其他算法的中位数均大于 2。以 R-squared、MAE 进行测量的结果亦类似（结果略），故而选定 CatBoost 作为建模算法。
Figure 3: 算法选择结果 基于10折交叉验证的RMSE结果分布比较
- 第三步为超参数调优。在确定算法之后，针对 CatBoost 的超参数展开调优工作，其目的在于获取最终的算法建模结果。依据通常的处理方式2，选取三个超参数进行调优，分别为：learning_rate、depth 以及 rsm。采用嵌套的三折交叉验证（外层迭代次数为 4），进而确定三个超参数的最优取值为：learning_rate=0.329、depth=9、rsm=0.493。
在获取超参数的最优取值之后，将其引入学习器以进行模型训练。随后，使用测试集进行预测评价，拟合效果甚佳（详见附图12)），所得 RMSE 为 2.151，R-squared 为 0.842。
在模型确定之后，接下来将进入对该模型的阐释环节。
4 模型解释
可解释机器学习（Interpretable Machine Learning）乃是机器学习的重要领域之一。其聚焦于如何促使机器学习模型的预测结果与决策过程更易于为人类所理解及阐释。传统的机器学习模型往往被视作“黑箱”，虽能做出精准预测，却难以阐明为何会做出如此预测。这致使人们在运用机器学习模型时，可能对其结果心存疑虑，尤其是在诸如医疗、金融以及法律等关键领域。可解释机器学习的目标在于开发方法与技术，以使得机器学习模型的预测结果及决策过程能够以可理解的形式呈现在人类用户面前。如此一来，用户便能更好地理解模型的工作原理，评估其可靠性与公正性，并在必要之时对其进行调整与改进。
在展开解释之前，首先针对本文所涉及的机器学习模型所确定的 11 个特征进行简单的描述统计分析。
4.1 特征描述与符号对应说明
实际投入模型的数据集乃是经过特征工程处理后的标准化特征，其均值为 0，标准差为 1。为能更为直观地了解特征的基本状况，在此呈现未经标准化的数据库基本描述结果（中位数、IQR、均值/标准差，详见表1）。
模型特征描述统计明确了所选取的 11 个特征。其中，第一组为人口学变量的滞后项，涵盖了第六次全国人口普查的0-14 岁人口比例 ，但未包含第五次全国人口普查数据。即 10 年前的0-14 岁人口比例 具有一定影响，而 20 年前的0-14 岁人口比例 几乎不具备任何预测意义，从一个侧面反映出人口转变速度较快。2016 年全面放开二孩政策对于农村户籍人口实际意义不大，在许多地方，农村户籍在 2016 年以前便可以生育二孩。这或许是第七次与第六次全国人口普查中0-14 岁人口比例分布变化不大的原因之一。同样，第六次全国人口普查的户规模具有影响，而第七次与第五次全国人口普查的户规模未被纳入特征之中。值得注意的是，男性第五次全国人口普查的受教育年限与女性第六次全国人口普查的受教育年限被选中，可理解为夫妻教育程度对生育行为具有重要影响，这也是众多已有研究已证明的影响因素（参见 Roser (2024) 的综述）。第五次与第六次全国人口普查的 0-14 岁人群即为第七次全国人口普查的育龄人口。当然，夫妻之间还存在年龄差距问题，一般而言，男性年龄大于女性年龄。第五次全国人口普查的性别比亦是被选中的影响因素之一。
第二组的主要特征体现为地理与交通方面的特征，具体涵盖海拔、经度、地形起伏度以及与所属地级市之间的距离。而其他诸如经济社会、灯光亮度、污染、土地利用、乡村数字化、普惠金融等要素均未被选取。在分组标志方面，仅有是否为贫困县这一要素入选，省、市特征编码均未纳入模型之中。
对特征的选择进行简单总结。除贫困县分组被选中之外，解释0-14 岁人口比例的特征主要分为两大类，即人口学特征与地理区位特征。此外，在当前的机器学习框架下，已有研究中常见的经济社会解释变量对于农村地区少子化现象（在包含已选中特征的情况下）均不具备预测能力。
表1提供了特征的符号、具体含义以及基本描述统计量。鉴于当前的机器学习程序不支持中文特征命名与标签，仅能带入符合规定的特征符号（通常为英文与数字的组合），部分图示仅能呈现特征符号，故而需对照表1以理解特征的确切含义。
| 特征符号 | 特征含义 | Median [Q1, Q3] (Mean,SD); n (%) |
|---|---|---|
| y | 七普14岁以下人群百分比 | 19.6 [15.6, 23.7] (19.8,5.4) |
| poorcounty | 贫困县 | NA |
| poorcounty | 非贫困县 | 1,087 (59%) |
| poorcounty | 贫困县 | 770 (41%) |
| cen_01 | 总人口数 | 323,031 [178,368, 540,572] (400,025,309,714) |
| fedu_10 | 女性六普平均受教育年限 | 7.86 [7.19, 8.37] (7.60,1.29) |
| height | 海拔 | 282 [66, 1,017] (718,1,002) |
| hsize_10 | 六普户规模 | 3.31 [3.01, 3.62] (3.35,0.54) |
| jv1 | 经度 | 112 [105, 117] (110,10) |
| jv5 | 与所属地级市距离 | 53 [32, 84] (64,48) |
| medu_00 | 男性五普平均受教育年限 | 7.72 [7.22, 8.03] (7.39,1.20) |
| qq1 | 地形起伏度 | 0.73 [0.18, 1.51] (1.16,1.32) |
| sexr_00 | 五普性别比 | 104.1 [100.5, 107.6] (104.3,5.8) |
| y_lag1 | 六普14岁以下人群百分比 | 18.4 [15.4, 22.2] (19.0,4.8) |
所有特征（包括目标）的相关性分析可参考图4。此处所使用的特征均为经过标准化处理后的变量。图4的下三角部分呈现相关系数，几乎所有特征的简单相关系数在统计上均具有显著异常性。对角线部分为特征的核密度（经标准化处理后的特征均值为 0，标准差为 1）；上三角部分为两两之间的散点图与 lowess 回归结果。图4暗示了机器学习方法的两大优势：其一，特征之间相关性普遍存在，且相关性基本显著，传统统计模型系数解释中“假定其他条件不变”的设定场景几乎难以成立；其二，非线性普遍存在。相关性与普遍存在的非线性特征相互叠加，使得即便包含交互项的传统建模也成为一种不切实际的过度简化。
Figure 4: 特征变量相关性呈现
4.2 基于置换方法的全局解释
机器学习的解释主要通过两种方式予以实现。其一为基于置换方法的解释；其二为基于 SHAP 方法的解释。就基于置换方法的解释而言，其主要用于评估特征对模型预测的重要程度。其基本理念是通过随机打乱（置换）某一特征的值，进而观察模型性能的变动情况。若在打乱特征值后，模型性能出现显著下降，则表明该特征对于模型的预测具有重要意义；反之，若模型性能未呈现明显变化，则说明该特征对于模型的预测并非十分重要。此方法的优势在于相对简洁直观，且无需对模型的内部结构有深入认知。
首先，运用置换方法对特征重要性与交互重要性展开分析。经计算所得的特征重要性及其交互项之间的重要性，直接通过网络图予以概括（参见图5）。图中的 Vimp（即 Variables importance 的含义，用于标识节点色域），节点色域结合节点大小用以说明特征重要性，颜色越深且尺寸越大的节点，代表单个特征的重要性越高。Vint（即 Variables interaction 的含义，用于标识连线色域），节点连线色域及粗细程度用以说明特征之间交互项的重要性。
（参见图@ref(fig:fig5)）。
Figure 5: 特征重要性网络结构
图5的最重要4个特征分别是六普0-14 岁人口比例 （即y_lag1，重要性均值为4.01），县质心点经度（即jv1，重要性为2.11），五普16岁以上男性受教育平均年限（即medu_00，重要性为1.786），六普家庭平均规模（即hsize_10，重要性为1.606）。其余特征重要性都小于1.4。
选取最重要的4个特征的ICE+PDP在图 6呈现。ICE（Individual Conditional Expectation）和 PDP（Partial Dependence Plot）是两种用于解释模型行为和理解特征影响的方法。PDP（Partial Dependence Plot，部分依赖图）：展示了一个或多个特征对模型预测结果的平均边际效应。通过固定其他特征的值，改变目标特征的值，并计算模型的平均预测输出，从而描绘出目标特征与模型预测输出之间的关系3。
Figure 6: 用于测试集的模型训练结果
图 6呈现出两个要点：第一，最重要特征的非线性影响。4个特征几乎都呈现非线性影响。最重要的六普0-14 岁人口比例在一个标准差以前都是一种正向的近乎线性的影响，但大于1个标准差以后的取值范围上（对照表1，拐点大致在六普0-14 岁人口比例为23.8%以上），拐点以上，本特征对于预测值的影响就几乎消失了（接近水平）。县质心经度来说，在115度附近出现一个影响的拐点，持续到120度附近重要性有明显下降。也就是山东/河南到浙江/江苏一段（桐乡/滨海/东台）是少子化的分界，过了分界区以后向沿海的区域，七普的0-14 岁人口比例比例明显降低。男性五普受教育年限则是在6.1年以后出现降低的拐点（大致相当于小学毕业）；六普的家庭规模则是在2.8-3.8人阶段出现一个上升爬坡，其他取值上几乎没有明显影响；第二，ICE不见得完全平行，即便在样本密集的取值区域（横轴的rug显示了样本点的实际分布）也不平行，这样暗示了交互项作用不能被忽视。
PDP假设特征之间相互独立，可能会掩盖特征之间的交互作用(Friedman, 2001)；ICE的非平行特点揭示了交互项不可忽视。下面转入对交互项重要性的简要分析。
交互项重要性计算结果在图7呈现。左图是每个特征的全部交互项的重要性加总排序的情况（也就是图5各个节点发出的连线合计），右图是最重要的特征六普0-14 岁人口比例全部交互项重要性（也就是图5y_lag1发出的连线）。一个基本的特点就是特征的自身重要性越高，交互项的重要性也就越大。
Figure 7: 交互项重要性：综合与配对
选取最重要的六普0-14 岁人口比例与五普男性平均受教育年限与县经度交互项计算结果在图8加以呈现。结果表明：六普0-14 岁人口比例越高预测七普0-14 岁人口比例也越高，但是这种正向关系强烈依赖于五普男性平均受教育年限与县经度。具体而言，五普男性受教育程度低会强化这种正向影响，教育程度高则会弱化六普对于七普的预测效果。县经度作用更加明显，经度越大的县，六普对于七普0-14 岁人口比例预测能力越是弱化。经度越小的县，预测能力越强。
Figure 8: 六普0-14 岁人口比例与五普男性平均受教育年限和县经度交互项特征
4.3 基于SHAP的全局解释
前面置换方法存在两个主要缺点：假设特征之间相互独立，在实际情况中这个假设往往不成立，可能导致不准确的解释。只能提供关于特征整体重要性的信息，难以解释单个样本的预测。下面转入SHAP（SHapley Additive exPlanations）方法解释进行补充论证。SHAP 值的基本思想来源于博弈论中的 Shapley 值，用于公平地分配每个特征对模型预测的贡献。主要优点体现在：能够为每个样本的每个特征提供一个具体的贡献值（SHAP 值），有助于解释单个预测。 具有较好的理论基础和数学解释。但是计算复杂度相对较高，特别是对于大型数据集和复杂模型。
本文使用核方法计算SHAP值的结果在图9呈现。SHAP 值表示每个特征对模型输出的贡献程度。对于给定的输入样本和预测模型，每个特征都有一个对应的 SHAP 值。 正的 SHAP 值表示该特征增加了模型的输出值，负的 SHAP 值表示该特征降低了模型的输出值。SHAP 值的大小反映了特征的影响力，绝对值越大，说明该特征对模型输出的影响越大。
图9左图是SHAP值的均值绝对值，特征重要性结果和置换方法几乎一致，六普0-14 岁人口比例的重要性远远超过其他特征。右图蜂群图是SHAP方法的优势所在，呈现了每个样本的特征贡献全景。例如六普0-14 岁人口比例的蜂群显示出一致的模式，特征的值越大，预测的值也就越大，而且呈现均匀的分布影响。县经度（jv1)则是特征值越大，则预测值越小，但影响点群聚在较小的SHAP正值附近，显得很不均匀。意味着县经度对模型的输出有一定的正向推动作用，但作用程度相对较小。其他特征的解释类似，几乎都有在较小的SHAP值附近的群聚，使得对于预测的重要性不高。
Figure 9: 特征重要性的SHAP值结果
六普0-14 岁人口比例与其他三个重要特征的交互影响 SHAP 值结果参见图10。左图表明：重要性高的影响点是六普0-14 岁人口比例高且经度相对较小（均值附近和低于均值）的一些县。中图表明：重要性高的影响点是六普0-14 岁人口比例且六普平均家庭规模也高的一些县。右图表明：在六普0-14 岁人口比例均值加上一个标准差的取值左侧，五普男性教育年限几乎没有什么影响。在六普0-14 岁人口比例均值加上一个标准差的取值右侧两个变量才开始产生交互作用，男性教育年限少且六普比例高这些县重要性比较高。
Figure 10: 三个重要特征依赖的SHAP值
SHAP交互作用的计算结果进一步表明，特征之间的依存普遍存在，而且非线性的特征与异质性特征明显。采用平均意义上的解释现实意义不大。反倒是针对具体地区的分类施策更有意义一些，例如教育资源养老资源规划布局等等。
4.4 局部解释
SHAP方法的一个优势就在于提供了局部解释和预测；本文选取两个具有典型性的县样本为例，局部解释参见图11。两个县的0-14 岁人口比例平均预测值（期望）都是19.8%，A县实际值为21.9%，B县实际值为12.0%。
Figure 11: 选取两个样本县的SHAP解释
瀑布图对A县的特征解释：从均值出发（最下方），两个不重要的特征增加了预测值，与所属地级市的距离（2.95，远高于均值）增加了预测值；不是贫困县（poorcounty=0）减少了预测值；海拔（0.198，略高于均值）降低了预测值（0.334个百分点）；逐渐向上，经度（-2.7）为内陆地区，使得预测值增加了1.12个百分点；六普0-14 岁人口比例（0.266）高于平均水平，使得预测值提高1.81个百分点；最终得到该县七普的0-14 岁人口比例实际值21.9%。
B县的瀑布图解释类似，B县更加靠近沿海（经度高）六普0-14 岁人口比例较低，还有地形起伏度较小（-0.367，低于均值），这几项主要作用使得最终B县的七普0-14 岁人口比例仅为12.9%。
常用的一种局部解释手段还有反事实论证（counter-factual）。例如针对上面的B县，怎样改变特征才能增加0-14 岁人口比例到某个合意水平。但仔细检视全部11个特征。这些特征几乎都是严格外生（地理区位等）或者前定的变量（六普0-14 岁人口比例/是否贫困县/教育年限等）。这些特征在现实中几乎无法改变，强行的设定所谓的反事实，就是荒谬。因此，“反事实论证”在本文场景下无意义。
5 结论与评述
本文采用机器学习方法对中国农村地区的七普0-14 岁人口比例进行了预测与解释，在包含了大量区域经济社会/地理区位/土地环境/人口学特征和省市异质性等特征以后，发现具有预测能力的特征包含11个变量，主要是前两次普查的人口学特征和地理区位变量。本文比较了8种机器学习算法的表现，确定了Catboost作为最终建模的选择，在进行了细致的超参数调优过程以后确定了模型。进过测试，模型的表现非常好。一般来说，R² 值大于 0.7 可能被认为是一个较好的拟合(Bischl et al., 2023)，本文的R² 值达到了0.842；并且与基准模型和传统的近邻回归/随机森林比较都有明显的改善。
本文采用置换方法和SHAP方法对模型的结果进行了全局和局部的解释，主要发现包括：其一，全局意义上，最重要的特征是六普0-14 岁人口比例，重要性大幅领先于其他所有特征；其次是五普的男性平均受教育年限/县经度/六普户规模这三个特征，其他特征的重要性相对比较低；其二，特征之间的非线性与交互效应普遍存在。线性设定与假定其他条件不变这些传统建模假设条件几乎无法成立。最重要的两个交互效应计算表明六普0-14 岁人口比例对七普0-14 岁人口比例有明显的正向影响，但是这种正向关系强烈依赖于五普男性平均受教育年限与县经度，教育年限低，影响更强；经度越小（靠近内陆），影响越强；其三，个体异质性普遍存在，很难简单使用平均趋势来概括全部的地区的人口少子化情况，SHAP值的分布体现了个体差异性极大，同时具有强烈的非线性特征。
本文的研究对于已有文献的拓展主要体现在两个方面。
第一方面是解释。在理论指导实证建模不够清晰与准确的前提下，采用机器学习方法来理解中国农村的少子化现状。已有研究提供的信息最重要的就是：人口少子化影响因素很多，是一个长期的复杂过程。机器学习适合于影响因素很多，且相互关联比较复杂的情形。本文解释部分表明滞后一期的少子化情况对于当期具有最重要的影响，传递了两个信息，一是人口学转变过程具有重要的路径依赖特点，现在的少子化状态依赖于过去；而且SHAP值个体分布也表明，无论过去的状态是高还是低，对于现在的少子化都有非常显著的影响作用；二是中国农村地区现在的少子化过程是快速的，20年前的五普少子化情况对于当下的预测能力就几乎可以被忽略掉了。靠近沿海和东北的高经度地区，少子化的转变更加迅速。育龄男女的受教育程度也是重要预测特征这也是被以往研究所反复证明过的。包含大量特征以后，表示经济和收入水平的因素，一个都没有被模型选中，这和已有研究是有差别的。一个理解是特征之间普遍存在的相关性（共线性），另一个理解则是系统的自适应，也就是说家庭生育决策根据经济状况/教育/社会保障/观念/居住/就业情况等方面进行了调整。因此在跨度为10年的普查数据解释上，这些因素的预测能力就被内化了。
第二方面是预测。在普遍存在非线性与相关性的情况下，获得个体的理解和预测更重要。比起：“假定其他条件不变，某个因素变化，将会使得因变量平均变化多少”的传统建模解释要现实得多。SHAP方法提供的针对单个县样本的局部解释是全局解释的必要补充，也是预测今后一段事情中国农村少子化的模型基础。特别是针对到具体的县来说。提供一个分县的，到2030年的0-14 岁人口比例预测模型api也是本文的一项重要产出。这是已有研究没有涉及的。
政策含义有两点：第一是在目前的状态下，农村人口的少子化趋势几乎无法改变，这些具有解释力的特征几乎都是无法短时期操控的变量；短期最好的方式是顺势而为。也就是社会对于农村少子化的自适应调节。包括教育资源/生育政策/产业转型/医疗与养老等的区域布局和优化；第二是尽快确立大数据大模型为底座的农村区域人口少子化预测模型，这样才能前瞻性的具备针对各个县的精准施策与制定规划。仅仅依靠单一学科的理解是很难获得相对可靠的预测结果的4，必须结合新的数据和新的工具。
本文包含的纳入特征挑选的变量局限于公开渠道可以获得的数据，一些重要的变量，例如迁移/就业情况等缺位；数据方面，使用历史数据可能无法完全涵盖未来的社会经济变化。模型方面，只研究了有限的机器学习算法和超参数，可能错过了更有效的模型或配置。另外，城市的人口变迁情况更加复杂，变化更大，本文没有纳入考虑，而且模型对于城市少子化的泛化也就存疑。这些是本文的主要不足。机器学习/深度学习技术的变化非常迅猛，人口学的预测与人口密度/灯光/交通等卫星图像技术相结合也是未来本领域的一个趋势。
6 附录
6.1 模型预测结果 附图12
Figure 12: 附图1 模型预测结果图示
6.2 缺失值概况 附图13
Figure 13: 附图2 全部变量缺失值概况
6.3 全部变量清单
y_lag1: 六普14岁以下人口百分比; y_lag2: 五普14岁以下人口百分比; cgroup: 县和县级市分组; region5: 地区分组; province: 省份; city: 地级市; v002: 行政区域土地面积(平方公里); v034: 地方财政一般预算收入(万元); v036: 地方财政一般预算支出(万元); v037: 城乡居民储蓄存款余额(万元); v077: 夜间灯光平均值; v079: 夜间灯光标准差; v082: 数字乡村指数; v083: 数字基础设施指数; cen_01: 人数合计; cen_06: 少数民族人口比重%; qq1: 地形起伏度; qq2: 距离港口的最近距离 (公里); jv1: 县经度; jv2: 县纬度; jv5: 与所属地级市距离; jv8: 与所属省会城市距离; pm25: pm2.5均值; maxpm: pm2.5最大值; minpm: pm2.5最小值; lw2: 路网总长度; lw3: 路网密度; cropland: 农田; forest: 森林; grassland: 草原; water: 水域; height: 海拔; ph01: 数字普惠金融总指数; oldrate_00: 五普65岁以上人口百分比; oldrate_10: 六普65岁以上人口百分比; oldrate_20: 七普65岁以上人口百分比; medu_00: 男性五普平均受教育年限; medu_10: 男性六普平均受教育年限; medu_20: 男性七普平均受教育年限; fedu_00: 女性五普平均受教育年限; fedu_10: 女性六普平均受教育年限; fedu_20: 女性七普平均受教育年限; marr_00: 五普有配偶人口比例; marr_10: 六普有配偶人口比例; marr_20: 七普有配偶人口比例; sexr_00: 五普性别比; sexr_10: 六普性别比; sexr_20: 七普性别比; hsize_00: 五普户规模; hsize_10: 六普户规模; hsize_20: 七普户规模; houser_00: 五普自有住房比例; houser_10: 六普自有住房比例; houser_20: 七普自有住房比例; urbanpopu: 城镇人口百分比; agripopu: 第一产业人员占比%; pergdp: 人均GDP（万元); gdpr1: 一产比例; gdpr2: 二产比例; gdpr3: 三产比例
程序的编译时间为：285.761秒。
参考文献
2011 年确定的国家级贫困县数量为 832 个，涉及全国中西部 21 个省区市. 经过多年的脱贫攻坚努力，2020 年 11 月 23 日，所有国家级贫困县全部脱贫摘帽↩︎
参阅官方说明，网址：https://github.com/catboost/catboost。2024年12月19日接入。↩︎
ICE 曲线为每个样本绘制了一条线，展示了当改变某个特征的值时，模型预测是如何变化的。当样本数量较多时，图形可能会变得杂乱，难以解读。通常，PDP 和 ICE 会结合使用，以更全面地理解模型中特征的作用和影响。↩︎
人口学界对于全面二孩政策效果的误判便是明证。本文的结果已经表明，除了人口学特征以外，其他的特征也会明显影响农村少子化的预测，这些特征究竟是地理学/经济学还是人口学，很难加以明确的区分。加之，系统的自适应过程起作用，简单的/平均意义上的干预效果评测容易产生误判，导致决策的失误。↩︎