Newer
Older
#'Calculates the Anomaly Correlation Coefficient for the ensemble mean of
#'each model and the corresponding references for each startdate and each
#'leadtime.\cr
#'The domain of interest can be specified by providing the list
#'of longitudes/latitudes (lon/lat) of the grid together with the corners
#'of the domain:\cr
#' lonlatbox = c(lonmin, lonmax, latmin, latmax).
#'@param var_exp Array of experimental anomalies with dimensions:
#' c(nexp, nsdates, nltimes, nlat, nlon) or
#' c(nexp, nsdates, nmembers, nltimes, nlat, nlon).
#'@param var_obs Array of observational anomalies, same dimensions as var_exp
#' except along the first dimension and the second if it corresponds to the
#' member dimension.
#'@param lon Array of longitudes of the var_exp/var_obs grids, optional.
#'@param lat Array of latitudes of the var_exp/var_obs grids, optional.
#'@param lonlatbox Domain to select: c(lonmin, lonmax, latmin, latmax),
#' optional.
#'@param conf TRUE/FALSE: confidence intervals and significance level
#' provided or not.
#'@param conftype "Parametric" provides a confidence interval for the ACC
#' computed by a Fisher transformation and a significance level for the ACC
#' from a one-sided student-T distribution. "Bootstrap" provides a confidence
#' interval for the ACC and MACC computed from bootstrapping on the members
#' with 100 drawings with replacement. To guarantee the statistical robustness
#' of the result, make sure that your experiments/oservations/startdates/
#' leadtimes always have the same number of members.
#'@param siglev The confidence level for the computed confidence intervals.
#' If \code{conf = TRUE}, array with dimensions:
#' c(nexp, nobs, nsdates, nleadtimes, 4)
#' The fifth dimension of length 4 corresponds to the lower limit of the
#' \code{siglev}\% confidence interval, the ACC, the upper limit of the
#' \code{siglev}\% confidence interval and the \code{siglev}\% significance
#' level. If \code{conf = FALSE}, the array of the Anomaly Correlation
#' Coefficient has dimensions c(nexp, nobs, nsdates, nleadtimes).
#'}
#'\item{MACC}{
#' The array of the Mean Anomaly Correlation Coefficient with dimensions
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#'
#'@examples
#'# See ?Load for explanations on the first part of this example.
#' \dontrun{
#'data_path <- system.file('sample_data', package = 's2dverification')
#'expA <- list(name = 'experiment', path = file.path(data_path,
#' 'model/$EXP_NAME$/$STORE_FREQ$_mean/$VAR_NAME$_3hourly',
#' '$VAR_NAME$_$START_DATE$.nc'))
#'obsX <- list(name = 'observation', path = file.path(data_path,
#' '$OBS_NAME$/$STORE_FREQ$_mean/$VAR_NAME$',
#' '$VAR_NAME$_$YEAR$$MONTH$.nc'))
#'
#'# Now we are ready to use Load().
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- Load('tos', list(expA), list(obsX), startDates,
#' leadtimemin = 1, leadtimemax = 4, output = 'lonlat',
#' latmin = 27, latmax = 48, lonmin = -12, lonmax = 40)
#' }
#' \dontshow{
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'),
#' c('observation'), startDates,
#' leadtimemin = 1,
#' leadtimemax = 4,
#' output = 'lonlat',
#' latmin = 27, latmax = 48,
#' lonmin = -12, lonmax = 40)
#' }
#'sampleData$mod <- Season(sampleData$mod, 4, 11, 12, 2)
#'sampleData$obs <- Season(sampleData$obs, 4, 11, 12, 2)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
#'acc <- ACC(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2))
#'@references Joliffe and Stephenson (2012). Forecast Verification: A
#' Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.
#'@keywords datagen
#' 0.1 - 2013-08 (V. Guemas) - Original code\cr
#' 1.0 - 2013-09 (N. Manubens) - Formatting to CRAN\cr
#' 1.1 - 2013-09 (C. Prodhomme) - optimization\cr
#' 1.2 - 2014-08 (V. Guemas) - Bug-fixes: handling of NA & selection of domain + Simplification of code\cr
#' 1.3.0 - 2014-08 (V. Guemas) - Boostrapping over members\cr
#' 1.3.1 - 2014-09 (C. Prodhomme) - Add comments and minor style changes\cr
#' 1.3.2 - 2015-02 (N. Manubens) - Fixed ACC documentation and examples
#'@importFrom abind abind
#'@importFrom stats qt qnorm quantile
ACC <- function(var_exp, var_obs, lon = NULL, lat = NULL,
lonlatbox = NULL, conf = TRUE, conftype = "parametric",
siglev = 0.95) {
Nicolau Manubens Gil
committed
dimsvar <- dim(var_exp)
if (length(dimsvar) == 5) {
checkfirst <- 2
}else if (length(dimsvar) == 6) {
checkfirst <- 3
nmembexp <- dimsvar[2]
nmembobs <- dim(var_obs)[2]
}else{
stop("var_exp & var_obs should have dimensions (nexp/nsobs, nsdates, nltimes, nlat, nlon)
or dimensions (nexp/nsobs, nmembers, nsdates, nltimes, nlat, nlon) ")
Nicolau Manubens Gil
committed
}
for (iind in checkfirst:length(dimsvar)) {
Virginie Guemas
committed
if (dim(var_obs)[iind] != dimsvar[iind]) {
stop("var_exp & var_obs must have same dimensions except the first one (number of experiments or number of observational datasets) ")
Nicolau Manubens Gil
committed
}
}
nexp <- dimsvar[1]
nobs <- dim(var_obs)[1]
nsdates <- dimsvar[checkfirst]
nltimes <- dimsvar[checkfirst+1]
nlat <- dimsvar[checkfirst+2]
nlon <- dimsvar[checkfirst+3]
Virginie Guemas
committed
if (is.null(lon) == FALSE & is.null(lat) == FALSE &
is.null(lonlatbox) == FALSE) {
Nicolau Manubens Gil
committed
for (jind in 1:2) {
while (lonlatbox[jind] < 0) {
lonlatbox[jind] <- lonlatbox[jind] + 360
}
while (lonlatbox[jind] > 360) {
lonlatbox[jind] <- lonlatbox[jind] - 360
}
}
indlon <- which((lon >= lonlatbox[1] & lon <= lonlatbox[2]) |
(lonlatbox[1] > lonlatbox[2] & (lon > lonlatbox[1] |
lon < lonlatbox[2])))
Nicolau Manubens Gil
committed
indlat <- which(lat >= lonlatbox[3] & lat <= lonlatbox[4])
} else {
indlon <- 1:nlon
indlat <- 1:nlat
}
if(conf == TRUE) {
ACC <- array(NA, dim = c(nexp, nobs, nsdates, nltimes, 4))
} else {
ACC <- array(NA, dim = c(nexp, nobs, nsdates, nltimes))
}
Virginie Guemas
committed
MACCaux <- array(0, dim = c(nexp, nobs, nsdates, nltimes, 3))
var_exp <- array(var_exp[,,,,indlat, indlon],
dim = c(nexp, nmembexp, nsdates, nltimes,
length(indlat), length(indlon)))
var_obs <- array(var_obs[,,,,indlat, indlon],
dim = c(nobs, nmembobs, nsdates, nltimes,
length(indlat), length(indlon)))
tmp01 <- Mean1Dim(var_exp,2)
tmp02 <- Mean1Dim(var_obs,2)
}else{
var_exp <- array(var_exp[,,,indlat, indlon],
dim = c(nexp, nsdates, nltimes,
length(indlat), length(indlon)))
var_obs <- array(var_obs[,,,indlat, indlon],
dim = c(nobs, nsdates, nltimes,
length(indlat), length(indlon)))
tmp01 <- var_exp
tmp02 <- var_obs
}
for( iobs in 1:nobs) {
for( iexp in 1:nexp) {
Virginie Guemas
committed
# tmp1 and tmp2 are splitted to handle NA before building tmp
tmp1 <- array(tmp01[iexp, , , , ], dim = c(1, nsdates, nltimes,
length(indlon) * length(indlat)))
tmp2 <- array(tmp02[iobs, , , , ], dim = c(1, nsdates, nltimes,
length(indlon) * length(indlat)))
Virginie Guemas
committed
# Variance(tmp1)should not take into account any point
Virginie Guemas
committed
# that is not available in tmp2 and therefore not accounted for
# in covariance(tmp1,tmp2) and vice-versa
Virginie Guemas
committed
tmp1[ is.na(tmp2) ] <- NA
tmp2[ is.na(tmp1) ] <- NA
tmp <- abind(tmp1, tmp2, along = 1)
top <- apply(tmp, c(2, 3), function(x)
sum(x[1, ]*x[2, ], na.rm = TRUE) )
bottom1 <- apply(tmp, c(2, 3), function(x)
sum(x[1, ]*x[1, ], na.rm = TRUE) )
bottom2 <- apply(tmp, c(2, 3), function(x)
sum(x[2, ]*x[2, ], na.rm = TRUE) )
bottom <- sqrt(bottom1 * bottom2 )
ACCaux <- top / bottom
#handle NA
tmpallNA <- which(is.na(bottom) | bottom == 0)
Virginie Guemas
committed
ACCaux[tmpallNA] <- NA
top[tmpallNA] = NA
bottom1[tmpallNA] = NA
bottom2[tmpallNA] = NA
Virginie Guemas
committed
#store the value to calculate the MACC
MACCaux[iexp, iobs, , , 1] <- top
MACCaux[iexp, iobs, , , 2] <- bottom1
MACCaux[iexp, iobs, , , 3] <- bottom2
Virginie Guemas
committed
if (conf == TRUE) {
ACC[iexp, iobs, , , 2] <- ACCaux
#calculate parametric confidence interval
Virginie Guemas
committed
eno <- Mean1Dim( Eno(tmp2, 4), 1)
t <- apply(eno, c(1, 2),
ACC[iexp, iobs, , , 4] <- apply(enot, c(1, 2), function(x)
sqrt((x[2] * x[2]) / ((x[2] * x[2]) + x[1] - 2)))
ACC[iexp, iobs, , , 1] <- apply(correno, c(1, 2), function(x)
tanh(atanh(x[1]) + qnorm(1 - (1 - siglev) / 2) / sqrt(x[2] - 3)))
ACC[iexp, iobs, , , 3] <- apply(correno, c(1, 2), function(x)
tanh(atanh(x[1]) + qnorm((1 - siglev) / 2) / sqrt(x[2] - 3)))
ACC[iexp, iobs, , ] <- ACCaux
Nicolau Manubens Gil
committed
}
Virginie Guemas
committed
Nicolau Manubens Gil
committed
}
}
# #na.rm should be TRUE to obtain a MACC even if a few
# #start dates are missing
Virginie Guemas
committed
topfinal <- apply(MACCaux, c(1, 2, 4), function(x)
Virginie Guemas
committed
bottomfinal <- apply(MACCaux, c(1, 2, 4), function(x)
sqrt(sum(x[, 2], na.rm = TRUE) * sum(x[, 3], na.rm = TRUE)))
#to avoid that some NA are called NaN or Inf
Virginie Guemas
committed
MACC <- topfinal / bottomfinal
MACC[tmpNA] <- NA
if (conf == TRUE & conftype == "bootstrap") {
if (length(dimsvar) != 6) {
stop("Var_exp and var_obs must have a member dimension")
}
ndraw <- 100
#create the matrix to store the random values
ACC_draw = array(dim=c(nexp,nobs,nsdates,nltimes,ndraw))
MACC_draw = array(dim=c(nexp,nobs,nltimes,ndraw))
#put the member dimension first
var_exp <- aperm(var_exp, c(2, 1, 3, 4, 5, 6))
var_obs <- aperm(var_obs, c(2, 1, 3, 4, 5, 6))
for (jdraw in 1:ndraw) {
#choose a randomly member index for each point of the matrix
indexp <- array(sample(nmembexp, size = (nexp*nmembexp*nsdates*nltimes),
replace = TRUE), dim = c(nmembexp, nexp, nsdates, nltimes,
length(indlat), length(indlon)) )
indobs <- array(sample(nmembobs, size = (nobs*nmembobs*nsdates*nltimes),
replace = TRUE), dim = c(nmembobs, nobs, nsdates, nltimes,
length(indlat), length(indlon)) )
#combine maxtrix of data and random index
varindexp <- abind(var_exp, indexp, along = 7 )
varindobs <- abind(var_obs, indobs, along = 7 )
#select randomly the members for each point of the matrix
varexpdraw <- aperm( array(
apply( varindexp, c(2, 3, 4, 5, 6), function(x) x[,1][x[,2]] ),
dim = c(nmembexp, nexp, nsdates, nltimes,
length(indlat), length(indlon))),
c(2, 1, 3, 4, 5, 6))
varobsdraw <- aperm( array(
apply( varindobs, c(2, 3, 4, 5, 6), function(x) x[,1][x[,2]] ),
dim = c(nmembobs, nobs, nsdates, nltimes,
length(indlat), length(indlon))),
c(2, 1, 3, 4, 5, 6))
#calculate the ACC of the randomized field
tmpACC <- ACC(varexpdraw, varobsdraw, conf = FALSE)
ACC_draw[,,,,jdraw] <- tmpACC$ACC
MACC_draw[,,,jdraw] <- tmpACC$MACC
#calculate the confidence interval
ACC[ , , , , 3] <- apply(ACC_draw, c(1, 2, 3, 4), function(x)
quantile(x, 1 - (1 - siglev) / 2, na.rm = TRUE))
Chloé Prodhomme
committed
ACC[ , , , , 1] <- apply(ACC_draw, c(1, 2, 3, 4), function(x)
MACC <- InsertDim(MACC, 4, 3)
MACC[ , , , 3] <- apply(MACC_draw, c(1, 2, 3), function(x)
quantile(x, 1 - (1 - siglev) / 2, na.rm = TRUE))
Chloé Prodhomme
committed
MACC[ , , , 1] <- apply(MACC_draw, c(1, 2, 3), function(x)
Nicolau Manubens Gil
committed
invisible(list(ACC = ACC, MACC = MACC))
}