r: уменьшение размерности с помощью PCA в растровом кирпиче - PullRequest
Новая
       16

r: уменьшение размерности с помощью PCA в растровом кирпиче

0 голосов
/ 18 марта 2019

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

[https://stats.stackexchange.com/questions/57467/how-to-perform-dimensionality-reduction-with-pca-in-r/57478#57478][1]

и

[https://stats.stackexchange.com/questions/229092/how-to-reverse-pca-and-reconstruct-original-variables-from-several-principal-com][2]

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

1 Ответ

1 голос
/ 19 марта 2019

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

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 путем подмножества ПК, но pca $ вращение имеет все переменные 

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

То, что вы действительно хотите, это plot(x)?

Вы используете плохой код для создания RasterBrick из 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

Сравните б и с 

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