章节概述

本章聚焦人口普查数据分析与地理人口特征分类。专题研究 1 使用美国 ACS 数据,在全国县级尺度上构建一个基础版的地理人口特征分类,展示如何从人口与社会经济变量出发识别区域类型。专题研究 2 进一步转向纽约市 census tract 尺度,在更细颗粒度上整合人口、社会经济、住房与通勤等多维变量,识别城市内部街区之间的社会空间差异。

方法层面,本章主要涉及人口普查数据获取变量重构与比例指标计算相关性检验主成分分析(PCA)变量筛选K-meansH-K-means 聚类分析,以及聚类结果的空间制图指数热力图类型解释等关键方法。

本章学习内容

学习导航

本章学习内容

欢迎开始 第 7 章 的学习。请点击 下方导航卡 进入相应小节:

💡 提示:学习完一个小节后,请再次点击 屏幕右下角的章节主页按钮回到本导航页

1.0 本章导论

1.0 本章导论

1.1 人口普查数据概论

基本概念

人口普查数据(census data)是指由国家统计机构特定时点,对全国或特定行政范围内的人口、住户及其社会经济特征进行系统登记与统计后形成的数据。其核心目标在于刻画一个国家或地区人口的规模结构分布居住特征,并为公共治理、资源配置与社会研究提供基础性证据。

人口普查数据通常具有三个突出特征。第一,它具有较强的覆盖性,统计对象通常面向全国人口或全体住户,而不是抽样意义上的局部样本。第二,它具有明确的空间对应性,往往可以匹配到不同层级的行政区或统计区。第三,它具有较强的多维性,包含丰富的人口与社会属性变量,因此特别适合用于区域差异分析与社会空间研究。

人口普查中常见的统计主题包括但不限于年龄(age)、性别(sex)、受教育程度(educational attainment)、职业(occupation)、住房条件(housing)、民族或族裔(ethnicity)、迁移状况(migration)以及家庭结构(household structure)等。

从研究应用看,人口普查数据的价值并不只在于提供人口总量等宏观信息,更在于其能够被进一步细化为一组可比较的结构性指标,以刻画不同地区的社会经济与居住特征。

例如,可基于年龄结构(age structure)、家庭类型(household type)、住房占有方式(housing tenure)、就业状态(employment status)、受教育水平(educational attainment)、通勤方式(commuting mode)等变量,构建反映地区社会结构差异的指标体系。

这样的处理方式既是后续地理人口统计特征分类(geodemographics)的重要基础,也为下一章节的综合指数构建(composite index)提供了关键的数据来源与变量支撑。

Figure 1. 简化版人口普查流程示意图

Figure 1. 简化版人口普查流程示意图


国际对比

不同国家的人口普查制度与数据发布口径并不完全一致,但总体上都由官方统计机构负责组织、汇总与发布。对于人口社会空间分析而言,较常见的案例包括美国英国中国。三者在数据获取方式统计时点空间单元体系公开粒度(即,可公开的最小地理单元尺度)上均存在一定差异,因此在使用时需要特别留意年份口径空间边界变量定义的一致性。

  • 美国(United States)
    美国人口普查数据主要由 U.S. Census Bureau 发布。研究中最常用的两类数据,一类是十年一次的人口普查(Decennial Census),主要提供较完整的人口计数与基础人口特征;另一类则是当前更为常用的 American Community Survey(ACS)。ACS 并非传统意义上的整体现时点普查,而是一项持续开展的滚动抽样调查,并按年发布估计结果

    ACS 常见产品包括 1-year estimates5-year estimates。二者的差异不仅体现在可发布的空间尺度上,也体现在估计所依据的时间窗口与统计稳定性上。前者基于 12 个月的数据汇总而成,时效性更强,更适合观察较新的社会经济变化;后者基于 60 个月的数据汇总而成,样本量更大,统计稳定性通常更强,也因此能够支持更小尺度地理单元的估计。但需要注意的是,5-year estimates 反映的是一个五年时期内的综合特征,而不是某一单独年份的即时状态,因此 1-year estimates5-year estimates 不宜直接进行横向比较。

    美国人口普查常用的地理统计单元由细到粗通常包括 block groupcensus tractcountystate。在社会经济结构分析与地理人口统计研究中,较常使用的尺度通常为 block groupcensus tract。官方平台通常支持表格检索、下载、地图浏览与 API 访问。

    官方入口:U.S. Census Bureaudata.census.govACS Data

  • 英国(United Kingdom)
    英国人口普查具有较长的制度延续性。就英格兰和威尔士而言,人口普查自 1801 年起通常每 10 年开展一次,1941 年因二战未进行。当前英国人口普查数据的发布主要可通过 Office for National Statistics(ONS) 获取;不过需要注意,英国并非所有地区都由同一统计机构负责,其中英格兰与威尔士主要由 ONS 统筹,北爱尔兰由 NISRA(Northern Ireland Statistics and Research Agency)负责,苏格兰则由 Scotland’s Census / National Records of Scotland负责。

    英国人口普查常用的地理统计单元由细到粗通常包括 Output Area(OA)Lower layer Super Output Area(LSOA)Middle layer Super Output Area(MSOA)local authority。这一分层体系使人口普查数据能够较方便地与边界文件结合,并支持不同尺度的社会空间分析。

    官方入口:ONS CensusCensus 2021 geographiesStatistical geographies

Figure 2. 英国2021人口普查问题节选图

Figure 2. 英国2021人口普查问题节选图

  • 中国(China)
    中国人口普查数据主要由国家统计局发布。新中国成立以来,我国已分别于 1953 年1964 年1982 年1990 年2000 年2010 年2020 年开展了七次全国人口普查;现行制度下,人口普查通常每 10 年开展一次。当前研究中最常引用的是第七次全国人口普查相关公报、主要数据资料与《中国人口普查年鉴-2020》。

    第七次全国人口普查的标准时点为 2020 年 11 月 1 日零时。其主要内容涵盖人口数量、结构、分布、城乡住房以及性别、年龄、民族、受教育程度、行业、职业、迁移流动、婚姻生育等信息。与前文美国的 ACS 或英国的小区域普查表不同,中国人口普查更常见的公开资料形式是公报、年鉴与汇总表。其中,国家统计局发布的《中国人口普查年鉴-2020》主要包括全部人口数据长表数据附录等部分,其中长表数据更多反映人口的各类结构性特征。

    从公开结果的组织层级来看,中国人口普查数据通常由细到粗涉及乡镇/街道县级地级省级全国等行政层级;但在公开可获得性上,研究中更常稳定使用的是县级及以上尺度的数据。与美国和英国相比,中国公开可直接获取的精细化小区域人口普查数据相对有限,尤其是在县级以下统一、标准化、可直接下载的空间单元层面,可得性通常较弱。

    官方入口:国家统计局·第七次全国人口普查主要数据中国人口普查年鉴-2020

Figure 3. 中国第六次人口普查短表示意图

Figure 3. 中国第六次人口普查短表示意图


常见的人口普查变量

人口普查数据的内容通常具有明显的层级组织特征。概括而言,可将其理解为从概念层领域层再到变量层的逐级展开。所谓概念(concepts),是指较高层次的社会人口主题,如人口结构(demographic structure)社会经济特征(socioeconomic characteristics)住房与居住状况(housing and living conditions);在每一个概念之下,又可以进一步细分为若干领域(domains),例如年龄、性别、教育、就业或住房占有方式;而每一个领域内部,还会包含更具体的变量(variables),如 0–4 岁、5–9 岁、10–14 岁等年龄分组,或自有住房、社会租赁、私人租赁等住房类别。

换言之,人口普查并不是简单提供若干零散变量,而是围绕若干主题,形成从宏观主题具体分类项的系统化统计框架。这种层级式组织方式,一方面有助于从整体上把握地区人口社会特征,另一方面也为后续的变量筛选、指标构建与聚类分析提供了基础。

注意: 这里的 concept 可以同义替换为 theme / topic 等
以英国 Census 2021 为例,ONS 将 England and Wales 的 topic summaries 组织为若干主题组。官方页面列出 9 个主题组

从分析实践看,不同国家的人口普查虽然在发布平台、主题名称与表格结构上存在差异,但其核心内容通常都围绕若干相对稳定的领域展开。例如,人口学领域常涉及年龄、性别、婚姻与家庭结构;社会经济领域常涉及教育、职业、就业状态与行业;住房领域常涉及住房类型、住房占有方式与拥挤程度;迁移与流动领域则可能涉及出生地、迁居、国籍或通勤方式等。

Figure 4. 人口普查数据层级概念

Figure 4. 人口普查数据层级概念

方法提示:从计数到可比较指标

人口普查原始表格通常以计数值(counts)为主,而不会直接提供适合聚类分析的比例指标(percentages)密度指标(densities)标准化指标(standardised measures)

因此,在后续分析中,研究者通常需要根据研究目标,对原始计数进一步处理。例如,可将某一类别人数除以总人口,得到结构比例;将人口数除以面积,得到人口密度;或在多变量比较中进一步进行标准化处理


1.2 地理人口特征(Geodemographics)

基本概念

“an analysis of people by where they live” - (Sleight, 1996)

地理人口特征(geodemographics)地理人口特征分类法(geodemographic classification)通常是指:基于小区域的人口、社会、经济与居住特征,将具有相似属性的地理单元归并为若干类型,从而形成一种邻里类型学(neighbourhood typology)区域分类体系(area classification)。它并不是对单个个体进行分类,而是对地区进行分类;其基本出发点在于,不同地区的人口社会结构、居住条件与生活方式往往存在相对稳定的空间差异。

底层逻辑:‘物以类聚,人以群分’ - 居住在同一邻里的人群往往具有类似的多维属性

地理人口特征分类 通常以小区域数据为基础,围绕地区的多维城市背景特征,选取人口、社会、经济、住房等多个变量,对空间单元进行聚类分析,并在结果基础上进一步开展类别解释与命名。其最终产出通常不是连续型数值,而是若干离散型组别。这些组别往往具有较强的解释性,其命名通常并非随意设定,而是依据该类地区在多维属性中最突出的结构特征加以概括,例如 “高龄社区”“多族裔城市租住区”“郊区家庭型社区” 等。

这类分类结果的意义,在于将原本维度较多、较分散的普查变量压缩为更易解释、也更适合比较的区域类型,从而为城市邻里差异社会空间分异区域类型识别提供一种多维度的分析视角。

注:小区域 往往是指街道、邻里等地理单元
(常见如美国的census tract及以下、英国的 LSOA及以下单元)

Figure 5. Esri Tapestry 分类 1B 人群肖像

Figure 5. Esri Tapestry 分类 1B 人群肖像


多维度方法

在社会空间研究中,单一变量方法往往只能揭示地区特征的某一个侧面。例如,仅用老年人口占比,只能描述年龄结构;仅用住房占有方式(housing tenure),也只能反映居住制度的一部分。即便在多元模型中,城市背景变量通常也仍然被拆解为若干相对独立的解释项,分别考察其边际影响。这样的做法有其明确用途,但对于理解城市邻里作为整体情境的内部结构而言,往往仍然是不充分的。

地理人口特征分类 所代表的,正是一种更偏向情境化(contextual)的分析思路。在这一思路下,邻里并不被视为若干彼此孤立变量的机械组合,而更像是一组相互交织、难以完全拆分的属性束。例如,财富、教育程度、住房成本、家庭结构与就业特征在现实城市空间中往往高度缠绕,其意义并不总能通过单独考察某一个变量而被充分理解。相较之下,地理人口特征分类 更关注的是地区在属性组合上的整体差异,而不是某一个单独变量上的增减变化。

这一多维方法的优势主要体现在两个方面。第一,它能够将人口结构社会经济特征住房状况就业特征迁移流动等多个维度结合起来,从整体上识别不同地区之间的相似性与差异性,因此所得结果更接近于一种综合社会画像,而不是单指标的局部刻画。第二,它能够将高维变量压缩为较易解释的类别结构,从而更便于制图展示、区域比较与结果沟通。

相关研究还指出(Singleton & Spielman, 2004),当小区域人口与社会经济变量本身存在较大估计误差(margin of error)时,情境化方法相较于逐变量分析可能具有额外优势。这一点在美国邻里尺度的 ACS 数据中尤为典型:当单一变量的测量并不十分精确时,将邻里理解为多属性组合并进行整体分类,有时比逐一解释单个变量更稳健。


发展脉络

从历史渊源看,geodemographics 的早期雏形可追溯到 19 世纪末的城市社会调查与社会制图实践。其中最具代表性的案例之一,是 Charles Booth 对伦敦贫困问题的系统调查及其著名的 London Poverty Maps。在 Booth 于 1886–1903 年开展的 Inquiry into Life and Labour in London 中,研究团队以街道为单位,对居民的收入与社会阶层进行分类,并通过不同颜色在地图上展示伦敦内部的贫困与富裕分布。这一工作虽然尚不属于现代意义上的 地理人口特征分类,但其核心思想已经体现出一个重要特征:即通过对小区域人口社会特征的归纳与制图,识别城市内部的社会空间差异。

Figure 6. Charles Booth 与 伦敦贫困地图

Figure 6. Charles Booth 与 伦敦贫困地图

不过,从更严格的学术脉络看,地理人口特征分类的正式理论起源通常并不追溯到 Booth 本人,而更多被认为来自 20 世纪 20–30 年代芝加哥学派(Chicago School)及其城市生态学(urban ecology / human ecology)传统。相关综述普遍指出,现代 geodemographics 的理论基础与芝加哥学派的人类生态学、其后发展的 social area analysisfactorial ecology 有更直接的关联。

芝加哥学派的核心关切,在于理解城市内部不同群体如何在空间上分化、聚集与竞争,并将城市视为一种具有相对秩序的社会生态系统。其代表性思想包括:不同社会群体会在城市内部形成具有差异性的居住区位,而这些差异并不是随机分布,而是与社会阶层、族裔结构、家庭生命周期及经济功能等因素密切相关。正是在这一意义上,城市生态学为后来的 geodemographic classification 提供了更明确的理论框架:即通过小区域社会人口特征的系统比较,识别城市内部可区分的邻里类型与社会空间结构。

地理人口特征学/分类 并不是脱离学术传统的独立技术,而是建立在长期的邻里分类与城市社会分异研究之上的一种现代化延伸。

现代意义上的 geodemographic systems 一般被认为是在 20 世纪 70 年代于英国和美国逐步确立的。其中,英国学者 Richard Webber 的工作具有关键奠基意义。Webber 在 1970 年代中后期推动了这一方法从社会区域分析走向更系统化、可操作化的邻里分类框架。此后,地理人口特征分类 很快进入市场细分消费者分析商业选址定向营销等商业场景,并在商业数据库与地理信息系统支持下迅速发展。 在随后的发展过程中,此方法逐渐从商业部门扩展到公共部门社会政策官方统计领域。

Figure 7. Singleton & Spielman 2014 英美两国Geodemographics简要发展历史梳理图

Figure 7. Singleton & Spielman 2014 英美两国Geodemographics简要发展历史梳理图

英国人口普查与公共地理人口特征分类

英国是 地理人口特征分类 发展最典型的国家之一。
英国国家统计局 ONS 及其前身机构长期基于人口普查开展小区域 area classifications。相关文献与项目资料显示,英国从 1971 年普查后的社会区域分类起,随后在 1981199120012011 以及 2021 等普查轮次后,均持续推进类似的区域分类工作。其中,近二十年来最常见的官方命名为 Output Area Classification(OAC);例如 2001 OAC2011 OAC,以及最新的 2021 residential-based area classifications


典型应用场景

地理人口特征分类 的应用领域较为广泛,既可用于商业市场分析,也可服务于公共政策、城市研究和社会科学建模。常见应用场景包括:

  • 商业与市场分析:用于消费者分类目标人群识别客户画像选址分析精准营销零售市场细分等。

  • 城市规划与社区研究:用于识别不同类型的城市邻里,分析居住分异邻里变化城市更新通勤结构职住关系等问题。

  • 公共政策与社会治理:用于支持资源配置服务对象识别社区分层干预空间不平等监测政策评估等。

  • 公共健康与健康地理:用于分析健康不平等疾病风险差异医疗服务可达性健康行为模式社区健康干预等。

  • 教育与社会机会分析:用于研究教育资源分布升学机会差异高等教育参与学生来源地特征等。

  • 住房、交通与基础设施:用于分析住房市场分异房地产需求结构出行需求差异交通服务配置公共服务可达性等。

  • 环境与可持续发展:用于研究环境暴露差异空气污染不平等绿地可达性能源消费模式环境公平等。

  • 官方统计与区域分类:用于构建标准化的小区域分类体系,例如英国国家统计局 ONS 推出的 Output Area Classification(OAC) 及其后续更新版本。

  • 社会科学建模与解释:可作为描述性分类变量社会经济地位代理变量回归模型控制变量多层模型分组变量小区域剥夺测度,用于解释不同社区背景下的社会行为、健康结果和空间不平等。

Figure 8. 地理人口特征分类分析潜在应用场景

Figure 8. 地理人口特征分类分析潜在应用场景

2.0 地理人口特征分类实战 1

2.0 专题研究 1: 全美县级分类

专题研究 1: 美国各县地理人口特征分类

专题导航

美国全国范围县级地理人口特征分类

“从人口普查数据到区域类型识别”

本节以美国 ACS 数据为例,介绍地理人口特征分类的基本流程,包括变量选择与获取数据重构与预处理K-means 聚类以及结果解释

点击下方卡片,快速跳转至相应模块:

💡 提示:学习完一个小节后,请再次点击 屏幕右下角的专题导航按钮回到本导航页

研究课题
2.1 研究课题与分析目标

实战教学

本节以美国人口普查局 (U.S. Census Bureau) 最新发布的 2024年 American Community Survey(ACS) 5-year estimates 为数据基础,结合 tidycensus 包,开展一个面向全美国县级单元的人口与社会经济地理人口特征分类实战。

研究重点是在县域尺度(county)上整合多个与人口结构教育收入就业等相关的指标,识别不同县之间在社会人口背景上的综合差异,并进一步构建具有解释性的地理人口特征分类(geodemographic classification)

与十年一次的人口普查相比,ACS 是一项持续开展的抽样调查,能够提供更为及时的社会经济统计信息;其中 5-year estimates 由于覆盖范围更广、空间粒度更细,因此特别适合用于地理人口特征分析。

Figure 9. ACS官方网站截图

Figure 9. ACS官方网站截图

第四章在“纽约市贫困率空间可视化”部分,已经利用 tidycensus 包演示过如何调用美国人口普查局发布的 ACS 5-year 数据。本节将在此前基础上,进一步将 tidycensus 用于全国尺度的人口普查数据整理,即从单一指标的空间可视化,转向面向地理人口特征分类的多变量数据构建。

复习tidycensus 是 R 环境中连接美国 Census API 的常用工具,能够将 Census 与 ACS 数据整理为适合 tidyverse 工作流的数据格式,并在需要时同步返回空间边界数据

数据来源说明
美国人口普查局 ACS 5-year Data:
https://www.census.gov/data/developers/data-sets/acs-5year.html
tidycensus 官方文档:
https://walker-data.com/tidycensus/
get_acs() 函数说明:
https://walker-data.com/tidycensus/reference/get_acs.html


示例课题

本节设置一个面向流程演示的示例研究:基于 tidycensus 获取美国 ACS 5-year estimates 数据,在全美国县级尺度上,围绕人口特征社会经济特征构建一个基础版的地理人口特征分类

作为教学示例,本节并不追求纳入尽可能全面的城市背景变量,而是有意识地选取一组较为常见、也较易理解的人口与社会经济指标,以展示地理人口特征分析的基本流程。住房、健康、迁移、行为等变量同样是相关研究中的常见维度,但本节暂不展开。

在空间单元上,本节选取县级(county)作为分析对象,主要出于教学实现与全国尺度演示的考虑。县级单元便于在全美范围内统一获取、整理与比较,也适合作为地理人口特征分类的入门示例。

在更完整的研究中,研究者往往会进一步使用更细颗粒度的小区域单元;因此,本节更适合被理解为一个方法流程示范,而不是完整意义上的精细化邻里分类实践。

本节的分析目标主要包括以下几个方面:

  • 学习如何使用 tidycensus 提取 ACS 5-year estimates 中与人口社会经济相关的变量,建立基础数据获取流程;
  • 完成原始普查数据的基本预处理,包括变量整理指标转换与分析数据框构建;
  • 按照常见地理人口特征分类的技术路线,逐步演示数据获取数据预处理聚类分析结果解释空间制图类别解读等核心环节;
  • 在全国尺度上,对美国本土县域进行基础聚类,并对所得类型进行初步解读。
数据选择与获取
2.2 变量选择与数据获取

变量选择

在正式提取数据之前,需要先明确本节示例研究所使用的变量框架。

考虑到本节主要服务于流程演示,而非构建一个高度严谨、完整的地理人口特征分类体系,因此这里暂时仅选择两个主题:一是人口结构,二是社会经济属性。在每个主题之下,再进一步细分为若干领域,并据此选取后续聚类分析所需的指标。

ACS 原始表格通常以计数值(count)为主,因此后续用于聚类的多数指标并不会直接以“百分比”形式下载,而是需要根据总量变量与分项变量进一步计算得到。例如,年龄占比、族裔占比、失业率或本科及以上比例,通常都需要在提取表格后再自行重构。tidycensus 提供了 get_acs() 用于提取数据,也可结合 load_variables() 先检索表号与变量编码。

本节示例研究的变量设计与 ACS 变量选取
主题 领域 变量 ACS变量代码
人口结构 年龄 0–17岁人口占比 (B01001_003 + B01001_004 + B01001_005 + B01001_006 + B01001_027 + B01001_028 + B01001_029 + B01001_030) / B01001_001
18–34岁人口占比 (B01001_007 + B01001_008 + B01001_009 + B01001_010 + B01001_011 + B01001_012 + B01001_031 + B01001_032 + B01001_033 + B01001_034 + B01001_035 + B01001_036) / B01001_001
35–64岁人口占比 (B01001_013 + B01001_014 + B01001_015 + B01001_016 + B01001_017 + B01001_018 + B01001_019 + B01001_037 + B01001_038 + B01001_039 + B01001_040 + B01001_041 + B01001_042 + B01001_043) / B01001_001
65岁及以上人口占比 (B01001_020 + B01001_021 + B01001_022 + B01001_023 + B01001_024 + B01001_025 + B01001_044 + B01001_045 + B01001_046 + B01001_047 + B01001_048 + B01001_049) / B01001_001
性别 女性人口占比 B01001_026 / B01001_001
族裔 白人人口占比 B02001_002 / B02001_001
黑人人口占比 B02001_003 / B02001_001
亚裔人口占比 B02001_005 / B02001_001
其他族裔人口占比 (B02001_004 + B02001_006 + B02001_007 + B02001_008) / B02001_001
社会经济属性 收入 家庭收入中位数 B19013_001
就业 失业率 B23025_005 / B23025_003
教育 本科及以上人口占比 (B15003_022 + B15003_023 + B15003_024 + B15003_025) / B15003_001
注:表中分母变量通常表示该 ACS 表对应的总体人口或基准人口。比例型变量需在下载原始 ACS 分项变量后进行汇总计算。

若在后续研究中希望进一步提高分类的丰富度,也可以继续纳入住房、健康、迁移、家庭结构或通勤等维度。


变量选择说明

本节将使用美国 ACS 2020–2024 5-year estimates 作为数据来源。

tidycensus 在实际操作中通常通过 get_acs() 提取表格数据,通过 load_variables() 查询表号与字段定义,因此本节后续也将沿用这一工作流。

本节主要涉及以下几张 ACS 常用详细表:

  • B01001
    Sex by Age,用于构造不同年龄段比例与性别结构;
  • B02001
    Race,用于构造主要族裔结构变量;
  • B19013
    Median Household Income,用于表示收入水平;
  • B23025
    Employment Status,用于构造失业率;
  • B15003
    Educational Attainment,用于构造本科及以上人口比例。

ACS 变量代码如何理解?

ACS 变量代码通常由表号变量后缀组成。表号表示变量所属的主题表,变量后缀则对应表格中的具体分类项。

B01001 为例,该表的主题是 Sex by Age。其中:

变量代码 变量含义
B01001_003 Male: Under 5 years
B01001_004 Male: 5 to 9 years
B01001_005 Male: 10 to 14 years
B01001_006 Male: 15 to 17 years

因此,若需要构造“0–17岁男性人口”,就可以将上述变量加总。类似地,构造比例型指标时,还需要选择合适的分母变量,例如 B01001_001 表示该表中的总人口。

(推荐查询工具:Census Reporter 适合快速浏览 ACS 表格结构;tidycensus 官方文档与变量搜索工具适合在 R 中进一步核对变量代码。)

实用查询工具

查询 ACS 变量时,可以优先使用两个工具:

  • Census Reporter:适合快速浏览 ACS 表格、变量层级和主题分类。网址:censusreporter.org
  • tidycensus 文档与变量搜索:适合在 R 中结合 load_variables() 查询和核对变量代码。网址:walker-data.com/tidycensus/articles/basic-usage.html

(建议先用关键词找到候选表号,再回到变量字典或官方说明中核对变量含义、分母口径和适用年份。)


重要

如果不知道 ACS 表号怎么办?

在使用 ACS 数据时,研究者往往先知道自己需要的概念变量,例如房屋类型教育水平通勤方式就业状况,但并不一定一开始就知道它们对应的 ACS 表号与变量代码。遇到这种情况时,可以按照以下顺序进行检索:

1. 优先在 R 中查询 ACS 变量字典

最推荐的方式是使用 tidycensus::load_variables() 读取 ACS 变量字典,然后通过关键词搜索相关变量。这样可以直接看到变量代码、变量名称和所属表格主题,便于后续下载与计算。

library(tidycensus)
library(dplyr)
library(stringr)

acs_vars <- load_variables(
  year = 2022,
  dataset = "acs5",
  cache = TRUE
)

acs_vars %>%
  filter(str_detect(str_to_lower(label), "units in structure")) %>%
  select(name, label, concept)

例如,如果想查找房屋类型,可以尝试搜索:

acs_vars %>%
  filter(
    str_detect(str_to_lower(label), "detached|attached|apartment|units") |
    str_detect(str_to_lower(concept), "units in structure")
  ) %>%
  select(name, label, concept)

这类检索通常可以帮助定位到 B25024,其表格主题为 Units in Structure,可用于构造独立式住宅、相连住宅和多户住宅等指标。


2. 使用表格主题关键词进行模糊搜索

如果不知道具体英文表达,可以先从较宽泛的主题词开始搜索。例如:

  • 房屋类型housingstructureunitsdetachedattached

  • 教育水平educationeducational attainmentbachelor

  • 就业状况employmentlabor forceunemployed

  • 职业结构occupationmanagementservice occupations

  • 通勤方式transportation to workmeans of transportation

  • 通勤时间travel time to work

  • 收入median household income

  • 贫困poverty status

acs_vars %>%
  filter(str_detect(str_to_lower(concept), "travel time to work")) %>%
  select(name, label, concept)

3. 查阅 tidycensus 官方变量浏览器或 Census 官方表格目录

如果 R 中搜索结果较多,可以进一步使用在线变量浏览器或官方 ACS 表格目录核对。常见做法是先在 R 中获得候选表号,例如 B25024B08301B15003,再回到官方说明中确认该表的主题、统计口径和变量定义。

(这一步尤其适合核对变量分母。例如教育变量通常以 25 岁及以上人口为基准,婚姻变量通常以 15 岁及以上人口为基准,就业变量通常以 16 岁及以上人口为基准。)


4. 使用 AI 辅助初步定位,再回到变量字典核实

推荐

在无法快速确定表号时,也可以先向 AI 描述研究需求,让其给出候选 ACS 表号。需要注意的是,AI 给出的结果应作为初步线索,最终仍需通过 load_variables() 或官方文档核对。

可以这样提问:

(我想从 ACS 5-year 数据中构造“房屋类型”变量,例如独立式住宅、相连住宅和多户住宅占比。请问可能对应哪张 ACS 表?具体变量代码和分母变量是什么?)

或者:

(我想从 ACS 数据中计算“本科及以上人口占比”。请问应使用哪张 ACS 表?分子变量和分母变量分别是什么?该变量的统计基准人口是什么?)

更推荐的提问方式是同时说明:

  • 研究尺度,例如 census tractcountyblock group
  • 数据类型,例如 ACS 5-year
  • 目标年份,例如 20222023
  • 需要构造的指标,例如比例、中位数或分类变量
  • 希望返回的内容,例如表号、变量代码、变量名称、分子变量和分母变量

数据获取 实战代码

环境配置部分,所需要的包,以及 acs年份等基础设置

# ==============================================================================
# 1. 环境构建与依赖项声明
# ==============================================================================
library(tidycensus)   # 获取美国人口普查 / ACS 数据
library(tidyverse)    # 数据整理与可视化
library(sf)           # 处理空间矢量数据
library(classInt)    # Fisher 断点
library(viridis)     # viridis 配色
library(ggspatial)   # 比例尺、指北针、经纬线
library(scales)      # 标签格式化
library(patchwork)   # 标签格式化

# ==============================================================================
# 2. Census API 认证设置
# ==============================================================================
# 注:
# - 首次使用 tidycensus 前,需先申请美国 Census API key
# - 申请地址:https://api.census.gov/data/key_signup.html
# - 首次配置时取消下行注释,并替换为自己的 key;运行一次后即可写入本地环境
# census_api_key("YOUR_CENSUS_API_KEY", install = TRUE, overwrite = TRUE)

# 若已完成 key 安装,建议重启 R 会话后再继续运行

# ==============================================================================
# 3. 参数设定:年份、调查类型与空间尺度
# ==============================================================================
acs_year   <- 2024       # 对应 ACS 2020–2024 5-year estimates 的发布年
acs_survey <- "acs5"     # 5-year estimates
geo_unit   <- "county"   # 县级单元(county)
# ==============================================================================
# 4. ACS 变量字典检索
# ==============================================================================
# 目的:
# - 预先读取 ACS 变量字典,便于核对变量名称与字段含义
# - 本步不直接参与分析,但有助于后续检查表号与代码是否正确

acs_vars_2024 <- load_variables(
  year = acs_year,
  dataset = acs_survey,
  cache = TRUE
)

# 快速查看本节涉及的核心表号
acs_vars_2024 %>%
  filter(str_detect(name, "^B01001|^B02001|^B19013|^B23025|^B15003")) %>%
  select(name, label, concept) %>%
  print(n = 5) # 在script里运行的时候可以输入 View()来看这个字典
# A tibble: 380 × 3
  name        label                                   concept                 
  <chr>       <chr>                                   <chr>                   
1 B01001A_001 Estimate!!Total:                        Sex by Age (White Alone)
2 B01001A_002 Estimate!!Total:!!Male:                 Sex by Age (White Alone)
3 B01001A_003 Estimate!!Total:!!Male:!!Under 5 years  Sex by Age (White Alone)
4 B01001A_004 Estimate!!Total:!!Male:!!5 to 9 years   Sex by Age (White Alone)
5 B01001A_005 Estimate!!Total:!!Male:!!10 to 14 years Sex by Age (White Alone)
# ℹ 375 more rows

数据提取

# ==============================================================================
# 5. 研究所需变量清单
# ==============================================================================
# 说明:
# - 本节先完成“原始变量提取”,暂不进行比例化、标准化与聚类
# - 因此,这里提取的是后续构造指标所需的底层变量
# - 其中:
#   * B01001  用于年龄结构与女性人口占比
#   * B02001  用于族裔结构
#   * B19013  用于家庭收入中位数
#   * B23025  用于失业率
#   * B15003  用于本科及以上人口占比

acs_var_dict <- c(
  # ------------------------------------------------------------
  # A. 年龄与性别(B01001: Sex by Age)
  # ------------------------------------------------------------
  total_pop      = "B01001_001",
  male_total     = "B01001_002",
  female_total   = "B01001_026",

  # 0–17 岁
  male_0_17_1    = "B01001_003",
  male_0_17_2    = "B01001_004",
  male_0_17_3    = "B01001_005",
  male_0_17_4    = "B01001_006",
  female_0_17_1  = "B01001_027",
  female_0_17_2  = "B01001_028",
  female_0_17_3  = "B01001_029",
  female_0_17_4  = "B01001_030",

  # 18–34 岁
  male_18_34_1   = "B01001_007",
  male_18_34_2   = "B01001_008",
  male_18_34_3   = "B01001_009",
  male_18_34_4   = "B01001_010",
  male_18_34_5   = "B01001_011",
  male_18_34_6   = "B01001_012",
  female_18_34_1 = "B01001_031",
  female_18_34_2 = "B01001_032",
  female_18_34_3 = "B01001_033",
  female_18_34_4 = "B01001_034",
  female_18_34_5 = "B01001_035",
  female_18_34_6 = "B01001_036",

  # 35–64 岁
  male_35_64_1   = "B01001_013",
  male_35_64_2   = "B01001_014",
  male_35_64_3   = "B01001_015",
  male_35_64_4   = "B01001_016",
  male_35_64_5   = "B01001_017",
  male_35_64_6   = "B01001_018",
  male_35_64_7   = "B01001_019",
  female_35_64_1 = "B01001_037",
  female_35_64_2 = "B01001_038",
  female_35_64_3 = "B01001_039",
  female_35_64_4 = "B01001_040",
  female_35_64_5 = "B01001_041",
  female_35_64_6 = "B01001_042",
  female_35_64_7 = "B01001_043",

  # 65 岁及以上
  male_65p_1     = "B01001_020",
  male_65p_2     = "B01001_021",
  male_65p_3     = "B01001_022",
  male_65p_4     = "B01001_023",
  male_65p_5     = "B01001_024",
  male_65p_6     = "B01001_025",
  female_65p_1   = "B01001_044",
  female_65p_2   = "B01001_045",
  female_65p_3   = "B01001_046",
  female_65p_4   = "B01001_047",
  female_65p_5   = "B01001_048",
  female_65p_6   = "B01001_049",

  # ------------------------------------------------------------
  # B. 族裔结构(B02001: Race)
  # ------------------------------------------------------------
  race_total     = "B02001_001",
  white_alone    = "B02001_002",
  black_alone    = "B02001_003",
  asian_alone    = "B02001_005",
  other_race_1   = "B02001_004",
  other_race_2   = "B02001_006",
  other_race_3   = "B02001_007",
  other_race_4   = "B02001_008",

  # ------------------------------------------------------------
  # C. 收入(B19013: Median Household Income)
  # ------------------------------------------------------------
  med_hh_income  = "B19013_001",

  # ------------------------------------------------------------
  # D. 就业(B23025: Employment Status)
  # ------------------------------------------------------------
  labor_force    = "B23025_003",
  unemployed     = "B23025_005",

  # ------------------------------------------------------------
  # E. 教育(B15003: Educational Attainment for Population 25 Years and Over)
  # ------------------------------------------------------------
  edu_total_25p  = "B15003_001",
  bachelor       = "B15003_022",
  master         = "B15003_023",
  professional   = "B15003_024",
  doctorate      = "B15003_025"
)

# 查看变量数量
length(acs_var_dict)
[1] 65

下面这一部分可能会因为网速波动中断卡顿, 可以多试几次

# ==============================================================================
# 6. 全国县级 ACS 原始数据提取(附带县级边界)
# ==============================================================================
# 说明:
# - 若本地已存在缓存文件,则直接读取
# - 若不存在,则调用 Census API 下载,并保存至 data/ 文件夹
# - geometry = TRUE:同步下载 county 边界,便于后续空间制图
# - output = "wide":每个变量展开为单独列,便于后续重构

dir.create("data", showWarnings = FALSE)

raw_file <- "data/acs_county_raw_2024.rds"

if (file.exists(raw_file)) {
  message("检测到本地缓存文件,直接读取原始 ACS 县级数据...")
  acs_county_raw <- readRDS(raw_file)
} else {
  message("本地无缓存文件,开始从 Census API 下载原始 ACS 县级数据...")
  
  acs_county_raw <- get_acs(
    geography = geo_unit,
    variables = acs_var_dict,
    year = acs_year,
    survey = acs_survey,
    geometry = TRUE,
    output = "wide"
  )
  
  saveRDS(acs_county_raw, raw_file)
  message("原始 ACS 县级数据已保存至:", raw_file)
}
Downloading: 1.8 kB     Downloading: 1.8 kB     Downloading: 1.8 kB     Downloading: 1.8 kB     
# ==============================================================================
# 7. 仅保留美国本土县域(CONUS)
# ==============================================================================
# 说明:
# - 去除阿拉斯加(02)、夏威夷(15)及海外属地(60, 66, 69, 72, 78)
# - 本节后续分析仅保留美国本土县级单元

conus_file <- "data/acs_county_conus_raw_2024.rds"

if (file.exists(conus_file)) {
  message("检测到本地缓存文件,直接读取美国本土县域数据...")
  acs_county_conus <- readRDS(conus_file)
} else {
  acs_county_conus <- acs_county_raw %>%
    filter(!str_sub(GEOID, 1, 2) %in% c("02", "15", "60", "66", "69", "72", "78"))
  
  saveRDS(acs_county_conus, conus_file)
  message("美国本土县域数据已保存至:", conus_file)
}

# ==============================================================================
# 8. 原始结果结构检验
# ==============================================================================
# print(acs_county_conus)

# glimpse(acs_county_conus)  # 因变量较多,此处不展开打印

数据展示制图

使用第四章空间制图方法,简单查看一下获取的数据(空间以及非空间数据)

# ==============================================================================
# 9. 美国本土县级家庭收入中位数可视化
# ==============================================================================

# ------------------------------------------------------------------------------
# 9.1 数据准备:提取家庭收入中位数并计算 Fisher 断点
# ------------------------------------------------------------------------------
map_income <- acs_county_conus %>%
  select(GEOID, NAME, med_hh_incomeE, geometry) %>%
  filter(!is.na(med_hh_incomeE))

# Fisher 分级数:这里设为 6 组
n_class <- 6

fisher_brks <- classIntervals(
  var = map_income$med_hh_incomeE,
  n = n_class,
  style = "fisher"
)

# 构造更易读的图例标签:统一换算为千美元(k USD)
brks <- fisher_brks$brks
income_labels <- paste0(
  "$", round(brks[-length(brks)] / 1000, 0), "k",
  " – ",
  "$", round(brks[-1] / 1000, 0), "k"
)

# 将连续变量切分为分级因子
# 注:此处保持“低值 → 高值”的自然顺序,不再反转因子
map_income <- map_income %>%
  mutate(
    income_class = cut(
      med_hh_incomeE,
      breaks = brks,
      include.lowest = TRUE,
      labels = income_labels
    )
  )

# ------------------------------------------------------------------------------
# 9.2 构建经纬线对象
# ------------------------------------------------------------------------------
graticule_sf <- st_graticule(map_income)

# ------------------------------------------------------------------------------
# 9.3 可视化:县级家庭收入中位数设色图
# ------------------------------------------------------------------------------
ggplot() +
  # 经纬线
  geom_sf(
    data = graticule_sf,
    color = "grey70",
    linewidth = 0.3,
    linetype = "dashed"
  ) +
  # 县级设色图层
  geom_sf(
    data = map_income,
    aes(fill = income_class),
    color = NA
  ) +
  # 比例尺
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    text_cex = 0.8
  ) +
  # 指北针
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    style = north_arrow_fancy_orienteering,
    height = unit(0.9, "cm"),
    width  = unit(0.9, "cm")
  ) +
  # 配色与图例
  scale_fill_viridis_d(
    option = "mako",
    direction = -1,
    name = "Median household income\n(1,000 USD)",
    guide = guide_legend(reverse = TRUE)
  ) +
  coord_sf(crs = st_crs(4326)) +
  labs(
    title = "Median Household Income across CONUS Counties",
    subtitle = "ACS 2020–2024 5-year estimates",
    caption = "Source: U.S. Census Bureau ACS via tidycensus"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 11),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 9),
    axis.title = element_blank()
  )

数据重构与预处理
2.3 数据重构与预处理

数据重构

前一节获取的是基于 ACS 原始表项提取的县级数据。这些数据虽然已经完成下载与空间匹配,但仍然属于底层统计字段,还不能直接用于地理人口特征分类分析。原因在于,当前字段中既包含总量变量,也包含不同年龄、族裔、教育与就业状态的分项计数,因此仍需进一步重构,才能形成适合比较的分析指标。

本节将依次完成以下工作:首先,根据原始 ACS 表项计算各类比例指标结构性指标;随后,将这些变量整理为适合聚类分析的分析数据框;最后,对不同量纲变量进行标准化处理,为下一节的聚类分析做好准备。

需要注意的是,原始 ACS 数据中既包含计数值,也包含收入中位数等连续指标。若不进行重构与标准化,直接进入聚类分析,结果往往更容易受变量量纲影响,而不一定反映真实的社会空间差异。


重构比例变量

# ==============================================================================
# 10. 县级分析变量重构
# ==============================================================================
# 说明:
# - 基于原始 ACS 表项,重构后续聚类分析所需的比例变量与结构性变量
# - 本节暂不进入聚类,仅完成变量计算与整理

acs_county_vars <- acs_county_conus %>%
  mutate(
    # --------------------------------------------------------------------------
    # A. 年龄结构
    # --------------------------------------------------------------------------
    age_0_17 = male_0_17_1E + male_0_17_2E + male_0_17_3E + male_0_17_4E +
               female_0_17_1E + female_0_17_2E + female_0_17_3E + female_0_17_4E,

    age_18_34 = male_18_34_1E + male_18_34_2E + male_18_34_3E +
                male_18_34_4E + male_18_34_5E + male_18_34_6E +
                female_18_34_1E + female_18_34_2E + female_18_34_3E +
                female_18_34_4E + female_18_34_5E + female_18_34_6E,

    age_35_64 = male_35_64_1E + male_35_64_2E + male_35_64_3E +
                male_35_64_4E + male_35_64_5E + male_35_64_6E + male_35_64_7E +
                female_35_64_1E + female_35_64_2E + female_35_64_3E +
                female_35_64_4E + female_35_64_5E + female_35_64_6E + female_35_64_7E,

    age_65p = male_65p_1E + male_65p_2E + male_65p_3E +
              male_65p_4E + male_65p_5E + male_65p_6E +
              female_65p_1E + female_65p_2E + female_65p_3E +
              female_65p_4E + female_65p_5E + female_65p_6E,

    pct_age_0_17  = age_0_17 / total_popE,
    pct_age_18_34 = age_18_34 / total_popE,
    pct_age_35_64 = age_35_64 / total_popE,
    pct_age_65p   = age_65p / total_popE,

    # --------------------------------------------------------------------------
    # B. 性别结构
    # --------------------------------------------------------------------------
    pct_female = female_totalE / total_popE,

    # --------------------------------------------------------------------------
    # C. 族裔结构
    # --------------------------------------------------------------------------
    other_race = other_race_1E + other_race_2E + other_race_3E + other_race_4E,

    pct_white = white_aloneE / race_totalE,
    pct_black = black_aloneE / race_totalE,
    pct_asian = asian_aloneE / race_totalE,
    pct_other_race = other_race / race_totalE,

    # --------------------------------------------------------------------------
    # D. 收入
    # --------------------------------------------------------------------------
    med_hh_income = med_hh_incomeE,

    # --------------------------------------------------------------------------
    # E. 就业
    # --------------------------------------------------------------------------
    unemployment_rate = unemployedE / labor_forceE,

    # --------------------------------------------------------------------------
    # F. 教育
    # --------------------------------------------------------------------------
    bachelor_plus = bachelorE + masterE + professionalE + doctorateE,
    pct_bachelor_plus = bachelor_plus / edu_total_25pE
  )
# ==============================================================================
# 11. 聚类分析用数据框整理
# ==============================================================================
# 说明:
# - 保留县级标识字段、最终分析变量与空间几何
# - 该对象将作为后续聚类分析与空间可视化的基础数据框

acs_county_analysis <- acs_county_vars %>%
  select(
    GEOID, NAME, geometry,
    pct_age_0_17, pct_age_18_34, pct_age_35_64, pct_age_65p,
    pct_female,
    pct_white, pct_black, pct_asian, pct_other_race,
    med_hh_income,
    unemployment_rate,
    pct_bachelor_plus
  ) %>%
  filter(
    if_all(
      c(
        pct_age_0_17, pct_age_18_34, pct_age_35_64, pct_age_65p,
        pct_female,
        pct_white, pct_black, pct_asian, pct_other_race,
        med_hh_income,
        unemployment_rate,
        pct_bachelor_plus
      ),
      ~ !is.na(.x) & is.finite(.x)
    )
  )

# 查看结果结构
head(acs_county_analysis)
Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -88.47323 ymin: 30.99708 xmax: -75.04894 ymax: 38.96187
Geodetic CRS:  NAD83
  GEOID                    NAME pct_age_0_17 pct_age_18_34 pct_age_35_64
1 01069 Houston County, Alabama    0.2290827     0.2033290     0.3807379
2 01023 Choctaw County, Alabama    0.1973344     0.1821486     0.3768982
3 01113 Russell County, Alabama    0.2431178     0.2232916     0.3809151
4 10005 Sussex County, Delaware    0.1809519     0.1565060     0.3603468
5 01071 Jackson County, Alabama    0.2079995     0.1901683     0.3929467
6 01089 Madison County, Alabama    0.2162044     0.2395605     0.3870965
  pct_age_65p pct_female pct_white  pct_black    pct_asian pct_other_race
1   0.1868504  0.5194378 0.6459035 0.27225818 0.0106990938     0.07113926
2   0.2436187  0.5216478 0.5674475 0.37552504 0.0002423263     0.05678514
3   0.1526755  0.5243577 0.4581966 0.45586710 0.0053731445     0.08056316
4   0.3021954  0.5148459 0.7414504 0.09942259 0.0126904149     0.14643659
5   0.2088855  0.5070024 0.8837201 0.03166645 0.0047688161     0.07984468
6   0.1571387  0.5074337 0.6302875 0.23907492 0.0255497661     0.10508777
  med_hh_income unemployment_rate pct_bachelor_plus
1         58626        0.04503522         0.2304748
2         47813        0.05034793         0.1576141
3         51197        0.06239683         0.1819799
4         81497        0.05330297         0.3298676
5         51908        0.04466211         0.1592226
6         86499        0.03605736         0.4642349
                        geometry
1 MULTIPOLYGON (((-85.71209 3...
2 MULTIPOLYGON (((-88.47323 3...
3 MULTIPOLYGON (((-85.4347 32...
4 MULTIPOLYGON (((-75.7226 38...
5 MULTIPOLYGON (((-86.15423 3...
6 MULTIPOLYGON (((-86.78955 3...
# ==============================================================================
# 13. 预处理结果可视化:重构变量分面地图(标准化前)
# ==============================================================================


# ------------------------------------------------------------------------------
# 13.1 选择示例变量
# ------------------------------------------------------------------------------
map_vars <- acs_county_analysis %>%
  select(
    GEOID, NAME, geometry,
    pct_age_65p,
    pct_bachelor_plus,
    unemployment_rate,
    med_hh_income
  )

# 构建经纬线对象
graticule_sf <- st_graticule(map_vars)

# ------------------------------------------------------------------------------
# 13.2 自定义绘图函数
# ------------------------------------------------------------------------------
# 说明:
# - 对每个变量分别应用 Fisher 断点
# - 各图使用不同 viridis 色盘,但统一保持“浅色 = 低值,深色 = 高值”
# - 通过 patchwork 拼接为 2 × 2 布局

plot_fisher_map <- function(data, var, title_text, legend_title, palette_option, digits = 2) {
  
  # 提取目标变量
  values <- st_drop_geometry(data)[[var]]
  
  # 计算 Fisher 断点
  brks <- classIntervals(values, n = 6, style = "fisher")$brks
  
  # 构造图例标签
  if (var == "med_hh_income") {
    labels <- paste0(
      "$", round(brks[-length(brks)] / 1000, 0), "k",
      " – ",
      "$", round(brks[-1] / 1000, 0), "k"
    )
  } else {
    labels <- paste0(
      percent(brks[-length(brks)], accuracy = 0.1),
      " – ",
      percent(brks[-1], accuracy = 0.1)
    )
  }
  
  # 分级
  plot_data <- data %>%
    mutate(
      class = cut(
        .data[[var]],
        breaks = brks,
        include.lowest = TRUE,
        labels = labels
      )
    )
  
  # 出图
  ggplot() +
    geom_sf(
      data = graticule_sf,
      color = "grey75",
      linewidth = 0.2,
      linetype = "dashed"
    ) +
    geom_sf(
      data = plot_data,
      aes(fill = class),
      color = NA
    ) +
    annotation_scale(
      location = "bl",
      width_hint = 0.22,
      text_cex = 0.65
    ) +
    annotation_north_arrow(
      location = "br",
      which_north = "true",
      style = north_arrow_fancy_orienteering,
      height = unit(0.75, "cm"),
      width  = unit(0.75, "cm")
    ) +
    scale_fill_viridis_d(
      option = palette_option,
      direction = -1,
      guide = guide_legend(reverse = TRUE),
      name = legend_title
    ) +
    coord_sf(crs = st_crs(4326)) +
    labs(title = title_text) +
    theme_minimal(base_size = 10.5) +
    theme(
      legend.position = "bottom",
      panel.grid.major = element_blank(),
      panel.background = element_rect(fill = "white", color = NA),
      plot.title = element_text(face = "bold", size = 11.5),
      legend.title = element_text(size = 9),
      legend.text = element_text(size = 8),
      axis.title = element_blank()
    )
}

# ------------------------------------------------------------------------------
# 13.3 分别绘制四幅地图
# ------------------------------------------------------------------------------
p_age65 <- plot_fisher_map(
  data = map_vars,
  var = "pct_age_65p",
  title_text = "Population aged 65+",
  legend_title = "Share of population",
  palette_option = "cividis"
)

p_edu <- plot_fisher_map(
  data = map_vars,
  var = "pct_bachelor_plus",
  title_text = "Bachelor's degree or above",
  legend_title = "Share of population 25+",
  palette_option = "viridis"
)

p_unemp <- plot_fisher_map(
  data = map_vars,
  var = "unemployment_rate",
  title_text = "Unemployment rate",
  legend_title = "Share of labour force",
  palette_option = "plasma"
)

p_income <- plot_fisher_map(
  data = map_vars,
  var = "med_hh_income",
  title_text = "Median household income",
  legend_title = "Income (1,000 USD)",
  palette_option = "magma"
)

# ------------------------------------------------------------------------------
# 13.4 拼接输出
# ------------------------------------------------------------------------------
(p_age65 | p_edu) /
  (p_unemp | p_income) +
  plot_annotation(
    title = "Selected reconstructed variables across CONUS counties",
    subtitle = "Fisher-classified maps before standardisation",
    caption = "Source: U.S. Census Bureau ACS via tidycensus"
  )


标准化处理

完成变量重构之后,可以看到当前分析变量并不处于同一量纲之下。例如,年龄、性别、族裔、失业率与本科及以上比例都属于比例变量,其取值通常落在 0 到 1 之间;而家庭收入中位数则以美元计量,数值范围明显更大。若直接将这些变量输入聚类模型,聚类结果往往更容易受数值尺度较大变量的影响,从而削弱不同变量在距离计算中的可比性。

因此,在进入聚类分析之前,通常需要先进行标准化处理,使各变量被转换到可比较的尺度上。这里采用最常见的 Z-score 标准化(均值为 0、标准差为 1),其基本形式为:

\[ Z = \frac{x - \mu}{\sigma} \]

其中,\(x\) 表示原始观测值,\(\mu\) 表示变量的均值,\(\sigma\) 表示变量的标准差。经过这一处理后,不同变量虽然原始单位不同,但都被转换到统一的标准尺度上,从而更适合后续的聚类分析

标准化并不只有一种方式。除 Z-score 之外,研究中也常见 Min-Max 归一化(将变量缩放到固定区间,如 0–1) 等方法。

# ==============================================================================
# 12. 聚类前标准化处理
# ==============================================================================
# 说明:
# - 对聚类输入变量进行 Z-score 标准化
# - 保留 GEOID 与 NAME 以便后续结果回接

cluster_vars <- acs_county_analysis %>%
  st_drop_geometry() %>%
  select(
    pct_age_0_17, pct_age_18_34, pct_age_35_64, pct_age_65p,
    pct_female,
    pct_white, pct_black, pct_asian, pct_other_race,
    med_hh_income,
    unemployment_rate,
    pct_bachelor_plus
  )

cluster_vars_scaled <- scale(cluster_vars) %>%
  as.data.frame()

acs_county_scaled <- acs_county_analysis %>%
  select(GEOID, NAME, geometry) %>%
  bind_cols(cluster_vars_scaled)

# 简要检查标准化结果
summary(cluster_vars_scaled)
  pct_age_0_17     pct_age_18_34     pct_age_35_64       pct_age_65p      
 Min.   :-6.0911   Min.   :-3.1919   Min.   :-8.79706   Min.   :-3.16917  
 1st Qu.:-0.5693   1st Qu.:-0.5448   1st Qu.:-0.45735   1st Qu.:-0.62858  
 Median : 0.0257   Median :-0.1243   Median : 0.08011   Median :-0.08968  
 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.5483   3rd Qu.: 0.3507   3rd Qu.: 0.58280   3rd Qu.: 0.48668  
 Max.   : 5.0618   Max.   : 9.4240   Max.   : 5.58812   Max.   :12.13915  
   pct_female         pct_white         pct_black          pct_asian        
 Min.   :-12.7630   Min.   :-3.8489   Min.   :-0.62038   Min.   :-0.532858  
 1st Qu.: -0.2353   1st Qu.:-0.6048   1st Qu.:-0.57268   1st Qu.:-0.422698  
 Median :  0.1596   Median : 0.3486   Median :-0.46261   Median :-0.286247  
 Mean   :  0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.000000  
 3rd Qu.:  0.4887   3rd Qu.: 0.8085   3rd Qu.: 0.07054   3rd Qu.:-0.002275  
 Max.   :  4.4495   Max.   : 1.3011   Max.   : 5.46941   Max.   :14.972561  
 pct_other_race    med_hh_income     unemployment_rate pct_bachelor_plus
 Min.   :-1.0771   Min.   :-2.5754   Min.   :-1.9704   Min.   :-2.3965  
 1st Qu.:-0.5967   1st Qu.:-0.6296   1st Qu.:-0.6149   1st Qu.:-0.6964  
 Median :-0.3547   Median :-0.1244   Median :-0.1048   Median :-0.2487  
 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
 3rd Qu.: 0.1887   3rd Qu.: 0.4420   3rd Qu.: 0.4329   3rd Qu.: 0.4687  
 Max.   : 6.4328   Max.   : 6.4345   Max.   :10.1963   Max.   : 5.4183  
聚类分析与类型识别
2.4 聚类分析导入

聚类分析

在完成数据重构标准化之后,当前数据已经具备进入聚类分析的基本条件。对于地理人口特征研究而言,聚类分析是其最核心的技术环节之一。其基本思路在于:将每一个地理单元视为一个由多维属性构成的观测对象,再依据这些属性之间的相似性或差异性,将具有相近社会人口背景的地区归并为若干类别。

聚类分析的目标并不是考察某一个变量与结果变量之间的线性关系,也不是验证某种因果机制,而是从多维属性组合出发,识别数据内部自然形成的群组结构。对于本示例研究来说,这意味着将美国本土各县视为一组由年龄结构、性别、族裔、收入、就业与教育等变量共同刻画的空间单元,并进一步识别哪些县在这些维度上更为相似,哪些县则表现出明显不同的社会人口特征。

地理人口特征分类的角度看,聚类分析的意义并不只是生成若干数学上的分组结果,更重要的是为后续的类别解释区域命名奠定基础。也就是说,聚类得到的每一类,并不只是“第 1 类”、“第 2类”这样的抽象标签,而是有可能进一步被解释为具有社会空间含义的区域类型。

聚类分析在地理人口特征研究中通常不是终点,而是连接数据结构识别类型解释的关键步骤。


K-means 聚类

在本节中,将首先采用 K-means 聚类(K-means clustering) 作为示例方法。K-means 是聚类分析中最常见、也最易理解的方法之一,其基本思想较为直观:通过不断迭代,将观测对象划分为若干组,并使组内差异尽可能小、组间差异尽可能大。对于当前这种已经完成标准化、且以数值变量为主的数据而言,K-means 能够提供一个较清晰的入门分析框架,便于展示 geodemographics 的基本技术流程。

当然,K-means 不是唯一可用的方法。相关研究与实践中,也常见层次聚类(hierarchical clustering)模糊聚类(fuzzy c-means)高斯混合模型(Gaussian mixture models)等不同方法。不同聚类算法在类别边界假设、对异常值的敏感性、类别形状的适应性以及结果解释方式上各有差异。

方法提示:聚类分析与非监督学习

机器学习的角度看,聚类分析属于非监督学习(unsupervised learning)的典型方法。与分类模型或回归模型不同,它不依赖预先给定的“正确答案”或标签变量,而是根据数据自身的结构特征,自动识别潜在的群组模式

这也是geodemographics与传统变量分析的一个重要区别:其研究重点并不在于预测某一个外部结果,而在于发现不同地区在多维属性组合上的内在相似性与差异性


本节任务

基于前述已经完成的预处理结果,接下来的分析将正式进入聚类分析阶段。本节将以 K-means 算法为核心,对标准化后的县域多维属性数据进行分组,从而识别美国本土不同县在社会人口背景上的潜在类型结构。

  • 第一,使用 K-means 对标准化后的变量进行聚类分析,并理解其基本实现逻辑。

  • 第二,结合聚类分析中的关键输入参数,尤其是聚类数量(number of clusters, k的设定,对不同分类结果进行初步比较与判断。

  • 第三,将聚类结果重新连接回原始空间数据,并通过面向离散类别的分层设色图展示其基本空间格局,为后续的类别解释与命名提供直观基础。


Step 1

# ==============================================================================
# 14. 聚类输入数据准备
# ==============================================================================
# 说明:
# - 本节使用前面已经完成标准化处理的数据
# - 聚类算法只接受数值型输入,因此这里先去除空间几何列与标识列

kmeans_input <- acs_county_scaled %>%
  st_drop_geometry() %>%
  select(-GEOID, -NAME)

# 检查输入数据结构
dim(kmeans_input) # 3108个地理单元 12个变量
[1] 3108   12
head(kmeans_input,3)
  pct_age_0_17 pct_age_18_34 pct_age_35_64 pct_age_65p pct_female  pct_white
1    0.3191397    0.02810222    0.18261087  -0.3553552  0.9461956 -0.6606159
2   -0.5692470   -0.43952906    0.04975571   0.7681317  1.0359355 -1.0952765
3    0.7118733    0.46884524    0.18874214  -1.0317004  1.1459706 -1.7005466
  pct_black  pct_asian pct_other_race med_hh_income unemployment_rate
1  1.305309 -0.1232341     -0.5019596    -0.5180874       -0.05625866
2  2.035717 -0.5235807     -0.6180074    -1.1286002        0.16954679
3  2.603978 -0.3271428     -0.4257710    -0.9375362        0.68165976
  pct_bachelor_plus
1        -0.1600036
2        -0.8670362
3        -0.6305927

Step 2 重要

确定聚类数量(k)

确定聚类数量的方法不止有以下示例中的两种;可以询问AI更多方法来确定k

# ==============================================================================
# 15. 聚类数量 k 的初步选择
# ==============================================================================
# 说明:
# - 这里使用两种常见方法辅助判断 k 的取值
#   1) Elbow method:观察组内平方和(WSS)的下降趋势
#   2) Silhouette method:观察平均轮廓系数的变化
# - 作为教学示例,这里先考察 2 到 10 类

library(factoextra) # 调用包

set.seed(1234)

p_elbow <- fviz_nbclust(
  x = kmeans_input,
  FUNcluster = kmeans,
  method = "wss",
  k.max = 10,
  nstart = 25
) +
  labs(
    title = "Elbow method",
    subtitle = "Within-cluster sum of squares"
  ) +
  theme_minimal(base_size = 11)

p_silhouette <- fviz_nbclust(
  x = kmeans_input,
  FUNcluster = kmeans,
  method = "silhouette",
  k.max = 10,
  nstart = 25
) +
  labs(
    title = "Silhouette method",
    subtitle = "Average silhouette width"
  ) +
  theme_minimal(base_size = 11)

p_elbow + p_silhouette

如何理解聚类数量(k)的选择?

K-means 聚类 中,聚类数量(number of clusters, k需要由研究者预先设定。一个常见做法是借助 Elbow method(肘部法) 进行初步判断。其基本思想是:随着 k 增大,组内平方和 (within-cluster sum of squares, WSS)通常会不断下降,因为类别越多,每一类内部自然更容易变得紧密。

在实际图形中,如果曲线在某一位置之前下降较快、之后下降幅度明显放缓,那么这个“由陡转缓”的位置通常可被理解为一个较有参考价值拐点。这意味着继续增加类别数虽然仍会改善组内紧密度,但改善幅度已经开始减弱。因此,研究中常将这一“肘部”附近的 k 视为较合理的候选值。

不过需要注意的是,肘部法通常只能提供一种启发式判断,并不总会给出一个非常清晰、唯一的答案。有些数据会呈现出明显的拐点,有些则可能只是整体平滑下降,此时就需要结合其他信息一并判断,例如轮廓系数类别空间分布类别规模是否过小,以及分类结果是否具有可解释性等。

地理人口特征分类研究中,k 的选择并不一定单纯追求“统计上最优”或“图上最明显”的那个值。

由于 geodemographics 本身具有较强的应用导向,研究者往往还需要结合实际问题、变量体系、空间尺度与结果解释的便利性来综合考虑。换言之,一个在统计指标上略逊一筹、但在社会空间解释上更清晰、在应用场景中更有意义的分类结果,有时反而更适合作为最终方案。


Step 3

聚类分析 (k = 4)

# ==============================================================================
# 16. 执行 K-means 聚类
# ==============================================================================
# 说明:
# - 这里先设定一个示例性的聚类数量 k=4 (根据上一步的图观察)
# - 后续可根据 Elbow 与 Silhouette 图形进一步调整
# - nstart = 500 用于提高结果稳定性 (真实情况下n的数值会更大)

set.seed(1234)

k_selected <- 4

kmeans_model <- kmeans(
  x = kmeans_input,
  centers = k_selected,
  nstart = 500
)

# 查看聚类结果概况
# kmeans_model
# table(kmeans_model$cluster) # 查看聚类里的地块数量 (和下面一致)

Step 4

聚类结果放回原是数据中,进行空间制图准备

# ==============================================================================
# 17. 将聚类结果回接到空间数据
# ==============================================================================
# 说明:
# - 将 K-means 的聚类编号写回县级空间数据
# - 后续即可进行空间制图与类别解释

acs_county_clustered <- acs_county_analysis %>%
  mutate(
    cluster = factor(kmeans_model$cluster)
  )

# 简要检查结果
acs_county_clustered %>%
  st_drop_geometry() %>%
  count(cluster)
  cluster    n
1       1  353
2       2  406
3       3  458
4       4 1891

Step 5

聚类结果空间制图

# ==============================================================================
# 18. 聚类结果空间制图
# ==============================================================================
# 说明:
# - 这里对离散型聚类结果进行分层设色
# - 使用定性配色方案展示不同类别
# - 同时加入经纬线、比例尺与指北针

library(RColorBrewer)

# 构建经纬线对象
graticule_sf <- st_graticule(acs_county_clustered)

ggplot() +
  # 经纬线
  geom_sf(
    data = graticule_sf,
    color = "grey75",
    linewidth = 0.2,
    linetype = "dashed"
  ) +
  # 聚类结果图层
  geom_sf(
    data = acs_county_clustered,
    aes(fill = cluster),
    color = NA
  ) +
  # 比例尺
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    text_cex = 0.75
  ) +
  # 指北针
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    style = north_arrow_fancy_orienteering,
    height = unit(0.85, "cm"),
    width  = unit(0.85, "cm")
  ) +
  # 定性配色
  scale_fill_brewer(
    palette = "Set1",
    name = "Cluster"
  ) +
  coord_sf(crs = st_crs(4326)) +
  labs(
    title = "K-means clustering of CONUS counties",
    subtitle = paste0("k = ", k_selected),
    caption = "Source: U.S. Census Bureau ACS via tidycensus"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 11),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 9),
    axis.title = element_blank()
  )

结果解释与命名
2.5 聚类解释与命名

结果解释

在完成 K-means聚类并对其空间分布进行初步可视化之后,地理人口特征分类的工作并未结束。相反,聚类结果的真正价值,往往体现在对各类别内部特征的进一步识别、比较与解释之中。

前一节得到的 1、2、3、4类,只是算法返回的离散编号,本身并不直接具有社会空间意义。若要将其转化为可理解、可交流、可应用的地理人口特征类型,还需要进一步分析每一类在多维属性上的相对特征,并据此完成类型命名类别描述

从 geodemographics 的角度看,聚类分析的目标并不是简单地将地区“分成几组”,而是识别出不同地理单元在人口结构社会经济属性及其他背景维度上的组合差异。因此,聚类结果解释的核心任务在于回答以下问题:每一类地区在哪些变量上表现得最为突出?这些特征是相对于全国平均水平而言偏高,还是偏低?不同类型之间最显著的差异体现在哪些维度上?只有在回答这些问题之后,聚类结果才真正从抽象的数学分组转化为具有解释力的区域类型。

为什么需要类型命名?

聚类算法生成的类别编号通常是任意的。例如,Cluster 1 并不天然意味着“最好”“最差”,也不代表某种固定的社会属性。若直接停留在数值编号层面,结果往往难以传达其社会空间含义,也不利于后续比较、制图与应用。因此,在地理人口特征分类研究中,研究者通常会进一步对各类别进行命名,使其从算法输出转换为具有解释性的地区类型

这种命名并不是随意的标签替换,而是建立在对类别内部多维特征进行系统比较的基础之上。通常而言,一个类别的名称往往会优先突出其最具辨识度的属性组合。例如,若某一类地区同时表现出较高的老年人口比例、较低的收入水平与相对较弱的教育背景,则其命名可能会围绕老龄化低收入教育水平偏低等关键词展开;若另一类地区在高收入、本科及以上人口比例较高、失业率较低等方面更为突出,则其命名可能更偏向高教育高收入优势型县域等表述。

在命名过程中,还需要特别注意语言的中性准确非污名化。应尽量避免使用带有歧义、歧视性、攻击性、道德评判色彩或可能强化刻板印象的词语,而更适合使用描述性结构性相对客观的表述。

因此,这里的类型命名并不是为了制造标签、强化偏见或挑起社会对立,而是为了更准确地概括不同地区在多维社会人口属性上的结构差异,使聚类结果更便于解释、交流与应用。


解释逻辑与类型画像

在实际分析中,对聚类结果的解释通常需要结合类别中心或各类别主要变量的均值水平,并将其与研究区域整体的平均水平进行比较。换言之,某一类的“突出特征”通常并不是绝对意义上的高或低,而是相对于研究区域总体更为明显。例如,当研究范围为全国时,可与全国平均水平比较;当研究范围为某一城市时,则更适合与城市整体平均水平比较。与此同时,还需要结合不同类别之间的差异以及其空间分布格局,进一步判断某一类最具代表性的多维属性组合。

聚类结果解释不是对每一个变量逐项罗列,而是从多个维度中提炼出最能代表该类别的核心特征,并进一步形成一种相对简洁、但具有解释力的类型画像。这类画像并非针对个体,而是针对地区类型本身,可理解为对某一类地理单元社会人口背景的概括性描述。

需要强调的是,这类“画像”只是对某一类别整体特征的概括,并不意味着其中每一个县都完全符合该类的所有属性。换言之,类型解释强调的是总体趋势而非个体层面的绝对判断。因此,在命名与描述时,更适合使用“相对较高”“相对突出”“更倾向于呈现……特征”或 more likely to 等较为稳健的表述,而不宜将某一类别简单理解为所有地区都具有同样特征。


本节任务

基于前一节已经得到的 4类聚类结果,本节将进一步进入聚类结果解释阶段,将围绕以下几个方面展开:

  • 结合多种可视化方法对聚类结果进行比较与解读;
  • 计算各类别在变量上的指数得分(Index Score),并将其与研究区域整体平均水平进行比较,以识别各类最突出的多维属性;
  • 使用热力图雷达图等方式展示不同类别相对于整体平均水平的特征差异,并辅助识别各类的核心画像;
  • 总结与反思部分(下一节):综合前述空间分布格局属性特征图表,对聚类类别进行简要描述,并尝试赋予更具解释性的类型名称。

为了更直观地解读每个类别的含义,本节引入指数得分(Index Score)的概念,其基本形式为:

\[ \text{Index} = \frac{\text{类别均值}}{\text{整体均值}} \times 100 \]

其中,Index = 100 表示与整体平均水平一致;Index > 100 表示高于平均水平;Index < 100 则表示低于平均水平。通过这一方式,可以更直观地判断某一类在不同变量上相对于整体背景的偏离方向与偏离程度。


实战代码 计算相对指数

# ==============================================================================
# 21. 计算各聚类相对于全国平均的指数得分(Index Score)
# ==============================================================================
# 说明:
# - Index Score = (聚类均值 / 全国平均) × 100
# - Index = 100 表示与全国平均一致
# - Index > 100 表示高于全国平均
# - Index < 100 表示低于全国平均

# ------------------------------------------------------------------------------
# 21.1 指定用于解释聚类的原始变量
# ------------------------------------------------------------------------------
profile_vars <- c(
  "pct_age_0_17",
  "pct_age_18_34",
  "pct_age_35_64",
  "pct_age_65p",
  "pct_female",
  "pct_white",
  "pct_black",
  "pct_asian",
  "pct_other_race",
  "med_hh_income",
  "unemployment_rate",
  "pct_bachelor_plus"
)

# ------------------------------------------------------------------------------
# 21.2 计算全国平均
# ------------------------------------------------------------------------------
overall_means <- acs_county_clustered %>%
  st_drop_geometry() %>%
  select(all_of(profile_vars)) %>%
  summarise(
    across(everything(), mean, na.rm = TRUE)
  )

# ------------------------------------------------------------------------------
# 21.3 计算各聚类均值
# ------------------------------------------------------------------------------
cluster_means <- acs_county_clustered %>%
  st_drop_geometry() %>%
  group_by(cluster) %>%
  summarise(
    across(all_of(profile_vars), mean, na.rm = TRUE),
    .groups = "drop"
  )

# ------------------------------------------------------------------------------
# 21.4 计算指数得分
# ------------------------------------------------------------------------------
index_scores <- cluster_means %>%
  mutate(
    across(
      all_of(profile_vars),
      ~ .x / overall_means[[cur_column()]] * 100
    )
  ) %>%
  mutate(
    across(where(is.numeric), ~ round(.x, 1))
  )

# 查看指数得分表
index_scores
# A tibble: 4 × 13
  cluster pct_age_0_17 pct_age_18_34 pct_age_35_64 pct_age_65p pct_female
  <fct>          <dbl>         <dbl>         <dbl>       <dbl>      <dbl>
1 1              119.          112.           93.8        79.5       98.2
2 2               99.8         105.          100          95        102. 
3 3               98.5         121.          100.         80.1      101. 
4 4               96.8          91.5         101.        110.        99.7
# ℹ 7 more variables: pct_white <dbl>, pct_black <dbl>, pct_asian <dbl>,
#   pct_other_race <dbl>, med_hh_income <dbl>, unemployment_rate <dbl>,
#   pct_bachelor_plus <dbl>
# ==============================================================================
# 22. 各聚类指数得分热力图(7组分级)
# ==============================================================================
# 说明:
# - x 轴:变量
# - y 轴:聚类类别
# - fill:Index Score 的分组区间
# - 100 为基准值,接近 100 的区间使用白色
# - 高于平均使用绿色系,低于平均使用红色系
# - 极端高值与极端低值分别单独设组,以避免压缩整体色阶

# ------------------------------------------------------------------------------
# 22.1 转为长表,并设置变量标签
# ------------------------------------------------------------------------------
variable_labels <- c(
  pct_age_0_17      = "Age 0–17",
  pct_age_18_34     = "Age 18–34",
  pct_age_35_64     = "Age 35–64",
  pct_age_65p       = "Age 65+",
  pct_female        = "Female",
  pct_white         = "White",
  pct_black         = "Black",
  pct_asian         = "Asian",
  pct_other_race    = "Other race",
  med_hh_income     = "Median income",
  unemployment_rate = "Unemployment",
  pct_bachelor_plus = "Bachelor+"
)

index_long <- index_scores %>%
  pivot_longer(
    cols = -cluster,
    names_to = "variable",
    values_to = "index"
  ) %>%
  mutate(
    variable = factor(variable, levels = profile_vars),

    # --------------------------------------------------------------------------
    # 22.2 按指数区间分组(7组)
    # --------------------------------------------------------------------------
    index_group = case_when(
      index <= 50                ~ "<= 50",
      index > 50  & index < 80   ~ "50–80",
      index >= 80 & index < 95   ~ "80–95",
      index >= 95 & index <= 105 ~ "95–105",
      index > 105 & index <= 120 ~ "105–120",
      index > 120 & index < 200  ~ "120–200",
      index >= 200               ~ ">= 200"
    ),
    index_group = factor(
      index_group,
      levels = c("<= 50", "50–80", "80–95", "95–105", "105–120", "120–200", ">= 200")
    )
  )

# ------------------------------------------------------------------------------
# 22.3 绘制热力图
# ------------------------------------------------------------------------------
ggplot(index_long, aes(x = variable, y = cluster, fill = index_group)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(
    aes(label = round(index, 1)),
    size = 3.6,
    color = "black"
  ) +
  scale_fill_manual(
    values = c(
      "<= 50"   = "#B2182B",
      "50–80"   = "#EF8A62",
      "80–95"   = "#FDDBC7",
      "95–105"  = "white",
      "105–120" = "#D9F0D3",
      "120–200" = "#7FBF7B",
      ">= 200"  = "#1B7837"
    ),
    name = "Index score"
  ) +
  scale_x_discrete(labels = variable_labels) +
  labs(
    title = "Cluster profiles relative to the national average",
    subtitle = "Grouped index score: 100 = national average",
    x = NULL,
    y = "Cluster"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1,
      size = 10
    ),
    axis.text.y = element_text(size = 11),
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 10.5),
    panel.grid = element_blank(),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 9)
  )+coord_fixed()


雷达图辅助解释 参考 · 可选

除热力图之外,雷达图(radar plot)也常用于展示不同类别在多维属性上的整体轮廓。相较于逐变量比较的表格或热力图,雷达图更强调各类别在多个指标上的相对形态差异,因此尤其适合用于辅助构建类型画像

这里继续采用前文计算得到的 index score 作为输入,其中 100 表示与全国平均水平一致。为避免极端高值或低值压缩整体图形尺度,后续将指数统一限制在一定范围内,即将 50 设为下限、250 设为上限,并将超出范围的值分别截断为对应边界值。通过这一处理,可以更清晰地比较不同聚类在多维属性上的相对突出方向与整体轮廓。

需要注意的是:当我们选取的变量特别多时候,建议还是使用热力图展示

# ==============================================================================
# 23. 雷达图输入数据准备
# ==============================================================================
# 说明:
# - 基于前面已经计算好的 index_scores
# - 将指数统一截断在 50 到 250 之间
# - 后续使用 fmsb::radarchart() 绘图

library(fmsb)

# ------------------------------------------------------------------------------
# 23.1 选择用于雷达图的变量
# ------------------------------------------------------------------------------
radar_vars <- c(
  "pct_age_0_17",
  "pct_age_18_34",
  "pct_age_35_64",
  "pct_age_65p",
  "pct_female",
  "pct_white",
  "pct_black",
  "pct_asian",
  "pct_other_race",
  "med_hh_income",
  "unemployment_rate",
  "pct_bachelor_plus"
)

# ------------------------------------------------------------------------------
# 23.2 构造变量标签
# ------------------------------------------------------------------------------
radar_labels <- c(
  pct_age_0_17      = "Age 0–17",
  pct_age_18_34     = "Age 18–34",
  pct_age_35_64     = "Age 35–64",
  pct_age_65p       = "Age 65+",
  pct_female        = "Female",
  pct_white         = "White",
  pct_black         = "Black",
  pct_asian         = "Asian",
  pct_other_race    = "Other race",
  med_hh_income     = "Income",
  unemployment_rate = "Unemployment",
  pct_bachelor_plus = "Bachelor+"
)

# ------------------------------------------------------------------------------
# 23.3 整理指数数据并进行截断
# ------------------------------------------------------------------------------
radar_index <- index_scores %>%
  select(cluster, all_of(radar_vars)) %>%
  mutate(
    across(
      all_of(radar_vars),
      ~ pmin(pmax(.x, 50), 250)
    )
  )

# 查看结果
radar_index
# A tibble: 4 × 13
  cluster pct_age_0_17 pct_age_18_34 pct_age_35_64 pct_age_65p pct_female
  <fct>          <dbl>         <dbl>         <dbl>       <dbl>      <dbl>
1 1              119.          112.           93.8        79.5       98.2
2 2               99.8         105.          100          95        102. 
3 3               98.5         121.          100.         80.1      101. 
4 4               96.8          91.5         101.        110.        99.7
# ℹ 7 more variables: pct_white <dbl>, pct_black <dbl>, pct_asian <dbl>,
#   pct_other_race <dbl>, med_hh_income <dbl>, unemployment_rate <dbl>,
#   pct_bachelor_plus <dbl>
# ==============================================================================
# 24. 各聚类指数轮廓雷达图
# ==============================================================================
# 说明:
# - 每个 cluster 单独绘制一张雷达图
# - 50 为下限,250 为上限
# - 100 对应全国平均水平,可作为解读基准

# ------------------------------------------------------------------------------
# 24.1 定义雷达图绘制函数
# ------------------------------------------------------------------------------
plot_radar_cluster <-  function(cluster_id, line_col) {
  
  # 提取单个 cluster 的数据
  df_cluster <- radar_index %>%
    filter(cluster == cluster_id) %>%
    select(all_of(radar_vars))
  
  # fmsb 要求前两行为 max 和 min
  df_plot <- rbind(
    rep(250, length(radar_vars)),
    rep(50,  length(radar_vars)),
    df_cluster
  )
  
  colnames(df_plot) <- radar_labels[radar_vars]
  
  # 作图
  radarchart(
    df_plot,
    axistype = 1,
    pcol = line_col,
    pfcol = scales::alpha(line_col, 0.25),
    plwd = 2,
    plty = 1,
    cglcol = "grey75",
    cglty = 1,
    cglwd = 0.8,
    axislabcol = "grey40",
    vlcex = 0.85,
    title = paste("Cluster", cluster_id)
  )
  
  # 添加基准说明
  mtext("Index baseline = 100", side = 1, line = 0.5, cex = 0.8, col = "grey35")
}

# ------------------------------------------------------------------------------
# 24.2 依次绘制四个聚类的小雷达图
# ------------------------------------------------------------------------------
par(mfrow = c(2, 2), mar = c(2, 2, 3, 2))

plot_radar_cluster("1", "#1b9e77")
plot_radar_cluster("2", "#d95f02")
plot_radar_cluster("3", "#7570b3")
plot_radar_cluster("4", "#e7298a")

结果说明

截至目前,本章已经完成聚类结果的生成,并通过空间分布图指数热力图雷达图等方式展示了不同类别在空间格局与内部属性上的基本差异。

不过,这一部分的重点仍然是结果呈现图形解读基础,而不是对各类别进行完整命名。关于不同聚类的简要描述类型命名,将放在下一节 3.0 总结与反思 中,结合其空间分布特征与内部属性表现进行统一讨论。

3.0 专题研究 1: 回顾与反思

3.0 专题研究 1:回顾与反思

重要

本章将对前述专题研究 1“美国全国范围县级地理人口特征分类”进行简要回顾。需要说明的是,受教学场景与篇幅所限,该示例更适合作为一个方法演示型案例,而非完整展开的正式科研项目。不过,其所呈现的地理人口特征分类(geodemographic classification)工作流,确实反映了人口普查数据分析、区域类型识别与社会空间研究中较为常见的一类分析路径。

具体而言,前文已经依次完成了变量选择与数据获取指标重构与标准化K-means 聚类以及空间可视化与结果展示等关键步骤。通过这一过程,原本分散的人口与社会经济变量被进一步整合为具有解释潜力的区域类型,也初步展示了美国本土各县在多维社会人口背景上的差异格局。

在此基础上,本节将不再重复方法实现过程,而是回到结果本身,对前面已经得到的空间分布图指数热力图雷达图进行综合回顾。重点在于:结合不同聚类在空间格局与内部属性上的相对特征,对各类别作出简要描述,并尝试赋予更具解释性的类型名称。


聚类描述与命名

基于前文的聚类空间分布图Index 热力图,美国各 county4 个聚类可作如下初步命名与简要描述:

  • Cluster 1:年轻多元族裔聚集区 红色
    这一类型主要分布在美国西南部边境地区,包括德克萨斯州南部、新墨西哥州、亚利桑那州,以及加利福尼亚州中部农业谷地和中西部部分城市周边。其最突出的特征是人口结构年轻少数族裔占比较高。0–17 岁人口指数达到 119.3,18–34 岁人口指数为 111.6,而 65 岁及以上人口指数仅为 79.5,显示出较强的年轻化特征。族裔结构上,白人人口占比低于全国平均水平,指数为 72.3,黑人人口占比也较低,指数为 49.6;而其他族裔指数高达 293.2。结合其空间分布,这一特征很可能主要反映了拉美裔人口的高度集聚。社会经济方面,该类地区本科及以上人口占比指数为 80.5,家庭收入中位数指数为 93.8,失业率指数则达到 124.4。整体来看,这类地区可理解为以拉美裔或其他少数族裔为主的年轻社区,具有较强的人口增长潜力,但也面临教育水平偏低、收入水平有限与就业压力较高等问题。

  • Cluster 2:南部非裔美国人传统社区 蓝色
    这一类型呈现出较强的地理连续性,主要集中在美国南部腹地,从德克萨斯州东部延伸至路易斯安那州、密西西比州、阿拉巴马州、佐治亚州和南卡罗来纳州。这一区域与美国历史上的 Black Belt 空间格局高度相关。其最核心的特征是黑人人口高度集中,黑人人口占比指数达到 444.8,显著高于全国平均水平;相应地,白人人口占比指数仅为 67.8。从年龄结构看,该类地区各年龄段指数大体接近全国平均水平,未表现出极端年轻化或老龄化特征。社会经济方面,这一类型的经济弱势较为突出,家庭收入中位数指数仅为 77.7,本科及以上人口占比指数为 76.2,失业率指数则高达 142.5,是四类中最高。整体来看,该类型反映了美国南部非裔美国人传统聚居区的延续性,也揭示出其在收入、教育和就业机会方面面临的结构性不利处境。

  • Cluster 3:高学历亚裔都市精英区 绿色
    这一类型在空间上多呈现点状或小片区分布,主要集中于美国主要大都市区及其富裕郊区,例如加利福尼亚沿海地区、旧金山湾区、洛杉矶、东北部城市走廊,以及西雅图、芝加哥等大城市周边。其最鲜明的特征是亚裔人口高度集聚高教育、高收入并存。亚裔人口占比指数达到 362.5,约为全国平均水平的 3.6 倍;18–34 岁人口指数为 121.3,显示出对年轻劳动力和专业人才的吸引力。社会经济方面,该类地区家庭收入中位数指数达到 136.9,本科及以上人口占比指数高达 170.1,均显著高于全国平均水平;失业率指数为 96.4,略低于全国平均。整体来看,这类地区更接近由高学历、高收入群体构成的都市型社区,往往与科技产业、金融服务、高等教育和知识密集型就业机会高度相关。

  • Cluster 4:白人主导的农村及老龄化社区 紫色
    这一类型覆盖面积最广,主要分布在美国中西部大平原、北部农村地区、阿巴拉契亚山脉部分地区,以及新英格兰的若干农村县域。其核心特征是白人人口占比较高族裔多样性较低人口老龄化趋势较明显。白人人口占比指数为 114.1,高于全国平均水平;黑人人口和亚裔人口指数分别仅为 33.946.3。年龄结构上,65 岁及以上人口指数为 109.7,而 18–34 岁人口指数为 91.5,显示出一定程度的老龄化和年轻人口外流特征。社会经济方面,该类地区家庭收入中位数指数为 97.0,本科及以上人口占比指数为 91.8,整体接近或略低于全国平均水平;失业率指数为 87.2,是四类中最低。整体来看,这类地区代表了美国广泛分布的农村与小城镇腹地,社会结构相对稳定,但在人口更新、高学历人才集聚和经济增长动力方面相对有限。


专题研究 1 反思

重要

任何地理人口特征分类流程都包含一定程度的变量筛选空间范围界定地理单元选择算法参数设定。前述案例虽然初步展示了美国本土县域在人口与社会经济属性上的类型差异,但从更严格的研究角度看,相关结果仍需结合研究目标、数据来源与方法选择进行审慎理解。


以下几个问题,可作为对本案例的进一步反思:

  • 变量选择的依据是什么?
    本案例仅选取了年龄、性别、族裔、收入、教育与就业等少量变量,主要服务于教学演示,因此更适合作为一个简化版的地理人口特征分类流程示例,而非完整的变量体系。在更正式的研究中,变量选择往往需要围绕研究目标展开,并结合既有文献、前期研究经验与数据可得性综合判断。

例如,若研究主题转向公共交通出行相关分类,则可能需要进一步考虑通勤方式、车辆拥有、人口密度、住房类型、职住关系与可达性等变量。由此也引出一个值得思考的问题:在面对大量可选变量时,哪些变量是真正必要的?哪些变量只是“看起来有用”,但未必与研究问题直接相关?

  • 空间范围与地理单元会影响分类结果吗?
    本案例采用的是美国本土县级单元,且研究范围覆盖全国。但若改为州内研究、都市区研究,或缩小到某一城市内部,分类结果很可能会发生变化。同样地,若将县级单元替换为更细的小区域单元,如 tractblock group 或其他邻里尺度,所识别出的类型也往往会更细致、更具局部差异。这提示出一个重要问题:地理人口特征分类并不是脱离空间尺度而独立存在的,地理范围与观测单元的改变,本身就可能重塑分类结构。因此,在实际研究中,如何界定研究区域、如何选择地理单元的颗粒度,通常需要同时考虑研究目的数据可用性数据获取成本结果解释尺度等因素。

  • K-means 的结果是否稳定?
    本案例使用 K-means 作为聚类方法,这一方法具有实现简洁、解释直观等优点,但它本身也包含若干值得反思的设定。例如,初始中心点的随机生成、nstart 的设置,以及最终聚类数 k 的选择,都可能影响结果的稳定性与可重复性。换言之,即便变量与研究区域不变,不同参数设定也可能得到略有差异的分类结果。由此可以进一步思考:当前案例中 k 的选择是否足够合理?若改用其他判别方式,结果会不会不同?而在一个强调解释性与应用性的 geodemographics 研究中,是否应只依据统计指标来决定最终分类方案?

  • 是否存在其他可替代的聚类方法?
    K-means 并不是唯一可用的方法。实际研究中,还常见层次聚类模糊聚类高斯混合模型等不同路径。不同方法在类别边界、异常值敏感性、类别形状假设以及结果解释方式上各有特点。因此,一个值得继续思考的问题是:对于不同的研究对象、不同的数据尺度与不同的变量体系,是否存在比 K-means 更合适的聚类方法?若更换算法,分类结果的整体轮廓是否仍能保持稳定?

  • 从方法演示到正式研究,还缺少什么?
    本案例较好展示了从数据获取变量重构、再到聚类分析结果解释的完整流程,但它仍然是一个以教学为导向的示例。若进入更正式的研究设计,通常还需要进一步考虑变量筛选的理论依据、不同尺度下结果的一致性、聚类方案的稳健性检验,以及分类结果在具体应用场景中的解释价值。也就是说,本案例所回答的更多是“如何做出一个基础分类”,而更深入的问题则是:这个分类为何成立是否稳健、以及能否服务于特定研究问题

延伸方法:什么是层次聚类(Hierarchical Clustering)?

K-means 之外,层次聚类(Hierarchical Clustering)也是 geodemographics 中非常常见的一类聚类方法。它的核心特点在于:不必一开始就固定类别中心,而是通过样本之间的相似性,逐步形成一个具有层级结构的聚类树 (dendrogram)。研究者随后可以根据树状结构,在不同高度切分出不同数量的类别。 oai_citation:1‡Seminar for Statistics

常见的层次聚类大致有两类:

  • 凝聚型(Agglomerative):从“每个样本各自成类”开始,逐步把更相似的样本或类别合并起来;
  • 分裂型(Divisive):从“所有样本先视为一类”开始,再逐步向下拆分。
    其中,在 R 的实际应用中,更常见的是凝聚型层次聚类oai_citation:2‡Seminar for Statistics

在 R 中,最常用的基础实现方式通常是:

  • dist():先计算样本之间的距离矩阵;
  • hclust():进行层次聚类;
  • cutree():将树切分为指定数量的类别;
  • plot()fviz_dend():绘制树状图。
    若想进一步做分裂型层次聚类,也可以使用 cluster 包中的 diana()oai_citation:3‡Seminar for Statistics

一个最基础的实现流程如下:

dist_mat <- dist(x, method = "euclidean")
hc <- hclust(dist_mat, method = "ward.D2")
plot(hc)
cluster_id <- cutree(hc, k = 4)

阅读与教程推荐

层次聚类的优势在于可以直接观察类别的层级关系,且不必像 K-means 那样强依赖初始中心点;但它在大样本下计算与可视化成本通常更高,因此在实际研究中常与 K-means、H-K-means 等方法配合使用。

4.0 地理人口特征分类实战 2

4.0 专题研究 2:纽约市tract分类

专题研究 2:纽约市内部地理人口特征分类

专题导航

纽约市内部地理人口特征分类

“从 ACS 变量到城市内部街区类型识别”

本专题以纽约市 census tract 尺度的 ACS 数据为基础,介绍城市内部地理人口特征分类的完整流程,包括变量选择与获取数据重构与预处理PCA 筛选H-K-means 聚类类型展示

点击下方卡片,快速跳转至相应模块:

💡 提示:学习完一个小节后,请再次点击屏幕右下角的专题导航按钮返回本导航页

研究课题
4.1 研究课题与分析目标

实战教学

本节将使用纽约市 census tract 尺度的 ACS 2024 5-year estimates 数据,开展一个面向城市内部社会空间分异地理人口特征分类实战。研究重点是在较细颗粒度的空间单元上,整合更多维度的人口、社会经济与居住背景变量,识别纽约市不同街区之间在社会人口结构上的潜在类型差异。

本节的方法流程将在前述基础上进一步细化。由于 tract 尺度下可纳入的变量更多,变量之间也更容易出现相关性,因此在正式聚类之前,将先进行相关性检验主成分分析(PCA),用于识别高相关变量、梳理变量结构,并完成必要的数据降维。随后,再结合 Kaiser rule(特征值大于 1) 初步筛选主成分,并参考各变量在主成分上的贡献情况,确定后续用于分类的变量。

在此基础上,本节将采用 H-K-means(Hierarchical K-means) 进行聚类分析,并结合多种方式对聚类数量k进行比较与判断。同时,还将对聚类结果进行基础质量检验,以评估分类结果的合理性与解释性。


示例课题

本节设置的研究课题为:基于纽约市 census tract 尺度的 ACS 2024 5-year estimates 数据,构建一个面向城市内部社会空间分异的地理人口特征分类框架,以识别纽约市不同街区在多维人口与社会经济属性上的潜在类型差异。

本节的分析目标主要包括以下几个方面:

  • 学习如何在纽约市 census tract 尺度上组织 ACS 多维变量,并构建适合城市内部分类的分析数据框;
  • 通过相关性检验PCA 对高维变量进行筛选与降维,减少变量冗余;
  • 结合 Kaiser rule变量贡献情况,确定后续用于聚类分析的变量;
  • 使用 H-K-means 完成聚类分析,并结合多种方法对聚类数量进行比较与判断;
  • 对聚类结果进行基础质量评估,并通过空间可视化展示纽约市内部不同类型街区的分布格局;
  • 在结果基础上,对不同类型进行简要解释,为后续城市内部社会空间分析提供分类框架。

本专题的几个核心方法

本专题将在地理人口特征分类流程中引入几个常见的统计与聚类工具

这些方法并不是彼此孤立的:相关性检验与 PCA 更多服务于变量筛选与降维,而 H-K-means 则用于后续的类型识别。三者结合,能够帮助构建一个更稳健的城市内部地理人口特征分类流程。

数据选择与获取
4.2 变量选择与数据获取

变量选择

在正式提取数据之前,需要先明确本节示例研究所使用的变量框架。与前一专题相比,本节将在纽约市 census tract 尺度下纳入更多维度的变量,以更细致地刻画城市内部的社会空间差异。考虑到本节仍以流程演示为主,而非构建一个最终定稿的完整分类体系,变量选择将遵循两个基本原则:一是尽量覆盖与邻里分异关系较强的核心维度;二是尽量优先选取 ACS 中较常用、解释性较强、后续便于重构的指标。

ACS 原始表格通常以计数值(count)为主,因此后续用于聚类的多数指标并不会直接以“百分比”形式下载,而是需要根据总量变量与分项变量进一步计算得到。例如,年龄占比、族裔占比、失业率、贫困人口比例或不同住房 tenure 的占比,通常都需要在提取表格后再自行重构。与此同时,像人口密度这类指标也往往需要结合人口总数与空间边界面积进一步计算。

变量选择

在正式提取数据之前,需要先明确本节示例研究所使用的变量框架。与前一专题相比,本节将在纽约市 census tract 尺度下纳入更多维度的变量,以更细致地刻画城市内部的社会空间差异。考虑到本节仍以流程演示为主,而非构建一个最终定稿的完整分类体系,变量选择将遵循两个基本原则:一是尽量覆盖与邻里分异关系较强的核心维度;二是尽量优先选取 ACS 中较常用、解释性较强、后续便于重构的指标。

ACS 原始表格通常以计数值(count)为主,因此后续用于聚类的多数指标并不会直接以“百分比”形式下载,而是需要根据总量变量与分项变量进一步计算得到。例如,年龄占比、族裔占比、失业率、贫困人口比例或不同住房占有方式的占比,通常都需要在提取表格后再自行重构。与此同时,像人口密度这类指标也往往需要结合人口总数与空间边界面积进一步计算。

本节示例研究的变量设计与 ACS 变量选取(纽约市 tract 尺度)
主题 领域 变量 分子变量 分母变量
人口结构 年龄 0–17岁人口占比 B01001_003, B01001_004, B01001_005, B01001_006, B01001_027, B01001_028, B01001_029, B01001_030 B01001_001
18–34岁人口占比 B01001_007, B01001_008, B01001_009, B01001_010, B01001_011, B01001_012, B01001_031, B01001_032, B01001_033, B01001_034, B01001_035, B01001_036 B01001_001
35–64岁人口占比 B01001_013, B01001_014, B01001_015, B01001_016, B01001_017, B01001_018, B01001_019, B01001_037, B01001_038, B01001_039, B01001_040, B01001_041, B01001_042, B01001_043 B01001_001
65岁及以上人口占比 B01001_020, B01001_021, B01001_022, B01001_023, B01001_024, B01001_025, B01001_044, B01001_045, B01001_046, B01001_047, B01001_048, B01001_049 B01001_001
性别 女性人口占比 B01001_026 B01001_001
族裔 非西班牙裔白人人口占比 B03002_003 B03002_001
非西班牙裔黑人人口占比 B03002_004 B03002_001
非西班牙裔亚裔人口占比 B03002_006 B03002_001
西班牙裔或拉丁裔人口占比 B03002_012 B03002_001
其他族裔人口占比 B03002_005, B03002_007, B03002_008, B03002_009, B03002_010, B03002_011 B03002_001
人口规模 人口密度 B01003_001 tract 面积
社会经济属性 婚姻状况 未婚人口占比 B12001_003, B12001_012 B12001_001
已婚或与配偶同居占比 B12001_004, B12001_013 B12001_001
分居人口占比 B12001_009, B12001_018 B12001_001
离异或丧偶人口占比 B12001_010, B12001_011, B12001_019, B12001_020 B12001_001
教育水平 低学历人口占比 B15003_002, B15003_003, B15003_004, B15003_005, B15003_006, B15003_007, B15003_008, B15003_009, B15003_010, B15003_011, B15003_012, B15003_013, B15003_014, B15003_015, B15003_016 B15003_001
中等学历人口占比 B15003_017, B15003_018, B15003_019, B15003_020, B15003_021 B15003_001
本科学历人口占比 B15003_022 B15003_001
研究生及以上学历人口占比 B15003_023, B15003_024, B15003_025 B15003_001
收入与贫困 家庭收入中位数 B19013_001 不适用
贫困人口比例 B17001_002 B17001_001
基尼系数 B19083_001 不适用
就业状况 失业率 B23025_005 B23025_003
非劳动力人口比例 B23025_007 B23025_001
管理与专业职业人口占比 C24010_003, C24010_039 C24010_001
服务业职业人口占比 C24010_019, C24010_055 C24010_001
销售与办公室职业人口占比 C24010_027, C24010_063 C24010_001
自然资源、建筑与维护职业人口占比 C24010_030, C24010_066 C24010_001
生产运输职业人口占比 C24010_034, C24010_070 C24010_001
住房与居住 住房占有方式 自有住房占比 B25003_002 B25003_001
租住房屋占比 B25003_003 B25003_001
租金 月租中位数 B25064_001 不适用
房屋类型 独立式住宅占比 B25024_002 B25024_001
联排或相连住宅占比 B25024_003 B25024_001
多户住宅占比 B25024_004, B25024_005, B25024_006, B25024_007, B25024_008, B25024_009 B25024_001
通勤与可达性 通勤方式 小汽车通勤占比 B08301_003, B08301_004, B08301_005 B08301_001
公共交通通勤占比 B08301_010 B08301_001
步行通勤占比 B08301_019 B08301_001
居家办公占比 B08301_021 B08301_001
通勤时间 小于5分钟通勤占比 B08303_002 B08303_001
5–14分钟通勤占比 B08303_003, B08303_004 B08303_001
15–29分钟通勤占比 B08303_005, B08303_006, B08303_007 B08303_001
30–59分钟通勤占比 B08303_008, B08303_009, B08303_010, B08303_011 B08303_001
60分钟及以上通勤占比 B08303_012, B08303_013 B08303_001
注:比例型指标一般使用分子变量汇总值除以分母变量。家庭收入中位数、月租中位数和基尼系数为 ACS 已提供的直接估计值。人口密度需将 B01003_001 与 tract 面积数据结合计算。教育水平变量以 25 岁及以上人口为统计基准;婚姻状况变量以 15 岁及以上人口为统计基准;就业变量以 16 岁及以上人口为统计基准;职业结构变量以 16 岁及以上就业人口为统计基准。

变量选择说明

本节将使用美国 ACS 2020–2024 5-year estimates 作为数据来源。由于本专题纳入的变量维度较前一专题更丰富,因此所涉及的 ACS 详细表也相应增加。

tidycensus 在实际操作中通常通过 get_acs() 提取表格数据,通过 load_variables() 查询表号与字段定义,因此本节后续也将沿用这一工作流。

本节主要涉及以下几张 ACS 常用详细表:

  • B01001
    Sex by Age,用于构造不同年龄段比例与性别结构;
  • B03002
    Hispanic or Latino Origin by Race,用于构造非西班牙裔白人、非西班牙裔黑人、非西班牙裔亚裔、西班牙裔或拉丁裔及其他族裔等结构变量;
  • B01003
    Total Population,用于提取总人口,并进一步计算人口密度;
  • B12001
    Sex by Marital Status by Age,用于构造未婚、已婚、离异或分居、丧偶等婚姻状况指标;
  • B15003
    Educational Attainment for the Population 25 Years and Over,用于构造不同教育层级的人口比例;
  • B19013
    Median Household Income in the Past 12 Months,用于表示家庭收入中位数;
  • B17001
    Poverty Status in the Past 12 Months by Sex by Age,用于构造贫困人口比例;
  • B19083
    Gini Index of Income Inequality,用于表示收入不平等程度;
  • B23025
    Employment Status,用于构造失业率;
  • C24010
    Sex by Occupation for the Civilian Employed Population 16 Years and Over,用于构造不同职业类型人口比例;
  • B25003
    Tenure,用于构造自有住房与租住房屋占比;
  • B25064
    Median Gross Rent,用于表示月租中位数;
  • B25024
    Units in Structure,用于构造独立式住宅、联排或相连住宅及多户住宅等房屋类型指标;
  • B08301
    Means of Transportation to Work,用于构造公共交通通勤与小汽车通勤比例;
  • B08303
    Travel Time to Work,用于表示通勤时间相关指标。

从变量组织方式看,这些表大致对应人口结构社会经济属性住房与居住以及通勤与可达性四个主题。后续分析中,并非所有变量都会原样进入聚类模型,而是还需结合相关性检验主成分分析与研究目标进一步筛选。


实战代码 环境配置

# ==============================================================================
# 1. 环境构建与依赖项声明
# ==============================================================================
library(tidycensus)   # 获取美国人口普查 / ACS 数据
library(tidyverse)    # 数据整理与可视化
library(sf)           # 处理空间矢量数据

# 可视化与地图
library(classInt)     # Fisher 断点
library(viridis)      # viridis 配色
library(ggspatial)    # 比例尺、指北针、经纬线
library(scales)       # 标签格式化
library(patchwork)    # 拼接多个图形

# 降维、相关性与聚类
library(corrplot)     # 相关系数矩阵可视化
library(FactoMineR)   # PCA 分析
library(factoextra)   # PCA 与聚类结果可视化
library(cluster)      # 聚类质量评估
library(NbClust)      # 辅助判断聚类数量
library(fpc)          # 聚类稳定性与质量检验

# 表格输出
library(knitr)
library(kableExtra)

# ==============================================================================
# 2. Census API 认证设置
# ==============================================================================
# 注:
# - 首次使用 tidycensus 前,需先申请美国 Census API key
# - 申请地址:https://api.census.gov/data/key_signup.html
# - 首次配置时取消下行注释,并替换为自己的 key;运行一次后即可写入本地环境
# census_api_key("YOUR_CENSUS_API_KEY", install = TRUE, overwrite = TRUE)

# 若已完成 key 安装,建议重启 R 会话后再继续运行

# ==============================================================================
# 3. 参数设定:年份、调查类型与空间尺度
# ==============================================================================
acs_year   <- 2024      # 对应 ACS 2020–2024 5-year estimates 的发布年
acs_survey <- "acs5"    # 5-year estimates
geo_unit   <- "tract"   # 研究尺度:census tract

# 纽约市五个 borough 对应的 county 名称
nyc_counties <- c("Bronx", "Kings", "New York", "Queens", "Richmond")

数据获取

# ==============================================================================
# 4. ACS 变量字典检索
# ==============================================================================
# 目的:
# - 预先读取 ACS 变量字典,便于核对变量名称与字段含义
# - 后续可结合关键词搜索,确认表号与字段定义

acs_vars_2024 <- load_variables(
  year = acs_year,
  dataset = acs_survey,
  cache = TRUE
)

# 快速查看本专题涉及的核心表号
acs_vars_2024 %>%
  filter(str_detect(name, "^B01001|^B03002|^B01003|^B12001|^B15003|^B19013|^B17001|^B19083|^B23025|^C24010|^B25003|^B25064|^B25024|^B08301|^B08303")) %>%
  select(name, label, concept) %>%
  print(n = 3)
# A tibble: 1,268 × 3
  name        label                                  concept                 
  <chr>       <chr>                                  <chr>                   
1 B01001A_001 Estimate!!Total:                       Sex by Age (White Alone)
2 B01001A_002 Estimate!!Total:!!Male:                Sex by Age (White Alone)
3 B01001A_003 Estimate!!Total:!!Male:!!Under 5 years Sex by Age (White Alone)
# ℹ 1,265 more rows
# ==============================================================================
# 5. 研究所需原始变量清单
# ==============================================================================
# 说明:
# - 本节先完成原始 ACS 表项提取
# - 后续再基于这些字段计算比例变量、密度变量与结构性指标
# - 这里优先保留各主题下便于后续重构的底层变量

acs_var_dict_nyc <- c(
  # ------------------------------------------------------------
  # A. 年龄与性别(B01001: Sex by Age)
  # ------------------------------------------------------------
  total_pop      = "B01001_001",
  male_total     = "B01001_002",
  female_total   = "B01001_026",

  # 0–17 岁
  male_0_17_1    = "B01001_003",
  male_0_17_2    = "B01001_004",
  male_0_17_3    = "B01001_005",
  male_0_17_4    = "B01001_006",
  female_0_17_1  = "B01001_027",
  female_0_17_2  = "B01001_028",
  female_0_17_3  = "B01001_029",
  female_0_17_4  = "B01001_030",

  # 18–34 岁
  male_18_34_1   = "B01001_007",
  male_18_34_2   = "B01001_008",
  male_18_34_3   = "B01001_009",
  male_18_34_4   = "B01001_010",
  male_18_34_5   = "B01001_011",
  male_18_34_6   = "B01001_012",
  female_18_34_1 = "B01001_031",
  female_18_34_2 = "B01001_032",
  female_18_34_3 = "B01001_033",
  female_18_34_4 = "B01001_034",
  female_18_34_5 = "B01001_035",
  female_18_34_6 = "B01001_036",

  # 35–64 岁
  male_35_64_1   = "B01001_013",
  male_35_64_2   = "B01001_014",
  male_35_64_3   = "B01001_015",
  male_35_64_4   = "B01001_016",
  male_35_64_5   = "B01001_017",
  male_35_64_6   = "B01001_018",
  male_35_64_7   = "B01001_019",
  female_35_64_1 = "B01001_037",
  female_35_64_2 = "B01001_038",
  female_35_64_3 = "B01001_039",
  female_35_64_4 = "B01001_040",
  female_35_64_5 = "B01001_041",
  female_35_64_6 = "B01001_042",
  female_35_64_7 = "B01001_043",

  # 65 岁及以上
  male_65p_1     = "B01001_020",
  male_65p_2     = "B01001_021",
  male_65p_3     = "B01001_022",
  male_65p_4     = "B01001_023",
  male_65p_5     = "B01001_024",
  male_65p_6     = "B01001_025",
  female_65p_1   = "B01001_044",
  female_65p_2   = "B01001_045",
  female_65p_3   = "B01001_046",
  female_65p_4   = "B01001_047",
  female_65p_5   = "B01001_048",
  female_65p_6   = "B01001_049",

  # ------------------------------------------------------------
  # B. 族裔结构(B03002: Hispanic or Latino Origin by Race)
  # ------------------------------------------------------------
  race_eth_total = "B03002_001",
  non_hisp_total = "B03002_002",
  nh_white       = "B03002_003",
  nh_black       = "B03002_004",
  nh_asian       = "B03002_006",
  hispanic_total = "B03002_012",
  other_eth_1    = "B03002_005",
  other_eth_2    = "B03002_007",
  other_eth_3    = "B03002_008",
  other_eth_4    = "B03002_009",

  # ------------------------------------------------------------
  # C. 总人口、收入与不平等
  # ------------------------------------------------------------
  total_pop_alt  = "B01003_001",
  med_hh_income  = "B19013_001",
  gini_index     = "B19083_001",

  # ------------------------------------------------------------
  # D. 婚姻状况(B12001: Sex by Marital Status by Age)
  # ------------------------------------------------------------
  marital_total        = "B12001_001",
  male_never_married   = "B12001_003",
  male_now_married     = "B12001_004",
  male_separated       = "B12001_007",
  male_widowed         = "B12001_009",
  male_divorced        = "B12001_010",
  female_never_married = "B12001_012",
  female_now_married   = "B12001_013",
  female_separated     = "B12001_016",
  female_widowed       = "B12001_018",
  female_divorced      = "B12001_019",

  # ------------------------------------------------------------
  # E. 教育水平(B15003: Educational Attainment for Population 25 Years and Over)
  # ------------------------------------------------------------
  edu_total_25p   = "B15003_001",
  no_school       = "B15003_002",
  nursery         = "B15003_003",
  kindergarten    = "B15003_004",
  grade_1         = "B15003_005",
  grade_2         = "B15003_006",
  grade_3         = "B15003_007",
  grade_4         = "B15003_008",
  grade_5         = "B15003_009",
  grade_6         = "B15003_010",
  grade_7         = "B15003_011",
  grade_8         = "B15003_012",
  grade_9         = "B15003_013",
  grade_10        = "B15003_014",
  grade_11        = "B15003_015",
  grade_12_no_dip = "B15003_016",
  hs_diploma      = "B15003_017",
  ged             = "B15003_018",
  some_college_1  = "B15003_019",
  some_college_2  = "B15003_020",
  associate       = "B15003_021",
  bachelor        = "B15003_022",
  master          = "B15003_023",
  professional    = "B15003_024",
  doctorate       = "B15003_025",

  # ------------------------------------------------------------
  # F. 贫困(B17001: Poverty Status in the Past 12 Months by Sex by Age)
  # ------------------------------------------------------------
  poverty_total        = "B17001_001",
  below_poverty_total  = "B17001_002",

  # ------------------------------------------------------------
  # G. 就业(B23025: Employment Status)
  # ------------------------------------------------------------
  labor_force    = "B23025_003",
  employed       = "B23025_004",
  unemployed     = "B23025_005",
  not_in_lf      = "B23025_007",

  # ------------------------------------------------------------
  # H. 职业结构(C24010: Sex by Occupation)
  # ------------------------------------------------------------
  occ_total            = "C24010_001",
  male_mgmt_prof       = "C24010_003",
  female_mgmt_prof     = "C24010_039",
  male_service         = "C24010_019",
  female_service       = "C24010_055",
  male_sales_office    = "C24010_027",
  female_sales_office  = "C24010_063",
  male_manual          = "C24010_030",
  female_manual        = "C24010_066",
  male_prod_trans      = "C24010_034",
  female_prod_trans    = "C24010_070",

  # ------------------------------------------------------------
  # I. 住房占有方式与租金(B25003 / B25064)
  # ------------------------------------------------------------
  owner_renter_total = "B25003_001",
  owner_occupied     = "B25003_002",
  renter_occupied    = "B25003_003",
  med_gross_rent     = "B25064_001",

  # ------------------------------------------------------------
  # J. 房屋类型(B25024: Units in Structure)
  # ------------------------------------------------------------
  units_total        = "B25024_001",
  detached_1unit     = "B25024_002",
  attached_1unit     = "B25024_003",
  units_2            = "B25024_004",
  units_3_4          = "B25024_005",
  units_5_9          = "B25024_006",
  units_10_19        = "B25024_007",
  units_20_49        = "B25024_008",
  units_50_plus      = "B25024_009",
  mobile_home        = "B25024_010",
  boat_rv_van        = "B25024_011",

  # ------------------------------------------------------------
  # K. 通勤方式(B08301: Means of Transportation to Work)
  # ------------------------------------------------------------
  commute_total      = "B08301_001",
  commute_car_total  = "B08301_002",
  commute_drove      = "B08301_003",
  commute_carpool    = "B08301_004",
  commute_transit    = "B08301_010",
  commute_bus        = "B08301_011",
  commute_subway     = "B08301_012",
  commute_rail       = "B08301_013",
  commute_taxi_ride  = "B08301_016",
  commute_bicycle    = "B08301_018",
  commute_walk       = "B08301_019",
  commute_home       = "B08301_021",

  # ------------------------------------------------------------
  # L. 通勤时间(B08303: Travel Time to Work)
  # ------------------------------------------------------------
  travel_time_total  = "B08303_001",
  tt_lt5             = "B08303_002",
  tt_5_9             = "B08303_003",
  tt_10_14           = "B08303_004",
  tt_15_19           = "B08303_005",
  tt_20_24           = "B08303_006",
  tt_25_29           = "B08303_007",
  tt_30_34           = "B08303_008",
  tt_35_39           = "B08303_009",
  tt_40_44           = "B08303_010",
  tt_45_59           = "B08303_011",
  tt_60_89           = "B08303_012",
  tt_90_plus         = "B08303_013"
)

# 查看变量数量
length(acs_var_dict_nyc)
[1] 155

下面这一部分可能会因为网速波动中断卡顿, 可以多试几次 注:由于下载数据量较多可能会失败,建议多试几次并存入本地文件中

# ==============================================================================
# 6. 纽约市 tract 尺度 ACS 原始数据提取(附带边界)
# ==============================================================================
# 说明:
# - 若本地已存在缓存文件,则直接读取
# - 若不存在,则调用 Census API 下载,并保存至 data/ 文件夹
# - geometry = TRUE:同步下载 tract 边界,便于后续空间制图
# - output = "wide":每个变量展开为单独列,便于后续重构

dir.create("data", showWarnings = FALSE)

nyc_raw_file <- "data/acs_nyc_tract_raw_2024.rds"

# 纽约市五个 borough 对应的 county 名称
nyc_counties <- c("Bronx", "Kings", "New York", "Queens", "Richmond")

if (file.exists(nyc_raw_file)) {
  message("检测到本地缓存文件,直接读取纽约市 tract 原始数据...")
  acs_nyc_raw <- readRDS(nyc_raw_file)
} else {
  message("本地无缓存文件,开始从 Census API 下载纽约市 tract 原始数据...")
  
  acs_nyc_raw <- get_acs(
    geography = geo_unit,
    variables = acs_var_dict_nyc,
    state = "NY",
    county = nyc_counties,
    year = acs_year,
    survey = acs_survey,
    geometry = TRUE,
    output = "wide"
  )
  
  saveRDS(acs_nyc_raw, nyc_raw_file)
  message("纽约市 tract 原始数据已保存至:", nyc_raw_file)
}

# 简要查看结果规模
dim(acs_nyc_raw)
[1] 2327  313

原始数据展示制图

使用第四章空间制图方法,简单查看一下获取的数据(空间以及非空间数据)

这里我们选择家庭收入中位数 + Gini系数作为展示出图

# ==============================================================================
# 9. 纽约市 tract 尺度原始变量可视化
# ==============================================================================

# ------------------------------------------------------------------------------
# 9.1 选择两个原始变量:家庭收入中位数 与 基尼系数
# ------------------------------------------------------------------------------
map_nyc_raw <- acs_nyc_raw %>%
  select(GEOID, NAME, med_hh_incomeE, gini_indexE, geometry) %>%
  filter(!is.na(med_hh_incomeE), !is.na(gini_indexE))

# ------------------------------------------------------------------------------
# 9.2 构建经纬线对象
# ------------------------------------------------------------------------------
graticule_sf <- st_graticule(map_nyc_raw)

# ------------------------------------------------------------------------------
# 9.3 自定义 Fisher 分级制图函数
# ------------------------------------------------------------------------------
plot_fisher_map <- function(data, var, title_text, legend_title, palette_option, is_income = FALSE) {
  
  # 提取变量
  values <- st_drop_geometry(data)[[var]]
  
  # Fisher 断点
  fisher_brks <- classIntervals(
    var = values,
    n = 6,
    style = "fisher"
  )
  
  brks <- fisher_brks$brks
  
  # 图例标签
  if (is_income) {
    labels <- paste0(
      "$", round(brks[-length(brks)] / 1000, 0), "k",
      " – ",
      "$", round(brks[-1] / 1000, 0), "k"
    )
  } else {
    labels <- paste0(
      round(brks[-length(brks)], 2),
      " – ",
      round(brks[-1], 2)
    )
  }
  
  # 分级
  plot_data <- data %>%
    mutate(
      class = cut(
        .data[[var]],
        breaks = brks,
        include.lowest = TRUE,
        labels = labels
      )
    )
  
  # 出图
  ggplot() +
    geom_sf(
      data = graticule_sf,
      color = "grey75",
      linewidth = 0.2,
      linetype = "dashed"
    ) +
    geom_sf(
      data = plot_data,
      aes(fill = class),
      color = NA
    ) +
    annotation_scale(
      location = "bl",
      width_hint = 0.25,
      text_cex = 0.75
    ) +
    annotation_north_arrow(
      location = "br",
      which_north = "true",
      style = north_arrow_fancy_orienteering,
      height = unit(0.85, "cm"),
      width  = unit(0.85, "cm")
    ) +
    scale_fill_viridis_d(
      option = palette_option,
      direction = -1,
      guide = guide_legend(reverse = TRUE),
      name = legend_title
    ) +
    coord_sf(crs = st_crs(4326)) +
    labs(title = title_text) +
    theme_minimal(base_size = 11) +
    theme(
      legend.position = "bottom",
      panel.grid.major = element_blank(),
      panel.background = element_rect(fill = "white", color = NA),
      plot.title = element_text(face = "bold", size = 12),
      legend.title = element_text(size = 9),
      legend.text = element_text(size = 8),
      axis.title = element_blank()
    )
}

# ------------------------------------------------------------------------------
# 9.4 分别绘制两个变量
# ------------------------------------------------------------------------------
p_income_nyc <- plot_fisher_map(
  data = map_nyc_raw,
  var = "med_hh_incomeE",
  title_text = "Median household income",
  legend_title = "Income (1,000 USD)",
  palette_option = "mako",
  is_income = TRUE
)

p_gini_nyc <- plot_fisher_map(
  data = map_nyc_raw,
  var = "gini_indexE",
  title_text = "Gini index",
  legend_title = "Gini index",
  palette_option = "plasma",
  is_income = FALSE
)

# ------------------------------------------------------------------------------
# 9.5 拼接输出
# ------------------------------------------------------------------------------
p_income_nyc | p_gini_nyc

什么是基尼系数(Gini index)?

基尼系数(Gini index)是衡量一个地区收入分配不平等程度的常用指标,取值通常介于 01 之间。其基本含义可以简单理解为:

  • 当数值越接近 0 时,表示该地区内部的收入分配越平均;
  • 当数值越接近 1 时,表示该地区内部的收入差距越大,不平等程度越高。

因此,基尼系数并不直接表示“这个地区是富裕还是贫困”,而是反映该地区内部不同人群之间的收入差距有多大。一个地区即使整体收入水平较高,也仍然可能具有较高的基尼系数;反之,一个整体收入并不高的地区,也可能因为内部收入差异较小而表现出较低的基尼系数。

结合上图可以看到,曼哈顿部分 tract 的基尼系数明显偏高,颜色也更深。这通常说明这些地区内部的收入分化更为显著,即高收入群体与低收入群体可能在同一空间范围内并存,整体呈现出更强的社会经济不均衡特征。换言之,曼哈顿中部分街区并不一定只是“收入高”,而更可能是高收入与低收入人群并置、内部差距较大的区域。

下一节将进行变量的重构以及预处理!

数据重构与预处理
4.3 数据重构与预处理

数据重构

前一节获取的是基于 ACS 原始表项提取的纽约市 tract 尺度数据。这些数据虽然已经完成下载,并具有明确的空间边界信息,但仍然主要属于底层统计字段,还不能直接用于后续的地理人口特征分类分析 (参考前一专题研究)

ACS 的多数表格以计数值(count)为主,例如某一年龄组人口数、某一婚姻状态人数、某一教育层级人数或某一通勤方式人数等。因此,在正式进入分析之前,通常需要先进行变量重构

所谓变量重构,是指根据研究目标,将原始表项进一步转换为更适合横向比较的分析指标。最常见的方式包括比例变量密度变量与少量可直接使用的连续变量。例如,年龄结构、婚姻状况、族裔构成、教育层级、住房占有方式与通勤方式等,通常都更适合表示为占比,而不是直接使用原始人数;家庭收入中位数、基尼系数与月租中位数等,则可以直接作为连续变量使用。

在这些转换中,比例变量最为常见,其基本形式可表示为:

\[ \text{比例} = \frac{\text{某一类别人数}}{\text{对应总体人数}} \]

例如,65 岁及以上人口占比可以由 65 岁及以上人口数除以总人口得到;租住房屋占比可以由租住房屋数除以全部住房单元数得到;公共交通通勤占比则可由公共交通通勤人数除以上班通勤总人数得到。通过这种方式,不同 tract 之间的结构差异可以在相同尺度上进行比较,而不会过度受到人口总量差异的影响。

除比例变量外,人口密度也是城市内部社会空间分析中较常见的指标。其基本形式为:

\[ \text{人口密度} = \frac{\text{总人口}}{\text{地理单元面积}} \]

在实际操作中,若面积以平方公里计算,则人口密度通常可表达为每平方公里人口数。这一指标并不是 ACS 直接提供的,而需要结合人口总数tract 边界面积进一步计算。因此,人口密度可以看作是一类由非空间属性空间边界信息共同派生出来的变量。

注意!计算面积的时候需要先查看空间数据的 CRS是否是投影坐标系!并且明确单位


实战代码 变量重构

# ==============================================================================
# 10. NYC tract 尺度分析变量重构
# ==============================================================================
# 说明:
# - 基于 ACS 原始表项,重构后续分析所需的比例变量、密度变量与连续变量
# - 需要安装 units包

# ------------------------------------------------------------------------------
# 10.1 先计算 tract 面积(平方公里)
# ------------------------------------------------------------------------------
# 注:
# - 原始数据坐标系为 NAD83
# - 面积计算时先投影到 NAD83 / New York Long Island(EPSG:6539)
# - 该坐标系以英尺为单位,因此需先转平方米,再换算为平方公里
acs_nyc_vars <- acs_nyc_raw %>%
  st_transform(6539) %>%
  mutate(
    tract_area_km2 = units::set_units(st_area(geometry), km^2) %>% as.numeric()
  ) %>%
  st_transform(4269) %>%
  mutate(
    # --------------------------------------------------------------------------
    # A. 年龄结构
    # --------------------------------------------------------------------------
    age_0_17 = male_0_17_1E + male_0_17_2E + male_0_17_3E + male_0_17_4E +
      female_0_17_1E + female_0_17_2E + female_0_17_3E + female_0_17_4E,
    
    age_18_34 = male_18_34_1E + male_18_34_2E + male_18_34_3E +
      male_18_34_4E + male_18_34_5E + male_18_34_6E +
      female_18_34_1E + female_18_34_2E + female_18_34_3E +
      female_18_34_4E + female_18_34_5E + female_18_34_6E,
    
    age_35_64 = male_35_64_1E + male_35_64_2E + male_35_64_3E +
      male_35_64_4E + male_35_64_5E + male_35_64_6E + male_35_64_7E +
      female_35_64_1E + female_35_64_2E + female_35_64_3E +
      female_35_64_4E + female_35_64_5E + female_35_64_6E + female_35_64_7E,
    
    age_65p = male_65p_1E + male_65p_2E + male_65p_3E +
      male_65p_4E + male_65p_5E + male_65p_6E +
      female_65p_1E + female_65p_2E + female_65p_3E +
      female_65p_4E + female_65p_5E + female_65p_6E,
    
    pct_age_0_17  = age_0_17 / total_popE,
    pct_age_18_34 = age_18_34 / total_popE,
    pct_age_35_64 = age_35_64 / total_popE,
    pct_age_65p   = age_65p / total_popE,
    
    # --------------------------------------------------------------------------
    # B. 性别结构
    # --------------------------------------------------------------------------
    pct_female = female_totalE / total_popE,
    
    # --------------------------------------------------------------------------
    # C. 族裔结构
    # --------------------------------------------------------------------------
    other_eth = other_eth_1E + other_eth_2E + other_eth_3E + other_eth_4E,
    
    pct_nh_white = nh_whiteE / race_eth_totalE,
    pct_nh_black = nh_blackE / race_eth_totalE,
    pct_nh_asian = nh_asianE / race_eth_totalE,
    pct_hispanic = hispanic_totalE / race_eth_totalE,
    pct_other_eth = other_eth / race_eth_totalE,
    
    # --------------------------------------------------------------------------
    # D. 人口密度
    # --------------------------------------------------------------------------
    pop_density = total_popE / tract_area_km2,
    
    # --------------------------------------------------------------------------
    # E. 婚姻状况
    # --------------------------------------------------------------------------
    never_married = male_never_marriedE + female_never_marriedE,
    married       = male_now_marriedE + female_now_marriedE,
    separated     = male_separatedE + female_separatedE,
    widowed       = male_widowedE + female_widowedE,
    divorced      = male_divorcedE + female_divorcedE,
    div_widowed   = divorced + widowed,
    
    pct_never_married = never_married / marital_totalE,
    pct_married       = married / marital_totalE,
    pct_separated     = separated / marital_totalE,
    pct_div_widowed   = div_widowed / marital_totalE,
    
    # --------------------------------------------------------------------------
    # F. 教育水平
    # --------------------------------------------------------------------------
    edu_low = no_schoolE + nurseryE + kindergartenE + grade_1E + grade_2E +
      grade_3E + grade_4E + grade_5E + grade_6E + grade_7E +
      grade_8E + grade_9E + grade_10E + grade_11E + grade_12_no_dipE,
    
    edu_secondary = hs_diplomaE + gedE,
    edu_bachelor  = bachelorE,
    edu_gradplus  = masterE + professionalE + doctorateE,
    
    pct_edu_low       = edu_low / edu_total_25pE,
    pct_edu_secondary = edu_secondary / edu_total_25pE,
    pct_edu_bachelor  = edu_bachelor / edu_total_25pE,
    pct_edu_gradplus  = edu_gradplus / edu_total_25pE,
    
    # --------------------------------------------------------------------------
    # G. 收入、贫困与不平等
    # --------------------------------------------------------------------------
    med_hh_income = med_hh_incomeE,
    gini_index    = gini_indexE,
    poverty_rate  = below_poverty_totalE / poverty_totalE,
    
    # --------------------------------------------------------------------------
    # H. 就业与职业结构
    # --------------------------------------------------------------------------
    unemployment_rate = unemployedE / labor_forceE,
    not_in_lf_rate    = not_in_lfE / (labor_forceE + not_in_lfE),
    
    mgmt_prof  = male_mgmt_profE + female_mgmt_profE,
    service    = male_serviceE + female_serviceE,
    sales_off  = male_sales_officeE + female_sales_officeE,
    manual     = male_manualE + female_manualE,
    prod_trans = male_prod_transE + female_prod_transE,
    
    pct_mgmt_prof  = mgmt_prof / occ_totalE,
    pct_service    = service / occ_totalE,
    pct_sales_off  = sales_off / occ_totalE,
    pct_manual     = manual / occ_totalE,
    pct_prod_trans = prod_trans / occ_totalE,
    
    # --------------------------------------------------------------------------
    # I. 住房与租金
    # --------------------------------------------------------------------------
    pct_owner_occ  = owner_occupiedE / owner_renter_totalE,
    pct_renter_occ = renter_occupiedE / owner_renter_totalE,
    med_gross_rent = med_gross_rentE,
    
    # --------------------------------------------------------------------------
    # J. 房屋类型
    # --------------------------------------------------------------------------
    multi_unit = units_2E + units_3_4E + units_5_9E +
      units_10_19E + units_20_49E + units_50_plusE,
    
    pct_detached_1unit = detached_1unitE / units_totalE,
    pct_attached_1unit = attached_1unitE / units_totalE,
    pct_multi_unit     = multi_unit / units_totalE,
    
    # --------------------------------------------------------------------------
    # K. 通勤方式
    # --------------------------------------------------------------------------
    pct_commute_car     = commute_car_totalE / commute_totalE,
    pct_commute_transit = commute_transitE / commute_totalE,
    pct_commute_walk    = commute_walkE / commute_totalE,
    pct_commute_home    = commute_homeE / commute_totalE,
    
    # --------------------------------------------------------------------------
    # L. 通勤时间
    # --------------------------------------------------------------------------
    tt_5_14   = tt_5_9E + tt_10_14E,
    tt_15_29  = tt_15_19E + tt_20_24E + tt_25_29E,
    tt_30_59  = tt_30_34E + tt_35_39E + tt_40_44E + tt_45_59E,
    tt_60plus = tt_60_89E + tt_90_plusE,
    
    pct_tt_lt5    = tt_lt5E / travel_time_totalE,
    pct_tt_5_14   = tt_5_14 / travel_time_totalE,
    pct_tt_15_29  = tt_15_29 / travel_time_totalE,
    pct_tt_30_59  = tt_30_59 / travel_time_totalE,
    pct_tt_60plus = tt_60plus / travel_time_totalE
  )
# ==============================================================================
# 11.  分析变量数据框整理
# ==============================================================================
acs_nyc_analysis <- acs_nyc_vars %>%
  select(
    GEOID, NAME, geometry,

    # 人口结构
    pct_age_0_17, pct_age_18_34, pct_age_35_64, pct_age_65p,
    pct_female,
    pct_nh_white, pct_nh_black, pct_nh_asian, pct_hispanic, pct_other_eth,
    pop_density,

    # 婚姻
    pct_never_married, pct_married, pct_separated, pct_div_widowed,

    # 教育
    pct_edu_low, pct_edu_secondary, pct_edu_bachelor, pct_edu_gradplus,

    # 收入与不平等
    med_hh_income, poverty_rate, gini_index,

    # 就业与职业
    unemployment_rate, not_in_lf_rate,
    pct_mgmt_prof, pct_service, pct_sales_off, pct_manual, pct_prod_trans,

    # 住房
    pct_owner_occ, pct_renter_occ, med_gross_rent,
    pct_detached_1unit, pct_attached_1unit, pct_multi_unit,

    # 通勤
    pct_commute_car, pct_commute_transit, pct_commute_walk, pct_commute_home,
    pct_tt_lt5, pct_tt_5_14, pct_tt_15_29, pct_tt_30_59, pct_tt_60plus
  ) %>%
  filter(
    if_all(
      where(is.numeric),
      ~ !is.na(.x) & is.finite(.x)
    )
  )

dim(acs_nyc_analysis) # 2170 census tract, 44+3个变量
[1] 2170   47

相关性检验

在进入 PCA 与后续聚类分析之前,通常可以先对变量之间的相关性作一个初步检查。一方面,相关系数矩阵有助于识别高度相关的变量组合,帮助发现潜在的变量冗余;另一方面,也能够帮助理解不同社会人口特征之间是否存在较明显的联动关系。这里先使用 Pearson 相关系数(Pearson correlation) 对当前分析变量进行初步检验,并通过相关性图与高相关变量对列表进行展示。

# ==============================================================================
# 12. 相关性检验:提取数值变量
# ==============================================================================
# 说明:
# - 这里只保留数值变量进行 Pearson 相关系数计算
# - 空间几何列与标识列不参与相关性检验

nyc_corr_vars <- acs_nyc_analysis %>%
  st_drop_geometry() %>%
  select(-GEOID, -NAME)

# ==============================================================================
# 13. 相关系数矩阵可视化
# ==============================================================================
# 说明:
# - 使用 Pearson 相关系数
# - 相关性图采用上三角布局,突出变量之间的线性关联模式

# 计算相关系数矩阵
corr_mat <- cor(
  nyc_corr_vars,
  use = "pairwise.complete.obs",
  method = "pearson"
)

# 绘制相关性图 可以使用 ?corrplot()查看可调参数
corrplot(
  corr_mat,
  method = "square",
  type = "lower",
  tl.col = "black",
  tl.cex = 0.68,
  tl.srt = 90,
  addCoef.col = NULL,
  col = colorRampPalette(c("#2166AC", "white", "#B2182B"))(250),
  diag = TRUE,
  outline = FALSE,
  cl.cex = 0.8,
  number.cex = 0.55,
  mar = c(0, 0, 1, 0),
  title = "Correlation matrix of reconstructed NYC tract variables"
)

# ==============================================================================
# 14.  提取相关性最高的变量对
# ==============================================================================
# 说明:
# - 提取 Pearson 相关系数绝对值最高的前 5 组变量
# - 仅保留变量对的唯一组合,避免重复
# - 这里同时展示相关系数原值与其绝对值,便于判断方向与强度

# ------------------------------------------------------------------------------
# 14.1  将相关系数矩阵转换为长表
# ------------------------------------------------------------------------------
top_corr_pairs <- as.data.frame(as.table(corr_mat)) %>%
  rename(
    var1 = Var1,
    var2 = Var2,
    correlation = Freq
  ) %>%
  filter(var1 != var2) %>%
  rowwise() %>%
  mutate(
    pair_id = paste(sort(c(as.character(var1), as.character(var2))), collapse = " | ")
  ) %>%
  ungroup() %>%
  distinct(pair_id, .keep_all = TRUE) %>%
  mutate(
    abs_correlation = abs(correlation)
  ) %>%
  arrange(desc(abs_correlation)) %>%
  slice_head(n = 5) %>%
  mutate(
    correlation = round(correlation, 3),
    abs_correlation = round(abs_correlation, 3)
  ) %>%
  select(
    变量1 = var1,
    变量2 = var2,
    相关系数 = correlation,
    相关系数绝对值 = abs_correlation
  )

# ------------------------------------------------------------------------------
# 14.2 以表格形式展示
# ------------------------------------------------------------------------------
kable(
  top_corr_pairs,
  align = "c",
  caption = "相关性最高的前 5 组变量"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
相关性最高的前 5 组变量
变量1 变量2 相关系数 相关系数绝对值
pct_renter_occ pct_owner_occ -1.000 1.000
pct_married pct_never_married -0.905 0.905
pct_multi_unit pct_detached_1unit -0.875 0.875
pct_mgmt_prof pct_edu_gradplus 0.843 0.843
pct_service pct_mgmt_prof -0.821 0.821

高相关变量初筛

在完成相关性检验之后,可以进一步对高度相关的变量作一个初步筛选。这样做的目的,并不是单纯为了减少变量数量,而是为了避免若干几乎表达同一结构维度的变量在后续 PCA 与聚类分析中被重复计入。

对于绝对相关系数特别高的变量对,本节优先保留解释性更强更符合纽约市邻里研究语境更常用于社会空间分析的变量,而将与其高度重复或互补的变量暂时移除。例如,自有住房占比与租住房屋占比本身就是一种近乎完全互补的关系,因此通常没有必要同时保留;类似地,独立式住宅占比与多户住宅占比、已婚人口占比与未婚人口占比,在城市内部分析中也往往体现出较强的结构对立。

pct_renter_occpct_owner_occ 相关系数为 -1; pct_marriedpct_never_married 相关系数约为 -0.91 等等;这里我们仅取变量对中的一个变量进行保留

注意:高相关性变量对中选取哪一个留下或剔除?可以询问AI;参考关键词 最小生成树!与 变量解释性

# ==============================================================================
# 15.  基于高相关变量对进行初步筛选
# ==============================================================================
# 说明:
# - 对绝对相关系数较高的变量对进行人工筛选
# - 优先保留解释性更强、在纽约市社会空间研究中更常用的变量
# - 这里先移除几组明显冗余或互补的变量:
#   1) pct_owner_occ      (与 pct_renter_occ 几乎完全互补)
#   2) pct_married        (与 pct_never_married 高度负相关)
#   3) pct_detached_1unit (与 pct_multi_unit 高度负相关)
# - 后续 PCA 将基于筛选后的变量继续运行

vars_to_remove_corr <- c(
  "pct_owner_occ",
  "pct_married",
  "pct_detached_1unit"
)

acs_nyc_analysis_pca <- acs_nyc_analysis %>%
  select(-all_of(vars_to_remove_corr))

# ------------------------------------------------------------------------------
# 15.1 查看被移除的变量
# ------------------------------------------------------------------------------
vars_to_remove_corr
[1] "pct_owner_occ"      "pct_married"        "pct_detached_1unit"
# ------------------------------------------------------------------------------
# 15.2 检查筛选后的变量数量
# ------------------------------------------------------------------------------
acs_nyc_analysis_pca %>%
  st_drop_geometry() %>%
  select(-GEOID, -NAME) %>%
  ncol()
[1] 41
# ------------------------------------------------------------------------------
# 15.3 查看筛选后仍保留的变量名称
# ------------------------------------------------------------------------------
 names(
  acs_nyc_analysis_pca %>%
    st_drop_geometry() %>%
    select(-GEOID, -NAME)
)
 [1] "pct_age_0_17"        "pct_age_18_34"       "pct_age_35_64"      
 [4] "pct_age_65p"         "pct_female"          "pct_nh_white"       
 [7] "pct_nh_black"        "pct_nh_asian"        "pct_hispanic"       
[10] "pct_other_eth"       "pop_density"         "pct_never_married"  
[13] "pct_separated"       "pct_div_widowed"     "pct_edu_low"        
[16] "pct_edu_secondary"   "pct_edu_bachelor"    "pct_edu_gradplus"   
[19] "med_hh_income"       "poverty_rate"        "gini_index"         
[22] "unemployment_rate"   "not_in_lf_rate"      "pct_mgmt_prof"      
[25] "pct_service"         "pct_sales_off"       "pct_manual"         
[28] "pct_prod_trans"      "pct_renter_occ"      "med_gross_rent"     
[31] "pct_attached_1unit"  "pct_multi_unit"      "pct_commute_car"    
[34] "pct_commute_transit" "pct_commute_walk"    "pct_commute_home"   
[37] "pct_tt_lt5"          "pct_tt_5_14"         "pct_tt_15_29"       
[40] "pct_tt_30_59"        "pct_tt_60plus"      

方法提示:什么是 Pearson 相关系数?

Pearson 相关系数(Pearson correlation coefficient)是衡量两个连续变量之间线性相关程度的常用指标,通常记为 \(r\),其取值范围在 -11 之间。其基本形式可写为:

\[ r = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2}\sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}} \]

其中,\(x_i\)\(y_i\) 表示两个变量在第 \(i\) 个观测上的取值,\(\bar{x}\)\(\bar{y}\) 分别表示它们的均值。

从解读上看,当 \(r\) 接近 1 时,表示两个变量之间存在较强的正线性相关;当 \(r\) 接近 -1 时,表示存在较强的负线性相关;当 \(r\) 接近 0 时,则说明两者之间的线性关系较弱。需要注意的是,Pearson 相关系数反映的是线性关系,而不是所有形式的关联。

在本节中,Pearson 相关系数的主要作用,是帮助识别高度相关的变量对,从而发现潜在的变量冗余,并为后续的 PCA 与聚类分析提供参考。例如,若两个变量之间高度正相关,则它们很可能在一定程度上表达了相近的信息。

不过,Pearson 相关系数也存在局限。它对异常值较敏感,并且默认关注的是线性关系;若变量之间呈现明显的非线性关系,Pearson 相关系数可能无法充分反映它们之间的关联程度。此外,当变量分布偏态较强或本身更适合看排序关系时,也可以考虑使用 Spearman 等级相关系数(Spearman correlation) 作为替代方法。

PCA数据筛选
4.4 PCA 分析导入

主成分分析 基础理论

在前一节完成变量重构相关性检验之后,可以发现纽约市 census tract 尺度下的分析变量已经明显增多,不同变量之间也存在不同程度的相关性。这种情况在人口普查数据分析中非常常见:一方面,多个变量可能共同反映相似的人口、社会经济维度;另一方面,若将所有变量不加筛选地直接纳入后续聚类分析,不仅会增加模型复杂度,也可能因为变量冗余而削弱分类结果的稳定性与解释性。

因此,在正式进入聚类分析之前,本节先引入主成分分析(Principal Component Analysis, PCA),作为一种常见的降维变量筛选辅助方法。这里使用 PCA 的目的,并不是要直接用主成分得分替代全部原始变量(虽然这种做法也较常见),而是借助它梳理变量之间的整体结构,识别哪些变量在相似维度上高度重叠,哪些变量则更能代表不同的社会空间特征。

本节的筛选思路大致如下:首先,对前一节整理好的分析变量进行标准化处理,并执行 PCA;其次,根据特征值解释方差比例初步判断哪些主成分值得保留;再次,结合 Kaiser rule(特征值大于 1)、碎石图以及各变量在主成分上的载荷贡献度,识别主要维度中更具代表性的变量;最后,在此基础上形成一个相对精简、但仍保留核心信息的变量集合,用于后续的 H-K-means 聚类分析。

需要注意的是,PCA 在这里更适合作为一种探索性工具。它能够帮助我们理解变量结构、辅助筛选变量,但并不意味着只要经过 PCA,后续变量选择就完全是“客观自动”的。最终保留哪些变量,仍需结合研究目标变量解释性城市社会空间分析的实际需要综合判断。换言之,PCA 提供的是一种结构化参考,而不是唯一答案。

方法提示:什么是 PCA?

主成分分析(PCA)是一种经典的降维方法(dimension reduction)。它的基本思想是:将原来多个彼此相关的变量,重新线性组合为少数几个彼此正交的新变量,即主成分(principal components)。这些主成分按照解释数据总体变异的能力依次排序,第一个主成分解释的变异最多,第二个主成分解释次之,依此类推。

从数学形式上看,第 \(k\) 个主成分可表示为:

\[ PC_k = a_{k1}x_1 + a_{k2}x_2 + \cdots + a_{kp}x_p \]

其中,\(x_1, x_2, \dots, x_p\) 表示原始变量,\(a_{k1}, a_{k2}, \dots, a_{kp}\) 表示对应主成分上的系数或载荷。不同主成分本质上代表了原始变量在某种综合方向上的主要变化模式。

在实际应用中,PCA 常用于:

  • 高维数据降维:在尽量保留主要信息的前提下,减少变量数量;
  • 变量结构识别:观察哪些变量在同一主成分上共同变化;
  • 变量筛选辅助:结合载荷与贡献度,识别更具代表性的变量;
  • 指数构建下一章内容
  • 可视化探索:通过主成分空间查看观测对象或变量的大致分布格局。

不过,PCA 也有其局限。它本质上是一种线性变换,更擅长提取线性结构;同时,主成分往往是多个变量的组合,因此在解释上有时不如原始变量那样直观。正因如此,在本专题中,PCA 的作用主要是辅助理解变量结构与支持变量筛选,而不是直接替代后续全部分析。


本节任务

基于前一节已经完成的变量整理与相关性检验,本节将正式进入 PCA 分析阶段,主要围绕以下几个任务展开:

  • 对当前纽约市 tract 尺度的多维分析变量进行标准化,并执行主成分分析
  • 结合特征值解释方差与碎石图,判断值得重点关注的主成分;
  • 查看各变量在主要主成分上的载荷贡献度,识别潜在的变量聚合结构;
  • 在统计结果与解释性并重的基础上,形成后续 H-K-means 聚类分析所使用的候选变量集合。

主成分分析 实战代码

# ==============================================================================
# 16. PCA 前的数据标准化
# ==============================================================================
# 说明:
# - PCA 对变量量纲较为敏感,因此需先进行标准化
# - 这里使用 Z-score 标准化,使各变量均值为 0、标准差为 1
# - 输入数据基于前一步已经完成高相关变量初筛的 acs_nyc_analysis_pca
# - 后续 PCA 与 H-K-means 都将基于该标准化数据展开

nyc_pca_vars <- acs_nyc_analysis_pca %>%
  st_drop_geometry() %>%
  select(-GEOID, -NAME)

nyc_pca_scaled <- scale(nyc_pca_vars) %>%
  as.data.frame()
# ==============================================================================
# 17. 执行 PCA 分析
# ==============================================================================
# 说明:
# - 这里使用 FactoMineR::PCA() 进行主成分分析
# - 因前一步已完成标准化,这里直接输入标准化后的数据
# - ncp 设为变量总数,避免默认仅保留前 5 个主成分
# - graph = FALSE 表示先不自动绘图,后续再分别提取结果并可视化

nyc_pca <- FactoMineR::PCA(
  nyc_pca_scaled,
  ncp = ncol(nyc_pca_scaled),
  graph = FALSE
)

# 查看 PCA 结果摘要
print(nyc_pca)
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 2170 individuals, described by 41 variables
*The results are available in the following objects:

   name               description                          
1  "$eig"             "eigenvalues"                        
2  "$var"             "results for the variables"          
3  "$var$coord"       "coord. for the variables"           
4  "$var$cor"         "correlations variables - dimensions"
5  "$var$cos2"        "cos2 for the variables"             
6  "$var$contrib"     "contributions of the variables"     
7  "$ind"             "results for the individuals"        
8  "$ind$coord"       "coord. for the individuals"         
9  "$ind$cos2"        "cos2 for the individuals"           
10 "$ind$contrib"     "contributions of the individuals"   
11 "$call"            "summary statistics"                 
12 "$call$centre"     "mean of the variables"              
13 "$call$ecart.type" "standard error of the variables"    
14 "$call$row.w"      "weights for the individuals"        
15 "$call$col.w"      "weights for the variables"          
# ==============================================================================
# 18. 提取 PCA 的特征值与解释方差
# ==============================================================================
# 说明:
# - 这里无需手动计算 eigenvalue、variance explained 等指标
# - factoextra::get_eigenvalue() 可直接提取:
#   1) eigenvalue:特征值
#   2) variance.percent:单个主成分的解释方差比例
#   3) cumulative.variance.percent:累计解释方差比例
# - 这些结果将为后续 scree plot 与主成分筛选提供依据
# - 这里进一步将结果转为 data.frame,便于后续筛选与可视化

pca_eigen <- factoextra::get_eigenvalue(nyc_pca) %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "PC")

# 查看 PCA 结果整理表
head(pca_eigen, 3)
     PC eigenvalue variance.percent cumulative.variance.percent
1 Dim.1  10.094740        24.621318                    24.62132
2 Dim.2   6.037202        14.724882                    39.34620
3 Dim.3   3.013068         7.348947                    46.69515

主成分筛选

在完成 PCA 之后,下一步通常需要判断哪些主成分值得保留。这里的判断并不是机械地“保留前两个主成分”,而是需要结合多个信息来源综合考虑,包括各主成分的特征值解释方差比例累计解释方差以及后续的碎石图(scree plot)

在常见做法中,Kaiser rule(保留特征值大于 1 的主成分)是一种较常用的初步判断标准;研究者也会结合碎石图中曲线由陡转缓的位置,识别解释能力开始明显减弱的拐点。后续主成分是否保留,还需要结合变量载荷与研究解释性进一步判断。

# ==============================================================================
# 19. PCA 碎石图:辅助筛选值得保留的主成分
# ==============================================================================
# 说明:
# - fviz_eig() 可直接基于 PCA 结果绘制 scree plot
# - 柱形表示各主成分的解释方差百分比
# - 红线表示累计解释方差变化趋势
# - 灰色虚线表示 Kaiser rule 的参考阈值
# - 由于本节基于标准化变量进行 PCA,因此:
#   Kaiser rule(特征值 > 1)可换算为解释方差百分比 > 100 / 变量数

kaiser_line <- 100 / ncol(nyc_pca_scaled)

factoextra::fviz_eig(
  nyc_pca,
  ncp = 20,
  addlabels = TRUE,
  ylim = c(0, max(pca_eigen %>% slice(1:20) %>% pull(variance.percent)) + 5),
  barfill = "steelblue",
  barcolor = "steelblue",
  linecolor = "red",
  main = "Scree plot of principal components",
  xlab = "Principal component",
  ylab = "Percentage of explained variance"
) +
  geom_hline(
    yintercept = kaiser_line,
    linetype = "dashed",
    color = "gray50",
    linewidth = 0.7
  ) +
  annotate(
    "text",
    x = 16,
    y = kaiser_line + 1.2,
    label = paste0("Kaiser rule (", round(kaiser_line, 2), "%)"),
    color = "gray50",
    size = 3.2
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 10)
  )

# ==============================================================================
# 20. 提取变量在 PCA 中的结果
# ==============================================================================
# 说明:
# - factoextra::get_pca_var() 可直接提取变量层面的 PCA 结果
# - 返回的结果主要包括:
#   1) coord   : 变量在主成分空间中的坐标
#   2) cor     : 变量与主成分之间的相关系数
#   3) cos2    : 变量在主成分空间中的表示质量
#   4) contrib : 变量对主成分的贡献率(%)
#
# 补充理解:
# - cos2 可理解为变量在某主成分上的“解释质量”,通常与坐标平方相关
# - contrib 可理解为变量对某主成分的重要性,其本质上是该变量 cos2
#   占该主成分总 cos2 的比例,再换算为百分比

pca_var <- factoextra::get_pca_var(nyc_pca)

# 查看对象结构
names(pca_var)
[1] "coord"   "cor"     "cos2"    "contrib"

变量贡献筛选

在通过 Kaiser rule 初步确定需要保留的主成分(PCs)之后,下一步可以进一步查看各变量在这些主成分上的贡献率(contribution)。贡献率越高,通常意味着该变量在对应主成分所代表的结构维度中越重要。

本节采用一种较为直观的初步筛选思路:先保留满足 Kaiser rule 的主成分,再计算每个变量在这些主成分上的平均贡献率。若某变量的平均贡献率高于全部变量的平均水平,则说明其在保留主成分中总体贡献相对更突出,可作为后续聚类分析的优先候选变量。

需要说明的是,这一规则更适合作为探索性筛选的第一步,而不是唯一标准。后续仍应结合变量解释性、变量冗余与研究目标进行综合判断。

# ==============================================================================
# 21.  根据 Kaiser rule 确定保留的主成分
# ==============================================================================
# 说明:
# - Kaiser rule:保留特征值 > 1 的主成分
# - 这里直接基于前面已经整理好的 pca_eigen 结果进行筛选

retained_pcs <- pca_eigen %>%
  filter(eigenvalue > 1) %>%
  pull(PC)

# 查看保留的主成分名称与数量
retained_pcs
 [1] "Dim.1"  "Dim.2"  "Dim.3"  "Dim.4"  "Dim.5"  "Dim.6"  "Dim.7"  "Dim.8" 
 [9] "Dim.9"  "Dim.10"
length(retained_pcs)
[1] 10

根据 Kaiser rule,这里仅保留 eigenvalue > 1 的主成分(PCs)

# ==============================================================================
# 22.  保留主成分上的变量贡献率可视化
# ==============================================================================
# 说明:
# - 展示保留主成分中各变量的 contribution
# - 若保留主成分较多,图会较长,这是正常现象
# - 这里用 factoextra 自带函数进行展示

retained_axes <- seq_along(retained_pcs)

factoextra::fviz_contrib(
  nyc_pca,
  choice = "var",
  axes = retained_axes,
  top = 40
) +
  labs(
    title = "Top variable contributions across retained principal components"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    axis.text.x = element_text(
      angle = 90,
      vjust = 0.5,
      hjust = 1
    )
  )

上图中我们主要关注红色虚线以上的变量名称,即贡献率高于平均值的变量。这些变量可视为后续聚类分析中优先保留的候选变量


可视化优化 可选

factoextra::fviz_contrib()提供了一个默认的可视化样式,即柱状图;下方代码展示的是对此柱状图进行适度美化,提供多种可视化的思路,如果遇到问题可以选择跳过

下方代码与上方图形表达的是同一思路,只是换用了热力图与棒棒糖图等方式进行展示;总体上仍然是只保留 contribution 高于平均值的变量

library(patchwork)
# ==============================================================================
# 22. 保留主成分上的变量贡献率与筛选结果
# ==============================================================================
# 说明:
# - 为保证与 fviz_contrib() 完全一致,这里直接使用 facto_summarize()
# - 右图棒棒糖图与变量筛选均基于同一份 fviz_rank 结果
# - 左图热力图仍展示各变量在各保留主成分上的 contribution 分布
# - 红色虚线表示平均 contribution 的参考阈值:100 / 变量总数

library(patchwork)

# ------------------------------------------------------------------------------
# 22.1 提取保留主成分
# ------------------------------------------------------------------------------
retained_idx <- which(pca_eigen$eigenvalue > 1)
retained_axes <- seq_along(retained_idx)
retained_pc_names <- colnames(pca_var$contrib[, retained_idx, drop = FALSE])

# ------------------------------------------------------------------------------
# 22.2 提取与 fviz_contrib() 完全一致的变量 contribution 排序结果
# ------------------------------------------------------------------------------
fviz_rank <- factoextra::facto_summarize(
  X = nyc_pca,
  element = "var",
  result = "contrib",
  axes = retained_axes
) %>%
  as_tibble() %>%
  rename(
    variable = name,
    contribution = contrib
  ) %>%
  arrange(desc(contribution))

# 平均 contribution 阈值(与 fviz_contrib 图中的虚线一致)
contrib_cutoff <- 100 / ncol(nyc_pca_scaled)

fviz_rank <- fviz_rank %>%
  mutate(
    selected = if_else(contribution >= contrib_cutoff, "Selected", "Removed")
  )

# 查看前 10 个变量,核对排序
fviz_rank %>% slice(1:10)
# A tibble: 10 × 3
   variable          contribution selected
   <fct>                    <dbl> <chr>   
 1 pct_age_18_34             3.13 Selected
 2 pct_mgmt_prof             3.12 Selected
 3 pct_age_0_17              3.03 Selected
 4 pct_never_married         2.90 Selected
 5 pct_age_65p               2.87 Selected
 6 pct_edu_gradplus          2.84 Selected
 7 pct_commute_car           2.82 Selected
 8 pct_renter_occ            2.80 Selected
 9 pct_multi_unit            2.79 Selected
10 pct_tt_30_59              2.78 Selected
# ------------------------------------------------------------------------------
# 22.3 整理逐 PC contribution 数据,用于热力图
# ------------------------------------------------------------------------------
contrib_long <- pca_var$contrib[, retained_idx, drop = FALSE] %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "variable") %>%
  pivot_longer(
    cols = -variable,
    names_to = "PC",
    values_to = "contribution_pc"
  )

# 为左右图统一变量顺序
var_order <- fviz_rank$variable

contrib_long <- contrib_long %>%
  mutate(
    variable = factor(variable, levels = rev(var_order)),
    PC = factor(
      PC,
      levels = retained_pc_names,
      labels = paste0("PC", seq_along(retained_pc_names))
    )
  )

fviz_rank <- fviz_rank %>%
  mutate(
    variable = factor(variable, levels = rev(var_order))
  )

# ------------------------------------------------------------------------------
# 22.4 左图:各保留主成分上的 contribution 热力图
# ------------------------------------------------------------------------------
p_heatmap <- ggplot(contrib_long, aes(x = PC, y = variable, fill = contribution_pc)) +
  geom_tile(color = "white", linewidth = 0.25) +
  scale_fill_gradient(
    low = "#f7fbff",
    high = "#2171b5",
    name = "Contribution (%)"
  ) +
  labs(
    title = "Contribution across retained principal components",
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 10.5) +
  theme(
    plot.title = element_text(face = "bold", size = 11.5),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    panel.grid = element_blank(),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 8)
  )

# ------------------------------------------------------------------------------
# 22.5 右图:与 fviz_contrib 一致的棒棒糖图
# ------------------------------------------------------------------------------
p_lollipop <- ggplot(fviz_rank, aes(x = contribution, y = variable, color = selected)) +
  geom_segment(
    aes(x = 0, xend = contribution, y = variable, yend = variable),
    linewidth = 0.5,
    color = "grey75"
  ) +
  geom_point(size = 2.8) +
  geom_vline(
    xintercept = contrib_cutoff,
    linetype = "dashed",
    linewidth = 0.8,
    color = "red"
  ) +
  scale_color_manual(
    values = c("Selected" = "#1b9e77", "Removed" = "#d95f02"),
    name = NULL
  ) +
  labs(
    title = "Mean contribution and variable selection",
    x = "Contribution (%)",
    y = NULL
  ) +
  theme_minimal(base_size = 10.5) +
  theme(
    plot.title = element_text(face = "bold", size = 11.5),
    panel.grid.major.y = element_blank(),
    legend.position = "top"
  )

# ------------------------------------------------------------------------------
# 22.6 拼接输出
# ------------------------------------------------------------------------------
p_heatmap + p_lollipop +
  plot_layout(widths = c(1.5, 1))


最终保留变量

根据上面的fviz_rank 表,选择contribution >= 平均值的变量;即所有Selected的标签

# ==============================================================================
# 23. 提取最终保留的候选变量
# ==============================================================================
# 说明:
# - 基于前一步 PCA contribution 筛选结果
# - 仅保留平均 contribution 高于阈值的变量
# - 这些变量将作为后续 H-K-means 聚类分析的候选输入变量

selected_vars_pca <- fviz_rank %>%
  filter(selected == "Selected") %>%
  pull(variable) %>%
  as.character()

# 查看保留变量名称
selected_vars_pca
 [1] "pct_age_18_34"       "pct_mgmt_prof"       "pct_age_0_17"       
 [4] "pct_never_married"   "pct_age_65p"         "pct_edu_gradplus"   
 [7] "pct_commute_car"     "pct_renter_occ"      "pct_multi_unit"     
[10] "pct_tt_30_59"        "pct_nh_white"        "pct_nh_black"       
[13] "med_hh_income"       "pct_nh_asian"        "pct_commute_walk"   
[16] "pct_edu_bachelor"    "poverty_rate"        "pct_commute_transit"
[19] "not_in_lf_rate"      "pct_tt_15_29"        "pct_edu_low"        
[22] "pct_service"         "pct_hispanic"        "med_gross_rent"     
# 查看保留变量数量
length(selected_vars_pca)
[1] 24

方法提醒:本节变量筛选具有教学演示性质

根据上面的结果,最终共有 24 个变量被保留下来,并带入后续聚类分析。

需要说明的是,这一部分主要服务于教学示范假想课题展示,因此这里采用了相对较为严格、也较为机械的 PCA contribution 筛选规则。这样做有助于清晰展示变量筛选流程,但也可能削弱某些领域(如住房、就业或通勤)在最终分类结果中的解释性。

在真实的地理人口特征分类研究中,研究者通常不会仅凭一次 PCA 结果就大幅删除变量,而是还需要综合考虑研究目标分类解释性实际应用场景以及不同利益相关者(stakeholders)的关注重点。换言之,统计筛选只是变量决策的一部分,而不是唯一依据。

聚类分析与类型识别
4.5 H-K-means 聚类分析导入

聚类分析 方法衔接

经过前一节的 PCA 分析与变量筛选之后,当前数据已经从原始的高维变量集合中提炼出一组更具代表性的候选变量。这一步并不意味着这些变量彼此完全独立,而是说明它们在前述主成分结构中表现出相对更高的贡献度,因此更适合作为后续地理人口特征分类的输入变量。

PCA 完成的是一个前置筛选过程,其结果是为下一步聚类分析提供一个相对精简、同时保留主要信息的变量基础。

在此基础上,本节将把前一节筛选出的变量正式带入 H-K-means(Hierarchical K-means) 聚类分析。其核心目标在于:基于纽约市 census tract 在多维社会人口、住房与通勤属性上的综合相似性,将这些空间单元划分为若干具有内部一致性、同时彼此之间存在明显差异的类型。与前一专题中较为基础的 K-means 示例不同,这里将进一步使用 H-K-means,以获得相对更稳健的初始分类结果,并为后续的空间可视化、类型解释与命名奠定基础。

本节的主要步骤包括:

  • 整理前述 PCA 保留下来的变量,构建最终聚类输入数据;
  • 结合多种方法对聚类数量 (k) 进行比较与判断;
  • 执行 H-K-means 聚类,并将结果重新连接回空间数据;
  • 结合空间分布图、类别均值与变量特征,对聚类结果进行初步解释。

方法提示:为什么使用 H-K-means?

相较于普通的 K-means,H-K-means 可以理解为一种结合了层次聚类(Hierarchical Clustering)K-means优势的聚类方式。普通 K-means 往往受到初始中心点选择的影响,不同随机起点有时会得到略有差异的结果;当变量较多、数据结构较复杂时,这种不稳定性会更加明显。

H-K-means 的基本思路是:先借助层次聚类提供更有结构性的初始分组,再在此基础上执行 K-means 迭代优化。这样做通常有助于减少对随机初始化的依赖,并获得相对更稳定、也更易解释的分类结果。

当然,H-K-means 也并非“完美方法”。它仍然需要研究者预先考虑聚类数量 k,而最终分类方案也仍需结合研究目标、变量体系与空间解释性综合判断。


实战代码 H-K-means聚类分析

实战代码 聚类输入准备

# ==============================================================================
# 24. H-K-means 聚类输入数据准备
# ==============================================================================
# 说明:
# - 这里使用前一节 PCA contribution 筛选后保留下来的变量
# - 输入数据继续使用标准化后的 nyc_pca_scaled
# - 后续 H-K-means 聚类将基于这些变量展开

nyc_cluster_input <- nyc_pca_scaled %>%
  select(all_of(selected_vars_pca))

# 查看输入数据规模
dim(nyc_cluster_input)
[1] 2170   24

聚类数量判断

在正式执行 H-K-means 之前,仍然需要先对聚类数量 (k) 作一个初步判断。与前一专题类似,这里不会只依赖单一标准,而是结合多种方法进行比较。较常见的思路包括:观察组内离差平方和变化趋势的肘部法(Elbow method),以及利用若干内部指标比较不同 k 下分类质量的变化。

需要注意的是,聚类数量的选择通常并不存在唯一的“标准答案”。图形指标与统计指标能够提供重要参考,但最终仍需结合分类结果的稳定性、类别规模与空间解释性综合判断。

# ==============================================================================
# 25. 使用 Elbow method 辅助判断聚类数量
# ==============================================================================
# 说明:
# - 这里先使用 WSS(组内平方和)观察拐点位置
# - 为保证结果可复现,设定随机种子
# - k.max 这里先考察 2 到 10 类

set.seed(1234)

p_elbow_nyc <- factoextra::fviz_nbclust(
  x = nyc_cluster_input,
  FUNcluster = factoextra::hkmeans,
  method = "wss",
  k.max = 10
) +
  labs(
    title = "Elbow method",
    subtitle = "Within-cluster sum of squares"
  ) +
  theme_minimal(base_size = 11)

p_elbow_nyc

这里我们发现 k = 3 / 6 / 9 时 折线图有明显的 elbow!

# ==============================================================================
# 26. 使用常见指标辅助比较不同聚类数量
# ==============================================================================
# 说明:
# - 相较于 NbClust(index = "alllong"),这里改用几种更常见、也更直观的方法:
#   1) Elbow method / WSS
#   2) Silhouette
#   3) Gap statistic
# - 这些方法可共同辅助判断 k 的合理范围
# - 其中 Gap statistic 计算相对更慢,因此这里设置 bootstrap 次数为 50

set.seed(1234)

# ------------------------------------------------------------------------------
# 26.1 Elbow method
# ------------------------------------------------------------------------------
p_elbow_nyc <- factoextra::fviz_nbclust(
  x = nyc_cluster_input,
  FUNcluster = factoextra::hkmeans,
  method = "wss",
  k.max = 10
) +
  labs(
    title = "Elbow method",
    subtitle = "Within-cluster sum of squares"
  ) +
  theme_minimal(base_size = 11)

# ------------------------------------------------------------------------------
# 26.2 Silhouette method
# ------------------------------------------------------------------------------
p_silhouette_nyc <- factoextra::fviz_nbclust(
  x = nyc_cluster_input,
  FUNcluster = factoextra::hkmeans,
  method = "silhouette",
  k.max = 10
) +
  labs(
    title = "Silhouette method",
    subtitle = "Average silhouette width"
  ) +
  theme_minimal(base_size = 11)

# ------------------------------------------------------------------------------
# 26.3 Gap statistic
# ------------------------------------------------------------------------------
gap_nyc <- cluster::clusGap(
  x = nyc_cluster_input,
  FUNcluster = function(x, k) {
    list(cluster = factoextra::hkmeans(x, k = k)$cluster)
  },
  K.max = 10,
  B = 50
)

p_gap_nyc <- factoextra::fviz_gap_stat(gap_nyc) +
  labs(
    title = "Gap statistic",
    subtitle = "Bootstrap comparison across different k"
  ) +
  theme_minimal(base_size = 11)

# ------------------------------------------------------------------------------
# 26.4 拼接输出
# ------------------------------------------------------------------------------
(p_elbow_nyc | p_silhouette_nyc) / p_gap_nyc

从上述结果看,不同指标对聚类数量的建议并不完全一致。Elbow method 在 k = 3 附近呈现出较明显的拐点,而 Silhouette method 更倾向于较小的聚类数,Gap statistic 则随着 k 增大持续上升。这样的结果说明,当前数据并不存在一个由所有指标共同支持的唯一最优解。
因此,本节不直接机械采用某一个指标给出的结果,而是将 k = 4k = 6 视为更值得重点比较的候选范围,并进一步结合类别规模、空间分布与解释性选择最终方案。作为教学展示,后续将先以 k = 6 为示例继续进行 H-K-means 聚类分析。

H-K-means 聚类 实战代码

# ==============================================================================
# 27. 执行 H-K-means 聚类
# ==============================================================================
# 说明:
# - 这里先以 k = 6 作为示例继续进行聚类分析
# - hkmeans() 先用层次聚类生成初始分组,再用 K-means 迭代优化
# - 输入数据为前面 PCA 筛选后保留变量构成的聚类输入数据

set.seed(1234)

k_selected_nyc <- 6

hk_nyc <- factoextra::hkmeans(
  x = nyc_cluster_input,
  k = k_selected_nyc
)

# 查看聚类对象结构
names(hk_nyc)
 [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
 [6] "betweenss"    "size"         "iter"         "ifault"       "data"        
[11] "hclust"      
# 查看各类别样本量 (确定有没有特别大,特别小的聚类产生)
table(hk_nyc$cluster)

  1   2   3   4   5   6 
436 395 309 418 341 271 

层次结构查看 可选

这一步可以跳过, 因为observation数量越大,运行速度越慢!

# ==============================================================================
# 28. 查看 H-K-means 的层次树状图
# ==============================================================================
# 说明:
# - hkmeans() 对象中包含层次聚类的结果
# - 可通过 fviz_dend() 查看其树状结构
# - 这里展示切分为 6 类时的树状图,用于辅助观察类别分裂关系

factoextra::fviz_dend(
  hk_nyc,
  k = k_selected_nyc,
  rect = TRUE,
  rect_border = "jco",
  rect_fill = TRUE,
  show_labels = FALSE,
  main = "Hierarchical tree of H-K-means clustering"
) +
  theme_minimal(base_size = 11)


聚类结果回接 重要

# ==============================================================================
# 29. 将 H-K-means 聚类结果回接到空间数据
# ==============================================================================
# 说明:
# - 将聚类编号写回原始 tract 空间数据
# - 后续即可进行空间制图与类型解释

acs_nyc_clustered <- acs_nyc_analysis %>%
  mutate(
    cluster = factor(hk_nyc$cluster)
  )

# 查看各类样本量
acs_nyc_clustered %>%
  st_drop_geometry() %>%
  count(cluster)
  cluster   n
1       1 436
2       2 395
3       3 309
4       4 418
5       5 341
6       6 271

聚类结果综合展示 空间分布 + 属性画像

在完成 H-K-means 聚类并将结果回接到空间数据之后,下一步需要同时从空间格局内部属性两个角度来理解分类结果。为此,下方代码将把两类图形放在同一组图中展示:左侧为纽约市 census tract 的聚类空间分布图,用于观察不同类别在城市内部的区位格局;右侧为各聚类相对于纽约市整体平均水平(Index = 100)的指数热力图,用于比较不同类别在年龄、族裔、教育、收入、住房与通勤等维度上的相对特征。

需要说明的是,下方代码框虽然较长,但其中的核心技巧其实在前文已经多次出现过,例如:聚类结果回接空间数据、计算整体均值与类别均值、构造 index score、整理长表以及绘制热力图与空间分布图等。因此,这里更多是对前述方法的综合应用,并不是全新的技术环节。

# ==============================================================================
# 30. NYC 聚类结果综合展示:空间分布图 + 指数热力图
# ==============================================================================
# 说明:
# - 左图:H-K-means 聚类结果的空间分布图
# - 右图:各 cluster 相对于 NYC 整体平均水平的指数热力图
# - Index = 100 表示与 NYC 整体平均一致
# - Index > 100 表示高于整体平均
# - Index < 100 表示低于整体平均
# - 为便于解释,这里对变量顺序进行人工整理,使同类变量尽量放在一起

library(patchwork)
library(RColorBrewer)

# ------------------------------------------------------------------------------
# 30.1 手动整理变量顺序
# ------------------------------------------------------------------------------
profile_vars_nyc <- c(
  # 年龄
  "pct_age_0_17",
  "pct_age_18_34",
  "pct_age_65p",
  
  # 族裔
  "pct_nh_white",
  "pct_nh_black",
  "pct_nh_asian",
  "pct_hispanic",
  
  # 教育与就业
  "pct_edu_low",
  "pct_edu_bachelor",
  "pct_edu_gradplus",
  "pct_mgmt_prof",
  "pct_service",
  "not_in_lf_rate",
  
  # 收入与贫困
  "med_hh_income",
  "poverty_rate",
  "med_gross_rent",
  
  # 婚姻与住房
  "pct_never_married",
  "pct_renter_occ",
  "pct_multi_unit",
  
  # 通勤
  "pct_commute_car",
  "pct_commute_transit",
  "pct_commute_walk",
  "pct_tt_15_29",
  "pct_tt_30_59"
)

# 只保留当前 PCA 最终进入聚类的变量,且按照上面顺序排列
profile_vars_nyc <- profile_vars_nyc[profile_vars_nyc %in% selected_vars_pca]

# ------------------------------------------------------------------------------
# 30.2 变量标签
# ------------------------------------------------------------------------------
variable_labels_nyc <- c(
  pct_age_0_17       = "Age 0–17",
  pct_age_18_34      = "Age 18–34",
  pct_age_65p        = "Age 65+",
  pct_nh_white       = "NH White",
  pct_nh_black       = "NH Black",
  pct_nh_asian       = "NH Asian",
  pct_hispanic       = "Hispanic",
  pct_edu_low        = "Low education",
  pct_edu_bachelor   = "Bachelor",
  pct_edu_gradplus   = "Graduate+",
  pct_mgmt_prof      = "Mgmt/prof",
  pct_service        = "Service",
  not_in_lf_rate     = "Not in labour force",
  med_hh_income      = "Median income",
  poverty_rate       = "Poverty",
  med_gross_rent     = "Median rent",
  pct_never_married  = "Never married",
  pct_renter_occ     = "Renter occupied",
  pct_multi_unit     = "Multi-unit housing",
  pct_commute_car    = "Commute by car",
  pct_commute_transit= "Commute by transit",
  pct_commute_walk   = "Commute by walk",
  pct_tt_15_29       = "Travel time 15–29",
  pct_tt_30_59       = "Travel time 30–59"
)

# ------------------------------------------------------------------------------
# 30.3 计算 NYC 整体平均值
# ------------------------------------------------------------------------------
overall_means_nyc <- acs_nyc_clustered %>%
  st_drop_geometry() %>%
  select(all_of(profile_vars_nyc)) %>%
  summarise(
    across(everything(), mean, na.rm = TRUE)
  )

# ------------------------------------------------------------------------------
# 30.4 计算各 cluster 均值
# ------------------------------------------------------------------------------
cluster_means_nyc <- acs_nyc_clustered %>%
  st_drop_geometry() %>%
  group_by(cluster) %>%
  summarise(
    across(all_of(profile_vars_nyc), mean, na.rm = TRUE),
    .groups = "drop"
  )

# ------------------------------------------------------------------------------
# 30.5 计算 index score(NYC mean = 100)
# ------------------------------------------------------------------------------
index_scores_nyc <- cluster_means_nyc %>%
  mutate(
    across(
      all_of(profile_vars_nyc),
      ~ .x / overall_means_nyc[[cur_column()]] * 100
    )
  ) %>%
  mutate(
    across(where(is.numeric), ~ round(.x, 1))
  )

# ------------------------------------------------------------------------------
# 30.6 转为热力图长表
# ------------------------------------------------------------------------------
index_long_nyc <- index_scores_nyc %>%
  pivot_longer(
    cols = -cluster,
    names_to = "variable",
    values_to = "index"
  ) %>%
  mutate(
    variable = factor(variable, levels = rev(profile_vars_nyc)),
    index_group = case_when(
      index <= 50                ~ "<= 50",
      index > 50  & index < 80   ~ "50–80",
      index >= 80 & index < 95   ~ "80–95",
      index >= 95 & index <= 105 ~ "95–105",
      index > 105 & index <= 120 ~ "105–120",
      index > 120 & index < 200  ~ "120–200",
      index >= 200               ~ ">= 200"
    ),
    index_group = factor(
      index_group,
      levels = c("<= 50", "50–80", "80–95", "95–105", "105–120", "120–200", ">= 200")
    )
  )

# ------------------------------------------------------------------------------
# 30.7 左图:聚类结果空间分布图
# ------------------------------------------------------------------------------
graticule_sf <- st_graticule(acs_nyc_clustered)

p_map_nyc <- ggplot() +
  geom_sf(
    data = graticule_sf,
    color = "grey80",
    linewidth = 0.2,
    linetype = "dashed"
  ) +
  geom_sf(
    data = acs_nyc_clustered,
    aes(fill = cluster),
    color = NA
  ) +
  annotation_scale(
    location = "bl",
    width_hint = 0.25,
    text_cex = 0.75
  ) +
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    style = north_arrow_fancy_orienteering,
    height = unit(0.8, "cm"),
    width  = unit(0.8, "cm")
  ) +
  scale_fill_brewer(
    palette = "Set2",
    name = "Cluster"
  ) +
  coord_sf(crs = st_crs(4326)) +
  labs(
    title = "H-K-means clustering of NYC census tracts",
    subtitle = paste0("k = ", k_selected_nyc)
  ) +
  theme_minimal(base_size = 11.5) +
  theme(
    panel.grid.major = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.title = element_text(face = "bold", size = 12.5),
    plot.subtitle = element_text(size = 10.5),
    legend.title = element_text(size = 9.5),
    legend.text = element_text(size = 9),
    axis.title = element_blank()
  )

# ------------------------------------------------------------------------------
# 30.8 右图:指数热力图(竖向)
# ------------------------------------------------------------------------------
p_heatmap_nyc <- ggplot(index_long_nyc, aes(x = cluster, y = variable, fill = index_group)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(
    aes(label = round(index, 1)),
    size = 3.1,
    color = "black"
  ) +
  scale_fill_manual(
    values = c(
      "<= 50"   = "#B2182B",
      "50–80"   = "#EF8A62",
      "80–95"   = "#FDDBC7",
      "95–105"  = "white",
      "105–120" = "#D9F0D3",
      "120–200" = "#7FBF7B",
      ">= 200"  = "#1B7837"
    ),
    name = "Index score"
  ) +
  scale_y_discrete(labels = function(x) variable_labels_nyc[x]) +
  labs(
    title = "Cluster profiles relative to the NYC average",
    subtitle = "Index score: 100 = NYC overall mean",
    x = "Cluster",
    y = NULL
  ) +
  theme_minimal(base_size = 10.5) +
  theme(
    plot.title = element_text(face = "bold", size = 12),
    plot.subtitle = element_text(size = 10),
    panel.grid = element_blank(),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 9),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 8)
  )

# ------------------------------------------------------------------------------
# 30.9 拼接输出
# ------------------------------------------------------------------------------
p_map_nyc + p_heatmap_nyc +
  plot_layout(widths = c(1.15, 1.15))

至此,我们已经完成了一个基于 PCA 变量筛选 的纽约市 geodemographics 基础构建流程,并展示了其空间分布格局与各聚类相对于纽约市整体平均水平的Index 特征。在 5.0 总结与反思 环节中,我们将进一步结合这两张图,对不同聚类进行简要描述,并尝试完成初步命名。

5.0 专题研究 2: 回顾与反思

5.0 专题研究 2:回顾与反思

重要

本节将对前述专题研究 2“纽约市内部地理人口特征分类”进行简要回顾。需要说明的是,受教学场景与篇幅所限,该示例更适合作为一个方法演示型案例,而非完整展开的正式科研项目。不过,其所呈现的地理人口特征分类(geodemographic classification)工作流,确实反映了城市内部社会空间研究、小区域人口普查分析与邻里类型识别中较为常见的一类分析路径。

与前一专题侧重全国尺度、县级单元的入门式分类不同,本专题进一步转向纽约市 census tract 尺度的城市内部分类实践。在这一尺度下,变量数量更多、结构更复杂,邻里之间的差异也更加细致。因此,前文在基础变量重构之外,进一步加入了相关性检验主成分分析(PCA)变量筛选以及 H-K-means 聚类分析等步骤,以构建一个更适合城市内部社会空间分异识别的分类流程。

具体而言,前文已经依次完成了变量选择与数据获取指标重构与预处理相关性检验与 PCA 筛选H-K-means 聚类以及空间可视化与指数热力图展示等关键环节。通过这一过程,原本较为分散的多维人口、社会经济、住房与通勤变量,被进一步整理为若干具有解释潜力的街区类型,也初步展示了纽约市内部在社会人口背景上的空间分异格局。

方法延伸:如果想构建一个具有特定目的的 geodemographics,该怎么做?

需要说明的是,无论是前面的专题研究 1,还是本节的专题研究 2,其所构建的地理人口特征分类都更接近一种通用型(general-purpose)分类。也就是说,这类分类更强调从人口、社会经济、住房与通勤等多个维度,对地区进行较为综合的画像,因此通常更适合用于公共部门官方统计或一般性的社会空间研究。英国 OAC(Output Area Classification) 就是典型案例。

但在实际应用中,也常常需要构建一种更具特定目的定制型分类(bespoke classification)。例如,若研究目标转向某类消费行为通勤偏好住房市场细分商业选址,那么分类流程通常需要作出相应调整。大致可以从以下几个方面理解:

  • 研究目标不同
    通用型分类更强调地区的综合社会画像;而定制型分类通常围绕某个明确问题展开,例如某类消费方式、特定服务需求或某类市场行为。

  • 变量体系不同
    定制型分类往往不会平均纳入所有领域,而是更强调与目标问题高度相关的变量。例如,若研究某类消费人群,则可能更关注年龄、收入、家庭结构、出行方式、数字接入或居住环境,而不一定需要完整保留所有普查变量。

  • 外部数据可能需要加入
    相较于前述主要基于 ACS 普查变量的流程,定制型分类往往还需要叠加交易数据消费数据出行行为数据平台用户数据或其他专题数据,以增强分类对特定问题的针对性。

  • 评价标准不同
    通用型分类通常更强调解释性与整体社会空间结构;而定制型分类往往还要考虑它是否真正有助于识别目标人群、支持决策,或提高某类场景下的应用价值。

  • 流程上仍有共通之处
    即便是定制型分类,其基本工作流通常仍与前文类似,包括变量选择数据重构标准化处理降维或筛选聚类分析结果解释等步骤。不同之处主要在于:变量是否围绕特定目标组织,以及最终分类是否服务于一个明确的应用场景。

因此,若想构建一个真正“有用”的定制型 geodemographics,关键往往不在于使用了多么复杂的算法,而在于:分类是否围绕明确问题展开,变量是否真正反映了该问题的核心维度,以及结果是否具有足够的解释性应用针对性


专题研究 2 反思

重要

任何面向城市内部的地理人口特征分类流程,都不只是一个“把变量放进聚类模型”的技术过程,而是同时涉及变量构建数据变换相关性处理降维策略聚类算法选择结果解释等多个环节。前述案例虽然展示了纽约市 census tract 尺度下从 ACS 数据到城市内部地理人口特征分类的完整流程,但从更严格的研究角度看,相关结果仍需结合研究目标、变量体系与方法选择进行审慎理解。


以下几个问题,可作为对本案例的进一步反思:

  • 变量是否需要做分布变换?
    本案例在进入 PCA 与聚类分析之前,主要进行了标准化处理,但并未进一步对偏态较强的变量作正态化变换,例如 logBox-Cox 或其他变换方式。在实际研究中,这一步往往是一个需要权衡的问题:一方面,某些变量如人口密度、收入或租金,本身可能存在明显右偏分布,适当变换有助于缓解极端值影响;但另一方面,若过度变换,也可能改变变量原本的社会空间含义,削弱结果的直观解释性。因此,一个值得思考的问题是:为了让变量更“适合模型”,是否值得牺牲部分原始分布特征?不同变量是否都应采用同一种变换方式?

  • 高度相关的变量到底要不要删除?
    在本案例中,我们先通过相关性检验识别高度相关的变量对,并对其中部分变量进行了初步删减。但在更正式的 地理人口特征分类 研究中,这一步其实并没有唯一标准。某些研究会根据变量解释性来选择保留哪一个变量;也有研究会借助最小生成树(minimum spanning tree)等方式,尽量减少变量之间的冗余连接;当然,也有研究选择在一定程度上保留高度相关变量,因为这些变量本身可能分别代表不同的实质性领域。由此也引出一个值得继续思考的问题:相关性高是否就一定意味着“必须删除”?如果要删,依据应更多来自统计结构,还是来自变量本身的社会学含义?

  • PCA 应该只用于变量筛选,还是直接用于聚类?
    本案例将 PCA 作为一种变量初筛结构梳理工具,而没有直接将主成分得分输入聚类模型。但在 geodemographics 的发展历程中,特别是较早期的一些分类研究中,直接以主成分得分进行聚类其实是一种非常常见的做法。其优点在于可以有效降维、减少共线性;但缺点则在于聚类对象不再是原始变量,而是若干综合主成分,结果解释有时会变得更抽象。随着数据规模、算力条件与算法方法的发展,学界对“是否必须先做降维”这一问题也一直存在讨论。因此,一个值得思考的问题是:在算力允许的情况下,直接输入更多原始变量,是否反而更能保留城市内部差异?还是说,适度降维仍然有助于提高分类的稳定性与可解释性?

  • H-K-means 之外,还有哪些聚类方法可选?
    本案例采用 H-K-means,是因为它在一定程度上缓解了普通 K-means 对随机初始化较为敏感的问题,并更适合当前这种变量较多、结构较复杂的教学示例。但它并不是唯一可行的方法。实际研究中,还常见层次聚类模糊聚类高斯混合模型谱聚类,甚至一些面向大规模数据的机器学习聚类方法。不同方法在类别边界、异常值敏感性、类别形状假设与结果解释方式上各有特点。因此,还可以进一步思考:对于城市内部这类小区域、多维属性的分类问题,是否存在比 H-K-means 更合适的方法?若更换算法,纽约市内部的分类轮廓是否仍能保持一致?

  • 聚类数量 k 应该如何确定?
    即便 H-K-means 在初始化方面较 K-means 更稳健,聚类数量 k 的选择仍然是一个普遍存在的问题。本案例中已经结合 Elbow methodSilhouetteGap statistic 等方法进行辅助判断,但不同指标给出的建议并不完全一致。这一现象本身就说明,k 的确定通常不是纯粹“客观自动”的,而需要研究者在统计合理性分类解释性之间做出取舍。换言之,一个在指标上略优、但类别过细且难以命名的方案,未必比一个稍微粗一些、但结构更清晰的方案更适合实际研究。

  • 城市内部地理人口特征分类的结果是否足够稳定?
    与全国县级分类相比,城市内部 tract 尺度的分类往往更容易受到变量选择、空间单元边界与算法参数的影响。在更正式的研究中,通常还需要进一步进行稳健性检验,例如比较不同 k、不同变量组合、不同算法下结果的一致性,或检查某些类别是否反复出现。由此也可以进一步追问:当前分类结果究竟反映的是纽约市内部相对稳定的社会空间结构,还是更多受到本次建模设定本身的影响?

  • 分类结果的解释性与应用性是否兼容?
    本案例已经展示了聚类的空间分布图与 Index 热力图,并尝试对各类街区作出命名与描述。但在真实应用中,一个地理人口特征分类是否“好”,并不只取决于它能否把地区分开,还取决于它是否足够易解释易沟通,以及是否能够真正服务于特定研究问题或政策场景。例如,若分类过于依赖统计筛选而忽略了某些重要领域,结果可能在数学上成立,却未必适合用于城市规划、社区治理或公共政策分析。因此,一个值得继续思考的问题是:分类结果究竟应优先服务于模型表现,还是优先服务于实际解释与应用?

延伸阅读:PCA 与变量筛选

如果希望进一步了解 PCA地理人口特征分类中的实际应用,以及其基本原理与实现过程,可结合以下两类材料进行阅读。

研究应用

该研究提出了一种基于 主成分分析 的自动化变量筛选框架,重点讨论了如何利用 PCA 辅助识别更具代表性的分类变量,对于理解本节中“先进行 PCA,再筛选变量”的思路具有较强参考价值。

Liu, Y., et al. (2019). A principal component analysis (PCA)-based framework for automated variable selection in geodemographic classification.
DOI: 10.1080/10095020.2019.1621549

基础入门

如果想进一步理解 PCA 的底层逻辑,也可以阅读这篇入门教程。该页面以 iris 数据为例,演示了 prcomp() 的基本输出、主成分得分的生成方式,以及如何从协方差矩阵进一步理解特征值特征向量的含义。对于初学者来说,这类材料更适合帮助建立对 PCA 数学结构与结果解释的直观认识。

Introduction to PCA. In A Machine Learning Compilation
教程链接:https://f0nzie.github.io/machine_learning_compilation/introduction-to-pca.html

建议阅读时重点关注以下几个问题:

  • PCA 在这里扮演的是降维工具,还是变量筛选工具
  • 特征值、解释方差与主成分载荷之间是什么关系?
  • 自动化筛选与研究解释性之间应如何平衡?
  • 在不同空间尺度与变量体系下,PCA 的筛选结果是否稳定?

6.0 本章练习

6.0 动手实战

提醒

专题章节中的 6.0 部分通常不再提供完整的可复现模板。请先认真学习本章两个专题研究中的代码与分析流程,再完成下列课堂练习。

课堂练习

本部分不再提供现成代码。请结合本章已经学习的人口普查数据获取变量重构相关性检验PCA 筛选聚类分析等方法,自主完成一个小型地理人口特征分类专题研究。

本章已经涉及ACS 数据提取小区域变量构建地理人口特征分类K-means / H-K-means 聚类以及空间可视化与类型解释等内容。本练习的目标,不再是单一步骤的复现,而是围绕一个明确研究区域,将这些方法组织为一套完整的 geodemographics 分析流程。


6.1 任务说明

请基于人口普查或人口调查数据,完成一次自主的地理人口特征分类分析。建议优先使用美国 ACS(American Community Survey) 数据,并通过 tidycensus 获取。你的研究区域可以自行选择,例如:

  • 全美国范围 (如 county level)
  • 某一个州 (如 California / Texas / New York)
  • 某一个城市或都市区 (如 Chicago、Los Angeles、Houston)

从分析尺度上,建议优先选择以下两类单元:

  1. county:适合较大范围、流程较简单的分析;
  2. census tract:适合城市内部或区域内部更细颗粒度的分类。

6.2 为什么本章练习建议优先使用 ACS?

本章练习建议优先使用美国 ACS 数据,主要原因并不在于美国数据本身“更重要”,而在于它与 R 环境之间具有较好的工具适配性。借助 tidycensus,研究者可以较方便地完成变量检索、表格下载与空间边界调用,因此更适合作为课程中的入门实践数据。

当然,其他国家的人口普查数据同样可以用于类似分析。例如:

  • 英国 Census:属于开放数据,空间颗粒度也较细,非常适合开展小区域社会空间研究;但在 R 中的数据获取与整理过程通常不如 tidycensus 直接。
  • 中国人口普查数据:理论上也可以用于类似分类研究,但前提是能够获取同一年、同尺度、且具有足够多维度的公开数据。相较而言,中国精细小区域普查数据的可得性通常较弱,因此课程练习中实施难度更高。

因此,本章练习更强调你是否掌握了地理人口特征分类的核心流程,而不是要求必须使用某一个特定国家的数据。若你能自行整理英国、中国或其他国家的多维小区域普查数据,同样可以作为练习对象。


6.3 建议选题方向

方向 A:通用型地理人口特征分类
例如:对某一州或某一城市内部进行综合型分类,纳入年龄、族裔、教育、收入、住房与通勤等变量,构建一个 general-purpose geodemographics。

方向 B:城市内部社会空间分异
例如:选择某一个大城市,在 census tract 尺度下识别不同街区类型,并分析其空间分布格局。

方向 C:具有特定目的的分类尝试
例如:围绕住房市场、通勤方式、老龄化社区或年轻专业人群聚居区,构建一个更具针对性的变量体系,并比较其与通用型分类的差异。

你也可以自行拟定题目,但应确保研究目标变量体系地理单元聚类方法之间保持一致。


6.4 硬性要求(Checklist)

你的 R Markdown 报告至少应包含以下内容:

  1. 研究区域与问题界定
    说明你选择的研究区域、空间尺度,以及你希望回答的核心问题。
    例如:你是希望识别某州内部的综合人口类型,还是希望分析某个城市内部的社会空间分异?

  2. 数据来源与变量体系
    说明你使用的数据来源 (如 ACS 2024 5-year estimates),并展示你选取了哪些变量。
    需说明这些变量分别来自哪些表格,属于哪些主题 (如人口结构、社会经济、住房、通勤等)

  3. 变量重构过程
    展示如何从原始普查表项重构比例变量、密度变量或连续变量。
    例如:年龄占比、租住房屋占比、贫困率、人口密度等。

  4. 基础预处理
    至少包含以下两项中的一项:

    • 相关性检验:识别高度相关变量;
    • PCA:用于辅助理解变量结构或初步筛选变量。

    并简要说明:为什么要做这一步?它对后续聚类分析有何帮助?

  5. 聚类分析
    至少使用 1 种聚类方法 完成分类,例如:

    • K-means
    • H-K-means
    • 层次聚类

    同时说明:你如何确定聚类数量 k?是否参考了 Elbow method、Silhouette 或其他方法?

  6. 空间可视化输出
    至少输出 1 张聚类结果地图,展示不同类别在研究区域中的空间分布。
    若条件允许,也建议进一步输出热力图画像图指数对比图

  7. 结果解释与命名
    至少对 2–3 个主要类别 进行简要解释,并尝试赋予具有描述性的类型名称。
    例如:

    • “高龄化郊区自有住房型社区”
    • “年轻高学历租住型核心城区”
    • “低收入少数族裔公共交通依赖型街区”
  8. 反思与讨论
    简要讨论你的分类结果可能受到哪些因素影响,例如:

    • 变量选择是否足够合理?
    • 是否需要删除高相关变量?
    • 是否需要做数据变换?
    • 若改变空间尺度或聚类方法,结果会不会不同?

6.5 建议

建议你的 R Markdown 报告尽量保持与本章专题研究相似的结构,即:

  • 研究课题与分析目标
  • 数据选择与获取
  • 数据重构与预处理
  • 聚类分析与类型识别
  • 结果解释与反思

重点不在于你是否得到了“完美”的分类结果,而在于你是否能够清晰展示:你为什么这样选变量、为什么这样处理数据、为什么这样做聚类,以及你如何解释最后得到的类别。