Styjun

Good day,

I am getting some difficulty when trying to calculate the variance from an inverse distance krig done in the gstat package. I would also like to run a cross-validation on an independent test set of variables, but I am not sure of how to do so in R with spatial data. Using the meuse dataset, this is what I attempted for calculating variance:

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

I can see why I am getting the error, but I am not sure how to proceed. I can do the cross-validation outside of R with a bit of trouble, but I really would like to keep all my working within R. Any suggestions would be most welcomed.

Kurt

To perform this kind of comparison, you need to use meuse and not meuse.grid as the newdata. Or even better, use krige.cv.

For example using the meuse dataset:

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

From this you can easily calculate the RMSE of the cross-validation. The automap package (disclaimer: which I wrote) contains a convienient function that can calculate a lot of these stats for you. Normally it only accepts the output of autoKrige.cv, but using a small hack you can still use it:

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