I am performing kriging on a large number of relatively modest data sets (144 data points on a 9 x 16 grid in each data set).
I was experimenting with different variogram models and other parameter tweaks to record errors, and noticed some data sets produced very large errors. Investigating these data sets revealed that in some cases the kriging predictions produce a very strange surface, and do not even pass through the original data (kriging predictions should have 0 error at locations in the source data).
There isn't any error message and the same data sets work perfectly well with alternate covariance models that appear to be identical fits to the empirical variogram.
Example to reproduce what I'm seeing:
# sessionInfo()
# R version 3.0.2 (2013-09-25)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
#
# locale:
# [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C
# [5] LC_TIME=English_United States.1252
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] rgl_0.93.991 geoR_1.7-4 MASS_7.3-29 sp_1.0-14
#
# loaded via a namespace (and not attached):
# [1] colorspace_1.2-4 dichromat_2.0-0 digest_0.6.4 ggplot2_0.9.3.1 grid_3.0.2 gtable_0.1.2 labeling_0.2 lattice_0.20-23
# [9] munsell_0.4.2 plyr_1.8 proto_0.3-10 RandomFields_3.0.5 RColorBrewer_1.0-5 reshape2_1.2.2 scales_0.2.3 splancs_2.01-34
# [17] stringr_0.6.2 tools_3.0.2
library(geoR)
library(rgl)
## This data came from an inclass example, but i'm seeing the same behavior on real data sets, though i don't really know if
## they are caused by the same circumstances.
dat <-data.frame(x1=c( 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 2),
x2 = c( 1, 2, 3, 4, 5, 1, 2, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 3),
Zi =c( 0.150599, 0.191034, 0.228643, 0.378379, 0.38171, 0.51886, 0.53333, 0.612982, 0.700318, 1.027502, 1.046067, 1.198059, 1.287741, 1.600192, 1.757163, 1.757786, 1.913516, 1.994275, 2.083261, 2.276455, 2.294944, 2.308114, 2.100058, 1.135489, 0.580503539))
## convert to geodata object
geodat <-as.geodata(dat)
## compute empirical variogram and plot
evariogram <-variog(geodat)
plot(evariogram)
## fit several models with default initial parameters.
exponentialfit <-variofit(evariogram, cov.model = "exponential")
linearfit <-variofit(evariogram, cov.model = "linear")
sphericalfit <-variofit(evariogram, cov.model = "spherical")
gaussianfit <-variofit(evariogram, cov.model = "gaussian")
## plot the fitted models
lines(exponentialfit, col = 'red')
lines(linearfit, col = 'black')
lines(sphericalfit, col = 'blue')
lines(gaussianfit, col = 'green')
## create krige.control objects from each model
ExpObj <-krige.control(obj.model = exponentialfit)
LinObj <-krige.control(obj.model = linearfit)
SpherObj <-krige.control(obj.model = sphericalfit)
GaussObj <-krige.control(obj.model = gaussianfit)
## predictions at a single known value
krigingG <-krige.conv(geodat, locations = c(2,3), krige = GaussObj)
krigingS <-krige.conv(geodat, locations = c(2,3), krige = SpherObj)
krigingL <-krige.conv(geodat, locations = c(2,3), krige = LinObj)
krigingE <-krige.conv(geodat, locations = c(2,3), krige = ExpObj)
## they are all identical (as they should be since the surface should pass exactly through the known values), except the prediction using
## the linear variogram model. It's not just the linear model, I have observed the same behavior in many different models in the actual
## data sets I'm using (NOAA/ESRL computer simulated environmental variables--data available in a sqaure grid of one degree resolution)
krigingG$predict
krigingS$predict
krigingL$predict
krigingE$predict
## some data to plot a smooth surface of the kriging predicitons
xgrid <-seq(1, 5, .05)
ygrid <-seq(1, 5, .05)
gridLoc <-expand.grid(xgrid, ygrid)
## make predictions along a fine grid to see a smooth prediction surface
krigingGS <-krige.conv(geodat, locations = gridLoc, krige = GaussObj)
krigingSS <-krige.conv(geodat, locations = gridLoc, krige = SpherObj)
krigingLS <-krige.conv(geodat, locations = gridLoc, krige = LinObj)
krigingES <-krige.conv(geodat, locations = gridLoc, krige = ExpObj)
## create a matrix of the prediction values (required for rgl surface3d function)
zlinear <-matrix(krigingLS$predict, nrow = length(xgrid), ncol = length(ygrid))
zgaussian <-matrix(krigingGS$predict, nrow = length(xgrid), ncol = length(ygrid))
zspherical <-matrix(krigingSS$predict, nrow = length(xgrid), ncol = length(ygrid))
zexponential <-matrix(krigingES$predict, nrow = length(xgrid), ncol = length(ygrid))
## plot the data
plot3d(x = dat$x1, y = dat$x2, z=dat$Zi, size = 5)
## plot the linear model prediction surface
surface3d(xgrid, ygrid, zlinear, back = "lines", col = "red", alpha = .4)
## plot the gaussian model prediction surface
surface3d(xgrid, ygrid, zgaussian, back = "lines", col = "blue", alpha = .4)
## plot the exponential model prediction surface
surface3d(xgrid, ygrid, zexponential, back = "lines", col = "green", alpha = .4)
## plot the spherical model prediction surface
surface3d(xgrid, ygrid, zspherical, back = "lines", col = "yellow", alpha = .4)
If you see the same things I see, then the Gaussian, spherical and exponential surfaces are nearly identical, but the linear surface is a spiky nightmare-inducing sort of thing.
Any insights? A bug? Unfortunate data? Or - say it isn't so - am I doing something silly?
