Rookie宇

V1

2022/07/12阅读:13主题:雁栖湖

R语言描述性时空统计模型

Part1描述性时空统计模型-时空克里格法(Spatio-Temporal Kriging)

本文的主要目标:

  • R语言实现时空克里格法(Spatio-Temporal Kriging)
  • 通用时空克里格法预测结果绘图
  • 通用时空克里格法预测标准差结果绘图

尽管时间是一维向前的,但回顾一段时期内目标观测的变化也是很有价值的。我们可以用空间-时间均值协方差函数来描述空间和时间的相互作用,而不必致力于建立一个动态表达相互作用的机械模型。

在本中我们使用"描述性"方法的时空模型。我们明确区分了数据代表测量所依据的现实世界过程的基本潜在过程( the underlying latent process that represents the real-world process upon which measurements were taken)。 也就是说,我们需要进行有条件地思考。我们考虑了:

  • 条件分布为高斯的数据模型
  • 条件分布为非高斯的数据模型

在这两种情况下,我们都以一个潜在的高斯时空过程为条件。

对于潜在的高斯时空过程,我们考虑用外生协变量(包括空间或时间位置的函数)来指定一阶(平均)结构,用时空协方差函数来指定二阶依赖。我们讨论了与静止性 stationarity、分离性separability和完全对称性full symmetry有关的此类模型的各种假设。

这类表示法非常适合于这样的问题:没有太多的观测数据或想要预测的时间和空间位置。通常我们关注并可以较轻松的知道依赖结构(可以用协方差函数表示),而不太在乎模型是否符合现实。在有大量数据集和/或大量预测位置的情况下,考虑使用基函数展开的二阶结构的随机效应表示,往往在计算上更有效率。基函数结构也使建模者不必开发有效的时空协方差函数,因为我们的条件基函数随机效应会引起有效的边际协方差函数。此外,我们简要地展示了这些使用基函数的时空混合模型与广义加性模型(GAM/GAMMs)的关系,这取决于基函数的选择和估计的方法。

在具有空间或时空随机效应的描述性模型中进行参数推断的一个潜在问题是混杂(confounding)的问题。传统上,这在空间和时空统计中并不是一个大问题,因为重点是预测。但是,随着这些方法在解释固定效应时越来越多地被用来解释依赖性,混杂问题得到了更多的关注。

无论是从边际方差模型的角度还是从条件基函数的角度,表征时空依赖性结构的最具挑战性的方面之一是对现实世界中发生的跨时空的相互作用进行建模的能力。 在这种情况下,基本过程通常最好由空间场(spatial fields)来描述,这些空间场根据管理时空变化的 "规则 "随时间演变。也就是说,它们代表了一个动态系统。

1使用 gstat 包实现时空克里格法(Spatio-Temporal Kriging)

在本实验中,我们将使用带有 gstat包的半变异函数(semivariogram)来完成时空通用克里格法(Spatio-Temporal universal kriging)的过程。我们重点关注 1993 年 7 月 NOAA 数据集 (Tmax) 中的最高温度数据。

加载所需要的包

library("sp")
library("spacetime"
library("ggplot2"
library("dplyr"
library("gstat"
library("RColorBrewer"
library("STRbook"
library("tidyr")

经验半变异函数(empirical semivariogram)

对于 1993 年 7 月最高温度数据集的 S-T 克里格法,我们需要将参数函数拟合到经验半变异函数(empirical semivariogram)中。经验半变异函数如下:

data("STObj3", package = "STRbook")
STObj4 <- STObj3[, "1993-07-01::1993-07-31"]

vv <- variogram(object = z ~ 1 + lat, # fixed effect component
                data = STObj4, # July data
                width = 80# spatial bin (80 km)
                cutoff = 1000# consider pts < 1000 km apart
                tlags = 0.01:6.01# 0 days to 6 days

gstat 提供了许多协方差函数模型 covariannce-function models;详细情况可以通过以下代码获得:

vignette("spatio-temporal-kriging")

使用gstat包里的函数 vgmST 构建时空半变异函数。参数 stModel = "separable" 用于定义可分离模型,而函数 vgm用于构建单独的半变异函数(一个用于空间,一个用于时间)。可以将几个参数传递给 vgm, 我们在下面使用的前四个参数,分别对应于部分测点(partial sill)、模型类型(model type)、范围(range)和块(nnugget)。提供给vgmST的参数sill定义了联合时空基。它们定义中使用的数字是提供给函数fit.StVariogram中用于拟合的优化例程的初始值,该函数将sepVgm拟合到vv中。这些初始值应该是合理的,例如:长度尺度 可以设置为横跨10%的空间/时间域的值,而变差/sill可以设置为与测量的总变差具有相似的数量级。

sepVgm <- vgmST(stModel = "separable",
                space = vgm(10"Exp"400, nugget = 0.1),
                time = vgm(10"Exp"1, nugget = 0.1),
                sill = 20)
sepVgm <- fit.StVariogram(vv, sepVgm)

可分离半变异函数(Separable (in Space and Time) Covariance Functions)

我们拟合的第二个模型具有可分离半变异的协方差函数。对于这个模型,函数 vgmST 将联合半变异函数作为参数,以及 (sill) 和比例因子 (stAni),用 in 。这个参数可以通过考虑数量级来初始设置——如果空间场在数百公里的尺度上演化,而时间演化的尺度在天级,那么初始值 stAni = 100 是合理的.

metricVgm <- vgmST(stModel = "metric",
                   joint = vgm(100"Exp"400, nugget = 0.1),
                   sill = 10,
                   stAni = 100)
metricVgm <- fit.StVariogram(vv, metricVgm)

我们可以通过检查拟合的均方误差来比较两个半变异函数的拟合。这些可以通过直接访问 fit.StVariogram 使用的优化器的最终函数值来找到。

metricMSE <- attr(metricVgm, "optim")$value 
sepMSE <- attr(sepVgm, "optim")$value

此处变量 metricMSE 为 2.1,而 sepMSE 为 1.4,表明在这种情况下可分离半变异函数比经验半变异函数更适合。可以使用标准绘图函数绘制拟合的半变异函数。

plot(vv, list(sepVgm, metricVgm), main = "Semi-variance")
图1.半变异函数
图1.半变异函数

拟合变异函数的等值线图显示在图1的底部面板中。相应的平稳 S-T 协方差函数是从 (4.15) 中获得的。

预测-使用通用S-T克里格法

接下来,基于已拟合的 S-T 协方差模型 S-T covariance models来使用 通用 S-T 克里格法 universal S-T kriging进行预测,并将纬度坐标视为协变量。

我们需要创建一个时空预测网格

  • 对于我们的空间网格,我们考虑了 100°W 和 80°W 之间的 20 个空间位置,以及 32°N 和 46°N 之间的 20 个空间位置。在下面的代码中,在转换为 SpatialPoints 时,我们确保预测网格的坐标参考系统 (CRS) 与观测值相同。
spat_pred_grid <- expand.grid(
  lon = seq(-100, -80, length = 20),
  lat = seq(3246, length = 20)) %>% 
  SpatialPoints(proj4string = CRS(proj4string(STObj3)))
  
gridded(spat_pred_grid) <- TRUE
  • 对于我们的时间网格,我们考虑 1993 年 7 月的六个等距天。
temp_pred_grid <- as.Date("1993-07-01") + seq(328, length = 6)
  • 然后我们可以结合 spat_pred_gridtemp_pred_grid 来构建一个 STF即我们的时空预测网格的对象
DE_pred <- STF(sp = spat_pred_grid, # spatial part 
               time = temp_pred_grid) # temporal part

由于 STObj4 中存在缺失值,我们首先需要将 STObj4 转换为 STSDFSTIDF,并删除记录缺失值的数据。为简单起见,我们考虑 STIDF(考虑到 STSDF 的速度大约是其两倍)。此外,为了显示 S-T 克里格法跨时间预测的能力,我们从数据集中省略了 1993 年 7 月 14 日的数据

STObj5 <- as(STObj4[, -14], "STIDF"# convert to STIDF 
STObj5 <- subset(STObj5, !is.na(STObj5$z)) # remove missing data

现在我们可以使用 STObj5 作为我们的数据调用 krigeST

pred_kriged <- krigeST(z ~ 1 + lat, # latitude trend
                       data = STObj5, # data set w/o 14 July
                       newdata = DE_pred, # prediction grid
                       modelList = sepVgm, # semivariogram
                       computeVar = TRUE# compute variances

要绘制预测结果和附带的预测标准误差,可以直接使用函数stplot

  • 调色 首先,我们使用函数brewer.pal和函数colorRampPalette定义我们的调色板(关于这些函数的详细内容,请参见帮助文件)。
color_pal <- rev(colorRampPalette(brewer.pal(11"Spectral"))(16))
  • 绘制预测结果 我们用包含结果的对象调用stplot函数。
stplot(pred_kriged,
       main = "Predictions (degrees Fahrenheit)"
       layout = c(32),
       col.regions = color_pal)
图2.预测结果
图2.预测结果

使用R软件包gstat对1993年7月6天(每隔5天)围成的方形纬度块内最高温度(华氏度)的预测,其中1993年7月14日的数据在原始数据集中被省略了。

  • 绘制预测标准误差 预测(克里格法)的标准误差也可以用类似的方式绘制出来。
pred_kriged$se <- sqrt(pred_kriged$var1.var) 
stplot(pred_kriged[, , "se"],
       main = "Prediction std. errors (degrees Fahrenheit)"
       layout = c(32),
       col.regions = color_pal)

图3.预测标准误差 使用R软件包gstat对1993年7月6天(每隔5天)围成的方形纬度块内最高温度(华氏度)的预测标准误差,其中1993年7月14日的数据在原始数据集中被省略了。

本文的时空克里格法对于小数据集来说是比较快速和容易实现的,但是随着数据集的增大,它开始变得难以实现,除非使用一些近似的方法。例如,函数krigeST允许人们使用参数nmax来确定进行预测时使用的最大观测值数量。该预测器不再是最优的,但在许多实际情况下,它已经足够接近于最优预测器。

2Reference

Wikle, C. K., Zammit-Mangion, A., & Cressie, N. A. C. (2019). Spatio-temporal statistics with R. CRC Press, Taylor & Francis Group. Shi, X., & Yeung, D.-Y. (2018). Machine Learning for Spatiotemporal Sequence Forecasting: A Survey. 关雪峰, 曾宇媚, Xuefeng, G., & Yumei, Z. (2018). 时空大数据背景下并行数据处理分析挖掘的进展及趋势. 地理科学进展, 37(10), 1314–1327. 邓敏蔡建南, & DENG Min, C. J. (n.d.). 多模态地理大数据时空分析方法. 地球信息科学学报, 22(1), 41–56.

分类:

人工智能

标签:

数据挖掘

作者介绍

Rookie宇
V1

公众号: R语言vs科研