Composite.R

Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) {

  if ( dim(var)[3] != length(occ) ) {
     stop("temporal dimension of var is not equal to length of occ.")
  }
  K         <- max(occ)
  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 (eno == TRUE) { 
    n_tot <- Eno(var, posdim = 3)
  } else {
    n_tot <- length(occ)
  }
  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(occ == k) + lag
    toberemoved = which(0 > indices | indices > dim(var)[3])

    if (length(toberemoved) > 0) {
        indices=indices[-toberemoved]
    }
    if (eno == TRUE) {
        n_k <- Eno(var[,, indices], posdim = 3)
    } else {
        n_k <- length(indices)
    }
    if (length(indices) == 1) {
        composite[,, k] <- var[,, indices] 
        warning(paste("Composite", k, "has length 1 and pvalue is NA."))
    }  else {
        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(fileout) == FALSE ) { 
    output <- list(composite = composite, pvalue = pvalue)   
    save(output, file = paste(fileout, '.sav', sep = ''))
  }
  
  invisible(list(composite = composite, pvalue = pvalue)) 
}