# Write a function that will return the cross validated error for one of the three
# snowload datasets we are comparing to.
snowload_trials <- function(tdata, message = FALSE){
# Pre-allocate a three dimensional array to store the errors.
errors <- matrix(0, nrow = nrow(tdata), ncol = 7)
### PRISM predictions
# Notice that we tune the predictions along the way for
# some parameters.
fun <- "prism"
preds <- try(crv_pred(tdata, formula = log(yr50) ~ ELEVATION,
fun = fun,
message = message,
distImp = 0.8,
minRad = c(30, 50, 75, 100),
wdistance = 2,
welevRange = list(lwr = c(100, 200, 300),
upr = c(500, 1000, 1500, 2000, 2500)),
welevation = c("ELEVATION", 1),
wbasin = c("HUC", 1:5),
bound = FALSE,
weights = FALSE,
tune_folds = 10,
tune_score = "MAE"), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# Record the prism CV errors
errors[, 1] <- tdata$yr50 - exp(preds)
### SKLM predictions
fun <- "rkriging"
preds <- try(crv_pred(tdata, formula = log(yr50) ~ ELEVATION,
fun = fun,
model = gstat::vgm("Sph"), sklm = TRUE,
debug.level = 0), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the sklm cv errors
errors[, 2] <- tdata$yr50 - exp(preds)
### UK predictions
fun <- "rkriging"
preds <- try(crv_pred(tdata, formula = log(yr50) ~ ELEVATION,
fun = fun,
model = gstat::vgm("Sph"), sklm = FALSE,
debug.level = 0), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the sklm cv errors
errors[, 3] <- tdata$yr50 - exp(preds)
### IDW_snow preds (as implemented in Idaho)
fun <- "idw_snow"
# Notice no log transform when doing NGSL predictions
preds <- try(crv_pred(tdata, formula = yr50 ~ ELEVATION,
fun = fun,
NGSL = TRUE, tlayer = 1219.2,
debug.level = 0), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the idw_snow cv errors
errors[, 4] <- tdata$yr50 - preds
### tri_snow preds
fun <- "tri_snow"
# Notice no log transform when doing NGSL predictions
preds <- try(crv_pred(tdata, formula = yr50 ~ ELEVATION,
fun = fun), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the tri_snow cv errors
errors[, 5] <- tdata$yr50 - preds
### lr_snow preds
fun <- "lm_snow"
preds <- try(crv_pred(tdata, formula = log(yr50) ~ ELEVATION,
fun = fun), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the lr_snow cv errors
errors[, 6] <- tdata$yr50 - exp(preds)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
### snlwf preds
fun <- "snlwf"
preds <- try(crv_pred(tdata, formula = ~ELEVATION,
fun = fun), silent = TRUE)
if(inherits(preds, "try-error")){
print(paste("Method", fun, "failed to predict..."))
preds <- rep(NA, nrow(tdata))
}
# record the snlwf cv errors
errors[, 7] <- tdata$yr50 - preds
return(errors)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.