使用光栅中的PCA进行降维

Question

基于此处的示例：

How to perform dimensionality reduction with PCA in R

和

How to reverse PCA and reconstruct original variables from several principal components?

我正在尝试在光栅砖（具有69层）上执行PCA，然后获得领先的PC，最后仅使用累积比例约为95％的PC来重建原始变量。

library(raster)
library(ncdf4)

ln <- "https://www.dropbox.com/s/d88iuvp9oio14zk/test.nc?dl=1" # ~400 kb size

### DOWNLOAD THE FILE
download.file(ln,
              destfile="test.nc",
              method="auto")

st <- brick("test.nc")
nlayers(st)

### DO THE PCA
pca <- prcomp(st[]) 
# to visualize pcs as rasters
x <- predict(st, pca, index=1:4)
spplot(x) # there are the first 4 PCs explaining most of the data. 

然后，由于我对这些变量的空间分布感兴趣，因此我试图从前4台PC重构原始变量：

### PCA DETAILS
summary(pca) # importance of components
plot (pca) # scree plot
loadings(pca) #eigens 

mu <- colMeans(as.matrix(st)) # get the column means to use after

#### REDUCTION
nComp <- 4
Xhat <- pca$x[,1:nComp] %*% t(pca$rotation[,1:nComp])
Xhat <- scale(Xhat, center = -mu, scale = FALSE)

我以为我只会得到前4台PC。但是，我最终得到的像以前一样是69：

### CHECK THE DIMENSIONS
dim(Xhat)

### THEN CREATE THE RASTER WITH THE PCs
coords <- coordinates(st[[1]]) # get the lon/lat

rst <- cbind(coords, Xhat) # bind the coordinates 
rst <- rasterFromXYZ(rst) # create the raster
plot(rst)

我在这里错过了什么？我不是PCA方面的专家，但最初的想法是使层数更少，能够解释原始数据中的模式。谢谢！

Answer 1

在此处提问时，请不要指向投寄箱上的文件，而应包含以下示例数据：

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
s <- stack(b, flip(b, "y"), setValues(raster(b), runif(ncell(b))))
names(s) <- paste0("var", 1:nlayers(s))    

pca <- prcomp(values(s)) 
x <- predict(s, pca, index=1:4)

您通过设置PC来创建Xhat，但是pca $ rotation具有所有变量

round(pca$rotation[,1:nComp],1)
      PC1  PC2  PC3  PC4
var1 -0.4  0.4 -0.4  0.4
var2 -0.4  0.4 -0.2  0.2
var3 -0.4  0.4  0.6 -0.6
var4 -0.4 -0.4 -0.4 -0.4
var5 -0.4 -0.4 -0.2 -0.2
var6 -0.4 -0.4  0.6  0.6
var7  0.0  0.0  0.0  0.0

这是有道理的，因为您说的是要“从前4台PC重构原始变量，因为我对这些PC的空间分布感兴趣”。所有变量都影响PC。

您可能真正想要的是plot(x)？

您正在使用不良代码从Xhat创建RasterBrick。您可以改为：

Xhat <- pca$x[,1:nComp] %*% t(pca$rotation[,1:nComp])
Xhat <- scale(Xhat, scale = FALSE)

b <- brick(s, values=FALSE)
b <- setValues(b, Xhat)
b
#class      : RasterBrick 
#dimensions : 77, 101, 7777, 7  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : memory
#names      :          var1,          var2,          var3,          var4,          var5,          var6,          var7 
#min values : -184.32344039, -184.48714657, -193.05823803, -184.32341010, -184.48718831, -193.05827512,   -0.01466663 
#max values :   73.33872354,   70.38724578,   63.48912986,   73.33875039,   70.38723822,   63.48906605,    0.01193009 

比较b和s

m <- cellStats(s, mean)
bb <- b + m
plot(bb, s)