Newer
Older
1
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
40
41
42
43
44
45
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
Composite <- function(var, occurrences, lag=0, tflag=FALSE, nameout=NULL) {
#
# Composites a 3-d field var(x,y,time) according to the indices of
# mode/cluster occurrences in time and computes the pvalues (t-test).
#
# Args :
# var : 3-dimensional tensor (x, y, time) containing the variable to composite.
# occurrences : 1-dimensional array for the occurrence time series of mode(s)/cluster(s)
# (*1) when one wants to composite all modes, e.g., all K=3 clusters then
# for example occurrences could look like: 1 1 2 3 2 3 1 3 3 2 3 2 2 3 2,
# (*2) otherwise for compositing only the 2nd mode or cluster of the above
# example occurrences should look like 0 0 1 0 1 0 0 0 0 1 0 1 1 0 1.
# lag : Lag time step (an integer), e.g., for lag=2 composite will use +2
# occurrences (i.e., shifted 2 time steps forward). Deafult is lag=0.
# tflag : For using the effective sample size (TRUE) or the total sample size
# (FALSE that is the default) for the number of degrees of freedom.
# nameout : Name of the .sav output file (NULL is the default).
#
# Return :
# $Composite : 3-d tensor (x, y, k) containing the composites k=1,..,K for case (*1)
# or only k=1 for any specific cluster, i.e., case (*2).
# $Pvalue : 3-d tensor (x, y, k) containing the pvalue of the
# composites obtained through a t-test that accounts for the serial
# dependence of the data with the same structure as Composite.
#
# History:
# 0.1 # 2014-08 (N.S. Fuckar, neven.fuckar@ic3.cat) # Original code
#
if ( dim(var)[3]!=length(occurrences) ) { stop("temporal dimension of var is not equal to length of occurrences") }
K <- max(occurrences)
composite <- array(dim = c(dim(var)[1:2], K))
tvalue <- array(dim = dim(var)[1:2])
dof <- array(dim = dim(var)[1:2])
pvalue <- array(dim = c(dim(var)[1:2], K))
if ( tflag==TRUE ) { n_tot <- Eno(var, posdim=3) }
else { n_tot <- length(occurrences) }
mean_tot <- Mean1Dim(var, posdim=3, narm=TRUE)
stdv_tot <- apply(var, c(1,2), sd, na.rm=TRUE)
for (k in 1:K) {
indices <- which(occurrences==k)+lag
toberemoved=which(0>indices|indices>dim(var)[3])
if ( length(toberemoved) > 0 ) { indices=indices[-toberemoved] }
if ( tflag==TRUE ) { n_k <- Eno(var[,,indices], posdim=3) }
else { n_k <- length(indices) }
composite[,,k] <- Mean1Dim(var[,,indices], posdim=3, narm=TRUE)
stdv_k <- apply(var[,,indices], c(1,2), sd, na.rm=TRUE)
tvalue[,] <- (mean_tot - composite[,,k])/sqrt(stdv_tot^2/n_tot + stdv_k^2/n_k)
dof[,] <- (stdv_tot^2/n_tot + stdv_k^2/n_k)^2/((stdv_tot^2/n_tot)^2/(n_tot - 1) + (stdv_k^2/n_k)^2/(n_k - 1))
pvalue[,,k] <- 2*pt(-abs(tvalue[,]), df=dof[,])
}
if ( is.null(nameout)==FALSE ) {
output <- list(Composite=composite, Pvalue=pvalue)
save(output,file=paste(nameout,'.sav',sep=''))
}
#
# Outputs
# ~~~~~~~~~
#
invisible(list(Composite = composite, Pvalue = pvalue))
}