1. 工具简介

1.1 开发背景

在人工智能深度介入医疗、地理与公共治理的今天，“模型为何如此预测”已变得与“预测是否准确”同等重要。可解释性模型正是回应这一追问的前沿工具。时空可解释性模型在此基础上更进一步，将影响因素、空间位置与时间维度纳入统一分析框架，回答“何地、何时、因何”的核心问题，兼具更高预测精度与模型透明度。

时空可解释分析工具“BSTVC桌面版”首次推出全界面化“零代码”建模功能，并依托三个全新真实健康案例，系统演示如何针对连续型（发病率、死亡率）、二值型（是否发病、是否死亡）和计数型（病例数、死亡数）三类典型结局变量，开展以下三类核心分析：

（1）局部时空可解释性：探测解释变量的时空异质性影响机制，刻画关键因素在具体时空单元的作用差异；

（2）全局时空可解释性：量化解释变量的时空贡献百分比，识别时空维度上的关键因素；

（3）时空动态预测：融合时空可解释先验知识，实现缺失值填补、未来预测与时空平滑。

BSTVC以此构建“局部刻画—全局识别—精准预测”的时空解释闭环，一站式赋能地理时空研究与决策。

相比于R版本，桌面版更加简洁易用，既能满足专业人士的深度需求，又降低了贝叶斯复杂建模的门槛，使更广泛的用户群体能够轻松应用先进的贝叶斯局域时空回归分析方法，解析和解释复杂的时空面板数据。适用于涉及地理时空数据分析的各类自然与人文科学领域，包括但不限于公共卫生、医学地理、环境健康、卫生经济和社会医学等学科。

1.2 特色功能

面向多种目标变量：支持三种主流目标变量类型：连续型（log-Gaussian）、二分类（logistic）和计数型（Poisson），满足不同分析场景的需求。

探测时空异质影响机制（局部时空可解释性）：通过拟合时空回归系数，揭示解释变量（X）与目标变量（Y）之间的局域时空差异，深入分析“因地制宜、因时制宜”的规律，探索时空异质性视角下的影响机制。

明确时空驱动因素（全局时空可解释性）：在识别时空异质影响机制的基础上，通过计算时空可解释百分比，明确关键驱动因素，为地理时空归因提供有力证据。

提升时空预测精度：考虑局域变量关系的时空异质性，显著提高模型拟合度和预测精度，用于时空缺失值填补、时空平滑和未来预测等。

贝叶斯模型评价：提供贝叶斯模型的全面评估，包括模型拟合度（DIC、WAIC）、复杂度（pd）与预测精度（LS）等指标，帮助用户全面了解模型性能。

丰富的可视化输出：同步提供多种时空可视化工具和代码，帮助用户直观理解模型结果，增强数据分析的可解释性，推动您的应用研究创新。

1.3 理论基础

贝叶斯时空变系数（BSTVC）模型是地理学领域新近发展的时空可解释基准方法，已配套开源R工具包（Song and Tang, 2025），提供了一个完整统一的“全地图”地理建模框架来精准捕捉具有时空差异的变量关系，旨在揭示多源解释变量对目标变量的时空异质影响机制，即时空非平稳效应。

时空异质耦合分析同时考虑了基于两个地理学定律的时空异质性和时空自相关，并在同一框架内结合时间和空间维度进行耦合分析。这种方法是当前地理数据分析中最全面、精准且具有最强证据力的分析理念，远超传统的非空间分析、单一空间或时间分析，以及时空分离分析方法。BSTVC作为时空异质耦合分析的前沿最新工具，其统计机理主要包括三类：

贝叶斯时空变系数（Bayesian Spatiotemporally Varying Coefficients，BSTVC）模型是一类基于贝叶斯统计内核的局域时空回归分析方法（Song et al.,2019;2020;2022），其显著优势在于，利用“全地图”单独建模框架对所有局部回归系数的时空变化进行统一拟合，从而能够精准捕捉解释变量对目标变量的时空异质性影响，即揭示时空非平稳性。BSTVC系列模型（目前包括STVI、STIVI、STVC、STIVC四类子模型）为揭示目标变量的复杂时空动态变化和时空影响机制提供了强有力的工具（Song et al., 2022）。

贝叶斯空间变系数（Bayesian Spatially Varying Coefficients，BSVC）模型是BSTVC模型的空间维度精简版，仅用于识别具有空间异质性的变量关系，即空间非平稳效应。其优势在于集成了BSTVC的“全地图”单独建模框架，保证了拟合的局域空间回归系数具有直接可比性。空间非平稳和时空非平稳是地理学第二定律的重要内涵。此外，无论是BSTVC还是BSVC，在拟合非平稳性的时候，均考虑了基于地理学第一定律的时空自相关和空间自相关特征。

时空方差分割指标（Spatiotemporal Variance Partitioning Index，STVPI）基于BSTVC/BSVC建模结果，通过量化并比较不同时空异质影响因素的可解释百分比（时空贡献度/时空相对重要性）来明确关键驱动因素（Song et al., 2022; Wan et al., 2022）。与当前主流的依赖绝对评价指标进行因子重要性排序的方法不同，STVPI是一种相对评价指标，因此能够为地理时空归因提供重要的先验依据（Wan et al.,2022）。

BSTVC方法优势：

“全地图”单独建模框架：该框架通过完整统一的贝叶斯层次建模机制，确保了局部时空回归系数的直接可比性，同时具备极强的拓展性，能够应对更复杂的实际应用挑战。正是由于该框架的设计，才使得计算相对时空贡献度（STVPI）成为可能。相比之下，在频率统计体系下，类似的分析通常采用“局部分开建模”的方式，即针对每个地图单元单独建模后再进行组合，可能引发一系列问题，例如不同小模型之间的可比性不足。

参数不确定性：无论是目标变量的局域预测，还是时空回归系数的局域拟合，BSTVC都能直接输出参数的不确定性评估结果，包括两种贝叶斯可信区间（50%和95%），以宽窄区间形式呈现。而在频率统计方法中，受限于统计机理，类似的分析无法评估不确定性。

缺失值友好：即使目标变量Y存在较多时空缺失值，或解释变量X存在少量缺失值，也不影响时空非平稳效应的识别。而在频率统计方法中，类似的分析通常不支持缺失值。

更多空间权重矩阵的支持：不仅支持基于距离的和k临近的空间权重矩阵，还支持10邻接矩阵来刻画数据中的空间自相关效应。而在频率统计方法中，类似的分析通常不支持10邻接矩阵。

1.4 参考文献

• Song, Chao, Yin, Hao, Shi, Xun, Xie, Mingyu, Yang, Shujuan, Zhou, Junmin, Wang, Xiuli, Tang, Zhangying, Yang, Yili, & Pan, Jay. (2022). Spatiotemporal disparities in regional public risk perception of COVID-19 using Bayesian Spatiotemporally Varying Coefficients (STVC) series models across Chinese cities. International Journal of Disaster Risk Reduction, 77, 103078.
• Wan, Qin, Tang, Zhangying, Pan, Jay, Xie, Mingyu, Wang, Shaobin, Yin, Hao, Li, Junmin, Liu, Xin, Yang, Yang, & Song, Chao. (2022). Spatiotemporal heterogeneity in associations of national population ageing with socioeconomic and environmental factors at the global scale. Journal of Cleaner Production, 373, 133781.
• Song, Chao, Shi, Xun, & Wang, Jinfeng. (2020). Spatiotemporally Varying Coefficients (STVC) model: a Bayesian local regression to detect spatial and temporal nonstationarity in variables relationships. Annals of GIS, 26(3), 277-291.
• Song, Chao, Shi, Xun, Bo, Yanchen, Wang, Jinfeng, Wang, Yong, & Huang, Dacang. (2019). Exploring Spatiotemporal Nonstationary Effects of Climate Factors on Hand, Foot, and Mouth Disease Using Bayesian Spatiotemporally Varying Coefficients (STVC) Model in Sichuan, China. Science of The Total Environment, 648, 550-560.
• Song, Chao, Tang, Xianteng. (2025). User’s Guide for the BSTVC R Package.

2. 安装与界面功能

详细操作步骤说明：

BSTVC主页面介绍：

该程序页面主要由左右两大部分组成。左侧为功能导航栏，用于切换软件的各个功能模块，自上而下可分为数据、建模、辅助功能以及系统信息四个板块；右侧为功能面板，用于展示当前模块的具体内容，并进行数据上传、参数设置、模型运行和结果查看等交互操作。

数据板块主要用于完成建模前的数据准备工作，包括“数据输入”和“数据检查”两个功能。“数据输入”用于导入分析所需的属性数据、空间地图数据等基础文件；“数据检查”用于检查输入数据的结构、字段和时空顺序是否满足建模要求，帮助在正式建模前发现并修正数据格式、空间单元匹配或时间顺序方面的问题。

建模板块是软件的核心分析区域，包括“BSTVC时空建模”和“BSVC空间建模”两个功能。“BSTVC时空建模”面向连续时序的时空面板数据，用于分析解释变量与响应变量之间关系在时间和空间上的变化特征，识别不同因素作用的时空异质性及其动态演变规律；“BSVC空间建模”则主要面向单一时点或空间截面数据，用于分析解释变量作用在不同空间单元之间的差异，揭示变量关系的空间异质性和局部驱动模式。

辅助功能板块用于支持数据整理和建模配置，包括“数据转换”“字段转换”和“自定义空间矩阵”等功能。“数据转换”用于将包含多个时间点的空间截面数据整理为软件建模所需的时空面板数据格式；“字段转换”用于对数据字段进行规范化处理，使字段名称和变量格式满足软件识别与建模要求；“自定义空间矩阵”用于导入或构建用户自定义的空间权重矩阵，以满足不同研究对象、邻接关系或空间距离设定下的建模需求。

系统信息板块主要包括“关于”功能，用于展示软件信息、联系信息、开发说明及相关链接等内容，方便用户了解软件来源、功能定位和更新信息。

3. Q&A/联系我们

常规情况说明：

• 文件上传限制：由于底层系统框架默认文件上传大小较小，而空间面板数据和地图文件通常体积较大，因此本系统已将文件上传上限设置为1GB。
• 计算复杂度：BSTVC和BSVC均属于贝叶斯局部建模方法。模型计算量会随着空间单元数量、时间点数量、解释变量数量、响应变量类型以及线程设置的增加而增加。因此，运行较大规模数据时需要计算机具备足够的内存和计算性能。
• 空间权重矩阵：在使用k近邻空间权重矩阵时，可能出现报错或无输出结果的情况，这可能与底层贝叶斯推断算法INLA的实现机制有关。用户如遇到该问题，可尝试更换空间权重矩阵类型或调整矩阵构建方式。
• 残差项说明：对于泊松型响应变量和逻辑回归型响应变量，当前版本暂未包含残差项。因此，在这两类模型中，STVPI暂时无法计算整体模型的解释百分比。该功能将在后续版本中进一步完善。
• 非平稳性假设：当前版本的BSTVC时空建模模型默认所有解释变量均具有时空非平稳性。后续版本计划支持更灵活的模型设定，使用户能够指定部分变量不具有时间非平稳性。

联系信息：

如用户在软件使用过程中遇到理论方法、模型理解或统计解释方面的问题，可联系宋超：chaosong.gis@gmail.com；如用户在软件安装、系统操作、数据导入、功能使用或运行报错等方面需要帮助，可联系唐先腾：tangxt.me@gmail.com。此外，用户也可通过以下网址：https://github.com/bayesianstvc/BSTVC-R与https://chaosong.blog/bayesian-stvc/，查看项目说明、版本更新，了解BSTVC模型背景、方法介绍及相关研究资料。另外，程序页面提供“医学地理信息与空间健康统计”官方公众号二维码。用户可关注该公众号，获取医学地理信息、空间健康统计方法以及BSTVC模型相关更新内容。最后，欢迎用户通过上述邮箱与开发者联系，反馈使用过程中遇到的问题、提出功能需求或参与程序改进。

4. 实证案例I：连续型目标变量

4.1 示例数据

本案例以全球健康调整寿命HALE（连续型结局变量）长时序时空面板数据为例进行演示，使用的数据来源于GBD、GIOVANNI、World Bank及空气污染相关数据库，收集了2000到2020年间177个国家和地区的社会经济、自然环境、健康行为和空气污染等因素的真实数据变量，具体数据介绍如表4.1所示。

注意：时空面板数据包含时间和空间两个维度，同一空间单元有多年重复观测；空间截面数据仅包含空间维度，所有观测来自同一时间点。本案例原始数据属于时空面板数据。

表4.1 案例Ⅰ示例数据

4.2 数据输入与验证

本教程以案例一示例数据为例进行数据导入，用户可以根据自己的数据进行相应操作（注：本案例已经为时空面板数据，若需要查看空间截面数据转换详见案例二）。

步骤说明：

在建模过程中，空间权重矩阵用于描述不同空间单元之间的邻近关系或空间联系强度。一般情况下，如果用户没有特殊研究需求，可直接使用系统默认的空间权重矩阵，即在建模页面的权重矩阵参数中不需要额外输入文件，系统会自动按照默认规则构建空间关系。

需要注意的是，如果用户使用的是点类型shp地图数据，系统默认的邻接型空间权重矩阵可能无法充分反映点之间的空间关系，甚至可能得到不理想的结果。此时，建议用户根据研究目的自定义空间权重矩阵，例如选择距离型矩阵、k近邻矩阵或其他合适的空间关系设定。

本软件提供了“自定义空间矩阵”辅助工具，如下图。该页面左侧为参数设置区，右侧为信息展示区。用户首先需要导入地图数据，导入时应选择完整的shapefile文件，而不能只导入shp文件本身。读取地图后，右侧会显示系统识别到的空间单元、几何类型、坐标系类型和范围等信息。

随后，用户可在左下方“空间矩阵参数”中选择需要构建的空间矩阵类型，并设置相应参数和权重标准化方式。该步骤需要用户具备一定的空间权重矩阵相关知识，以便根据研究对象和空间关系特征选择合适的矩阵类型。点击“构建空间矩阵”后，系统会生成对应的空间权重矩阵，并在右侧进行空间关系可视化展示。构建完成后，用户应下载并保存生成的rds文件，后续建模时即可将该文件作为自定义空间权重矩阵导入使用。

【备注】在BSTVC支持的时空面板数据中，变量Y可以包含缺失值，变量X也可容忍少量缺失值。在用户导入自己的数据前，需注意：①缺失值请使用“NA”表示，便于识别；②检查表格数据与地图数据中的空间唯一标识码的字段类型，确保两个表中的该字段的名称和类型一致。

4.3 数据检查

4.3.1 空间单元顺序检查

建模表格数据中的每个时间节点下空间单元的顺序需要与建模地图数据中的shapefile文件中空间单元的顺序完全一致，不一致会造成最终结果不正确，但模型运行并不会报错的情况。例如，地图TestMap_Continuous中代表国家与地区单元的ISO_A3字段的顺序为FJI，TZA，ESH，CAN，USA，则TestData_Continuous.csv数据中的ISO_A3字段也必须为此排序。BSTVC建模系统中提供一键式自动比对数据表与地图中的空间单元唯一值字段的功能，完成空间序列的自动化重排与对齐。

步骤说明：

4.3.2 字段类型检查

在建模前，用户需要特别注意字段类型问题。当原始数据中存在缺失值，并使用“NA”表示时，部分CSV或Excel文件在读取过程中可能会将原本应为数值型的变量识别为字符型字段。如果该字段后续被作为解释变量或响应变量参与建模，可能导致模型运行时出现指向不明确或难以判断原因的报错。因此，在建模前建议用户先检查关键变量的字段类型是否正确。

针对这一情况，软件也提供了“字段转换”辅助工具，用户无需再到其他软件中单独处理数据。使用时，用户只需导入原始数据文件并点击“读取字段转换数据”，随后在“需要转换类型的字段”中选择一个或多个需要处理的变量，再在“目标字段类型”中选择转换后的字段类型，例如数值型numeric，最后点击“开始字段转换”即可。

转换结果会显示在右侧“字段转换结果”面板中。右侧上半部分用于展示原始数据表格的字段属性，下半部分用于展示转换后的新表格字段属性，方便用户对比转换前后的变化。确认无误后，用户可通过下方下载按钮将转换后的数据保存到本地。

4.4 BSTVC时空建模

构建BSTVC时空模型所需的参数主要有9个，具体参数描述如下：

• 数据来源：如果上一步数据检查显示顺序已经一致，请选择“原始数据”；如果检查提示已完成自动重排，请选择“检查后数据”！
• 响应变量Y：模型中的因变量，即用户主要研究和尝试解释的核心目标变量。例如，本案例中的健康调整预期寿命（HALE）。
• 解释变量X：模型中的自变量，即可能对响应变量Y产生影响的因素。本平台支持选择一个或多个变量。例如，本案例中的高体重指数、当前卫生支出占GDP比重等。
• 响应类型：响应变量Y的数据分布形态。本平台提供的选项包括连续型continuous、二分类binary以及计数型count，需要根据响应变量Y的实际数据特征来选择。
• 时间字段Time：数据集中代表时间维度的列。它帮助模型识别数据的多时间点属性，从而分析时间序列上的变化。
• 空间字段Space：数据集中代表空间维度的列。需确保选择的字段与建模表格数据中的列名完全一致，且建模表格数据与建模地图数据中应同时具有名为此字符串的字段。例如，在本案例中，空间字段Space为177个国家及地区的唯一标识ISO_A3字段。
• 解释变量标准化：决定是否对解释变量X进行无量纲化处理。我们强烈推荐选择“是”，因为这不仅能确保所有指标在同一尺度上进行比较，有助于衡量各个变量对模型的贡献大小，还能加速模型的计算过程。
• 空间权重矩阵：默认为QUEEN规则下的二分式（具有邻接关系的赋值为1，否则赋值为0）空间权重矩阵。用户可输入自行构建好的空间权重矩阵文件，且支持k临近、反距离、固定距离等多类型的空间权重矩阵；若用户未输入，则本平台将构建默认矩阵。
• 线程数threads：设定分配给该模型运算的CPU线程数量，默认为6。由于时空模型的参数估计通常计算量庞大，在硬件支持的前提下，合理设置多线程可有效缩短运行时间。

步骤说明：

4.5 模型结果

BSTVC函数输出的结果共包括6个部分，具体输出部分描述如下表

表4.2 BSTVC模型输出说明

BSTVC的输出结果中的每一部分都是一个sheet，用户可以直接下载并保存以及用于后续的结果可视化。

4.6 BSTVC模型输出及其可视化

在本文档中，我们将详细提供一系列的绘制模型输出结果的步骤和代码示例，以便您能够轻松地在R环境中创建各种图表。我们的目标是让您不仅能够理解模型的输出，还能够通过图表的形式，将这些复杂的数据以更直观、更吸引人的方式呈现出来。

本文档在此部分展示的所有图表和可视化结果均基于示例数据构建的模型的输出。以下绘图代码及图件供参考。

4.6.1 基础R包配置与数据读入

 # ========================================================= # 1. 基础R包配置与数据选取读入 # ========================================================= # 加载所需的全部 R 包 library(readxl) library(openxlsx) library(ggplot2) library(dplyr) library(ggthemes) library(tidyr) library(sf) cat("请在弹窗中选HALE_all_result_tables.xlsx文件 (包含 model.evaluation, STVPI, time.coefficients 等 sheet)...\n") excel_file <- file.choose() cat("请在弹窗中选择TestMap_Continuous.shp地图矢量文件 (.shp)...\n") map_data <- st_read(file.choose()) # 提前读取 STVPI 数据以检测所有影响因素名称 (用于后续所有模块自动化) STVPI.data <- read.xlsx(excel_file, sheet = "STVPI") colnames(STVPI.data) <- c("id", "STVPI_2.5", "STVPI_25", "STVPI_mean", "STVPI_75", "STVPI_97.5", "effect_name") STVPI.data$effect_name <- trimws(STVPI.data$effect_name) # 定义全局背景变量并自动剥离前缀 overall_vars <- c("space.time.all", "residual", "space.all", "time.all") all_raw_names <- unique(STVPI.data$effect_name) factor_only_names <- all_raw_names[!all_raw_names %in% overall_vars] detected_factors <- unique(gsub("space\\.time\\.|space\\.|time\\.", "", factor_only_names)) detected_factors <- detected_factors[detected_factors != ""] cat("\n======================================================\n") cat("【系统自动诊断】成功从您的表格中检测到以下", length(detected_factors), "个影响因素：\n") print(detected_factors) cat("======================================================\n\n") 

4.6.2 模型的贝叶斯评价

模型的贝叶斯评价包括模型名称（包含建模起始时间）、DIC、LS、WAIC、pd2等指标。其中，DIC、WAIC表示模型拟合度，eff（DIC）和pd2（WAIC）表示模型复杂度，LS表示预测精度。具体而言，DIC和WAIC值越小表示模型拟合度越好，eff和pd2值越小表示模型复杂度越低，LS值越接近0表示预测精度越高。

 # ========================================================= # 2. 模型的贝叶斯评价 # ========================================================= cat(">>> 正在处理: [1/4] 模型的贝叶斯评价...\n") # 读取 Excel 文件中指定的 Sheet model.eval <- read_excel(path = excel_file, sheet = "model.evaluation") cat("\n--- 模型评价结果 ---\n") print(model.eval) # 转置数据框进行指标垂直展示 cat("\n--- 转置后的模型评价结果 ---\n") model.eval_t <- t(model.eval) print(model.eval_t) 

结果如下：

4.6.3 时间回归系数（时间非平稳）

BSTVC模型不仅可以计算时间回归系数来刻画变量关系的时间异质性，还可以估算出每个回归系数值的贝叶斯宽（95%）、窄（50%）可信区间，用于直接评价结果的不确定性。

绘制方案一：

 # ========================================================= # 3. 时间回归系数（时间非平稳）的两种绘制方式 # ========================================================= cat("\n>>> 正在处理: [2/4] 时间回归系数 (两种绘制方式)...\n") Time.Coef_raw <- read_excel(excel_file, sheet = "time.coefficients") # --- 方式一：时间回归系数演化趋势图 (分面、带置信区间) --- Time.Coef_1 <- Time.Coef_raw %>% filter(explain_variable %in% detected_factors) %>% mutate(Coefficients = round(mean, 2), Year = time_index + 1999) p_time_1 <- ggplot(data = Time.Coef_1, aes(x = Year, y = Coefficients, fill = explain_variable)) + geom_ribbon(aes(ymin = `0.025quant`, ymax = `0.975quant`, fill = explain_variable), alpha = 0.1) + geom_ribbon(aes(ymin = `0.25quant`, ymax = `0.75quant`, fill = explain_variable), alpha = 0.5) + geom_line(aes(x = Year, y = Coefficients, colour = explain_variable), alpha = 0.8) + facet_wrap(~explain_variable, scale = "free_y", ncol = 3) + theme_few() + scale_x_continuous(breaks = seq(2000, 2020, by = 2)) + labs(title = "各因素的时间演变趋势", x = "年份", y = "时间系数 (TCs)") + theme( legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 10), strip.text = element_text(size = 12, face = "bold"), plot.title = element_text(hjust = 0.5, face = "bold", size = 15, margin = margin(b = 15)) ) # 在 R 中实时显示图表 print(p_time_1) # 保存本地 ggsave(filename = "Time_Coefficients_Method1_Faceted.png", plot = p_time_1, width = 12, height = 8, dpi = 300) 

绘制方案二（为避免对时间变化形式作出线性假设，本例采用LOESS局部加权回归对年度时间系数进行平滑展示）：

 # --- 方式二：时间不平稳线平滑图 (Loess集中拟合) --- start_year <- 2000 Time.Coef_2 <- Time.Coef_raw %>% filter(explain_variable %in% detected_factors) %>% mutate(Year = start_year + time_index - 1) p_time_2 <- ggplot(data = Time.Coef_2, aes(x = Year, y = mean, color = explain_variable)) + geom_smooth(method = "loess", se = FALSE, span = 1, linewidth = 1.2) + labs(x = "Year", y = "Time_Coefficients", color = "Variables") + scale_x_continuous(breaks = seq(min(Time.Coef_2$Year, na.rm = TRUE), max(Time.Coef_2$Year, na.rm = TRUE), by = 5)) + theme_classic() + theme( legend.position = "right", legend.title = element_text(face = "bold", size = 11), legend.text = element_text(size = 10), axis.text = element_text(size = 11, color = "black"), axis.title = element_text(size = 13, face = "bold"), panel.border = element_rect(color = "black", fill = NA, linewidth = 1), plot.title = element_blank(), aspect.ratio = 1.5 ) # 在 R 中实时显示图表 print(p_time_2) # 保存本地 ggsave(filename = "Time_Coefficients_Method2_Smooth.png", plot = p_time_2, width = 4, height = 6, dpi = 600, bg = "white") 

4.6.4 空间回归系数（空间非平稳）

与时间回归系数类似，BSTVC模型不仅可以计算空间回归系数表示空间异质的变量关系，还可以估算其不确定性，即贝叶斯宽（95%）、窄（50%）可信区间。

提取出模型输出中的空间回归系数数据框后，建议在ArcGIS或ArcGIS Pro等专业地图制图软件中进行绘制。此处基于R语言，简单绘制出基础的空间回归系数地图展示结果，仅供样本制图参考。

 # ========================================================= # 4. 空间回归系数（空间非平稳）的绘制 # ========================================================= cat(">>> 正在处理: [3/4] 空间回归系数地图...\n") space.coef <- read.xlsx(excel_file, sheet = "space.coefficients") Space.Coef.merge <- merge(map_data, space.coef, by = "ISO_A3", all.x = TRUE) # 自动生成所有需要提取的列名 expected_space_cols <- paste0(detected_factors, "_mean") actual_space_cols <- intersect(expected_space_cols, colnames(Space.Coef.merge)) Space.coef.panel <- Space.Coef.merge %>% pivot_longer( cols = all_of(actual_space_cols), names_to = "variable", values_to = "SCs" ) %>% mutate(variable = gsub("_mean", "", variable)) %>% select(ISO_A3, variable, SCs, geometry) p_space <- ggplot() + geom_sf(data = Space.coef.panel, aes(fill = SCs), color = "white", size = 0.1) + facet_wrap(~variable, shrink = FALSE, drop = FALSE, ncol = 3) + scale_fill_viridis_c(option = "viridis", name = "空间回归系数 (SCs)") + theme_few() + labs(title = "HALE影响因素的空间分布 (2000-2020)", x = "", y = "") + theme( legend.position = "bottom", legend.key.width = unit(2, "cm"), axis.text = element_blank(), axis.ticks = element_blank(), strip.text = element_text(size = 12, face = "bold"), plot.title = element_text(hjust = 0.5, face = "bold", size = 15, margin = margin(b = 15)) ) # 在 R 中实时显示图表 print(p_space) # 保存本地 ggsave(filename = "Space_Coefficients_Map_Auto.png", plot = p_space, width = 14, height = 10, dpi = 300) 

4.6.5 时空贡献度评价

时空方差分割指标（Spatiotemporal Variance Partitioning Index，STVPI）旨在量化并比较不同时空异质影响因素的可解释百分比，通过计算时空贡献度来明确关键驱动因素。基于BSTVC建模结果，STVPI进一步将每个解释变量的总方差分解为时间和空间两个独立的组分，即时间非平稳随机效应和空间非平稳随机效应，进而揭示数据的时空变异来源。

此工具的优点有：①不仅能识别解释因子的时空总体贡献，还能分别识别时间和空间维度的贡献；②相比只能识别绝对贡献排序的主流方法（如随机森林、SHAP等），STVPI识别的是相对贡献（可解释百分比），为地理时空归为提供重要证据基础；③直接评价其不确定性（贝叶斯宽窄可信区间）

 # ========================================================= # 5. 时空贡献度评价 (STVPI) # ========================================================= cat(">>> 正在处理: [4/4] 时空贡献度评价 (STVPI)...\n") # 自动生成匹配的完整名称列表 space_vars <- paste0("space.", detected_factors) time_vars <- paste0("time.", detected_factors) st_vars <- paste0("space.time.", detected_factors) my_vars <- c(overall_vars, space_vars, time_vars, st_vars) plot.data <- STVPI.data %>% filter(effect_name %in% my_vars) %>% mutate(Group = case_when( effect_name %in% overall_vars ~ "Overall", effect_name %in% c("space.all", "time.all") ~ "Space & Time", effect_name %in% space_vars ~ "space", effect_name %in% time_vars ~ "time", effect_name %in% st_vars ~ "space - time" )) %>% filter(!is.na(Group) & Group != "NA") # 逆转展示顺序 plot.data$effect_name <- factor(plot.data$effect_name, levels = rev(my_vars)) p_stvpi <- ggplot(plot.data, aes(x = STVPI_mean, y = effect_name, group = Group)) + geom_errorbar(aes(xmin = STVPI_2.5, xmax = STVPI_97.5, colour = effect_name), width = 0, linewidth = 1.2, alpha = 0.3) + geom_errorbar(aes(xmin = STVPI_25, xmax = STVPI_75, colour = effect_name), width = 0, linewidth = 1.2, alpha = 0.6) + geom_point(aes(colour = effect_name), shape = 21, fill = "white", size = 2, stroke = 1) + # 自动标值层：修改为百分数形式（保留1位小数） geom_text(aes(label = sprintf("%.1f%%", STVPI_mean * 100)), vjust = -1.5, size = 3.2, color = "black", fontface = "bold", show.legend = FALSE) + labs(title = "各因素的时空贡献度 (STVPI) 及具体数值", x = "", y = "") + facet_wrap(~Group, scale = "free", ncol = 3) + theme_few() + scale_x_continuous(labels = scales::percent_format(), breaks = scales::breaks_pretty(n = 4)) + theme( legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 10), panel.spacing = unit(1.5, "lines"), strip.text = element_text(size = 12, face = "bold"), plot.title = element_text(hjust = 0.5, face = "bold", size = 15, margin = margin(b = 15)) ) # 在 R 中实时显示图表 print(p_stvpi) # 保存本地 ggsave(filename = "STVPI_Auto_Labeled.png", plot = p_stvpi, width = 12, height = 8, dpi = 300) # ========================================================= # 完成提示 # ========================================================= cat("\n====== 所有绘图任务已全自动匹配并成功完成！高清图片已全部保存至以下文件夹：======\n") cat(getwd(), "\n") 

4.6.6 目标变量的时空预测

对于数据中目标变量Y的所有缺失值或非缺失值，局部预测结果中都会输出预测值，以及预测值的贝叶斯宽窄可信区间，用于直接评价不确定性。在实证案例Ⅰ中，采用了12个解释变量，分别是Y3、HB1、HB2、HB6、SE3、SE5、SE6、SE7、SE9、NE2、NE3、PM2.5，其他参数不变，进行上述BSTVC建模操作，下载全部表格保存所有结果，保存后命名为HALE12_all_result.xlsx，最终预测结果和散点图如下：

 # --------------------------------------------------------------------- # 1. 加载所有必要的 R 包 # --------------------------------------------------------------------- library(readxl) # 用于读取 Excel 预测数据 library(ggplot2) # 核心绘图包 library(dplyr) # 核心数据清洗与处理包 library(ggthemes) # 提供图表主题 (theme_few) library(sf) # 核心空间数据处理包 library(viridis) # 提供学术级连续型色带 # --------------------------------------------------------------------- # 2. 数据读取与预处理 # --------------------------------------------------------------------- cat("\n👉 步骤 1/2: 请在弹出的窗口中选择HALE12_all_result_tables.xlsx\n") excel_file <- file.choose() predict_df <- read_excel(excel_file, sheet = "local.prediction") # 查看输出字段结构，确保包含 'y' 和 'predict_y' cat("\n--- Excel 数据结构预览 ---\n") str(predict_df) cat("\n👉 步骤 2/2: 请在弹出的窗口中选择TestMap_Continuous.shp\n") shp_file <- file.choose() world_map <- st_read(shp_file) 

其中，对于变量Y所有类型（连续型、计数型与二值型）的局部预测结果，输出数据框中Year与ISO_A3为时间字段与空间字段，y字段代表目标变量Y的原始值，fill_y字段代表对目标变量Y中缺失值进行预测填补后的完整字段，predict_y代表为所有单元的模型预测值。

【备注】针对连续型的变量Y，由于模型设置的是log-Gaussian先验分布，目标变量Y会自动进行对数变换，故我们对输出结果进行了再变换，并设置“predict_”前缀。因此，(predict_0.025quant,predict_0.975quant)即为预测值的宽可信区间，(predict_0.25quant,predict_0.75quant)即为预测值的窄可信区间；

此外，基于局部预测的结果，可进行主流的预测精度（如散点图、R2、RMSE等）的评价。下面以连续型建模结果为例，我们给出绘制散点图、计算R2、RMSE以及绘制预测后的Y的时空分面可视化地图的代码：

 # --------------------------------------------------------------------- # 3. 预测精度评估：散点图与指标计算 # --------------------------------------------------------------------- cat("\n--- 正在处理数据并绘制精度评估散点图 ---\n") # 数据清洗：剔除带有缺失值 (NA) 的行 predict_normal_clean <- predict_df %>% filter(!is.na(y) & !is.na(predict_y)) # 计算评估指标：决定系数 (R²) 和 均方根误差 (RMSE) R2_normal <- 1 - ((sum((predict_normal_clean$predict_y - predict_normal_clean$y)^2)) / (sum((mean(predict_normal_clean$y) - predict_normal_clean$y)^2))) RMSE_normal <- sqrt(mean((predict_normal_clean$predict_y - predict_normal_clean$y)^2)) # 控制台打印结果 cat(paste("\n计算结果 -> R²:", round(R2_normal, 4), " | RMSE:", round(RMSE_normal, 4), "\n")) # 构造要显示在图表右下角的文本内容 text_normal <- paste0("R² = ", round(R2_normal, 4), "\nRMSE = ", round(RMSE_normal, 4)) # 绘制高清晰度散点图 p_scatter <- ggplot(predict_normal_clean, aes(x = y, y = predict_y)) + geom_abline(slope = 1, intercept = 0, color = "#e74c3c", linewidth = 1, linetype = "dashed") + geom_smooth(method = "lm", se = FALSE, color = "#2980b9", linewidth = 1) + geom_point(aes(fill = y), size = 3, shape = 21, color = "gray40", alpha = 0.8, stroke = 0.5) + scale_fill_gradient(low = "#88d8db", high = "#b8aeeb") + annotate("text", x = Inf, y = -Inf, label = text_normal, hjust = 1.1, vjust = -0.5, family = "serif", fontface = "bold", size = 5) + labs(x = "Measured Values", y = "Predicted Values", title = "Scatter Plot of Global Predictions") + theme_few() + theme( text = element_text(family = "serif", face = "bold"), axis.text = element_text(family = "serif", size = 12, face = "bold"), plot.title = element_text(hjust = 0.5, size = 15, margin = margin(b = 15)), legend.position = "none" ) # 在 R 环境中先行展示散点图 cat("\n📊 正在 R 绘图窗口中展示散点图...\n") print(p_scatter) 

 # --------------------------------------------------------------------- # 4. 时空分面地图：关键时间节点可视化 (2000-2020) # --------------------------------------------------------------------- cat("\n--- 正在生成时空分面地图 ---\n") # 数据拼接与年份筛选 map_data <- world_map %>% left_join(predict_df, by = "ISO_A3") %>% filter(!is.na(Year)) %>% filter(Year %in% c(2000, 2005, 2010, 2015, 2020)) # 绘制时空分面地图 p_map <- ggplot(data = map_data) + geom_sf(aes(fill = predict_y), color = "white", linewidth = 0.1) + coord_sf(expand = FALSE) + facet_wrap(~ Year, ncol = 2) + scale_fill_gradientn( colors = c("#313695", "#abd9e9", "#ffffbf", "#fdae61", "#a50026"), na.value = "gray90" ) + labs( title = "Spatio-Temporal Distribution of Predicted Values", fill = "Predicted Y\n" ) + theme_void() + theme( text = element_text(family = "serif", face = "bold"), plot.title = element_text(hjust = 0.5, size = 22, margin = margin(t = 10, b = 25)), strip.text = element_text(size = 18, margin = margin(t = 10, b = 15)), panel.spacing.x = unit(1.5, "lines"), panel.spacing.y = unit(2.0, "lines"), legend.position = "bottom", legend.title = element_text(size = 14, vjust = 0.8), legend.text = element_text(size = 12), legend.key.width = unit(4, "cm"), legend.key.height = unit(0.5, "cm"), plot.background = element_rect(fill = "white", color = NA), panel.background = element_rect(fill = "white", color = NA), plot.margin = margin(t = 20, r = 20, b = 20, l = 20) ) # 【关键步骤】在 R 环境中先行展示空间分面地图 cat("\n🗺️.r内展示空间有限，图片混乱是正常的，建议在保存目录中查看.\n") print(p_map) # --------------------------------------------------------------------- # 5. 图表高清导出保存 # --------------------------------------------------------------------- cat("\n--- 正在导出高分辨率图片至工作目录 ---\n") # 导出散点图 ggsave("Scatter_Plot_HighRes.png", plot = p_scatter, width = 8, height = 6, dpi = 600) # 导出时空分面地图 ggsave("Spatio_Temporal_Map_5Years.png", plot = p_map, width = 16, height = 18, dpi = 600, limitsize = FALSE) cat("\n✅ 所有处理完毕！图片已保存到你的当前工作目录 (Working Directory) 下。\n") 

5. 实证案例II: 二分类目标变量

5.1 示例数据

本案例以中国省级层面的新冠肺炎疫情发生情况（二分类结局变量）时空面板数据为例进行演示，使用的数据整合自国家卫生健康委员会、多源遥感卫星监测以及百度人口迁徙大数据平台，收集了2020年全年中国大陆31个省份（自治区、直辖市）地区的人口流动、社会经济与自然环境等因素的真实数据变量，具体数据介绍如表5.1所示。

表5.1 案例Ⅱ示例数据

5.2 数据输入与验证

本教程以案例Ⅱ示例数据为例进行数据导入，用户可以根据自己的数据进行相应操作。

步骤说明：

数据转换：本案例原始数据为含多个时间信息的空间截面数据，故需先进行数据转换。而这里的“截面数据”和“时空面板数据”本质上包含相同的信息，区别主要在于数据的组织方式不同，如下图所示。截面数据格式是将每个地区作为一行，将不同时间点的变量值展开为不同列。例如，地区A一行中同时包含变量1在T1、T2的取值，也包含变量2、变量3、变量4在T1、T2的取值。这种格式更像是“宽表”，适合用户从原始统计表或多期截面数据中整理数据。时空面板数据格式则是将“时间ID”和“地区ID”共同作为行标识，每一行表示某一地区在某一时间点的观测值。例如，T1-地区A是一行，T2-地区A又是另一行，变量1至变量4分别作为列。这种格式更像是“长表”，也是软件进行BSTVC时空建模时使用的标准输入格式。需要注意的是，地区A、B、C、D的顺序原则上应与地图文件中的地区顺序保持一致；如果顺序不一致也没关系，本软件提供数据检查和自动转换模块，可帮助用户检查并调整数据顺序。

【备注】在BSTVC支持的时空面板数据中，变量Y可以包含缺失值，变量X也可容忍少量缺失值。在用户导入自己的数据前，需注意：①缺失值请使用“NA”表示，便于识别；②检查表格数据与地图数据中的空间唯一标识码的字段类型，确保两个表中的该字段的名称和类型一致。

5.3 数据检查

5.3.1 空间单元顺序检查

建模表格数据中的每个时间节点下空间单元的顺序需要与建模地图数据中的shapefile文件中空间单元的顺序完全一致，不一致会造成最终结果不正确，但模型运行并不会报错的情况。BSTVC时空可解释平台提供一键式自动比对数据表与地图中的空间单元唯一值字段的功能，完成空间序列的自动化重排与对齐。

步骤说明：

5.3.2 字段类型检查

在建模前，用户需要特别注意字段类型问题。具体原因、步骤及说明详见4.3.2小节。

5.4 BSTVC建模

构建BSTVC时空模型所需的参数主要有9个，具体参数描述详见4.4小节BSTVC时空建模。

操作步骤如下：

5.5 模型结果

BSTVC函数输出的结果共包括6个部分，具体输出部分描述同4.5小节。

5.6 BSTVC模型输出及其可视化

在本文档中，我们将详细提供一系列的绘制模型输出结果的步骤和代码示例，以便您能够轻松地在R环境中创建各种图表。我们的目标是让您不仅能够理解模型的输出，还能够通过图表的形式，将这些复杂的数据以更直观、更吸引人的方式呈现出来。

本文档在此部分展示的所有图表和可视化结果均基于示例数据构建的模型的输出。以下绘图代码及图件供参考。

5.6.1 模型的贝叶斯评价

模型的贝叶斯评价包括模型名称（包含建模起始时间）、DIC、LS、WAIC、pd2等指标。其中，DIC、WAIC表示模型拟合度，eff（DIC）和pd2（WAIC）表示模型复杂度，LS表示预测精度。具体而言，DIC和WAIC值越小表示模型拟合度越好，eff和pd2值越小表示模型复杂度越低，LS值越接近0表示预测精度越高。

 # 2 模型的贝叶斯评价 --------------------------------------------------------------- # 1. 安装并加载必要的R包（如果还没有安装，请先取消下一行的注释进行安装） # install.packages("readxl") library(readxl) # 2. 设置你的文件路径 file_path <- "/Users/xxx/Binary/BSTVC_all_result_tables.xlsx" # 3. 读取 Excel 文件中的 Sheet1 (即 model.evaluation) # 假设数据在第一个Sheet中，也可以用 sheet = "Sheet1" 指定 model.eval <- read_excel(file_path, sheet = 1) # 4. 指标水平展示（原始格式） cat("\n## 指标水平展示\n") print(model.eval) # 5. 指标垂直展示 cat("\n## 指标垂直展示）\n") # 因首列含中文故使用t()转置，此时所有数据会变成字符矩阵 model.eval_t <- t(model.eval) print(model.eval_t) 

• 结果如下：

5.6.2 时间回归系数（时间非平稳）

BSTVC模型不仅可以计算时间回归系数来刻画变量关系的时间异质性，还可以估算出每个回归系数值的贝叶斯宽（95%）、窄（50%）可信区间，用于直接评价结果的不确定性。

 # 时间非平稳----时间回归系数的时间序列展示形式 rm(list=ls())
gc() library(ggplot2) library(dplyr) library(readxl) library(showtext) showtext_auto() # --- 1. 读取数据--- excel_file <- "/Users/xxx/BSTVC-covid/Binary/BSTVC_all_result_tables.xlsx" Time.Coef_raw <- as.data.frame(read_excel(excel_file, sheet = "time.coefficients")) # --- 2. 数据清洗与格式转换 --- Time.Coef <- Time.Coef_raw %>% mutate(Coefficients = round(mean, 2), Month = time_index) names(Time.Coef)[names(Time.Coef) == "0.025quant"] <- "q025" names(Time.Coef)[names(Time.Coef) == "0.975quant"] <- "q975" names(Time.Coef)[names(Time.Coef) == "0.25quant"] <- "q25" names(Time.Coef)[names(Time.Coef) == "0.75quant"] <- "q75" Time.Coef$q025 <- as.numeric(Time.Coef$q025) Time.Coef$q975 <- as.numeric(Time.Coef$q975) Time.Coef$q25 <- as.numeric(Time.Coef$q25) Time.Coef$q75 <- as.numeric(Time.Coef$q75) Time.Coef$explain_variable <- gsub("time\\.", "", Time.Coef$explain_variable) # --- 3. 绘制时间系数图 --- p_time <- ggplot(data = Time.Coef, aes(x = Month, y = Coefficients, fill = explain_variable)) + geom_ribbon(aes(ymin = q025, ymax = q975, fill = explain_variable), alpha = 0.2) + geom_ribbon(aes(ymin = q25, ymax = q75, fill = explain_variable), alpha = 0.5) + geom_line(aes(x = Month, y = Coefficients, colour = explain_variable), linewidth = 1) + facet_wrap(~explain_variable, scales = "free_y", ncol = 3) + theme_bw() + scale_x_continuous(breaks = 1:12) + labs(title = "各因素影响的时间演变趋势 (2020年月度)", x = "月份", y = "时间系数 (Log-Odds)") + theme( legend.position = "none", # 图表四周的字体配置 plot.title = element_text(hjust = 0.5, face = "bold", size = 26, margin = ggplot2::margin(b = 15)), # 1. 分面标题（变量名） strip.background = element_rect(fill = "grey90", color = "black", linewidth = 0.8), strip.text = element_text(size = 24, face = "bold", color = "black"), # 2. 坐标轴刻度数字 axis.text = element_text(size = 24, color = "black"), # 3. 坐标轴标题 axis.title = element_text(size = 24, face = "bold", margin = ggplot2::margin(t = 10, r = 10)), panel.border = element_rect(color = "black", fill = NA, linewidth = 1), panel.grid = element_blank(), text = element_text(family = "wqy-microhei") ) # ========================================================= # --- 4. 导出 --- # ========================================================= ggsave(filename = "~/Desktop/Time_Coefficients_6Vars_.png", plot = p_time, width = 14, height = 8, dpi = 300) 

5.6.3 空间回归系数（空间非平稳）

与时间回归系数类似，BSTVC模型不仅可以计算空间回归系数表示空间异质的变量关系，还可以 估算其不确定性，即贝叶斯宽（95%）、窄（50%）可信区间。

提取出模型输出中的空间回归系数数据框后，建议在ArcGIS或ArcGIS Pro等专业地图制图软件中进行绘制。此处基于R语言，简单绘制出基础的空间回归系数 地图展示结果，仅供样本制图参考。

 # 4 空间非平稳 ----------------------------------------------------------------- # 如果未安装 patchwork，请先取消下方注释进行安装 # install.packages("patchwork") library(patchwork) # ========================================================= # 提取所有需要画图的变量名 variables_to_plot <- unique(Space.coef.panel$variable) # 使用 lapply 循环为每个变量单独画图 plot_list <- lapply(variables_to_plot, function(var_name) { # 提取当前变量的数据 current_data <- Space.coef.panel %>% filter(variable == var_name) # 单独绘图 p <- ggplot(current_data) + geom_sf(aes(fill = SCs), color = "gray60", linewidth = 0.2) + scale_fill_gradient2( low = "#4575b4", mid = "white", high = "#d73027", midpoint = 0, name = "SCs" ) + theme_bw() + # 将变量名作为每个小图的标题 labs(title = var_name) + # 统一美化设置 theme( legend.position = "bottom", legend.key.width = unit(1.2, "cm"), legend.title = element_text(face = "bold", size = 32), legend.text = element_text(size = 30), plot.title = element_text(hjust = 0.5, face = "bold", size = 38, margin = ggplot2::margin(b = 10)), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), panel.border = element_rect(color = "black", fill = NA, linewidth = 1), text = element_text(family = "wqy-microhei") ) return(p)
}) # 将 6 张图按 2行3列 拼接，并添加总标题 p_final <- wrap_plots(plot_list, ncol = 3) + plot_annotation( title = "各影响因子空间效应的非平稳性分布", theme = theme( plot.title = element_text(hjust = 0.5, face = "bold", size = 34, family = "wqy-microhei", margin = ggplot2::margin(b = 20)) ) ) # 导出图片 ggsave( filename = "~/Desktop/Space_Coefficients_Map_Patchwork.png", plot = p_final, width = 16, # 稍微加宽一点画板以容纳 6 个独立图例 height = 10, dpi = 300 ) 

5.6.4 时空贡献度评价

时空方差分割指标（Spatiotemporal Variance Partitioning Index，STVPI）旨在量化并比较不同时空异质影响因素的可解释百分比，通过计算时空贡献度来明确关键驱动因素。基于BSTVC建模结果，STVPI进一步将每个解释变量的总方差分解为时间和空间两个独立的组分，即时间非平稳随机效应和空间非平稳随机效应，进而揭示数据的时空变异来源。

此工具的优点有：①不仅能识别解释因子的时空总体贡献，还能分别识别时间和空间维度的贡献；②相比只能识别绝对贡献排序的主流方法（如随机森林、SHAP等），STVPI识别的是相对贡献 （可解释百分比），为地理时空归为提供重要证据基础；③直接评价其不确定性（贝叶斯宽窄可信区间）

 # 5 STVPI时空方差分割指标 ------------------------------------------------------------- rm(list=ls())
gc() library(ggplot2) library(dplyr) library(readxl) library(showtext) showtext_auto() # 1. 设置结果表格路径 my_path <- "/Users/xxx/BSTVC_all_result_tables(1).xlsx" STVPI.data <- as.data.frame(read_excel(my_path, sheet = "STVPI")) STVPI.data[,2:6] <- lapply(STVPI.data[,2:6], as.numeric) colnames(STVPI.data) <- c("id", "STVPI_2.5", "STVPI_25", "STVPI_mean", "STVPI_75", "STVPI_97.5", "effect_name") STVPI.data$effect_name <- trimws(STVPI.data$effect_name) # 提取变量名 time_raw <- unique(STVPI.data$effect_name[grepl("^time\\.", STVPI.data$effect_name) & STVPI.data$effect_name != "time.all"]) my_factors <- unique(gsub("^time\\.", "", time_raw)) my_factors <- my_factors[!is.na(my_factors) & my_factors != ""] # 划分四个组别 st_all_vars <- c("space.all", "time.all") space_vars <- paste0("space.", my_factors) time_vars <- paste0("time.", my_factors) st_vars <- paste0("space.time.", my_factors) # 2. 过滤并构建 4 个学术组别的数据 plot.data <- STVPI.data %>% filter(effect_name %in% c(st_all_vars, space_vars, time_vars, st_vars)) %>% distinct(effect_name, .keep_all = TRUE) %>% mutate(Group = case_when( effect_name %in% st_all_vars ~ "1. 时空总贡献度 (Space & Time Totals)", effect_name %in% space_vars ~ "2. 空间贡献 (Space)", effect_name %in% time_vars ~ "3. 时间贡献 (Time)", effect_name %in% st_vars ~ "4. 时空贡献 (Space-Time)" )) # 控制 Y 轴因子的上下排列顺序 plot.data$effect_name <- factor(plot.data$effect_name, levels = rev(c(st_all_vars, space_vars, time_vars, st_vars))) # 3. 绘制STVPI图 p_combined <- ggplot(plot.data, aes(x = STVPI_mean, y = effect_name)) + # 95% 置信区间 geom_errorbar(aes(xmin = STVPI_2.5, xmax = STVPI_97.5, colour = Group), width = 0, linewidth = 1.2, alpha = 0.5) + # 50% 置信区间 geom_errorbar(aes(xmin = STVPI_25, xmax = STVPI_75, colour = Group), width = 0, linewidth = 2.5, alpha = 0.8) + # 样式设置 geom_point(aes(colour = Group), shape = 124, size = 21, stroke = 0.9) + geom_text(aes(label = sprintf("%.2f%%", STVPI_mean * 100)), vjust = -1.6, hjust = 0.5, size = 10.6, fontface = "bold", family = "wqy-microhei") + # 统一布局 facet_wrap(~Group, scales = "free_y", ncol = 2) + theme_bw() + scale_x_continuous(labels = scales::percent_format()) + labs(title = "模型时空方差贡献度(STVPI)", x = "STVPI 贡献率百分比", y = "") + theme( legend.position = "none", # 顶部小标题 strip.background = element_rect(fill = "grey90", color = "black", linewidth = 1), strip.text = element_text(family = "wqy-microhei", size = 41, face = "bold", color = "black"), # 设置字号样式 plot.title = element_text(hjust = 0.5, face = "bold", size = 41, family = "wqy-microhei", margin = ggplot2::margin(b = 20)), axis.text.x = element_text(size = 41, color = "black"), axis.text.y = element_text(size = 41, color = "black"), axis.title.x = element_text(size = 41, face = "bold", margin = ggplot2::margin(t = 15)), panel.border = element_rect(color = "black", fill = NA, linewidth = 1), panel.grid.major = element_line(color = "grey95"), panel.grid.minor = element_blank(), # 增加画板四周的内边距 plot.margin = ggplot2::margin(t = 20, r = 20, b = 20, l = 20), text = element_text(family = "wqy-microhei") ) # 4. 导出图片 ggsave(filename = "~/Desktop/STVPI_Combined_4Facets_LongLine.png", plot = p_combined, width = 16, height = 13, dpi = 300) 

5.6.5 目标变量的时空预测

对于数据中目标变量Y的所有缺失值或非缺失值，局部预测结果中都会输出预测值，以及预测值的贝叶斯宽窄可信区间，用于直接评价不确定性。

其中，对于变量Y所有类型（连续型、计数型与二值型）的局部预测结果，y字段代表目标变量Y的原始值，fill_y字段代表对目标变量Y中缺失值进行预测填补后的完整字段，predict_y代表为所有单元的模型预测值。

1）观察读入数据结构

2）时空预测可视化

本案例选择绘制ROC曲线来动态展示灵敏度与误报率的权衡，并提供全局量化指标（AUC值），AUC值越大，说明预测效果越好。

 ## 绘制 Logistic-BSTVC: ROC 曲线 --- rm(list=ls())
gc() library(ggplot2) library(dplyr) library(readxl) library(showtext) library(pROC) showtext_auto() # 1. 文件路径 my_path <- "/Users/xx/BSTVC_all_result_tables(1).xlsx" predict_data <- as.data.frame(read_excel(my_path, sheet = "local.prediction")) # 提取没有缺失值的真实验证集部分 predict_clean <- predict_data %>% filter(!is.na(y)) # 2. 生成 ROC 对象并获取 AUC 值 roc_obj <- roc(predict_clean$y, predict_clean$predict_y, quiet = TRUE) auc_value <- as.numeric(auc(roc_obj)) cat("\n==========================================\n") cat(sprintf("模型预测精度 (AUC值): %.4f\n", auc_value)) cat("评价标准: 0.7~0.8(良好), 0.8~0.9(优秀), >0.9(极佳)\n") cat("==========================================\n") # 3. 提取用于 ggplot2 绘图的坐标点 roc_df <- data.frame( FPR = 1 - roc_obj$specificities, # 假阳性率 (1 - 特异度) TPR = roc_obj$sensitivities # 真阳性率 (灵敏度) ) %>% arrange(FPR, TPR) # 确保 (0,0) 始终是连线的第一个点 # 4. 绘制ROC 图 p_roc <- ggplot(roc_df, aes(x = FPR, y = TPR)) + # 绘制对角线和ROC曲线 geom_abline(slope = 1, intercept = 0, color = "grey50", linewidth = 1.5, linetype = "dashed") + geom_path(color = "#b2182b", linewidth = 1) + # 标注 AUC 值 annotate("text", x = 0.70, y = 0.25, label = sprintf("AUC = %.3f", auc_value), size = 24, fontface = "bold", color = "black", family = "wqy-microhei") + theme_bw() + # 为四周留出 2% 的物理边距 scale_x_continuous(limits = c(0, 1), expand = expansion(add = 0.02)) + scale_y_continuous(limits = c(0, 1), expand = expansion(add = 0.02)) + # 设置X 和 Y 轴比例 coord_fixed() + labs( title = "Logistic-BSTVC 模型预测精度评估 (ROC 曲线)", x = "假阳性率 (1 - 特异度, FPR)", y = "真阳性率 (灵敏度, TPR)" ) + theme( plot.title = element_text(hjust = 0.5, face = "bold", size = 46, margin = ggplot2::margin(b = 25)), axis.text = element_text(size = 42, color = "black"), axis.title.x = element_text(size = 44, face = "bold", margin = ggplot2::margin(t = 20)), axis.title.y = element_text(size = 44, face = "bold", margin = ggplot2::margin(r = 20)), # 外部边框 panel.border = element_rect(color = "black", fill = NA, linewidth = 2), panel.grid.major = element_line(color = "grey90"), panel.grid.minor = element_blank(), plot.margin = ggplot2::margin(t = 20, r = 20, b = 20, l = 20), text = element_text(family = "wqy-microhei") ) # 5. 导出图片 ggsave(filename = "~/Desktop/Model_Evaluation_ROC_NoShade_Final.png", plot = p_roc, width = 10, height = 10, dpi = 300) 

3）预测后的Y的时空分面可视化地图

 ## 预测后的 Y 的时空分面可视化地图 )# --- 清理环境与加载包 ---
rm(list=ls())
gc() library(sf) library(ggplot2) library(dplyr) library(readxl) library(showtext) library(ggthemes) showtext_auto() # 1. 设置文件路径 excel_file <- "/Users/xxx/BSTVC_all_result_tables(1).xlsx" shp_file <- "/Users/xxx/中国省级地图_审图号GS(2024)0650号.shp" # 2. 读取地图数据与预测结果 china_map <- st_read(shp_file) predict_data <- as.data.frame(read_excel(excel_file, sheet = "local.prediction")) # 3. 数据合并 predict_map <- china_map %>% mutate(gb = as.character(gb)) %>% left_join( predict_data %>% mutate(gb = as.character(gb)), by = "gb" ) %>% mutate(Month_CN = factor(Month, levels = 1:12, labels = c("一月", "二月", "三月", "四月", "五月", "六月", "七月", "八月", "九月", "十月", "十一月", "十二月"))) # 4. 绘制时空分面可视化地图 p7 <- ggplot(predict_map) + geom_sf(aes(fill = predict_y), color = "white", linewidth = 0.25) + facet_wrap(~ Month_CN, ncol = 4) + scale_fill_gradient( low = "#e0f3db", high = "#e34a33", na.value = "grey90", name = "预测概率" ) + labs( title = "各省份疫情发生概率的时空演变预测 (2020年)", x = "", y = "" ) + theme_few() + theme( text = element_text(family = "wqy-microhei"), plot.title = element_text(hjust = 0.5, face = "bold", size = 110, margin = ggplot2::margin(b = 50)), axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank(), panel.grid.major = element_blank(), panel.border = element_rect(color = "grey35", fill = NA, linewidth = 0.6), strip.background = element_rect(fill = "grey90", color = "grey35", linewidth = 0.6), strip.text = element_text(size = 70, face = "bold"), legend.title = element_text(size = 70, face = "bold", margin = ggplot2::margin(b = 20)), legend.text = element_text(size = 60), legend.position = "right", legend.key.height = unit(5, "cm"), legend.key.width = unit(2, "cm") ) # 5. 导出高分辨率图片 ggsave( filename = "~/Desktop/p7_predict_time_map_CN_NoLatLon_5xFont.png", plot = p7, width = 40, height = 25, dpi = 300, limitsize = FALSE ) 

6. 实证案例III: 计数型目标变量

6.1 示例数据

本案例以美国洲级层面的急性乙肝监测病例（泊松型结局变量）时空面板数据为例进行演示，使用的数据美国CDC、美国社区调查（ACS）及美国国家环境信息中心（NOAA），收集了2012到2023年间美国阿巴拉契亚地区5个州（肯塔基州KY，北卡罗来纳州NC，田纳西州TN，弗吉尼亚州VA，西弗吉尼亚州WV）地区的社会经济与自然环境等因素的真实数据变量，具体数据介绍如表6.1所示，本示例数据无缺失。

表6.1 案例Ⅲ示例数据

6.2 数据输入与验证

本教程以案例Ⅲ示例数据为例进行数据导入，用户可以根据自己的数据进行相应操作。

步骤说明：

数据转换：本案例原始数据为含多个时间信息的空间截面数据，故需先进行数据转换。关于截面数据与时空面板数据的示例、示意图与解释说明，可见5.2小节。在左侧导航栏选择数据转换模块，点击读取转换数据，导入原始数据

【备注】在BSTVC支持的时空面板数据中，变量Y可以包含缺失值，变量X也可容忍少量缺失值。在用户导入自己的数据前，需注意：①缺失值请使用“NA”表示，便于识别；②检查表格数据与地图数据中的空间唯一标识码的字段类型，确保两个表中的该字段的名称和类型一致。

6.3 数据检查

6.3.1 空间单元顺序检查

建模表格数据中的每个时间截面下空间单元的顺序需要与建模地图数据中的shapefile文件中空间单元的顺序完全一致，不一致会造成最终结果不正确，但模型运行并不会报错的情况。例如，地图us_a5中代表国家与地区单元的GEOID字段的顺序为21,37,47则BSTVC_panel_converted.xlsx数据中的GEOID字段也必须为此排序。

BSTVC时空可解释平台提供一键式自动比对数据表与地图中的空间单元唯一值字段的功能，完成空间序列的自动化重排与对齐

步骤说明：

6.3.2 字段类型检查

在建模前，用户需要特别注意字段类型问题。具体原因、步骤及说明详见4.3.2小节。

6.4 BSTVC建模

构建BSTVC时空模型所需的参数主要有9个，具体参数描述详见4.4小节BSTVC时空建模。

6.5 模型结果

BSTVC函数输出的结果共包括6个部分，具体输出部分描述同4.5小节。

6.6 BSTVC模型输出及其可视化

在本文档中，我们将详细提供一系列的绘制模型输出结果的步骤和代码示例，以便您能够轻松地在R环境中创建各种图表。我们的目标是让您不仅能够理解模型的输出，还能够通过图表的形式，将这些复杂的数据以更直观、更吸引人的方式呈现出来。

本文档在此部分展示的所有图表和可视化结果均基于示例数据构建的模型的输出。以下绘图代码及图件供参考

6.6.1 基础R包配置与数据读入

 #BSTVC 函数的输出及其可视化 # 0 R包配置 ---------------------------------------------------------- library(readxl) # 读取 Excel 文件，如 .xls 和 .xlsx library(sf) # 处理空间矢量数据，如 shp 文件、坐标系和空间分析 library(stringr) # 字符串处理，如提取、替换、匹配文本 library(ggplot2) # 绘图核心包，用于构建各类统计图形 library(dplyr) # 数据整理，如筛选、分组、汇总和变量变换 library(tidyr) # 数据整形，如宽表转长表、长表转宽表 library(ggthemes) # 提供额外的 ggplot2 主题，如 theme_few() library(forestplot) # 绘制森林图，常用于展示回归系数和置信区间 library(ggbeeswarm) # 绘制蜂群图，展示分组散点且减少点重叠 library(scales) # 坐标轴、图例和数值格式化，如百分比、颜色渐变 # 1 数据读入 ----------------------------------------------------------- setwd("C:/Users/XXXX ") #实际存储工作目录 # 读取.xlsx文件中的各类sheet表结果，并存为对应数据框 model.eval <- read_excel("hbv_all_result.xlsx", sheet = "model.evaluation") predict <- read_excel("hbv_all_result.xlsx", sheet = "local.prediction") time.coefficients <- read_excel("hbv_all_result.xlsx", sheet = "time.coefficients") space.coefficients<- read_excel("hbv_all_result.xlsx", sheet = "space.coefficients") stvpi<- read_excel("hbv_all_result.xlsx", sheet = "STVPI") #读取地图数据 us_a5<- st_read("C:/Users/XXX/us_a5") 

6.6.2 模型的贝叶斯评价

模型的贝叶斯评价包括模型名称（包含建模起始时间）、DIC、LS、WAIC、pd2等指标。其中，DIC、WAIC表示模型拟合度，eff（DIC）和pd2（WAIC）表示模型复杂度，LS表示预测精度。具体而言，DIC和WAIC值越小表示模型拟合度越好，eff和pd2值越小表示模型复杂度越低，LS值越接近0表示预测精度越高。

 # 2 模型的贝叶斯评价 --------------------------------------------------------------- ##指标水平展示 print(model.eval) ##指标垂直展示，因首列含中文故需使用t()转置 model.eval_t <- t(model.eval) print(model.eval_t) 

• 结果如下：

6.6.3 时间回归系数（时间非平稳）

BSTVC模型不仅可以计算时间回归系数来刻画变量关系的时间异质性，还可以估算出每个回归系数值的贝叶斯宽（95%）、窄（50%）可信区间，用于直接评价结果的不确定性。

1）初始数据处理

 # 3 时间非平稳 ----------------------------------------------------------------- # 提取模型输出中的时间回归系数，添加系数字段与可信区间字段 Time.Coef <- time.coefficients %>% mutate( Coefficients = round(mean, 2), `95%CI` = paste0("(", round(`0.025quant`, 2), " - ", round(`0.975quant`, 2), ")") ) # 将时间字段转换为年份（可做可不做） Time.Coef <- Time.Coef %>% mutate(Year = case_when( time_index == 1 ~ 2012, time_index == 2 ~ 2013, time_index == 3 ~ 2014, time_index == 4 ~ 2015, time_index == 5~ 2016, time_index == 6~ 2017, time_index == 7~ 2018, time_index == 8~ 2019, time_index == 9~ 2020, time_index == 10~ 2021, time_index == 11~ 2022, time_index == 12~ 2023, )) 

2）第一种制图-时间系数不确定性带图

 ## 时间回归系数的时间序列展示形式（第一种制图-时间系数不确定性带图） p3_1 <- ggplot(data=Time.Coef,aes(x=Year,y= Coefficients ,fill=explain_variable))+ geom_ribbon(aes(ymin=`0.025quant`,ymax=`0.975quant`,fill=explain_variable),alpha=0.1)+ geom_ribbon(aes(ymin=`0.25quant`,ymax=`0.75quant`,fill=explain_variable),alpha=0.5)+ geom_line(aes(x=Year,y= Coefficients,colour=explain_variable),alpha=0.8)+ facet_wrap(~explain_variable,scale = "free_y",axes = "all_x")+ theme_few()+ labs( title = NULL, x = "年份", y = "时间系数 (TCs)", fill = "解释因子", colour = "解释因子" ) p3_1 #保存图片结果 ggsave( "time_coefficients_1.png", p3_1, width = 12, height = 8, dpi = 300 ) 

3）第二种制图-多变量折线图

为避免对时间变化形式作出线性假设，本例采用LOESS局部加权回归对年度时间系数进行平滑展示，并保留原始估计点以反映真实年度波动。

 ## 时间回归系数（第二类制图：原始点 + LOESS 平滑曲线 + 95%置信区间） p3_2 <- ggplot( Time.Coef, aes(x = Year, y = Coefficients, colour = explain_variable, fill = explain_variable) ) + geom_point( shape = 21, size = 2.2, alpha = 0.75, color = "white", stroke = 0.4 ) + geom_smooth( se = TRUE, method = "loess", linewidth = 1.2, span = 0.8, alpha = 0.18 ) + scale_x_continuous( breaks = sort(unique(Time.Coef$Year)) ) + theme_few() + labs( title = NULL, x = "年份", y = "时间系数 (TCs)", fill = "解释因子", colour = "解释因子" ) + theme( legend.title = element_text(face = "bold"), axis.title = element_text(face = "bold") ) p3_2 #保存图片结果 ggsave( "time_coefficients_2.png", p3_2, width = 12, height = 8, dpi = 300 ) 

6.6.4 空间回归系数（空间非平稳）

与时间回归系数类似，BSTVC模型不仅可以计算空间回归系数表示空间异质的变量关系，还可以 估算其不确定性，即贝叶斯宽（95%）、窄（50%）可信区间。

提取出模型输出中的空间回归系数数据框后，建议在ArcGIS或ArcGIS Pro等专业地图制图软件中进行绘制。此处基于R语言，简单绘制出基础的空间回归系数地图展示结果，仅供样本制图参考。

 # 4 空间非平稳 ----------------------------------------------------------------- # 统一 GEOID 格式 us_a5 <- us_a5 %>% mutate(GEOID = str_pad(as.character(GEOID), width = 2, pad = "0")) space.coefficients <- space.coefficients %>% mutate(GEOID = str_pad(as.character(GEOID), width = 2, pad = "0")) # 合并地图和空间系数 Space.Coef.merge <- left_join(us_a5, space.coefficients, by = "GEOID") # 需要绘图的空间系数字段,务必确保一一对应! coef_vars <- paste0( c( "PD", "POV", "UNINS", "PRECIP"
),
"_mean" ) # 转成长表 Space.coef.panel <- Space.Coef.merge %>% pivot_longer( cols = all_of(coef_vars), names_to = "variable", values_to = "SCs" ) %>% mutate( variable = gsub("_mean$", "", variable), SCs = as.numeric(SCs) ) # 检查 variable 是否存在 names(Space.coef.panel) # 绘图 p4 <- ggplot(Space.coef.panel) + geom_sf(aes(fill = SCs)) + facet_wrap(~ variable, shrink = FALSE, drop = FALSE) + scale_fill_gradient( low = "#d9f1e6", high = "#4e5180", name = "空间系数 (SCs)" ) + theme_few() p4 #保存图片结果 ggsave( "space_coefficients.png", p4, width = 12, height = 8, dpi = 300 ) 

6.6.5 时空贡献度评价

时空方差分割指标（Spatiotemporal Variance Partitioning Index，STVPI）旨在量化并比较不同时空异质影响因素的可解释百分比，通过计算时空贡献度来明确关键驱动因素。基于BSTVC建模结果，STVPI进一步将每个解释变量的总方差分解为时间和空间两个独立的组分，即时间非平稳随机效应和空间非平稳随机效应，进而揭示数据的时空变异来源。

此工具的优点有：①不仅能识别解释因子的时空总体贡献，还能分别识别时间和空间维度的贡献；②相比只能识别绝对贡献排序的主流方法（如随机森林、SHAP等），STVPI识别的是相对贡献 （可解释百分比），为地理时空归为提供重要证据基础；③直接评价其不确定性（贝叶斯宽窄可信区间）

 # 5 STVPI时空方差分割指标 ------------------------------------------------------------- random.effects <- c( "space.all", "time.all", "space.PD", "space.POV", "space.PRECIP", "space.UNINS", "time.PD", "time.POV", "time.PRECIP", "time.UNINS", "space.time.PD", "space.time.POV", "space.time.PRECIP", "space.time.UNINS" ) plot.data <- as.data.frame( stvpi[stvpi$`random effects` %in% random.effects, ] ) %>% mutate( Group = case_when( `random effects` %in% c("space.all", "time.all") ~ "Space & Time", grepl("^space\\.time\\.", `random effects`) ~ "space - time", grepl("^space\\.", `random effects`) ~ "space", grepl("^time\\.", `random effects`) ~ "time", TRUE ~ NA_character_
) ) %>% filter(!is.na(Group)) %>% mutate( Group = factor( Group, levels = c("space", "time", "space - time", "Space & Time") ), label = scales::percent(STVPI_mean, accuracy = 0.01) ) plot.data <- plot.data %>% group_by(Group) %>% mutate( y_order = factor( `random effects`, levels = rev(`random effects`) ) ) %>% ungroup() p5 <- ggplot(plot.data, aes(x = STVPI_mean, y = y_order)) + geom_errorbarh( aes( xmin = `STVPI_2.5%`, xmax = `STVPI_97.5%`, colour = `random effects` ), height = 0, linewidth = 1.8, alpha = 0.3 ) + geom_errorbarh( aes( xmin = `STVPI_25%`, xmax = `STVPI_75%`, colour = `random effects` ), height = 0, linewidth = 1.8, alpha = 0.6 ) + geom_point( aes(colour = `random effects`), shape = 21, fill = "white", size = 1.8, alpha = 1 ) + geom_text( aes(label = label), hjust = -0.15, size = 3, show.legend = FALSE ) + facet_wrap( ~ Group, scales = "free", ncol = 3 ) + scale_x_continuous( labels = scales::percent_format(accuracy = 1), expand = expansion(mult = c(0.05, 0.25)) ) + labs( title = "", x = "STVPI", y = "" ) + theme_few() + theme( legend.position = "none", strip.text = element_text(size = 12), axis.text.y = element_text(size = 10) ) p5 #保存图片结果 ggsave( "STVPI.png", p5, width = 12, height = 8, dpi = 300 ) 

6.6.6 目标变量的时空预测

对于数据中目标变量Y的所有缺失值或非缺失值，局部预测结果中都会输出预测值，以及预测值的贝叶斯宽窄可信区间，用于直接评价不确定性。

其中，对于变量Y所有类型（连续型、计数型与二值型）的局部预测结果，y字段代表目标变量Y的原始值，fill_y字段代表对目标变量Y中缺失值进行预测填补后的完整字段，predict_y代表为所有单元的模型预测值。

1）观察读入数据结构

2）评估指标计算与散点图可视化

对于计数型目标变量Y，若模型预测结果已回到原始计数尺度，则可基于观测值与预测均值绘制散点图，并计算MAE、Mean Poisson Deviance等预测精度指标。

①MAE（平均绝对误差）：表示模型平均预测错了多少个单位

②Mean Poisson Deviance（平均偏差）：模型预测的计数分布，与真实计数之间平均偏离多少

③Deviance R²（基于泊松偏差的伪R²）：越接近1，说明模型比“只用平均值预测”的空模型越好。解释为相比空模型，当前模型解释或降低了约82%的泊松偏差。

注：尽管本案例没有缺失数据，此处只是做示例处理若有缺失数据的情况

 predict_clean <- predict %>% filter(!is.na(y), !is.na(predict_y)) %>% mutate( y = as.numeric(y), predict_y = pmax(as.numeric(predict_y), 1e-8) ) ## 泊松型预测评价指标 poisson_deviance <- function(y, mu) { 2 * sum(ifelse(y == 0, 0, y * log(y / mu)) - (y - mu)) } MAE <- mean(abs(predict_clean$y - predict_clean$predict_y)) model_deviance <- poisson_deviance( y = predict_clean$y, mu = predict_clean$predict_y ) Mean_Poisson_Deviance <- model_deviance / nrow(predict_clean) null_mu <- mean(predict_clean$y) null_deviance <- poisson_deviance( y = predict_clean$y, mu = rep(null_mu, nrow(predict_clean)) ) Deviance_R2 <- 1 - model_deviance / null_deviance label_text <- paste0( "MAE = ", round(MAE, 4), "\nMean Poisson dev. = ", round(Mean_Poisson_Deviance, 4) ) ## 打印结果 print(paste("MAE:", round(MAE, 4))) print(paste("Mean Poisson Deviance:", round(Mean_Poisson_Deviance, 4))) print(paste("Deviance R²:", round(Deviance_R2, 4))) 

 ## 绘制预测散点图 p6 <- ggplot(predict_clean, aes(x = y, y = predict_y)) + geom_abline( slope = 1, intercept = 0, color = "red", linewidth = 1, linetype = "dashed" ) + geom_point(aes(fill = y), size = 2.5, shape = 22, color = "white") + annotate(
"text", x = Inf, y = -Inf, label = label_text, hjust = 1.05, vjust = -0.35, family = "serif", fontface = "bold", size = 4.5 ) + scale_fill_gradient(low = "#88d8db", high = "#b8aeeb") + labs( x = "Measured Values", y = "Predicted Values" ) + theme_few() + theme( text = element_text(family = "serif", face = "bold"), axis.text = element_text(family = "serif", size = 12, face = "bold"), axis.title = element_text(size = 12), legend.position = "none" ) p6 #保存图片结果 ggsave( "predict_point.png", p6, width = 12, height = 8, dpi = 300 ) 

3）预测后的Y的时空分面可视化地图

 ## 预测后的 Y 的时空分面可视化地图 #观察数据结构 便于绘制 names(predict)
class(predict)
head(predict)
names(us_a5)
class(us_a5)
head(us_a5) predict_map <- us_a5 %>% mutate(GEOID = as.character(GEOID)) %>% left_join( predict %>% mutate(GEOID = as.character(GEOID)), by = "GEOID" ) ## 预测后的 Y 的时空分面可视化地图，本例设置为4列3行 p7 <- ggplot(predict_map) + geom_sf(aes(fill = predict_y), color = "white", linewidth = 0.25) + facet_wrap(~ YEAR, ncol = 4) + scale_fill_gradient( low = "#88d8db", high = "#b8aeeb", na.value = "grey90", name = "预测值" ) + scale_x_continuous( breaks = seq(-90, -75, by = 5), labels = function(x) paste0(abs(x), "°W") ) + scale_y_continuous( breaks = seq(34, 41, by = 2), labels = function(y) paste0(y, "°N") ) + coord_sf( xlim = c(-91, -74), ylim = c(33.5, 41), expand = FALSE ) + labs( title = NULL, x = "经度", y = "纬度" ) + theme_few() + theme( text = element_text(family = "serif", face = "bold"), axis.text = element_text(size = 7, color = "black"), axis.title = element_text(size = 10, face = "bold"), axis.ticks = element_line(color = "grey40", linewidth = 0.25), panel.grid.major = element_line(color = "grey85", linewidth = 0.25), panel.border = element_rect(color = "grey35", fill = NA, linewidth = 0.4), strip.text = element_text(size = 10, face = "bold"), legend.title = element_text(size = 10, face = "bold"), legend.text = element_text(size = 8), legend.position = "right" ) p7 ggsave( filename = "p7_predict_time_map.png", plot = p7, width = 12, height = 8, dpi = 300 )