Блог

На основе приведенных здесь примеров:

Как выполнить уменьшение размерности с помощью PCA в R

и

Как отменить PCA и восстановить исходные переменные из нескольких основных компонентов?

Я пытаюсь выполнить PCA для растрового блока (с 69 слоями), затем получить ведущие ПК и, наконец, восстановить исходные переменные, используя только ПК с совокупной долей около ~ 95%.

 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 компьютеров, поскольку меня интересует их пространственное распределение:

 ### 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 шт. Тем не менее, я заканчиваю с 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, но первоначальная идея заключалась в том, чтобы иметь меньшее количество слоев, способных объяснить закономерности в исходных данных.
Спасибо!

Задавая вопрос здесь, пожалуйста, не указывайте на файл в dropbox, а приводите некоторые примеры данных, подобные этому:

 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)
  

Вы создаете Xhat путем подмножества PCS, но 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 компьютеров, поскольку меня интересует их пространственное распределение». Все переменные вносят свой вклад в PCS.

Возможно, вам действительно нужно plot(x) ?

Вы используете плохой код для создания растрового кирпича из Xhat . Вы можете сделать это вместо:

 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)