简介:《河南省区域经济发展差异的统计分析》是一份面向经济学、区域经济与统计类专业学生及研究者的分析类文档,聚焦改革开放以来河南省内各区域发展不平衡问题,围绕产业结构、基础设施、人力资源、教育投入和城镇化率等维度梳理差异现状,并采用主成分分析法,从科技投入、固定资产投资、教育投入、人口素质、城镇化率、外贸与就业等变量中提取关键影响因子,衡量各因素的作用程度,进而提出优化资源配置、强化科技创新、改善基础设施、提升教育水平、促进城镇化健康发展与扩大对外开放等对策建议,适合课程论文写作、区域经济课题研究及统计指标方法参考。压缩包内为1个doc文档,约590KB,可直接查阅与引用。目前已有108人浏览学习,可作为区域经济差异测度与主成分分析应用的参考范本。
1. 河南省区域经济发展差异的统计分析:一份 .doc 需求背后的数据管线
需求方往往只丢过来一句话——"河南省区域经济发展差异的统计分析",交付物是一份 .doc。真正耗时间的不是写字,而是把 18 个统计单元、十几年的面板数据对齐口径:同一年的人均 GDP,统计年鉴和统计公报能差出小几百元,换一次常住人口来源,基尼系数就跳一下。这篇内容讲的是一条能跑的链路:先定统计单元和指标口径,用 pandas 把年鉴宽表拗成面板长表,再用 Python 复现极差、变异系数、基尼系数和泰尔指数的组间组内分解,最后补一层空间自相关,看差异是散着还是聚着。适合手里握着统计年鉴、需要把结论落到表格和图上的人,也适合想拿一份真实面板练手空间计量的开发者。
2. 统计单元与指标口径:河南区域差异分析的数据地基
2.1 18 个统计单元和县域两套口径怎么选
河南省的统计口径里,常见的区域分析单元有两层:一层是 18 个统计单元,即 17 个地级市加上济源示范区;另一层是县域口径,包含市辖区、县级市和县,数量在百位级别。两套口径各有各的用处——市级口径适合做跨地市比较和空间自相关,县域口径适合看"市域内部是不是也在分化"。
我一般会先做市级面板跑通全流程,再决定要不要下钻到县域。原因是县域口径的年鉴数据缺失更严重,不少县的规上工业增加值、R&D 经费支出会有整年空档,插值补出来的数放进泰尔指数分解里,组内差异会被稀释,结论容易失真。
提示:济源示范区在部分年份的年鉴里单列,在部分年份的省级汇总里并入其他口径。做面板前先确认这一列在目标年份区间内是否连续,不连续就单独标注,不要默默插值。
行政区划调整是另一个坑。个别县市在观察期内发生过隶属关系或统计归属的变化,表现为某个市的常住人口或 GDP 出现一次台阶式跳变。处理办法是在面板里加一列caliber_note,把断点年份记下来,做趋势图时用竖线标注,做增长率时把断点那一年的增速置为缺失。
2.2 指标选择:人均 GDP 之外还该放什么
只用人均 GDP 一个指标做区域差异分析,结论很容易被质疑"单一维度"。实测下来,一个能撑住 .doc 的指标体系至少包含产出、结构、收入、财政四类。
| 指标 | 口径 | 单位 | 常见坑 |
|---|---|---|---|
| 人均 GDP(现价) | GDP / 常住人口 | 元/人 | 常住人口口径前后有调整 |
| 人均 GDP(不变价) | 用 GDP 平减指数折算 | 元/人 | 市级平减指数缺失,只能用全省 |
| 三次产业结构 | 各产业增加值占比 | % | 部分年份分行业数据有缺 |
| 城镇化率 | 城镇常住人口 / 常住人口 | % | 与户籍城镇化率混用 |
| 人均可支配收入 | 全体居民口径 | 元/人 | 2013 年前按城乡分列,不能直接续接 |
| 人均一般公共预算收入 | 一般公共预算收入 / 常住人口 | 元/人 | 与财政总收入混淆 |
人均可支配收入这一项要特别注意口径断裂:早期年份城镇和农村分开统计,全体居民口径是后来才有的。跨断点做时间序列,要么只用断点之后的年份,要么把城镇、农村两组分别分析,不要用插值硬接。
2.3 用 pandas 把年鉴宽表拗成面板长表
年鉴导出的表格通常是"一行一个地区、一列一个年份"的宽表,先把它转成city_code, year, indicator, value的长表,后续所有计算都基于这一张表。
import pandas as pd import numpy as np # 年鉴导出的宽表:一行一个统计单元,列名形如 "2015年" raw = pd.read_excel("henan_yearbook_wide.xlsx", sheet_name="人均GDP") raw = raw.rename(columns={"地区": "city_name"}) long = raw.melt(id_vars="city_name", var_name="year", value_name="value") # 列名里抽出四位年份,非年份列(备注、数据来源)会被抽成 NaN 后丢弃 long["year"] = pd.to_numeric( long["year"].astype(str).str.extract(r"(\d{4})")[0], errors="coerce" ) long = long.dropna(subset=["year"]).copy() long["year"] = long["year"].astype(int) # 年鉴里的 "—"、"…" 先转成 NaN,再统一处理 long["value"] = pd.to_numeric(long["value"], errors="coerce") long["city_code"] = long["city_name"].map(CITY_CODE) # 行政区划代码对照表 long["indicator"] = "gdp_pc_current" panel = long[["city_code", "city_name", "year", "indicator", "value"]] panel = panel.sort_values(["city_code", "year"]).reset_index(drop=True) assert panel["value"].notna().mean() > 0.95, "缺失率过高,先回去核对年鉴页码"逻辑上有三处值得说明。str.extract(r"(\d{4})")解决的是列名不统一的问题,年鉴里会出现"2015年""2015 年""2015"几种写法,用正则一把抓比手工改列名稳。errors="coerce"是把所有非数值内容统一变成NaN,避免后面算均值时字符串参与运算。最后的assert是关键,缺失率超过 5% 时通常说明表格抓错了页或者列名映射出错,这时候停下来核对,比后面带着一行NaN跑完全流程要省事。
2.4 可比价换算和缺失值到底怎么补
跨年份比较人均 GDP,必须处理价格因素。市级的 GDP 平减指数一般查不到,常见做法是用全省的平减指数统一折算,换来的是"相对差距"可比较,代价是各市自己的价格结构差异被抹平。这一步要在 .doc 的数据说明里写清楚。
# deflator.csv: year, deflator(以 2015 年为 1.0) deflator = pd.read_csv("henan_deflator.csv") panel = panel.merge(deflator, on="year", how="left") # 现价 -> 2015 年不变价 panel["gdp_pc_real"] = panel["gdp_pc_current"] / panel["deflator"] # 只对中间年份做线性插值,首尾年份的缺口不补,留缺失 panel = panel.sort_values(["city_code", "year"]) panel["gdp_pc_real"] = ( panel.groupby("city_code")["gdp_pc_real"] .transform(lambda s: s.interpolate(method="linear", limit_area="inside")) )limit_area="inside"这个参数很实用:它只填中间的缺口,序列开头和结尾的缺失保持原样。理由是首尾缺失往往意味着这个统计单元在观察期两端不在口径内,插值会造出不存在的趋势。中间年份的单年缺失则多半是漏报或未发布,线性插值的影响可控。
3. 差异测度:极差、变异系数、基尼系数与泰尔指数分解
3.1 极值比和变异系数:先给出最直观的一个数
极值比(最高单元人均 GDP / 最低单元)是最好解释的指标,一句话就能说完,也最容易被追问。它的毛病是对极端值极其敏感,郑州和济源之间的比值变化,可能只是某一年的一个单元口径调整,跟整体分化程度关系不大。
变异系数用标准差除以均值,无量纲,可以跨指标比较——比如产业结构差异和收入差异谁更大。但它假设的是正态分布,区域经济数据通常右偏,变异系数会系统性低估真实离散程度。我一般把这两个数当作"开场白",真正下结论还是看下面两个。
3.2 基尼系数:人口加权与不加权的差别
区域基尼系数和收入分配里的基尼系数不是一回事。区域基尼系数衡量的是"各统计单位人均产出"的分布,必须用人口加权,否则一个人口几百万的地级市和一个几十万人的单元在曲线里权重相同,结果会明显偏高。
| 测度 | 量纲 | 可分解 | 极端值敏感度 | 适合场景 |
|---|---|---|---|---|
| 极值比 | 倍数 | 否 | 极敏感 | 对外一句话结论 |
| 变异系数 | 无量纲 | 否 | 较敏感 | 跨指标横向比较 |
| 基尼系数 | 0~1 | 否 | 中等 | 与收入分配研究对齐 |
| 泰尔指数 T | 非负 | 组间 + 组内 | 对高值敏感 | 区域板块分解 |
| 泰尔指数 L | 非负 | 组间 + 组内 | 对低值敏感 | 关注落后地区 |
3.3 泰尔指数分解:组间差异到底占多少
泰尔指数 T 的价值不在它本身的大小,而在可加分解。把 18 个统计单元按豫北、豫西、豫中、豫东、豫南这类板块分组后,总差异可以写成组间差异加上各组的组内差异按产出份额加权求和。这个分解能直接回答"D 型分化还是板块分化"——如果组间贡献率长期在三成以下,说明差异主要来自板块内部,跨板块的平衡政策作用有限;如果组间贡献率持续上行,那板块之间的落差才是主要矛盾。
泰尔指数 T 的表达式里,份额比取对数,高值单元贡献被放大;泰尔指数 L 反过来对低值单元更敏感。两个指数一起报,能避免只报一个时被质疑结论与指数选择绑定。
3.4 一次算完:measure.py 的完整实现
import numpy as np def gini(y, w=None): """人口加权的区域基尼系数。y 为人均产出,w 为人口。""" y = np.asarray(y, float) w = np.ones_like(y) if w is None else np.asarray(w, float) order = np.argsort(y) # 必须按人均产出升序排 y, w = y[order], w[order] cum_w = np.cumsum(w) / w.sum() cum_y = np.cumsum(y * w) / (y * w).sum() return float(1 - 2 * np.trapz(cum_y, cum_w)) # 洛伦兹曲线下方面积 def theil_t(y, w=None): y = np.asarray(y, float) w = np.ones_like(y) if w is None else np.asarray(w, float) s_y, s_w = y / y.sum(), w / w.sum() return float(np.sum(s_y * np.log(s_y / s_w))) def theil_decompose(y, w, group): """组间 + 组内分解。group 为板块标签数组。""" y, w, g = np.asarray(y, float), np.asarray(w, float), np.asarray(group) s_y, s_w = y / y.sum(), w / w.sum() total = float(np.sum(s_y * np.log(s_y / s_w))) between, within, detail = 0.0, 0.0, {} for gk in np.unique(g): m = (g == gk) sy, sw = s_y[m].sum(), s_w[m].sum() between += sy * np.log(sy / sw) # 把组当成一个"超级单元" t_g = float(np.sum(s_y[m] / sy * np.log((s_y[m] / sy) / (s_w[m] / sw)))) within += sy * t_g # 组内指数按产出份额加权 detail[gk] = {"theil": t_g, "share_y": float(sy), "share_p": float(sw)} return {"total": total, "between": between, "within": within, "between_ratio": between / total, "detail": detail}三个函数的参数含义是一致的:y传人均不变价 GDP 乘人口得到的产出总量,w传常住人口。注意y用总量而不是人均值——泰尔指数的份额比在总量和人均值上结果相同,但基尼系数的加权逻辑必须建立在总量上,否则权重会被重复计算一次。
rows = [] for yr, sub in panel.groupby("year"): sub = sub.dropna(subset=["gdp_real", "pop"]) r = theil_decompose(sub["gdp_real"], sub["pop"], sub["region"]) rows.append({ "year": yr, "gini": gini(sub["gdp_real"] / sub["pop"], sub["pop"]), "theil": r["total"], "between_ratio": r["between_ratio"], }) summary = pd.DataFrame(rows).set_index("year")跑完之后先看between_ratio的时间序列,再看gini。两者同时上行,说明板块间落差在扩大;只有基尼上行而组间贡献率平稳,说明是板块内部头部单元拉开的差距。这两种情形的应对思路完全不同,这也是为什么单报一个基尼系数不够用。
4. 空间视角:Moran's I 与 LISA 揭示的河南经济集聚
4.1 空间权重矩阵:邻接、距离衰减与 K 近邻怎么选
常规的差异测度默认各单元彼此独立,可地理上相邻的市往往同涨同落。空间权重矩阵就是把"相邻"这件事写成一个 n×n 的矩阵。三种常见构造:
一是邻接矩阵,共用边界记为 1,其余为 0,符合"行政边界相邻"的直觉,但河南西部山地单元邻居少,东部平原单元邻居多,度数差异大会带来偏误。二是距离衰减矩阵,权重取1/d^2或exp(-d),能刻画远距离影响的衰减,但要选一个截断距离,主观性强。三是 K 近邻矩阵,每个单元固定取最近的 K 个邻居,度数一致,适合单元规模悬殊的场景。
我一般主用 K 近邻(K 取 4 到 6),把邻接矩阵的结果作为稳健性对照。两种矩阵都做行标准化,让每行权重和为 1,避免邻居多的单元把 Moran's I 拉高。
4.2 全局 Moran's I 的手写实现
依赖 libpysal 当然更省事,但版本变动和geopandas的编译问题经常卡住环境。用经纬度坐标手写一遍,二十行代码就够,也更容易在 .doc 的方法部分讲清楚。
import numpy as np def knn_weight(coords, k=5): """coords: (n, 2) 经纬度或投影坐标。返回行标准化权重矩阵。""" coords = np.asarray(coords, float) n = len(coords) W = np.zeros((n, n)) for i in range(n): d = np.sqrt(((coords - coords[i]) ** 2).sum(axis=1)) idx = np.argsort(d)[1:k + 1] # 第 0 个是自己,跳过 W[i, idx] = 1.0 return W / W.sum(axis=1, keepdims=True) def morans_i(x, W, permutations=999, seed=42): x = np.asarray(x, float) z = x - x.mean() Ws = W.sum() I = (len(x) / Ws) * float(z @ W @ z) / float(z @ z) rng = np.random.default_rng(seed) sim = np.empty(permutations) for p in range(permutations): zp = rng.permutation(z) # 随机置换打破空间对应关系 sim[p] = (len(x) / Ws) * float(zp @ W @ zp) / float(zp @ zp) pval = (np.sum(np.abs(sim) >= abs(I)) + 1) / (permutations + 1) return {"I": I, "pseudo_p": float(pval), "sim_mean": sim.mean(), "sim_sd": sim.std()}seed固定是为了让 .doc 里的 p 值可复现,评审问起来能当场重跑。置换检验用双侧比较np.abs(sim) >= abs(I),因为 Moran's I 既可能显著为正(集聚)也可能显著为负(离散),单侧会漏掉一半情况。p 值分母加 1、分子加 1 是标准的无偏估计写法,避免出现 p = 0。
4.3 局部 LISA 与 999 次置换检验
全局 Moran's I 只给出一个总体判断,真要落到"哪些市是热点",得算局部 Moran's I,再逐单元做置换检验。
def local_morans(x, W, permutations=999, seed=42): x = np.asarray(x, float) z = x - x.mean() m2 = float((z ** 2).mean()) Ii = z * (W @ z) / m2 # 每个单元一个局部指数 rng = np.random.default_rng(seed) pvals = np.ones(len(x)) for i in range(len(x)): others = np.delete(z, i) cnt = 0 for _ in range(permutations): zi = rng.choice(others) # 从其余单元的取值里随机抽 if abs(zi * (W[i] @ z) / m2) >= abs(Ii[i]): cnt += 1 pvals[i] = (cnt + 1) / (permutations + 1) return Ii, pvals逐单元置换时从"其余单元的 z 值"里抽样,而不是对全部 z 打乱,是因为当前位置的取值参与了自身指数的计算,用全体打乱会低估显著性。999 次是个折中,跑 18 个单元大概几秒,p 值分辨率到 0.001 级别,够用。
4.4 结果怎么读:HH、LL、HL、LH 四类象限
把局部指数和它的空间滞后项画成散点图,四个象限对应四类单元。
| 象限 | 局部指数符号 | 邻居取值 | 含义 |
|---|---|---|---|
| HH | 正 | 高 | 高值被高值包围,核心增长区 |
| LL | 正 | 低 | 低值被低值包围,连片欠发达区 |
| HL | 负 | 低 | 高值孤岛,辐射未外溢 |
| LH | 负 | 高 | 低值塌陷,被高值包围 |
要强调的是:只有通过置换检验(比如 p < 0.05)的单元才值得在 .doc 里点名。不显著的 HL 单元很可能是随机波动,写进结论里经不起复核。另外 HH 和 LL 同时增多、HL 和 LH 减少,意味着空间极化在加深;反过来,HL 增多说明出现了新的增长极但还没带动周边。
5. 稳健性检验与 .doc 交付:从统计量到能站住的结论
5.1 四个方向的敏感性检验
一套差异测度跑出来的结论,最容易被问的就是"换个指标还成立吗"。我通常固定做四组对照:把人均 GDP 换成人均可支配收入重跑一遍基尼和泰尔;把全省平减指数换成按各市产业结构加权的平减指数;把济源示范区剔除后重算;把空间权重从 K 近邻换成邻接矩阵。四组里至少三组方向一致,结论才写进正文,不一致的放进附表并注明差异来源。
from itertools import product results = {} for metric, drop_jy, wmode in product(["gdp_pc_real", "income_pc_real"], [False, True], ["knn", "queen"]): sub = panel[panel["year"] == 2022] if drop_jy: sub = sub[sub["city_name"] != "济源示范区"] W = knn_weight(sub[["lon", "lat"]].values, k=5) if wmode == "knn" else queen_weight(sub) results[(metric, drop_jy, wmode)] = morans_i(sub[metric].values, W)["I"] pd.Series(results).unstack()四组对照跑完不到一分钟,换来的是方法部分一段可以写实的稳健性说明。
5.2 出图前先把中文字体配好
matplotlib 默认字体画不出中文,方框一出图就废。开头统一设置一次,全局生效。
import matplotlib matplotlib.rcParams["font.sans-serif"] = ["Noto Sans CJK SC", "SimHei", "Microsoft YaHei"] matplotlib.rcParams["axes.unicode_minus"] = False # 负号不显示成方块第二行经常被漏掉,导致 Moran 散点图里负的局部指数轴标签变成方框。字体列表按顺序回退,Linux 上装fonts-noto-cjk,Windows 上用SimHei兜底。
5.3 用 python-docx 把表格和结论写回 .doc
既然交付物是 .doc,与其手工复制粘贴,不如把汇总表直接写进 Word,改一次数据重跑一次。
from docx import Document doc = Document() doc.add_heading("河南省区域经济发展差异统计分析", level=0) doc.add_paragraph("数据来源:河南统计年鉴及各地市统计公报,GDP 按 2015 年不变价折算。") table = doc.add_table(rows=1, cols=4) table.style = "Table Grid" for cell, name in zip(table.rows[0].cells, ["年份", "基尼系数", "泰尔指数", "组间贡献率"]): cell.text = name for yr, row in summary.iterrows(): # summary 来自第 3 章的循环结果 cells = table.add_row().cells cells[0].text = str(yr) cells[1].text = f"{row['gini']:.4f}" cells[2].text = f"{row['theil']:.4f}" cells[3].text = f"{row['between_ratio']:.2%}" doc.save("河南省区域经济发展差异的统计分析.docx")Table Grid这个样式名不能写错,写成TableGrid会直接抛KeyError。数字格式化统一用四位小数,贡献率用百分号,避免文档里同一列出现不同精度的数。真要把 .doc 再转成老版二进制格式,用 LibreOffice 命令行批量转:soffice --headless --convert-to doc --outdir ./out ./*.docx,比在 Word 里手动另存为省事,也方便塞进 CI。
最后一个能省不少时间的技巧:把summary的列名和表头映射写成字典,指标口径一变只改字典不动循环体,重跑时整份 .doc 的表格、图注和口径说明会一起更新。
本文还有配套的精品资源,点击获取