Nothing
#' Verify how well NASA-POWER Data Represent Observed Data
#'
#' Calculates scalar measures of accuracy.
#'
#' @param obs_est
#' A 2-column matrix. The reference or observed and the estimated or predicted
#' data. See \code{ObsEst} object as an example.
#' @param conf.int
#' A character variable (\code{Yes} or \code{No}) defining if the function must
#' calculate confidence intervals. Default is \dQuote{Yes}.
#' @param sig.level
#' A numeric variable (between 0.90 and 0.95) defining the significance level
#' for parameter the confidence intervals. Default is 0.95.
#'
#' @return
#' An object of class \code{PowerSDI.Accuracy}, a \code{list}, which contains:
#'
#' \itemize{
#' \item Absolute mean error (AME),
#' \item square root of the mean squared error (RMSE),
#' \item Willmott's indices of agreement:
#' \itemize{
#' \item original (dorig),
#' \item modified (dmod) and
#' \item refined (dref)
#' }
#' \item Pearson determination coefficient (R2), and
#' \item if \code{conf.int = "Yes"}, confidence intervals.
#' }
#'
#' @examples
#' a <- Accuracy(obs_est = ObsEst, conf.int = "No")
#' a
#'
#' # A generic plotting method is also supplied
#' plot(a)
#'
#' @export
Accuracy <- function(obs_est,
conf.int = "Yes",
sig.level = 0.95) {
# convert conf.int to lowercase for easy checking
conf.int <- tolower(conf.int)
check.confint(conf.int, sig.level)
o <- obs_est[, 1]
p <- obs_est[, 2]
N <- length(o)
if (any(is.na(o)) || (any(is.na(p)))) {
stop("Missing data are not allowed.",
"Please, remove them from the input file object, `obs_est`.",
call. = FALSE)
}
min.length <- length(as.matrix(obs_est)) / 2
if (min.length < 10) {
stop(
"You must have at least 10 pairs of ObsVsEst data.
Please, provide a longer period.",
call. = FALSE
)
}
databoot <- matrix(NA, N, 2)
Nboots <- 10000
dorigboot <- matrix(NA, Nboots, 1)
drefboot <- matrix(NA, Nboots, 1)
dmodboot <- matrix(NA, Nboots, 1)
AMEboot <- matrix(NA, Nboots, 1)
RMSEboot <- matrix(NA, Nboots, 1)
RQuadboot <- matrix(NA, Nboots, 1)
# AME
AME <- sum((abs(p - o))) / N
# RMSE
RMSE <- sqrt(sum(((p - o) ^ 2)) / N)
# Original Willmott index
Num.orig <- sum((o - p) ^ 2)
Den.orig <- sum((abs(p - mean(o)) + abs(o - mean(o))) ^ 2)
dorig <- 1 - (Num.orig / Den.orig)
# Modified Willmott index
Num.mod <- sum(abs(o - p))
Den.mod <- sum(abs(p - mean(o)) + abs(o - mean(o)))
dmod <- 1 - (Num.mod / Den.mod)
# Refined Willmott's index
if (abs(sum(p - o)) <= 2 * sum(abs(o - mean(o)))) {
dref <- 1 - ((sum(abs(p - o))) / (2 * sum(abs(o - mean(
o
)))))
} else {
(dref <- ((2 * sum(abs(
o - mean(o)
))) / sum(abs(p - o))) - 1)
}
# Rquad
RQuad <- (cor(o, p, method = "pearson"))^2
if (conf.int == "yes") {
message("Just a sec")
# Bootstraping
sig.level1 <- sig.level / 2
for (i in 1:Nboots) {
databoot <- obs_est[sample(nrow(obs_est), replace = TRUE), ]
oboot <- as.matrix(databoot[, 1])
pboot <- as.matrix(databoot[, 2])
# AME
AMEboot[i, 1] <- sum((abs(pboot - oboot))) / N
RMSEboot[i, 1] <- sqrt(sum(((pboot - oboot) ^ 2)) / N)
# Original Willmott index
Numboot <- sum((oboot - pboot) ^ 2)
Denboot <-
sum((abs(pboot - mean(oboot)) + abs(oboot - mean(oboot))) ^ 2)
dorigboot[i, 1] <- 1 - ((Numboot) / Denboot)
# Modified Willmott index
Numboot.mod <- sum(abs(oboot - pboot))
Denboot.mod <-
sum(abs(pboot - mean(oboot)) + abs(oboot - mean(oboot)))
dmodboot[i, 1] <- 1 - (Numboot.mod / Denboot.mod)
# Refined Willmott index
if (abs(sum(pboot - oboot)) <= 2 * sum(abs(oboot - mean(oboot)))) {
drefboot[i, 1] <-
1 - ((sum(abs(
pboot - oboot
))) / (2 * sum(abs(
oboot - mean(oboot)
))))
} else {
(drefboot[i, 1] <-
((2 * sum(
abs(oboot - mean(oboot))
)) / sum(abs(
pboot - oboot
))) - 1)
}
# Rquad
RQuadboot[i, 1] <-
(cor(oboot, pboot, method = "pearson"))^2
}
# Defining confidence intervals
AME_CIinf <-
quantile(AMEboot, probs = sig.level1, na.rm = TRUE)
AME_CIsup <-
quantile(AMEboot,
probs = (1 - sig.level1),
na.rm = TRUE)
RMSE_CIinf <-
quantile(RMSEboot, probs = sig.level1, na.rm = TRUE)
RMSE_CIsup <-
quantile(RMSEboot,
probs = (1 - sig.level1),
na.rm = TRUE)
dorig_CIinf <-
quantile(dorigboot, probs = sig.level1, na.rm = TRUE)
dorig_CIsup <-
quantile(dorigboot,
probs = (1 - sig.level1),
na.rm = TRUE)
dmod_CIinf <-
quantile(dmodboot, probs = sig.level1, na.rm = TRUE)
dmod_CIsup <-
quantile(dmodboot,
probs = (1 - sig.level1),
na.rm = TRUE)
dref_CIinf <-
quantile(drefboot, probs = sig.level1, na.rm = TRUE)
dref_CIsup <-
quantile(drefboot,
probs = (1 - sig.level1),
na.rm = TRUE)
RQuad_CIinf <-
quantile(RQuadboot, probs = sig.level1, na.rm = TRUE)
RQuad_CIsup <-
quantile(RQuadboot,
probs = (1 - sig.level1),
na.rm = TRUE)
ModelAccuracy <- cbind(
AME_CIinf,
AME,
AME_CIsup,
RMSE_CIinf,
RMSE,
RMSE_CIsup,
dorig_CIinf,
dorig,
dorig_CIsup,
dmod_CIinf,
dmod,
dmod_CIsup,
dref_CIinf,
dref,
dref_CIsup,
RQuad_CIinf,
RQuad,
RQuad_CIsup
)
colnames(ModelAccuracy) <- c(
"AMECIinf",
"AME",
"AMECIsup",
"RMSECIinf",
"RMSE",
"RMSECIsup",
"dorig_CIinf",
"dorig",
"dorigCIsup",
"dmodCIinf",
"dmod",
"dmodCIsup",
"dref_CIinf",
"dref",
"dref_CIsup",
"RQuad_CIinf",
"RQuad",
"RQuad_CIsup"
)
row.names(ModelAccuracy) <- c("")
} else {
ModelAccuracy <- cbind(AME, RMSE, dorig, dmod, dref, RQuad)
colnames(ModelAccuracy) <-
c("AME", "RMSE", "dorig", "dmod", "dref", "RQuad")
row.names(ModelAccuracy) <- c("")
}
out <-
list(o, p, ModelAccuracy, class = "PowerSDI.Accuracy")
names(out) <- c("o", "p", "ModelAccuracy")
class(out) <- union("PowerSDI.Accuracy", class(out))
return(out)
}
#' Prints PowerSDI.Accuracy Objects
#'
#' Custom \code{print()} method for \code{PowerSDI.Accuracy} objects.
#'
#' @param x a \code{PowerSDI.Accuracy} object
#' @param digits The number of digits to be used after the decimal when
#' displaying accuracy values.
#' @param ... ignored
#' @return Nothing. Side-effect: pretty prints a \code{PowerSDI.Accuracy} object
#' in the \R console.
#' @export
#'
print.PowerSDI.Accuracy <- function(x,
digits = max(3L, getOption("digits") - 3L),
...) {
print(x$ModelAccuracy)
invisible(x)
}
#' Plots PowerSDI.Accuracy Objects
#'
#' Custom \code{plot()} method for \code{PowerSDI.Accuracy} objects.
#'
#' @param x a `PowerSDI.Accuracy` object
#' @param ... Other parameters as passed to \code{plot()}
#' @return No return value, called for side effects. Using this will display a
#' scatter plot of reference ETP data (x-axis) and estimated ETP (y-axis) in
#' the active \R session.
#' @export
#'
plot.PowerSDI.Accuracy <- function(x, ...) {
plot(
x = x$o,
y = x$p,
main = "",
xlab = "Reference",
ylab = "Estimated"
)
}
check.confint <- function(conf.int, sig.level) {
# fail early
if (conf.int != "yes" && conf.int != "no") {
stop(
"`conf.int` should be set to either 'Yes' or 'No'.\n
If 'Yes', the `sig.level` must be defined (a value between 0.9 and 0.95;
default is 0.95).",
call. = FALSE
)
}
if (conf.int == "yes") {
check.sig.level(sig.level)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.