简介《河南省区域经济发展差异的统计分析》是一份面向经济学、区域经济与统计类专业学生及研究者的分析类文档聚焦改革开放以来河南省内各区域发展不平衡问题围绕产业结构、基础设施、人力资源、教育投入和城镇化率等维度梳理差异现状并采用主成分分析法从科技投入、固定资产投资、教育投入、人口素质、城镇化率、外贸与就业等变量中提取关键影响因子衡量各因素的作用程度进而提出优化资源配置、强化科技创新、改善基础设施、提升教育水平、促进城镇化健康发展与扩大对外开放等对策建议适合课程论文写作、区域经济课题研究及统计指标方法参考。压缩包内为1个doc文档约590KB可直接查阅与引用。目前已有108人浏览学习可作为区域经济差异测度与主成分分析应用的参考范本。1. 河南省区域经济发展差异的统计分析一份 .doc 需求背后的数据管线需求方往往只丢过来一句话——河南省区域经济发展差异的统计分析交付物是一份 .doc。真正耗时间的不是写字而是把 18 个统计单元、十几年的面板数据对齐口径同一年的人均 GDP统计年鉴和统计公报能差出小几百元换一次常住人口来源基尼系数就跳一下。这篇内容讲的是一条能跑的链路先定统计单元和指标口径用 pandas 把年鉴宽表拗成面板长表再用 Python 复现极差、变异系数、基尼系数和泰尔指数的组间组内分解最后补一层空间自相关看差异是散着还是聚着。适合手里握着统计年鉴、需要把结论落到表格和图上的人也适合想拿一份真实面板练手空间计量的开发者。2. 统计单元与指标口径河南区域差异分析的数据地基2.1 18 个统计单元和县域两套口径怎么选河南省的统计口径里常见的区域分析单元有两层一层是 18 个统计单元即 17 个地级市加上济源示范区另一层是县域口径包含市辖区、县级市和县数量在百位级别。两套口径各有各的用处——市级口径适合做跨地市比较和空间自相关县域口径适合看市域内部是不是也在分化。我一般会先做市级面板跑通全流程再决定要不要下钻到县域。原因是县域口径的年鉴数据缺失更严重不少县的规上工业增加值、RD 经费支出会有整年空档插值补出来的数放进泰尔指数分解里组内差异会被稀释结论容易失真。提示济源示范区在部分年份的年鉴里单列在部分年份的省级汇总里并入其他口径。做面板前先确认这一列在目标年份区间内是否连续不连续就单独标注不要默默插值。行政区划调整是另一个坑。个别县市在观察期内发生过隶属关系或统计归属的变化表现为某个市的常住人口或 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_varscity_name, var_nameyear, value_namevalue) # 列名里抽出四位年份非年份列备注、数据来源会被抽成 NaN 后丢弃 long[year] pd.to_numeric( long[year].astype(str).str.extract(r(\d{4}))[0], errorscoerce ) long long.dropna(subset[year]).copy() long[year] long[year].astype(int) # 年鉴里的 —、… 先转成 NaN再统一处理 long[value] pd.to_numeric(long[value], errorscoerce) 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(dropTrue) assert panel[value].notna().mean() 0.95, 缺失率过高先回去核对年鉴页码逻辑上有三处值得说明。str.extract(r(\d{4}))解决的是列名不统一的问题年鉴里会出现2015年2015 年2015几种写法用正则一把抓比手工改列名稳。errorscoerce是把所有非数值内容统一变成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, onyear, howleft) # 现价 - 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(methodlinear, limit_areainside)) )limit_areainside这个参数很实用它只填中间的缺口序列开头和结尾的缺失保持原样。理由是首尾缺失往往意味着这个统计单元在观察期两端不在口径内插值会造出不存在的趋势。中间年份的单年缺失则多半是漏报或未发布线性插值的影响可控。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, wNone): 人口加权的区域基尼系数。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, wNone): 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. 空间视角Morans I 与 LISA 揭示的河南经济集聚4.1 空间权重矩阵邻接、距离衰减与 K 近邻怎么选常规的差异测度默认各单元彼此独立可地理上相邻的市往往同涨同落。空间权重矩阵就是把相邻这件事写成一个 n×n 的矩阵。三种常见构造一是邻接矩阵共用边界记为 1其余为 0符合行政边界相邻的直觉但河南西部山地单元邻居少东部平原单元邻居多度数差异大会带来偏误。二是距离衰减矩阵权重取1/d^2或exp(-d)能刻画远距离影响的衰减但要选一个截断距离主观性强。三是 K 近邻矩阵每个单元固定取最近的 K 个邻居度数一致适合单元规模悬殊的场景。我一般主用 K 近邻K 取 4 到 6把邻接矩阵的结果作为稳健性对照。两种矩阵都做行标准化让每行权重和为 1避免邻居多的单元把 Morans I 拉高。4.2 全局 Morans I 的手写实现依赖 libpysal 当然更省事但版本变动和geopandas的编译问题经常卡住环境。用经纬度坐标手写一遍二十行代码就够也更容易在 .doc 的方法部分讲清楚。import numpy as np def knn_weight(coords, k5): 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(axis1)) idx np.argsort(d)[1:k 1] # 第 0 个是自己跳过 W[i, idx] 1.0 return W / W.sum(axis1, keepdimsTrue) def morans_i(x, W, permutations999, seed42): 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)因为 Morans I 既可能显著为正集聚也可能显著为负离散单侧会漏掉一半情况。p 值分母加 1、分子加 1 是标准的无偏估计写法避免出现 p 0。4.3 局部 LISA 与 999 次置换检验全局 Morans I 只给出一个总体判断真要落到哪些市是热点得算局部 Morans I再逐单元做置换检验。def local_morans(x, W, permutations999, seed42): 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, k5) 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-cjkWindows 上用SimHei兜底。5.3 用 python-docx 把表格和结论写回 .doc既然交付物是 .doc与其手工复制粘贴不如把汇总表直接写进 Word改一次数据重跑一次。from docx import Document doc Document() doc.add_heading(河南省区域经济发展差异统计分析, level0) doc.add_paragraph(数据来源河南统计年鉴及各地市统计公报GDP 按 2015 年不变价折算。) table doc.add_table(rows1, cols4) 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 的表格、图注和口径说明会一起更新。本文还有配套的精品资源点击获取
