WSM stress data A to E; bug fix in weighted_rayleigh_test(); cran sub…

…mission fixes
tobiste · May 15, 2024 · a0517f5 · a0517f5
1 parent 662045f
commit a0517f5
Show file tree

Hide file tree

Showing 19 changed files with 107 additions and 133 deletions.
diff --git a/R/interpolation.R b/R/interpolation.R
@@ -592,7 +592,7 @@ kernel_dispersion <- function(x,
           if (stat == "nchisq") {
             y <- norm_chisq(datas$azi[ids_R], prd = datas$prd[ids_R], datas$unc[ids_R])
           } else if (stat == "rayleigh") {
-            y <- weighted_rayleigh(datas$azi[ids_R], prd = datas$prd[ids_R], unc = datas$unc[ids_R], ...)$statistic
+            y <- weighted_rayleigh(datas$azi[ids_R], mu = datas$prd[ids_R], w = 1 / datas$unc[ids_R], ...)$statistic
           } else {
             y <- circular_dispersion(datas$azi[ids_R], y = datas$prd[ids_R], w = 1 / datas$unc[ids_R], ...)
           }

diff --git a/R/model_shmax.R b/R/model_shmax.R
@@ -118,7 +118,8 @@ deviation_norm <- function(x, y = NULL) {
     ny <- rep(y, nx)
   }
   d <- (x %% 180) - (y %% 180)
-  ifelse(d < 90, d, 180 - d)
+  d <- ifelse(d < 90, d, 180 - d)
+  abs(d)
 }
 
 

diff --git a/R/plotting.R b/R/plotting.R
@@ -741,7 +741,7 @@ quick_plot <- function(
 
   nchisq <- norm_chisq(azi, prd, unc)
   suppressMessages(
-    rt <- weighted_rayleigh(azi, prd = prd, unc = unc)
+    rt <- weighted_rayleigh(azi, mu = prd, w = 1 / unc)
   )
   azi.PoR.mean <- circular_mean(azi, 1 / unc)
   azi.PoR.sd <- circular_sd(azi, 1 / unc)

diff --git a/R/roll_statistics.R b/R/roll_statistics.R
@@ -17,7 +17,7 @@
 #' range. If `TRUE` (default), then the subset of indexes that are in range
 #' are passed to `FUN`. A numeric argument to partial can be used to determine
 #' the minimal window size for partial computations. See below for more details.
-#' @inheritDotParams zoo::rollapply -data
+#' @param ... optional arguments passed to [zoo::rollapply()]
 #'
 #' @returns numeric vector  with the results of the rolling function.
 #'
@@ -134,9 +134,7 @@ roll_normchisq <- function(obs, prd, unc = NULL,
                            partial = TRUE,
                            fill = NA,
                            ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(obs)
-  }
+  if (is.null(width)) width <- optimal_rollwidth(obs)
 
   zoo::rollapply(
     cbind(obs, prd, unc),
@@ -158,12 +156,10 @@ roll_rayleigh <- function(obs, prd, unc = NULL,
                           partial = TRUE,
                           fill = NA,
                           ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(obs)
-  }
+  if (is.null(width)) width <- optimal_rollwidth(obs)
 
   zoo::rollapply(
-    cbind(obs, prd, unc),
+    cbind(obs, prd, 1 / unc),
     width = width,
     FUN = function(x) {
       suppressMessages(weighted_rayleigh(x[, 1], x[, 2], x[, 3]))$statistic
@@ -182,23 +178,17 @@ roll_dispersion <- function(x, y, w = NULL, w.y = NULL,
                             partial = TRUE,
                             fill = NA,
                             ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(x)
-  }
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-  if (is.null(w.y)) {
-    w.y <- rep(1, length(x))
-  }
+  if (is.null(width)) width <- optimal_rollwidth(x)
+  if (is.null(w)) w <- rep(1, length(x))
+  if (is.null(w.y)) w.y <- rep(1, length(x))
 
   zoo::rollapply(
     cbind(x, y, w, w.y),
     width = width,
     FUN = function(x) {
       suppressMessages(
         circular_dispersion(x[, 1], x[, 2], x[, 3], x[, 4], norm = TRUE)
-        )
+      )
     },
     by.column = by.column,
     partial = partial,
@@ -213,9 +203,7 @@ roll_confidence <- function(x, conf.level = .95, w = NULL, axial = TRUE,
                             width = NULL, by.column = FALSE, partial = TRUE,
                             fill = NA,
                             ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(x)
-  }
+  if (is.null(width)) width <- optimal_rollwidth(x)
 
   zoo::rollapply(
     cbind(x, rep(conf.level, length(x)), w),
@@ -234,15 +222,9 @@ roll_confidence <- function(x, conf.level = .95, w = NULL, axial = TRUE,
 #' @export
 roll_dispersion_CI <- function(x, y, w = NULL, w.y = NULL, R, conf.level = .95,
                                width = NULL, by.column = FALSE, partial = TRUE, fill = NA, ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(x)
-  }
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-  if (is.null(w.y)) {
-    w.y <- rep(1, length(x))
-  }
+  if (is.null(width)) width <- optimal_rollwidth(x)
+  if (is.null(w)) w <- rep(1, length(x))
+  if (is.null(w.y)) w.y <- rep(1, length(x))
 
   zoo::rollapply(
     cbind(x, y, w, w.y),
@@ -261,15 +243,9 @@ roll_dispersion_CI <- function(x, y, w = NULL, w.y = NULL, R, conf.level = .95,
 #' @export
 roll_dispersion_sde <- function(x, y, w = NULL, w.y = NULL, R, conf.level = .95,
                                 width = NULL, by.column = FALSE, partial = TRUE, fill = NA, ...) {
-  if (is.null(width)) {
-    width <- optimal_rollwidth(x)
-  }
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-  if (is.null(w.y)) {
-    w.y <- rep(1, length(x))
-  }
+  if (is.null(width)) width <- optimal_rollwidth(x)
+  if (is.null(w)) w <- rep(1, length(x))
+  if (is.null(w.y)) w.y <- rep(1, length(x))
 
   zoo::rollapply(
     cbind(x, y, w, w.y),
@@ -278,8 +254,9 @@ roll_dispersion_sde <- function(x, y, w = NULL, w.y = NULL, R, conf.level = .95,
       suppressMessages(
         circular_dispersion_boot(
           x[, 1], x[, 2], x[, 3], x[, 4],
-          R = R, conf.level = conf.level, ...)$sde
-        )
+          R = R, conf.level = conf.level, ...
+        )$sde
+      )
     },
     by.column = by.column,
     partial = partial,
@@ -351,9 +328,7 @@ distroll_circstats <- function(x, distance, FUN, width = NULL, min_n = 2,
                                sort = TRUE, ...) {
   align <- match.arg(align)
   stopifnot(length(x) == length(distance))
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
+  if (is.null(w)) w <- rep(1, length(x))
 
   if (is.null(width)) {
     width <- seq(from = min(distance, na.rm = TRUE), to = max(distance, na.rm = TRUE), length.out = 10) |>
@@ -394,10 +369,7 @@ distroll_confidence <- function(x, distance, w = NULL, width = NULL, min_n = 2,
                                 sort = TRUE, ...) {
   align <- match.arg(align)
   stopifnot(length(x) == length(distance))
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-
+  if (is.null(w)) w <- rep(1, length(x))
   if (is.null(width)) {
     width <- seq(min(distance, na.rm = TRUE), max(distance, na.rm = TRUE), length.out = 10) |>
       diff() |>
@@ -438,15 +410,9 @@ distroll_dispersion <- function(x, y, w = NULL, w.y = NULL, distance,
                                 sort = TRUE, ...) {
   align <- match.arg(align)
   stopifnot(length(x) == length(distance))
-  if (length(y) == 1) {
-    y <- rep(y, length(x))
-  }
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-  if (is.null(w.y)) {
-    w.y <- rep(1, length(x))
-  }
+  if (length(y) == 1) y <- rep(y, length(x))
+  if (is.null(w)) w <- rep(1, length(x))
+  if (is.null(w.y)) w.y <- rep(1, length(x))
 
   if (is.null(width)) {
     width <- seq(min(distance, na.rm = TRUE), max(distance, na.rm = TRUE), length.out = 10) |>
@@ -488,15 +454,9 @@ distroll_dispersion_sde <- function(x, y, w = NULL, w.y = NULL, distance,
                                     sort = TRUE, ...) {
   align <- match.arg(align)
   stopifnot(length(x) == length(distance))
-  if (length(y) == 1) {
-    y <- rep(y, length(x))
-  }
-  if (is.null(w)) {
-    w <- rep(1, length(x))
-  }
-  if (is.null(w.y)) {
-    w.y <- rep(1, length(x))
-  }
+  if (length(y) == 1) y <- rep(y, length(x))
+  if (is.null(w)) w <- rep(1, length(x))
+  if (is.null(w.y)) w.y <- rep(1, length(x))
 
   if (is.null(width)) {
     width <- seq(min(distance, na.rm = TRUE), max(distance, na.rm = TRUE), length.out = 10) |>

diff --git a/R/stat_tests.R b/R/stat_tests.R
@@ -170,7 +170,6 @@ rayleigh_test <- function(x, mu = NULL, axial = TRUE) {
     1
   }
 
-
   if (is.null(mu)) {
     x <- (na.omit(x) * f) %% 360
     n <- length(x)
@@ -258,9 +257,9 @@ rayleigh_p_value2 <- function(K, n) {
 #' Weighted version of the Rayleigh test (or V0-test) for uniformity against a
 #' distribution with a priori expected von Mises concentration.
 #' @param x numeric vector. Values in degrees
-#' @param unc numeric. The standard deviations of `x`. If `NULL`, the non-weighted
-#' Rayleigh test is performed.
-#' @param prd The a priori expected direction (in degrees) for the alternative
+#' @param w numeric vector weights of length `length(x)`. If `NULL`, the
+#' non-weighted Rayleigh test is performed.
+#' @param mu The *a priori* expected direction (in degrees) for the alternative
 #' hypothesis.
 #' @param axial logical. Whether the data are axial, i.e. \eqn{\pi}-periodical
 #' (`TRUE`, the default) or directional, i.e. \eqn{2 \pi}-periodical (`FALSE`).
@@ -295,40 +294,36 @@ rayleigh_p_value2 <- function(K, n) {
 #' tibet.por <- PoR_shmax(tibet, PoR.tib, "in")
 #'
 #' # GOF test:
-#' weighted_rayleigh(tibet.por$azi.PoR, prd = 90, unc = tibet$unc)
-#' weighted_rayleigh(ice.por$azi.PoR, prd = 0, unc = iceland$unc)
-#' weighted_rayleigh(sa.por$azi.PoR, prd = 135, unc = san_andreas$unc)
-weighted_rayleigh <- function(x, prd = NULL, unc, axial = TRUE) {
-  if (is.null(unc)) {
-    rayleigh_test(x, mu = prd, axial = axial)
+#' weighted_rayleigh(tibet.por$azi.PoR, mu = 90, w = 1 / tibet$unc)
+#' weighted_rayleigh(ice.por$azi.PoR, mu = 0, w = 1 / iceland$unc)
+#' weighted_rayleigh(sa.por$azi.PoR, mu = 135, w = 1 / san_andreas$unc)
+weighted_rayleigh <- function(x, mu = NULL, w = NULL, axial = TRUE) {
+  if (is.null(w)) {
+    rayleigh_test(x, mu = mu, axial = axial)
   } else {
-    data <- cbind(x = x, unc = unc)
+    data <- cbind(x = x, w = w)
     data <- data[stats::complete.cases(data), ] # remove NA values
 
 
-    unc <- data[, "unc"]
-    w <- 1 / unc
+    w <- data[, "w"]
     Z <- sum(w)
-    n <- length(unc)
+    n <- length(w)
 
-    if (is.null(prd)) {
-      prd <- circular_mean(x, w, axial, na.rm = FALSE)
+    if (is.null(mu)) {
+      mu <- circular_mean(x, w, axial, na.rm = FALSE)
     }
 
-    d <- data[, "x"] - prd
+    d <- deviation_norm(data[, "x"], mu)
 
-    f <- 1
-    if (axial) {
-      f <- 2
-    }
+    f <- ifelse(axial, 2, 1)
     cosd <- cosd(f * d)
     wcosd <- w * cosd
 
-    md <- 1
+    md <- 2
     # if(norm){
     #   md <- 2
     # }
-    wmd <- md * w # = w * (1 - cos(pi)) = w* (1 - (-1))
+    wmd <- md * w # = w * (1 - cos(pi)) = w * (1 - (-1))
 
     C <- (sum(wcosd) / sum(wmd))
 
@@ -389,11 +384,8 @@ kuiper_test <- function(x, alpha = 0, axial = TRUE) {
     c(0.15, 0.1, 0.05, 0.025, 0.01),
     c(1.537, 1.62, 1.747, 1.862, 2.001)
   )
-  f <- if (axial) {
-    2
-  } else {
-    1
-  }
+  f <- ifelse(axial, 2, 1)
+
   x <- (na.omit(x) * f) %% 360
   u <- sort(deg2rad(x) %% (2 * pi)) / (2 * pi)
   n <- length(x)
@@ -440,11 +432,14 @@ kuiper_test <- function(x, alpha = 0, axial = TRUE) {
 #' Watson's test statistic is a rotation-invariant Cramer - von Mises test
 #'
 #' @param x numeric vector. Values in degrees
-#' @param alpha Significance level of the test. Valid levels are `0.01`, `0.05`, and `0.1`.
-#' This argument may be omitted (`NULL`, the default), in which case, a range for the p-value will be returned.
+#' @param alpha Significance level of the test. Valid levels are `0.01`, `0.05`,
+#' and `0.1`.
+#' This argument may be omitted (`NULL`, the default), in which case, a range
+#' for the p-value will be returned.
 #' @param axial logical. Whether the data are axial, i.e. \eqn{\pi}-periodical
 #' (`TRUE`, the default) or circular, i.e. \eqn{2 \pi}-periodical (`FALSE`).
-#' @param mu (optional) The specified or known mean direction (in degrees) in alternative hypothesis
+#' @param mu (optional) The specified mean direction (in degrees) in alternative
+#'  hypothesis
 #' @param dist Distribution to test for. The default, `"uniform"`, is the
 #' uniform distribution. `"vonmises"` tests the von Mises distribution.
 #'
@@ -455,7 +450,8 @@ kuiper_test <- function(x, alpha = 0, axial = TRUE) {
 #' If `statistic > p.value`, the null hypothesis is rejected.
 #' If not, randomness (uniform distribution) cannot be excluded.
 #'
-#' @references Mardia and Jupp (2000). Directional Statistics. John Wiley and Sons.
+#' @references Mardia and Jupp (2000). Directional Statistics. John Wiley and
+#' Sons.
 #'
 #' @export
 #'
@@ -720,7 +716,7 @@ A1inv <- function(x) {
 #' @param bias logical parameter determining whether a bias correction is used
 #' in the computation of the MLE. Default for bias is `FALSE` for no bias
 #' correction.
-#' @param ... optional paramters passed to `circular_mean()`
+#' @param ... optional parameters passed to `circular_mean()`
 #'
 #' @returns numeric.
 #' @export

diff --git a/R/various.R b/R/various.R
@@ -103,7 +103,7 @@ stress_analysis <- function(x, PoR, type = c("none", "in", "out", "right", "left
   disp <- circular_dispersion(res$azi.PoR, prd, 1 / x$unc)
   conf <- confidence_angle(res$azi.PoR, w = 1 / x$unc)
   nchisq <- norm_chisq(res$azi.PoR, prd, unc = x$unc)
-  rayleigh <- weighted_rayleigh(res$azi.PoR, prd, unc = x$unc)
+  rayleigh <- weighted_rayleigh(res$azi.PoR, prd, w = 1 / x$unc)
 
   if (plot) {
     PoR_map(x, PoR, pb, type = type, deviation = TRUE)
@@ -157,7 +157,6 @@ line_azimuth <- function(x) {
     sf::st_coordinates()
 
   n <- nrow(mat)
-
   a <- numeric()
   for (i in 1:(n - 1)) {
     a[i] <- get_azimuth(mat[i, 2], mat[i, 1], mat[i + 1, 2], mat[i + 1, 1])