Cross-validation of TPS and kriging outputs from the fields package in R

Question

I am looking for existing package or functions that will allow me to cross-validate multiple outputs of Thin plate spline (TPS) and Kriging from the fields package.

Does anyone know of any available package or function in R that would do this?

Answer 1

If you are just after a metric of performance, this is a fairly straight forward type of analysis to specify. 1) specify a model, 2) predict the model(s) at the data, 3) apply an accuracy/performance metric based on observed vs. predicted. We can step through the process thus (note; this is a dummy example so, ignore the REML errors):

First lets specify Kriging and Thin Plate Spline models.

library(fields)
dat <- data.frame(ozone[["x"]], y=ozone[["y"]])
dat <- rbind( dat, data.frame(East.West=jitter(dat$East.West), 
              North.South=jitter(dat$North.South),
              y=jitter(dat$y)))

krig.fit <- Krig(dat[,1:2], dat$y, theta=20)  
tps.fit <- Tps(dat[,1:2], dat$y)

We then predict the models at the data to create a data.frame of the observed data and predictions for each model.

obs.pred <- data.frame(y=dat$y, krig=predict(krig.fit),
                       tps=predict(tps.fit) )

We can then write a function for the desired evaluation metric, in this case Root Mean Squared Error (RMSE), and apply it to the observed vs predicted data.

# RMSE function 
rmse <- function(x, y) sqrt(mean((x - y)^2, na.rm=TRUE))

# Kriging RMSE
rmse(obs.pred$krig, obs.pred$y) 

# Spline RMSE
rmse(obs.pred$tps, obs.pred$y)

If you wanted to implement a data withhold type approach you could use sample in a for loop to accumulate error against iteratively withheld data. This has the advantage of assessing the sensitivity of your models to the data.

tps.rmse <- vector()    
n <- 5
  for(i in 1:100) {
    s <- sample(1:dim(dat)[1],n)
    obs <- dat[s,]
    m.dat <- dat[-s,]
    md <- Tps(m.dat[,1:2], m.dat$y) 
    pred <- predict(md)[s]
    tps.rmse <- append(tps.rmse, rmse(pred,obs$y))
}
summary( tps.rmse )