Merge pull request #1515 from exploratory-io/fix/issue-36055

Add probability distribution option to One-Sample Proportion Test and Two-Sample t-Test

Author: kei51e

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: exploratory
 Type: Package
 Title: R package for Exploratory
-Version: 15.1.11
+Version: 15.1.12
 Date: 2026-05-30
 Authors@R: c(person("Hideaki", "Hayashi", email = "hideaki@exploratory.io", role = c("aut", "cre")), person("Hide", "Kojima", email = "hide@exploratory.io", role = c("aut")), person("Kan", "Nishida", email = "kan@exploratory.io", role = c("aut")), person("Kei", "Saito", email = "kei@exploratory.io", role = c("aut")), person("Yosuke", "Yasuda", email = "double.y.919.quick@gmail.com", role = c("aut")))
 URL: https://github.com/exploratory-io/exploratory_func

R/test_wrapper.R

@@ -1294,6 +1294,33 @@ tidy.ttest_exploratory <- function(x, type="model", conf_level=0.95) {
     ret <- generate_ttest_density_data(x$statistic, p.value=x$p.value, x$parameter, sig_level=x$test_sig_level, alternative=x$alternative)
     ret
   }
+  else if (type == "prob_dist_diff") {
+    if ("error" %in% class(x)) {
+      ret <- tibble::tibble()
+      return(ret)
+    }
+    # Data for the probability distribution (line) chart on the DIFFERENCE scale. This is
+    # the sampling distribution of the difference of means under H0 "no difference", i.e.
+    # the t-value distribution relocated/rescaled onto the difference axis via
+    # difference = 0 + t * stderr. We reuse the t-value generator and remap its x column
+    # so the shading semantics (critical region, sig.level based) stay identical to the
+    # t-value tab. The marked statistic point lands at t * stderr = the observed
+    # difference, and the curve peaks at 0 (the hypothesized difference).
+    # estimate is c(mean of group1, mean of group2) for an unpaired test, or a single
+    # value (mean of the differences) for a paired test; stderr = difference / t.
+    est <- unname(x$estimate)
+    difference <- if (length(est) >= 2) est[1] - est[2] else est[1]
+    t_stat <- unname(x$statistic)
+    se_val <- if (!is.na(t_stat) && t_stat != 0) abs(difference / t_stat) else NA_real_
+    if (is.na(se_val) || se_val == 0) {
+      ret <- tibble::tibble(x = numeric(0), y = numeric(0)) # degenerate input -> chart no-ops
+    } else {
+      ret <- generate_ttest_density_data(t = t_stat, p.value = x$p.value, df = unname(x$parameter),
+                                         sig_level = x$test_sig_level, alternative = x$alternative) %>%
+        dplyr::mutate(x = x * se_val)
+    }
+    ret
+  }
   else { # type == "data"
     if ("error" %in% class(x)) {
       ret <- tibble::tibble()
@@ -3429,7 +3456,7 @@ exp_one_sample_prop_test <- function(df, var, p = 0.5, alternative = "two.sided"
 
       model <- list(htest = res, x = x, n = n, p = p, observed_prop = observed_prop,
                     success_value = "TRUE", alternative = alternative,
-                    method_used = method_used, conf_level = conf_level, z = z,
+                    method_used = method_used, conf_level = conf_level, z = z, se = se,
                     cohens_h = cohens_h, power = power, var_col = var_col,
                     sig.level = sig.level)
       class(model) <- c("one_sample_prop_test_exploratory", class(model))
@@ -3488,6 +3515,22 @@ tidy.one_sample_prop_test_exploratory <- function(x, type = "model") {
     # statistic (x$z, the same value shown in the summary table) on it.
     generate_norm_density_data(x$z, p.value = x$htest$p.value, mu = 0, sigma = 1,
                                sig_level = x$sig.level, alternative = x$alternative)
+  } else if (type == "prob_dist_prop") {
+    # Data for the probability distribution (line) chart on the PROPORTION scale.
+    # This is the sampling distribution of the sample proportion under H0
+    # "the proportion equals the benchmark": a normal distribution centered at the
+    # benchmark proportion (x$p) with standard deviation equal to the standard error
+    # sqrt(p0 * (1 - p0) / n). We mark the observed proportion (x$observed_prop, the
+    # same value shown in the summary table) on it. The shaded critical region matches
+    # the Z value chart (sig.level based), so both probability distribution charts read
+    # the same way.
+    if (is.na(x$se) || x$se <= 0) {
+      tibble::tibble(x = numeric(0), y = numeric(0)) # degenerate input -> chart no-ops
+    } else {
+      generate_norm_density_data(x$observed_prop, p.value = x$htest$p.value, mu = x$p, sigma = x$se,
+                                 sig_level = x$sig.level, alternative = x$alternative) %>%
+        dplyr::filter(x >= 0 & x <= 1)
+    }
   } else { # type == "data"
     # Return the original data so that the data-level charts (Error Bar Plot,
     # Data Distribution) can map to the target column. The chart templates

tests/testthat/test_one_sample_prop_test.R

@@ -360,3 +360,60 @@ test_that("tidy prob_dist critical region follows the alternative direction", {
   ts <- tidy(exp_one_sample_prop_test(df, outcome, p = 0.4, alternative = "two.sided", method = "approximate")$model[[1]], type = "prob_dist")
   expect_true(any(ts$x[which(ts$critical)] > 0) && any(ts$x[which(ts$critical)] < 0))
 })
+
+# --- tidy(type = "prob_dist_prop") ---
+
+test_that("tidy prob_dist_prop marks the observed proportion on the proportion scale", {
+  df <- data.frame(outcome = c(rep(TRUE, 50), rep(FALSE, 50)))
+  model <- exp_one_sample_prop_test(df, outcome, p = 0.4, method = "approximate")$model[[1]]
+  pd <- tidy(model, type = "prob_dist_prop")
+  expect_true(all(c("x", "y", "statistic", "critical", "p.value") %in% names(pd)))
+  # SE = sqrt(p0*(1-p0)/n); the marked point sits at the observed proportion (NOT z).
+  se <- sqrt(0.4 * 0.6 / 100)
+  marked <- pd[which(pd$statistic), ]
+  expect_equal(nrow(marked), 1)
+  expect_equal(marked$x, 0.5)
+  expect_equal(marked$p.value, model$htest$p.value)
+  # Normal centered at the benchmark proportion p0 with sd = SE.
+  expect_equal(unique(pd$mean), 0.4)
+  expect_equal(unique(pd$sd), se)
+  expect_true(max(pd$y, na.rm = TRUE) <= dnorm(0.4, mean = 0.4, sd = se) + 1e-9)
+})
+
+test_that("tidy prob_dist_prop always uses the proportion-scale normal even for the exact method", {
+  df <- data.frame(outcome = c(rep(TRUE, 12), rep(FALSE, 88)))
+  model <- exp_one_sample_prop_test(df, outcome, p = 0.1, method = "exact")$model[[1]]
+  expect_equal(model$method_used, "Exact Binomial Test")
+  pd <- tidy(model, type = "prob_dist_prop")
+  se <- sqrt(0.1 * 0.9 / 100)
+  marked <- pd[which(pd$statistic), ]
+  expect_equal(nrow(marked), 1)
+  expect_equal(marked$x, 0.12)
+  expect_equal(unique(pd$mean), 0.1)
+  expect_equal(unique(pd$sd), se)
+})
+
+test_that("tidy prob_dist_prop critical region follows the alternative direction", {
+  df <- data.frame(outcome = c(rep(TRUE, 50), rep(FALSE, 50)))
+  p0 <- 0.4
+  # greater -> rejection region is the upper tail only (above the benchmark proportion).
+  gt <- tidy(exp_one_sample_prop_test(df, outcome, p = p0, alternative = "greater", method = "approximate")$model[[1]], type = "prob_dist_prop")
+  expect_true(all(gt$x[which(gt$critical)] > p0))
+  # less -> rejection region is the lower tail only.
+  lt <- tidy(exp_one_sample_prop_test(df, outcome, p = p0, alternative = "less", method = "approximate")$model[[1]], type = "prob_dist_prop")
+  expect_true(all(lt$x[which(lt$critical)] < p0))
+  # two.sided -> both tails relative to the benchmark proportion.
+  ts <- tidy(exp_one_sample_prop_test(df, outcome, p = p0, alternative = "two.sided", method = "approximate")$model[[1]], type = "prob_dist_prop")
+  expect_true(any(ts$x[which(ts$critical)] > p0) && any(ts$x[which(ts$critical)] < p0))
+})
+
+test_that("tidy prob_dist_prop returns an empty tibble when the SE is 0/NA", {
+  # n*p0 produces a valid SE normally; build a degenerate model object directly to
+  # exercise the guard (the benchmark SE is 0 only for a degenerate p0).
+  model <- structure(
+    list(htest = list(p.value = NA_real_), p = 0.4, observed_prop = 0.5,
+         se = 0, sig.level = 0.05, alternative = "two.sided"),
+    class = c("one_sample_prop_test_exploratory", "list"))
+  pd <- tidy(model, type = "prob_dist_prop")
+  expect_equal(nrow(pd), 0)
+})

tests/testthat/test_test_wrapper_1.R

@@ -141,6 +141,9 @@ test_that("test exp_ttest_aggregated", {
   ret <- model_df %>% tidy_rowwise(model, type="data_summary")
   ret <- model_df %>% tidy_rowwise(model, type="prob_dist")
   expect_true("p.value" %in% colnames(ret))
+  # Difference-scale probability distribution (shares the ttest_exploratory tidy branch).
+  ret_diff <- model_df %>% tidy_rowwise(model, type="prob_dist_diff")
+  expect_true(all(c("x", "y", "statistic", "critical", "p.value") %in% colnames(ret_diff)))
 })
 
 test_that("test two sample t-test with column name", {

tests/testthat/test_test_wrapper_2.R

@@ -231,3 +231,38 @@ test_that("test exp_wilcox with group-level error", {
   ret <- model_df %>% tidy_rowwise(model, type='prob_dist')
   expect_equal(nrow(ret), 0)
 })
+
+# --- tidy(type = "prob_dist_diff") for two-sample t-test ---
+
+test_that("exp_ttest tidy prob_dist_diff marks the observed difference on the difference scale", {
+  # Two clearly separated groups so the mean difference is non-trivial.
+  test_df <- data.frame(
+    cat = c(rep("cat1", 20), rep("cat2", 20)),
+    val = c(seq(1, 20), seq(4, 23))
+  )
+  model_df <- test_df %>% exp_ttest(val, cat)
+  pd <- model_df %>% tidy_rowwise(model, type = "prob_dist_diff")
+  expect_true(all(c("x", "y", "statistic", "critical", "p.value") %in% names(pd)))
+  diff <- (model_df %>% tidy_rowwise(model, type = "model"))$Difference
+  marked <- pd[which(pd$statistic), ]
+  expect_equal(nrow(marked), 1)
+  # The marked point sits at the observed mean difference, NOT the standardized t.
+  expect_equal(marked$x, diff, tolerance = 1e-6)
+  # The difference-scale distribution is centered at 0 (H0: no difference): the density
+  # peak (the curve rows, excluding the marked statistic point) is at x = 0.
+  curve <- pd[!is.na(pd$y) & is.na(pd$statistic), ]
+  expect_equal(curve$x[which.max(curve$y)], 0, tolerance = 1e-6)
+})
+
+test_that("exp_ttest tidy prob_dist_diff critical region follows the alternative direction", {
+  test_df <- data.frame(
+    cat = c(rep("cat1", 20), rep("cat2", 20)),
+    val = c(seq(1, 20), seq(4, 23))
+  )
+  gt <- test_df %>% exp_ttest(val, cat, alternative = "greater") %>% tidy_rowwise(model, type = "prob_dist_diff")
+  expect_true(all(gt$x[which(gt$critical)] > 0))
+  lt <- test_df %>% exp_ttest(val, cat, alternative = "less") %>% tidy_rowwise(model, type = "prob_dist_diff")
+  expect_true(all(lt$x[which(lt$critical)] < 0))
+  ts <- test_df %>% exp_ttest(val, cat, alternative = "two.sided") %>% tidy_rowwise(model, type = "prob_dist_diff")
+  expect_true(any(ts$x[which(ts$critical)] > 0) && any(ts$x[which(ts$critical)] < 0))
+})