时间序列数据的主成分分析(PCA):空间和时间模式

Question

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

我prcomp在R中使用函数,我想知道我怎么能做同样的事情.换句话说,如何从prcomp函数的结果中提取"空间模式"和"时间模式" ？谢谢.

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

"时间模式"解释了所有网格中时间序列的主要时间变化,它由PCA 的主成分(PC,多个时间序列)表示.在R中,它是prcomp(data)$x[,'PC1']最重要的PC,PC1.

"空间模式"解释了PC依赖于某些变量(在您的情况下是地理位置)的强大程度,并且它由每个主要组件的负载表示.例如,对于PC1,它是prcomp(data)$rotation[,'PC1'].

以下是使用您的数据构建R中时空数据的PCA并显示时间变化和空间异质性的示例.

首先,必须将数据转换为具有变量(空间网格)和观察(yyyy-mm)的数据框架.

加载和转换数据:

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)

如上所示,所提供的链接给出的网格化数据包括代表加拿大降雨量(某些类型的指数？)的值.

主成分分析:

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

它表明前两台PC解释了降雨量数据的10.5%和9.2%的变化.

我提取前两台PC的负载和PC时间序列本身:"空间模式"(负载),显示趋势强度(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)

"时间模式",即前两个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存在严重的季节性和年度变化(当然,夏季降雨量更多,看看PC的弱尾).

您可以通过延迟数据或通过回归消除年度趋势来改善PCA,如文献所示.但这已经超出了我们的主题.