diff --git a/DESCRIPTION b/DESCRIPTION
index 45e7ea5fb..2b0918124 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: NNS
 Type: Package
 Title: Nonlinear Nonparametric Statistics
 Version: 10.9.2
-Date: 2024-08-29
+Date: 2024-09-01
 Authors@R: c(
   person("Fred", "Viole", role=c("aut","cre"), email="ovvo.financial.systems@gmail.com"),
   person("Roberto", "Spadim", role=c("ctb"))
diff --git a/NNS_10.9.2.tar.gz b/NNS_10.9.2.tar.gz
index 5a004121f..b0fe3a1fe 100644
Binary files a/NNS_10.9.2.tar.gz and b/NNS_10.9.2.tar.gz differ
diff --git a/NNS_10.9.2.zip b/NNS_10.9.2.zip
index 7564de1e9..c9fc9f810 100644
Binary files a/NNS_10.9.2.zip and b/NNS_10.9.2.zip differ
diff --git a/R/ANOVA.R b/R/ANOVA.R
index b5dbd3814..0a0fdb394 100644
--- a/R/ANOVA.R
+++ b/R/ANOVA.R
@@ -5,6 +5,7 @@
 #' @param control a numeric vector, matrix or data frame, or list if unequal vector lengths.
 #' @param treatment \code{NULL} (default) a numeric vector, matrix or data frame.
 #' @param means.only logical; \code{FALSE} (default) test whether difference in sample means only is zero.
+#' @param medians logical; \code{FALSE} (default) test whether difference in sample medians only is zero.  Requires \code{means.only = TRUE}. 
 #' @param confidence.interval numeric [0, 1]; The confidence interval surrounding the \code{control} mean, defaults to \code{(confidence.interval = 0.95)}.
 #' @param tails options: ("Left", "Right", "Both").  \code{tails = "Both"}(Default) Selects the tail of the distribution to determine effect size.
 #' @param pairwise logical; \code{FALSE} (default) Returns pairwise certainty tests when set to \code{pairwise = TRUE}.
@@ -37,6 +38,9 @@
 #' ### Two variable analysis with no control variable
 #' A <- cbind(x, y)
 #' NNS.ANOVA(A)
+#' 
+#' ### Medians test
+#' NNS.ANOVA(A, means.only = TRUE, medians = TRUE)
 #'
 #' ### Multiple variable analysis with no control variable
 #' set.seed(123)
@@ -56,6 +60,7 @@ NNS.ANOVA <- function(
     control,
     treatment,
     means.only = FALSE,
+    medians = FALSE,
     confidence.interval = 0.95,
     tails = "Both",
     pairwise = FALSE,
@@ -82,7 +87,7 @@ NNS.ANOVA <- function(
       
       nns.certainties <- sapply(
         1:ncol(control_matrix), 
-        function(g) NNS.ANOVA.bin(control_matrix[,g], treatment_matrix[,g], means.only = means.only, plot = FALSE)$Certainty
+        function(g) NNS.ANOVA.bin(control_matrix[,g], treatment_matrix[,g], means.only = means.only, medians = medians, plot = FALSE)$Certainty
       )
       
       cer_lower_CI <- LPM.VaR(.025, 1, nns.certainties[-1])
@@ -105,6 +110,7 @@ NNS.ANOVA <- function(
             control, 
             treatment,
             means.only = means.only,
+            medians = medians,
             confidence.interval = confidence.interval, 
             plot = plot, 
             tails = tails, 
@@ -121,6 +127,7 @@ NNS.ANOVA <- function(
           control, 
           treatment, 
           means.only = means.only,
+          medians = medians,
           confidence.interval = confidence.interval, 
           plot = plot, 
           tails = tails
@@ -146,15 +153,21 @@ NNS.ANOVA <- function(
     A <- control
   }
   
-  mean.of.means <- mean(colMeans(A, na.rm = T))
+  if(medians) mean.of.means <- mean(apply(A, 2, function(i) median(i, na.rm = TRUE))) else mean.of.means <- mean(colMeans(A, na.rm = T))
  
   if(!pairwise){
     #Continuous CDF for each variable from Mean of Means
-    LPM_ratio        <- sapply(1:n, function(b) LPM.ratio(1, mean.of.means, na.omit(unlist(A[ , b]))))
+    if(medians){
+      LPM_ratio <- sapply(1:n, function(b) LPM.ratio(0, mean.of.means, na.omit(unlist(A[ , b]))))
+    } else {
+      LPM_ratio <- sapply(1:n, function(b) LPM.ratio(1, mean.of.means, na.omit(unlist(A[ , b]))))
+    }
+    
     lower.25.target  <- mean(sapply(1:n, function(i) LPM.VaR(.25,  1, na.omit(unlist(A[,i])))))
     upper.25.target  <- mean(sapply(1:n, function(i) UPM.VaR(.25,  1, na.omit(unlist(A[,i])))))
     lower.125.target <- mean(sapply(1:n, function(i) LPM.VaR(.125, 1, na.omit(unlist(A[,i])))))
     upper.125.target <- mean(sapply(1:n, function(i) UPM.VaR(.125, 1, na.omit(unlist(A[,i])))))
+
     raw.certainties <- list(n - 1)
     for(i in 1:(n - 1)){
       raw.certainties[[i]] <- sapply(
@@ -163,6 +176,7 @@ NNS.ANOVA <- function(
           na.omit(unlist(A[ , i])),
           na.omit(unlist(A[ , b])),
           means.only = means.only,
+          medians = medians,
           mean.of.means = mean.of.means,
           upper.25.target = upper.25.target,
           lower.25.target = lower.25.target,
@@ -188,7 +202,7 @@ NNS.ANOVA <- function(
       )
       #For ANOVA Visualization
       abline(v = mean.of.means, col = "red", lwd = 4)
-      mtext("Grand Mean", side = 3,col = "red", at = mean.of.means)
+      if(medians) mtext("Grand Median", side = 3,col = "red", at = mean.of.means) else mtext("Grand Mean", side = 3,col = "red", at = mean.of.means)
     }
     return(c("Certainty" = NNS.ANOVA.rho))
   }
@@ -197,7 +211,7 @@ NNS.ANOVA <- function(
   for(i in 1:(n - 1)){
     raw.certainties[[i]] <- sapply(
       (i + 1) : n, 
-      function(b) NNS.ANOVA.bin(na.omit(unlist(A[ , i])), na.omit(unlist(A[ , b])), means.only = means.only, plot = FALSE)$Certainty
+      function(b) NNS.ANOVA.bin(na.omit(unlist(A[ , i])), na.omit(unlist(A[ , b])), means.only = means.only, medians = medians, plot = FALSE)$Certainty
     )
   }
   
@@ -217,7 +231,7 @@ NNS.ANOVA <- function(
       col = c('steelblue', rainbow(n - 1))
     )
     abline(v = mean.of.means, col = "red", lwd = 4)
-    mtext("Grand Mean", side = 3,col = "red", at = mean.of.means)
+    if(medians) mtext("Grand Median", side = 3,col = "red", at = mean.of.means) else mtext("Grand Mean", side = 3,col = "red", at = mean.of.means)
   }
   return(certainties)
 }
diff --git a/R/Binary_ANOVA.R b/R/Binary_ANOVA.R
index a951af44f..6b873b11e 100644
--- a/R/Binary_ANOVA.R
+++ b/R/Binary_ANOVA.R
@@ -1,5 +1,6 @@
 NNS.ANOVA.bin <- function(control, treatment,
                          means.only = FALSE,
+                         medians = FALSE,
                          mean.of.means = NULL,
                          upper.25.target = NULL,
                          lower.25.target = NULL,
@@ -8,7 +9,7 @@ NNS.ANOVA.bin <- function(control, treatment,
                          confidence.interval = NULL, tails = NULL, plot = TRUE, par = NULL){
 
   if(is.null(upper.25.target) && is.null(lower.25.target)){
-        mean.of.means <- mean(c(mean(control), mean(treatment)))
+        if(medians) mean.of.means <- mean(c(median(control), median(treatment))) else mean.of.means <- mean(c(mean(control), mean(treatment)))
         upper.25.target <- mean(c(UPM.VaR(.25, 1, control), UPM.VaR(.25, 1, treatment)))
         lower.25.target <- mean(c(LPM.VaR(.25, 1, control), LPM.VaR(.25, 1, treatment)))
         upper.125.target <- mean(c(UPM.VaR(.125, 1, control), UPM.VaR(.125, 1, treatment)))
@@ -18,9 +19,14 @@ NNS.ANOVA.bin <- function(control, treatment,
 
 
   #Continuous CDF for each variable from Mean of Means
-        LPM_ratio.1 <- LPM.ratio(1, mean.of.means, control)
-        LPM_ratio.2 <- LPM.ratio(1, mean.of.means, treatment)
-
+        if(medians){
+          LPM_ratio.1 <- LPM.ratio(0, mean.of.means, control)
+          LPM_ratio.2 <- LPM.ratio(0, mean.of.means, treatment)
+        } else {
+          LPM_ratio.1 <- LPM.ratio(1, mean.of.means, control)
+          LPM_ratio.2 <- LPM.ratio(1, mean.of.means, treatment)  
+        }
+        
         Upper_25_ratio.1 <- UPM.ratio(1, upper.25.target, control)
         Upper_25_ratio.2 <- UPM.ratio(1, upper.25.target, treatment)
         Upper_25_ratio <- mean(c(Upper_25_ratio.1, Upper_25_ratio.2))
@@ -59,13 +65,9 @@ NNS.ANOVA.bin <- function(control, treatment,
         }
         
     
-        
-        
-        
         pop.adjustment <- ((length(control) + length(treatment) - 2) / (length(control)  + length(treatment) )) ^ 2
 
   #Graphs
-
         if(plot){
             if(is.null(par)) original.par <- par(no.readonly = TRUE) else original.par <- par
 
@@ -73,17 +75,26 @@ NNS.ANOVA.bin <- function(control, treatment,
 
             #For ANOVA Visualization
             abline(v = mean.of.means, col = "red", lwd = 4)
-            mtext("Grand Mean", side = 3, col = "red", at = mean.of.means)
+            if(medians) mtext("Grand Median", side = 3, col = "red", at = mean.of.means) else mtext("Grand Mean", side = 3, col = "red", at = mean.of.means)
         }
 
     if(is.null(confidence.interval)){
 
+      if(medians){
+        return(list("Control Median" = median(control),
+                    "Treatment Median" = median(treatment),
+                    "Grand Median" = mean.of.means,
+                    "Control CDF" = LPM_ratio.1,
+                    "Treatment CDF" = LPM_ratio.2,
+                    "Certainty" = min(1, NNS.ANOVA.rho * pop.adjustment)))
+      } else {
           return(list("Control Mean" = mean(control),
                 "Treatment Mean" = mean(treatment),
                 "Grand Mean" = mean.of.means,
                 "Control CDF" = LPM_ratio.1,
                 "Treatment CDF" = LPM_ratio.2,
                 "Certainty" = min(1, NNS.ANOVA.rho * pop.adjustment)))
+      }
     } else {
 
         #Upper end of CDF confidence interval for control mean
@@ -93,15 +104,18 @@ NNS.ANOVA.bin <- function(control, treatment,
         
             # Resample control means
             y_p <- replicate(1000, sample.int(length(control), replace = TRUE))
-            control_means <- Rfast::colmeans(matrix(control[y_p], ncol = ncol(y_p), byrow = T))
-
+            if(medians){
+              control_means <- apply(matrix(control[y_p], ncol = ncol(y_p), byrow = T), 1, function(i) median(i, na.rm = T))
+            } else {
+              control_means <- apply(matrix(control[y_p], ncol = ncol(y_p), byrow = T), 1, function(i) mean(i, na.rm = T))
+            }
             a <- UPM.VaR(CI, 0, control_means)
             b <- mean(control_means)
             
             if(plot){
                 if(tails == "both" | tails == "right"){
                     abline(v = max(a, b), col = "green", lwd = 2, lty = 3)
-                    text(max(a, b), pos = 4, 0.5, paste0("<--- ", "ctl mu+ ", CI * 100, "%" ), col = "green")
+                    text(max(a, b), pos = 4, 0.5, paste0("<--- ", ifelse(medians, "ctl med+ " , "ctl mu+ "), CI * 100, "%" ), col = "green")
                 }
             }
 
@@ -112,39 +126,48 @@ NNS.ANOVA.bin <- function(control, treatment,
             if(plot){
                 if(tails == "both" | tails == "left"){
                     abline(v = min(c, d), col = "blue", lwd = 2, lty = 3)
-                    text(min(c, d), pos = 2, 0.5, paste0("ctl mu- ", paste0(CI * 100, "% --->")) , col = "blue")
+                    text(min(c, d), pos = 2, 0.5, paste0(ifelse(medians, "ctl med- ", "ctl mu- "), paste0(CI * 100, "% --->")) , col = "blue")
                 }
 
                 par(original.par)
             }
 
             #Effect Size Lower Bound
-            if(tails == "both") Lower.Bound.Effect <- mean(treatment) - max(a, b)
-            
-            if(tails == "left") Lower.Bound.Effect <- mean(treatment) - max(c, d)
-            
-            if(tails == "right") Lower.Bound.Effect <- mean(treatment) - max(a, b)
+            if(tails == "both") if(medians) Lower.Bound.Effect <- median(treatment) - max(a, b) else Lower.Bound.Effect <- mean(treatment) - max(a, b)
+            if(tails == "left") if(medians) Lower.Bound.Effect <- median(treatment) - max(c, d) else Lower.Bound.Effect <- mean(treatment) - max(c, d)
+            if(tails == "right") if(medians) Lower.Bound.Effect <- mean(treatment) - max(a, b) else Lower.Bound.Effect <- mean(treatment) - max(a, b)
             
 
 
             #Effect Size Upper Bound
-            if(tails == "both") Upper.Bound.Effect <- mean(treatment) - min(c, d)
-            
-            if(tails == "left") Upper.Bound.Effect <- mean(treatment) - min(c, d)
-            
-            if(tails == "right") Upper.Bound.Effect <- mean(treatment) - min(a, b)
+            if(tails == "both") if(medians) Upper.Bound.Effect <- median(treatment) - min(c, d) else Upper.Bound.Effect <- mean(treatment) - min(c, d)
+            if(tails == "left") if(medians) Upper.Bound.Effect <- median(treatment) - min(c, d) else Upper.Bound.Effect <- mean(treatment) - min(c, d)
+            if(tails == "right") if(medians) Upper.Bound.Effect <- median(treatment) - min(a, b) else Upper.Bound.Effect <- mean(treatment) - min(a, b)
             
 
 
 
         #Certainty Statistic and Effect Size Given Confidence Interval
-        return(list("Control Mean" = mean(control),
-                    "Treatment Mean" = mean(treatment),
-                    "Grand Mean" = mean.of.means,
-                    "Control CDF" = LPM_ratio.1,
-                    "Treatment CDF" = LPM_ratio.2,
-                    "Certainty" = min(1, NNS.ANOVA.rho * pop.adjustment),
-                    "Lower Bound Effect" = Lower.Bound.Effect,
-                    "Upper Bound Effect" = Upper.Bound.Effect))
+            if(medians){
+              return(list("Control Median" = median(control),
+                          "Treatment Median" = median(treatment),
+                          "Grand Median" = mean.of.means,
+                          "Control CDF" = LPM_ratio.1,
+                          "Treatment CDF" = LPM_ratio.2,
+                          "Certainty" = min(1, NNS.ANOVA.rho * pop.adjustment),
+                          "Lower Bound Effect" = Lower.Bound.Effect,
+                          "Upper Bound Effect" = Upper.Bound.Effect))
+            } else {
+              return(list("Control Mean" = mean(control),
+                          "Treatment Mean" = mean(treatment),
+                          "Grand Mean" = mean.of.means,
+                          "Control CDF" = LPM_ratio.1,
+                          "Treatment CDF" = LPM_ratio.2,
+                          "Certainty" = min(1, NNS.ANOVA.rho * pop.adjustment),
+                          "Lower Bound Effect" = Lower.Bound.Effect,
+                          "Upper Bound Effect" = Upper.Bound.Effect))      
+              
+            }
+        
   }
 }
diff --git a/man/NNS.ANOVA.Rd b/man/NNS.ANOVA.Rd
index 9dae0c637..0fa912581 100644
--- a/man/NNS.ANOVA.Rd
+++ b/man/NNS.ANOVA.Rd
@@ -8,6 +8,7 @@ NNS.ANOVA(
   control,
   treatment,
   means.only = FALSE,
+  medians = FALSE,
   confidence.interval = 0.95,
   tails = "Both",
   pairwise = FALSE,
@@ -22,6 +23,8 @@ NNS.ANOVA(
 
 \item{means.only}{logical; \code{FALSE} (default) test whether difference in sample means only is zero.}
 
+\item{medians}{logical; \code{FALSE} (default) test whether difference in sample medians only is zero.  Requires \code{means.only = TRUE}.}
+
 \item{confidence.interval}{numeric [0, 1]; The confidence interval surrounding the \code{control} mean, defaults to \code{(confidence.interval = 0.95)}.}
 
 \item{tails}{options: ("Left", "Right", "Both").  \code{tails = "Both"}(Default) Selects the tail of the distribution to determine effect size.}
@@ -59,6 +62,9 @@ NNS.ANOVA(control = x, treatment = y)
 A <- cbind(x, y)
 NNS.ANOVA(A)
 
+### Medians test
+NNS.ANOVA(A, means.only = TRUE, medians = TRUE)
+
 ### Multiple variable analysis with no control variable
 set.seed(123)
 x <- rnorm(100) ; y <- rnorm(100) ; z <- rnorm(100)