R/Composite.R

@@ -17,6 +17,7 @@
 #'  use +2 occurrences (i.e., shifted 2 time steps forward). Default is lag = 0.
 #'@param eno For using the effective sample size (TRUE) or the total sample 
 #'  size (FALSE that is the default) for the number of degrees of freedom.
+#'@param K numeric value indicating the maximum number of composites. By default (NULL), it takes the maximum value provided in occ.
 #'@param fileout Name of the .sav output file (NULL is the default).
 #'
 #'@return 
@@ -67,18 +68,26 @@
 #'occ2[c(3, 9, 15, 21)] <- 1
 #'
 #'filled.contour(Composite(var=f1, occ=occ2)$composite[,,1])
+#'
+#'Example with one missing composite (#3) in occ:
+#'data <- 1 : (4 * 5 * 6)
+#'dim(data) <- c(lon = 4, lat = 5, case = 6)
+#'occ <- c(1, 1, 2, 2, 3, 3) 
+#'res <- Composite(data, occ, K = 4)
 #'@importFrom stats sd pt
 #'@export
-Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) {
+Composite <- function(var, occ, lag = 0, eno = FALSE, K = NULL, 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))
+  if (is.null(K)) {
+     K <- max(occ)
+  }
+  composite <- array(dim = c(dim(var)[1:2], composite = K))
   tvalue    <- array(dim = dim(var)[1:2])
   dof       <- array(dim = dim(var)[1:2])
-  pvalue    <- array(dim = c(dim(var)[1:2], K))
+  pvalue    <- array(dim = c(dim(var)[1:2], composite = K))
 
   if (eno == TRUE) { 
     n_tot <- Eno(var, posdim = 3)
@@ -89,34 +98,34 @@ Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) {
   stdv_tot <- apply(var, c(1, 2), sd, na.rm = TRUE) 
 
   for (k in 1 : K) {
+    if (length(which(occ == k)) >= 1) {
+      indices <- which(occ == k) + lag
+      toberemoved <-  which(0 > indices | indices > dim(var)[3])
 
-    indices <- which(occ == k) + lag
-    toberemoved <-  which(0 > indices | indices > dim(var)[3])
-
-    if (length(toberemoved) > 0) {
+      if (length(toberemoved) > 0) {
         indices <- indices[-toberemoved]
-    }
-    if (eno == TRUE) {
+      }
+      if (eno == TRUE) {
         n_k <- Eno(var[, , indices], posdim = 3)
-    } else {
+      } else {
         n_k <- length(indices)
-    }
-    if (length(indices) == 1) {
+      }
+      if (length(indices) == 1) {
         composite[, , k] <- var[, , indices] 
         warning(paste("Composite", k, "has length 1 and pvalue is NA."))
-    }  else {
+      }  else {
         composite[, , k] <- Mean1Dim(var[, , indices], posdim = 3, narm = TRUE)
-    }
-    stdv_k <- apply(var[, , indices], c(1, 2), sd, na.rm = 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)
-  }  
-
+      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 = ''))

tests/testthat/test-Composite.R

@@ -11,7 +11,7 @@ test_that("Sanity checks", {
 
   var <- array(rep(c(1, 3, 2, 1, 2), 8), dim = c(x = 2, y = 4, time = 5))
   occ <- c(1, 2, 2, 2, 1)
-  output <- c(x = 2, y = 4, 2)  #dim(asd$composite)
+  output <- c(x = 2, y = 4, composite = 2)  #dim(asd$composite)
   expect_equal(
     dim(Composite(var, occ)$composite),
     output
@@ -29,6 +29,20 @@ test_that("Sanity checks", {
     Composite(var, occ)$composite[, , 2],
     output
   )
+  expect_equal(
+    dim(Composite(var, occ, K = 2)$composite),
+    c(x = 2, y = 4, composite = 2)
+  )
+  expect_equal(
+    Composite(var, occ, K = 3),
+    Composite(var, occ)
+  )
+  expect_equal(
+    Composite(var, occ, K = 4)$composite[1, 1, ],
+    c(1.5, 1.5, 2.0, NA),
+    tolerance = 0.01
+  )
+
 
 })