A quick demo of manual back transformation of geoR kriging predictions
library(geoR)
data(s100)
s100$data <- exp(s100$data)
reml <- likfit(s100, ini = c(0.5, 0.5), fix.nug = TRUE, lik.met = "REML",
lambda = 0.1, fix.lambda = F)
summary(reml)
s100$data <- BCtransform(s100$data, reml$lambda)$data
# Replace the lambda value in the model object. It is located elsewhere in
# the model object but this is the location the xvalid function looks for.
reml$parameters.summary["lambda", 2] <- 1
# Using geoR
xvalid <- xvalid(s100, model = reml)
# A function to return the LOOCV results for either a gstat::krige.cv or
# geoR::xvalid.
LOOCV <- function(Xvalid) {
ret <- list()
# This function now handles the results from either gstat or geoR.
if (class(Xvalid) == "data.frame") {
ret$RMSE <- sqrt(sum((Xvalid$observed - Xvalid$var1.pred)^2)/length(Xvalid$observed))
ret$ME <- sum((Xvalid$observed - Xvalid$var1.pred))/length(Xvalid$observed)
ret$thetaMean <- mean(((Xvalid$observed - Xvalid$var1.pred)^2)/Xvalid$var1.var)
ret$thetaMedian <- median(((Xvalid$observed - Xvalid$var1.pred)^2)/Xvalid$var1.var)
} else {
ret$RMSE <- sqrt(sum((Xvalid$data - Xvalid$predicted)^2)/length(Xvalid$data))
ret$ME <- sum((Xvalid$data - Xvalid$predicted))/length(Xvalid$data)
ret$thetaMean <- mean(((Xvalid$data - Xvalid$predicted)^2)/Xvalid$krige.var)
ret$thetaMedian <- median(((Xvalid$data - Xvalid$predicted)^2)/Xvalid$krige.var)
}
return(ret)
}
LOOCV(xvalid)
# Back transform the predicted values.
backTransformed <- backtransform.moments(lambda = reml$lambda, mean = xvalid$predicted,
variance = xvalid$krige.var)
# then you can insert the backtransformed values back into the object
xvalid$predicted <- backTransformed$mean
xvalid$krige.var <- backTransformed$variance
plot(xvalid) # This would make more sense with actual kriging rather then cross validation.