R中`gstat`包中的交叉验证和方差计算

Question

美好的一天，

在尝试计算gstat包中的反距离krig的方差时，我遇到了一些困难。我还想在一个独立的变量测试集上运行交叉验证，但我不确定如何在R中使用空间数据。使用meuse数据集，这是我尝试计算方差的原因：

data(meuse); coordinates(meuse) <- ~x+y 

#randomly sample to get training and test data for later cross-validation
set.seed = (123)

sub1 <- nrow(meuse@data); len1 <- ceiling(sub1*2/3)

m.train <- meuse
m.train@data <- meuse@data[1:len1,]
m.train@coords <- meuse@coords[1:len1,]

m.test <- meuse
m.test@data <- meuse@data[(len1+1):sub1,]
m.test@coords  <- meuse@coords[(len1+1):sub1,]

 ## load grids:
data(meuse.grid); coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE; fullgrid(meuse.grid) <- TRUE

zinc.id <- krige(zinc~1, m.train, meuse.grid) ## inverse distance weighting

# --- My attempt at calculation of variance
rmse.id <- sqrt(mean((meuse.test@data$zinc - zinc.id@data$var1.pred)^2))

Warning message:
In meuse.test@data$z - zinc.id@data$var1.pred :
longer object length is not a multiple of shorter object length

我可以看到为什么我收到错误，但我不知道如何继续。我可以在R之外进行交叉验证，但有点麻烦，但我真的想让我在R工作。任何建议都会受到欢迎。

库尔特

Answer 1

要进行这种比较，你需要使用meuse而不是meuse.grid作为newdata。甚至更好，使用krige.cv。

例如，使用meuse数据集：

kr_cv = krige.cv(log(zinc)~1, meuse, vgm(.59, "Sph", 874, .04))
kr_cv[1:5,]
       coordinates var1.pred  var1.var observed    residual      zscore fold
1 (181072, 333611)  6.784729 0.1681011 6.929517  0.14478795  0.35314023    1
2 (181025, 333558)  6.777372 0.1635077 7.039660  0.26228828  0.64864901    2
3 (181165, 333537)  6.294508 0.1723531 6.461468  0.16696067  0.40216530    3
4 (181298, 333484)  6.033072 0.2191244 5.549076 -0.48399603 -1.03394256    4
5 (181307, 333330)  5.576879 0.1643513 5.594711  0.01783242  0.04398694    5

由此您可以轻松计算交叉验证的RMSE。 automap包（免责声明：我写的）包含一个方便的功能，可以为你计算很多这些统计数据。通常它只接受autoKrige.cv的输出，但是使用小hack你仍然可以使用它：

library(automap)
compare.cv(list(krige.cv_output = kr_cv))
            krige.cv_output
mean_error        0.0003146
me_mean           5.345e-05
MAE                  0.2898
MSE                  0.1515
MSNE                 0.8607
cor_obspred          0.8416
cor_predres         0.05449
RMSE                 0.3892
RMSE_sd              0.5391
URMSE                0.3892
iqr                  0.3949