R/homogeneity.R

In eFCM: Exponential Factor Copula Model

#' Homogeneity Test
#'
#' Modified versions of HW, and AD tests
#'
#' @param x Numeric vector of exceedances.
#' @param cod Vector of region/station codes, same length as x.
#' @param Nsim Number of simulations (default: 500).
#'
#' @return Named numeric vector with standardized HW test statistics (H1, H2) and regional tau3.
#' @noRd
#' @importFrom nsRFA Lmoments par.kappa rand.kappa regionalLmoments
#' @srrstats {G1.4a} *Implements a novel combination of regional homogeneity statistics*
fcmHW.tests <- function(x, cod, Nsim = 500) {
  if (length(x) != length(cod)) {
    stop("x and cod must have the same length")
  }
  fac <- factor(cod)
  ni <- tapply(x, fac, length)
  k <- nlevels(fac)
  Lm <- sapply(split(x, fac), Lmoments)
  rLm <- regionalLmoments(x, fac)
  ti <- as.numeric(Lm[3, ])
  t3i <- as.numeric(Lm[4, ])
  t4i <- as.numeric(Lm[5, ])
  lambda1Reg <- as.numeric(rLm[1])
  lambda2Reg <- as.numeric(rLm[2])
  tauReg <- as.numeric(rLm[3])
  tau3Reg <- as.numeric(rLm[4])
  tau4Reg <- as.numeric(rLm[5])
  V1 <- (sum(ni * (ti - tauReg)^2) / sum(ni))^0.5
  V2 <- sum(ni * ((ti - tauReg)^2 + (t3i - tau3Reg)^2)^0.5) / sum(ni)
  parkappa <- par.kappa(lambda1Reg, lambda2Reg, tau3Reg, tau4Reg)
  xi <- parkappa$xi
  alfa <- parkappa$alfa
  kappa <- parkappa$k
  hacca <- parkappa$h
  V1s <- rep(NA, Nsim)
  V2s <- rep(NA, Nsim)
  for (i in 1:Nsim) {
    ti.sim <- rep(NA, k)
    t3i.sim <- rep(NA, k)
    for (j in 1:k) {
      campione <- rand.kappa(ni[j], xi, alfa, kappa, hacca)
      campione.ad <- campione / mean(campione)
      lmom <- Lmoments(campione.ad)
      ti.sim[j] <- lmom[3]
      t3i.sim[j] <- lmom[4]
    }
    tauReg.sim <- sum(ni * ti.sim) / sum(ni)
    tau3Reg.sim <- sum(ni * t3i.sim) / sum(ni)
    V1s[i] <- (sum(ni * (ti.sim - tauReg.sim)^2) / sum(ni))^0.5
    V2s[i] <- sum(ni * ((ti.sim - tauReg.sim)^2 + (t3i.sim - tau3Reg.sim)^2)^0.5) / sum(ni)
  }
  muV1 <- mean(V1s)
  stdV1 <- sd(V1s)
  muV2 <- mean(V2s)
  stdV2 <- sd(V2s)
  H1 <- (V1 - muV1) / stdV1
  H2 <- (V2 - muV2) / stdV2
  output <- c(H1, H2, tau3Reg.sim)
  names(output) <- c("H1", "H2", "t3R")
  return(output)
}
#' @noRd
#' @importFrom nsRFA ksampleA2 nonparboot
fcmADbootstrap.test <- function(x, cod, Nsim = 500, index = 2, alpha = 0.05) {
  if (length(x) != length(cod)) {
    stop("x and cod must have the same length")
  }
  fac <- factor(cod)
  ni <- tapply(x, fac, length)
  k <- nlevels(fac)
  N <- sum(ni)
  if (index == 1) {
    indexflood <- function(x) {
      m <- mean(x)
      return(m)
    }
  } else if (index == 2) {
    indexflood <- function(x) {
      m <- median(x)
      return(m)
    }
  }
  med <- tapply(x, fac, indexflood)
  regione.adim <- x / unsplit(med, fac)
  A2kN <- ksampleA2(regione.adim, fac)
  A2kNs <- rep(NA, Nsim)
  for (i in 1:Nsim) {
    regione.simul <- nonparboot(regione.adim)
    med.simul <- tapply(regione.simul, fac, indexflood)
    regione.simul.adim <- regione.simul / unsplit(med.simul, fac)
    A2kNs[i] <- ksampleA2(regione.simul.adim, fac)
  }
  ecdfA2kNs <- ecdf(A2kNs)
  probabilita <- ecdfA2kNs(A2kN)
  output <- c(A2kN, probabilita)
  names(output) <- c("A2kN", "P")
  return(output)
}
#'
#' @noRd
#' @importFrom stats quantile
assess.hom <- function(data, cols, pr = 0.9, alpha = 0.05, rm.zeros = FALSE, DK.test = FALSE, which.test = c(1, 2)) {
  N.s0 <- data[, cols]
  if (rm.zeros) {
    N.s0[N.s0 == 0] <- NA
  }
  z <- NULL
  for (j in 1:length(cols)) {
    y <- N.s0[, j]
    y <- y[!is.na(y)]
    if (length(y) > 0) {
      thres <- quantile(y, pr)
      y <- y[y > thres]
      temp <- cbind(y, j)
      z <- rbind(z, temp)
    }
  }
  x <- z[, 1]
  cod <- z[, 2]
  if (all(which.test == c(1, 2))) {
    test <- tryCatch(fcmHW.tests(x, cod), error = function(e) e)
    if (!inherits(test, "error")) {
      if (test[3] <= 0.23) {
        value <- test[1]
        if (value < 1) {
          code <- 1
        } else {
          code <- 0
        } # 'acceptably homogeneous' if HW < 1
      } else {
        test2 <- fcmADbootstrap.test(x, cod, Nsim = 500, alpha = alpha)
        if (test2[2] < (1 - alpha)) {
          code <- 1
        } else {
          code <- 0
        }
      }
    } else {
      test2 <- fcmADbootstrap.test(x, cod, Nsim = 500, alpha = alpha)
      if (test2[2] < (1 - alpha)) {
        code <- 1
      } else {
        code <- 0
      }
    }
  } else {
    if (which.test == 1) {
      test <- fcmHW.tests(x, cod)
      value <- test[1]
      if (value < 1) {
        code <- 1
      } else {
        code <- 0
      }
    }
    if (which.test == 2) {
      test2 <- fcmADbootstrap.test(x, cod, Nsim = 500, alpha = alpha)
      if (test2[2] < (1 - alpha)) {
        code <- 1
      } else {
        code <- 0
      }
    }
  }
  list(code = code, idcols = unique(cod))
}