5

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?

  • Well, what type of cross validation? Are you wanting to assess the observed vs predicted (eg., RMSE), perform an n-fold cross-validation or obtain a measure of uncertainty akin to a conditional simulation? Please be specific and detailed in your questions as well as indicating what you have already tried. – Jeffrey Evans Apr 10 '15 at 17:03
5

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 )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.