时间序列数据的主成分分析（PCA）：空间和时间模式

Question

假设我有 100 个站从 1951 年到 1980 年的年降水量数据。在一些论文中，我发现人们将 PCA 应用于时间序列，然后绘制空间载荷图（值从 -1 到 1），并且还绘制了 PC 的时间序列。例如https://publicaciones.unirioja.es/ojs/index.php/cig/article/view/2931/2696中的图6就是PC的空间分布

我在 R 中使用函数 prcomp，我想知道如何做同样的事情。换句话说，如何从 prcomp 函数的结果中提取 "spatial pattern" 和 "temporal pattern"？谢谢。

set.seed(1234)
rainfall = sample(x=100:1000,size = 100*30,replace = T)
rainfall=matrix(rainfall,nrow=100)
colnames(rainfall)=1951:1980
PCA = prcomp(rainfall,retx=T)

或者https://1drv.ms/u/s!AnVl_zW00EHegxAprS4s7PDaYQVr

处有真实数据

Answer 1

"Temporal pattern"解释了时间序列在所有网格中占主导地位的时间变化，用PCA的主成分（PCs，多个时间序列）表示.在 R 中，它是 prcomp(data)$x[,'PC1'] 最重要的 PC，PC1。

"Spatial pattern" 解释了 PC 对某些变量（在您的情况下是地理位置）的依赖程度，它由每个主要成分的负载表示。例如，对于 PC1，它是 prcomp(data)$rotation[,'PC1'].

这里是一个示例，使用您的数据在 R 中构建时空数据的 PCA 并显示时间变化和空间异质性。

首先，必须将数据转换为具有变量（空间网格）和观测值（yyyy-mm）的data.frame。

加载和转换数据：

load('spei03_df.rdata') 
str(spei03_df) # the time dimension is saved as names (in yyyy-mm format) in the list
lat <- spei03_df$lat # latitude of each values of data
lon <- spei03_df$lon # longitude
rainfall <- spei03_df 
rainfall$lat <- NULL
rainfall$lon <- NULL
date <- names(rainfall)
rainfall <- t(as.data.frame(rainfall)) # columns are where the values belong, rows are the times

要了解数据，请在地图上绘制 1950 年 1 月的数据：

library(mapdata)
library(ggplot2) # for map drawing

drawing <- function(data, map, lonlim = c(-180,180), latlim = c(-90,90)) {
  major.label.x = c("180", "150W", "120W", "90W", "60W", "30W", "0", 
                    "30E", "60E", "90E", "120E", "150E", "180")
  major.breaks.x <- seq(-180,180,by = 30)
  minor.breaks.x <- seq(-180,180,by = 10)

  major.label.y = c("90S","60S","30S","0","30N","60N","90N")
  major.breaks.y <- seq(-90,90,by = 30)
  minor.breaks.y <- seq(-90,90,by = 10)
  panel.expand <- c(0,0)

  drawing <- ggplot() + 
    geom_path(aes(x = long, y = lat, group = group), data = map) + 
    geom_tile(data = data, aes(x = lon, y = lat, fill = val), alpha = 0.3, height = 2) + 
    scale_fill_gradient(low = 'white', high = 'red') + 
    scale_x_continuous(breaks = major.breaks.x, minor_breaks = minor.breaks.x, labels = major.label.x, 
                       expand = panel.expand,limits = lonlim) +
    scale_y_continuous(breaks = major.breaks.y, minor_breaks = minor.breaks.y, labels = major.label.y,
                       expand = panel.expand, limits = latlim) +
    theme(panel.grid = element_blank(), panel.background = element_blank(),
          panel.border = element_rect(fill = NA, color = 'black'),
          axis.ticks.length = unit(3,"mm"),
          axis.title = element_text(size = 0),
          legend.key.height = unit(1.5,"cm"))

  return(drawing)
}

map.global <- fortify(map(fill=TRUE, plot=FALSE))
dat <- data.frame(lon = lon, lat = lat, val = rainfall["1950-01",])
sample_plot <- drawing(dat, map.global, lonlim = c(-180,180), c(-90,90))
ggsave("sample_plot.png", sample_plot,width = 6,height=4,units = "in",dpi = 600)

如上所示，link提供的网格数据包括代表加拿大降雨量（某种指数？）的值。

主成分分析：

PCArainfall <- prcomp(rainfall, scale = TRUE)
summaryPCArainfall <- summary(PCArainfall)
summaryPCArainfall$importance[,c(1,2)]

这表明前两个 PC 解释了降雨数据中 10.5% 和 9.2% 的方差。

我提取了前两台 PC 的负载和 PC 时间序列本身： "spatial pattern"（载荷），显示趋势强度的空间异质性（PC1 和 PC2）。

loading.PC1 <- data.frame(lon=lon,lat=lat,val=PCArainfall$rotation[,'PC1'])
loading.PC2 <- data.frame(lon=lon,lat=lat,val=PCArainfall$rotation[,'PC2'])

drawing.loadingPC1 <- drawing(loading.PC1,map.global, lonlim = c(-180,-30), latlim = c(40,90)) + ggtitle("PC1")
drawing.loadingPC2 <- drawing(loading.PC2,map.global, lonlim = c(-180,-30), latlim = c(40,90)) + ggtitle("PC2")

ggsave("loading_PC1.png",drawing.loadingPC1,width = 6,height=4,units = "in",dpi = 600)
ggsave("loading_PC2.png",drawing.loadingPC2,width = 6,height=4,units = "in",dpi = 600)

"temporal pattern"，前两个 PC 时间序列，显示数据的主要时间趋势

library(xts)
PC1 <- ts(PCArainfall$x[,'PC1'],start=c(1950,1),end=c(2014,12),frequency = 12)
PC2 <- ts(PCArainfall$x[,'PC2'],start=c(1950,1),end=c(2014,12),frequency = 12)
png("PC-ts.png",width = 6,height = 4,res = 600,units = "in")
plot(as.xts(PC1),major.format = "%Y-%b", type = 'l', ylim = c(-100, 100), main = "PC") # the black one is PC1
lines(as.xts(PC2),col='blue',type="l") # the blue one is PC2
dev.off()

但是，此示例绝不是您数据的最佳 PCA，因为 PC1 和 PC2 存在严重的季节性和年度变化（当然，夏季下雨更多，并查看个人电脑）。

您可以像文献中建议的那样，通过对数据进行去季节化或通过回归消除年度趋势来改进 PCA。但这已经超出了我们的主题。

时间序列数据的主成分分析（PCA）：空间和时间模式

Principal component analysis (PCA) of time series data: spatial and temporal pattern

r

spatial

temporal

pca