Newer
Older
ACC <- function(var_exp, var_obs, lon = NULL, lat = NULL,
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
lonlatbox = NULL, conf = TRUE) {
# Matrix var_exp & var_obs should have dimensions (nexp/nobs, nsdates,
# nltimes, nlat, nlon).
# ACC computes the Anomaly Correlation Coefficient for each jexp in
# 1:nexp and each jobs in 1:nobs which gives nexp x nobs ACC for each
# startdate and each leadtime.
# A domain can be selected by providing the list of longitudes/latitudes
# (lon/lat) of the grid together with the corner of the domain:
# lonlatbox = c(lonmin, lonmax, latmin, latmax)
#
# Args:
# var_exp: Matrix of experimental data.
# var_obs: Matrix of observational data, same dimensions as var_exp except
# along the first dimension.
# lon: Array of longitudes of the var_exp/var_obs grids, optional.
# lat: Array of latitudes of the var_exp/var_obs grids, optional.
# lonlatbox: Domain to select c(lonmin, lonmax, latmin, latmax), optional.
# conf: TRUE/FALSE if conf equal TRUE:
# The confidence interval is computed by a Fisher transformation.
# The significance level relies on a one-sided student-T distribution.
# If conf equal FALSE: no confidence interval is calculated
#
# Returns:
# $ACC: if conf equal TRUE:
# Matrix with c(nexp, nobs, nsdates, nleadtimes, 4) dimensions.
# The fourth dimension of length 4 corresponds to the lower limit of
# the 95% confidence interval, the computed ACC, the upper limit of
# the 95% confidence interval and the 95% significance level given by
# a one-sided T-test.
# if conf equal FALSE:
# ACC with c(nexp, nobs, nleadtimes) dimensions.
# $MACC: Mean Anomaly Correlation Coefficient with c(nexp, nobs, nleadtimes)
# dimensions.
#
# History:
# 1.0 # 2013-08 (V. Guemas, vguemas@ic3.cat) # Original code
# 1.1 # 2014-05 (C. Prodhomme, chloe.prodhomme@ic3.cat) # optimisation
Nicolau Manubens Gil
committed
dimsvar <- dim(var_exp)
stop("var_exp & var_obs should have dimensions (nexp/nsobs, nsdates,
nltimes, nlat, nlon) ")
Nicolau Manubens Gil
committed
}
for (iind in 2:length(dimsvar)) {
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[2]
nltimes <- dimsvar[3]
nlat <- dimsvar[4]
nlon <- dimsvar[5]
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
}
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
if(conf == TRUE) {
ACC <- array(NA, dim = c(nexp, nobs, nsdates, nltimes, 4))
} else {
ACC <- array(NA, dim = c(nexp, nobs, nsdates, nltimes))
}
MACCaux <- array(0, dim = c(nexp, nobs, nsdates,nltimes,3))
MACC <- array(0, dim = c(nexp, nobs,nltimes))
for( iobs in 1:nobs) {
for( iexp in 1:nexp) {
tmp <- abind(
array(var_exp[iexp, , , , ],
dim = c(1, nsdates, nltimes,
length(indlon) * length(indlat))),
array(var_obs[iobs, , , , ],
dim = c(1, nsdates, nltimes,
length(indlon) * length(indlat))), along = 1)
tmpallNA <- apply(tmp, c(1, 2, 3), function(x)
isTRUE(
all.equal(x,
array(NA,
length(indlon) * length(indlat)))))
#calc of correlation
#correlation parametric test
# #(to be compatible with the previous ACC function)
ACC[iexp, iobs, , , 2] <- apply(tmp, c(2, 3),
function(x) cor(x[1, ], x[2, ],
use = "na.or.complete"))
eno <- Eno(tmp[2, , , ], 3)
t <- apply( eno, c(1, 2),
function(x) qt(0.95, x - 2))
enot = abind(eno, t, along = 3)
ACC[iexp, iobs, , , 4] <- apply(enot, c(1, 2), function(x)
sqrt((x[2] * x[2]) /
((x[2] * x[2]) + x[1] - 2)))
correno = abind(array(ACC[iexp, iobs, , , 2],
dim = c(nsdates, nltimes)), eno, along = 3)
ACC[iexp, iobs, , , 1] <- apply(correno, c(1, 2),
function(x)
tanh(atanh(x[1]) +
qnorm(0.975) / sqrt(x[2] - 3)))
ACC[iexp, iobs, , , 3] <- apply(correno, c(1, 2), function(x)
tanh(atanh(x[1]) +
qnorm(0.025) / sqrt(x[2] - 3)))
} else {
ACC[iexp, iobs, , ] <- apply(tmp, c(2, 3),
function(x)
cor(x[1, ], x[2, ], use = "na.or.complete"))
Nicolau Manubens Gil
committed
}
MACCaux[iexp, iobs, , , 1] <- apply(tmp, c(2, 3),
function(x)
sum(x[1, ] * x[2, ], na.rm = TRUE))
MACCaux[iexp, iobs, , , 2] <- apply(tmp[1, , , ], c(1, 2),
function(x)
sum(x^2, na.rm = TRUE))
MACCaux[iexp, iobs, , , 3] <- apply(tmp[2, , , ], c(1, 2),
function(x) sum(x^2, na.rm = TRUE))
MACCaux1 <- MACCaux[iexp, iobs, , , 1]
MACCaux2 <- MACCaux[iexp, iobs, , , 2]
MACCaux3 <- MACCaux[iexp, iobs, , , 3]
MACCaux1[tmpallNA[1, , ]] = NA
MACCaux2[tmpallNA[1, , ]] = NA
MACCaux3[tmpallNA[1, , ]] = NA
MACCaux1[tmpallNA[2, , ]] = NA
MACCaux2[tmpallNA[2, , ]] = NA
MACCaux3[tmpallNA[2, , ]] = NA
MACCaux[iexp, iobs, , , 1] = MACCaux1
MACCaux[iexp, iobs, , , 2] = MACCaux2
MACCaux[iexp, iobs, , , 3] = MACCaux3
Nicolau Manubens Gil
committed
}
}
# if one of the date if a not a number the result will be not a number
MACC<- apply(MACCaux, c(1, 2, 4), function(x) sum(x[, 1], na.rm = FALSE) /
sqrt(sum(x[, 2], na.rm = FALSE) * sum(x[, 3], na.rm = FALSE)))
Nicolau Manubens Gil
committed
invisible(list(ACC = ACC, MACC = MACC))