From c444fa6b47d993ad76705d47988ca8b2638e8bae Mon Sep 17 00:00:00 2001 From: Sebastian Jentschke Date: Sun, 16 Nov 2025 23:16:27 +1100 Subject: [PATCH 1/5] Adapt to changes in the handling of the table function (changes in R 4.6) --- R/descriptives.b.R | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/R/descriptives.b.R b/R/descriptives.b.R index b5efbfb2..45f9c2cb 100644 --- a/R/descriptives.b.R +++ b/R/descriptives.b.R @@ -990,7 +990,11 @@ descriptivesClass <- R6::R6Class( next() table <- tables$get(var) - freq <- freqs[[var]] + if (hasName(freqs, var)) { + freq <- freqs[[var]] + } else { + freq <- freqs + } tableVars <- c(var, splitBy) allLevels <- lapply(jmvcore::select(self$data, tableVars), levels) @@ -1000,7 +1004,7 @@ descriptivesClass <- R6::R6Class( cumsum <- 0 for (row in seq_len(nrow(grid))) { - counts <- do.call("[", c(list(freq), grid[row, ])) + counts <- as.numeric(do.call("[", c(list(freq), grid[row, ]))) cumsum <- cumsum + counts pc <- counts / n cumpc <- cumsum / n From 6040d6239c11ce3cde1bc7a7dce350329fee0fa3 Mon Sep 17 00:00:00 2001 From: Sebastian Jentschke Date: Fri, 16 Jan 2026 13:07:59 +0100 Subject: [PATCH 2/5] Replaced / re-implemented Durbin-Conover-test (removes dependency on PMCMR, which is deprecated) --- DESCRIPTION | 1 - R/00jmv.R | 7 ------- R/anovarmnp.b.R | 32 ++++++++++++++++++++++++-------- R/anovarmnp.h.R | 1 - jamovi/00refs.yaml | 8 -------- jamovi/anovarmnp.r.yaml | 1 - 6 files changed, 24 insertions(+), 26 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index dc1784e9..c52a5df2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,7 +25,6 @@ Imports: car (>= 3.0.0), multcomp, ggplot2 (>= 2.2.1), - PMCMR, emmeans (>= 1.4.2), vcd, vcdExtra, diff --git a/R/00jmv.R b/R/00jmv.R index 1cdf44b9..a7999905 100644 --- a/R/00jmv.R +++ b/R/00jmv.R @@ -66,13 +66,6 @@ `title`="mvnormtest: Normality test for multivariate variables", `publisher`="[R package]. Retrieved from https://cran.r-project.org/package=mvnormtest", `url`="https://cran.r-project.org/package=mvnormtest"), - `PMCMR`=list( - `type`="software", - `author`="Pohlert, T.", - `year`=2021, - `title`="PMCMR: Calculate Pairwise Multiple Comparisons of Mean Rank Sums", - `publisher`="[R package]. Retrieved from https://cran.r-project.org/package=PMCMR", - `url`="https://cran.r-project.org/package=PMCMR"), `ROCR`=list( `type`="software", `author`="Sing, T., Sander, O., Beerenwinkel, N., & Lengauer, T., Unterthiner, T., & Ernst, F. G. M.", diff --git a/R/anovarmnp.b.R b/R/anovarmnp.b.R index c4d38a5a..30d9b67b 100644 --- a/R/anovarmnp.b.R +++ b/R/anovarmnp.b.R @@ -31,16 +31,32 @@ anovaRMNPClass <- R6::R6Class( if (self$options$get('pairs')) { table <- self$results$get('comp') - result <- PMCMR::posthoc.durbin.test(mat, p.adjust='none') - n <- length(measureNames)-1 + # adapted from https://www.nist.gov/ -> search "Friedman Test" and + # https://en.wikipedia.org/wiki/Durbin_test + # NB: assumes that rows with NAs are omitted + # b denotes blocks (participant / unit of observation) + b <- nrow(mat) + # k denotes treatments (conditions / measures) + k <- ncol(mat) + + # assign ranks per block and calculate the rank sum per treatment + Rij <- apply(mat, 1, rank) + Rj <- apply(Rij, 1, sum) + + # use the formulas on the NIST web page to calculate the t-value + # for the rank sum, which is then used to determine the p-value + A1 <- sum(Rij ^ 2) + C1 <- (b * k * (k + 1) ^ 2) / 4 + T1 <- (k - 1) / (A1 - C1) * (sum(Rj ^ 2) - b * C1) + DF <- (b - 1) * (k - 1) + D <- sqrt(((2 * (A1 - C1) * b) / DF) * (1 - T1 / (b * (k - 1)))) + rowNo <- 1 - for (j in 1:n) { - for (k in j:n) { - table$setRow(rowNo=rowNo, list( - stat=result$statistic[k,j], - p=result$p.value[k,j] - )) + for (i in 1:(k - 1)) { + for (j in (i + 1):k) { + T <- abs(Rj[i] - Rj[j]) / D + table$setRow(rowNo = rowNo, list(stat = T, p = 2 * (1 - pt(T, DF)))) rowNo <- rowNo + 1 } } diff --git a/R/anovarmnp.h.R b/R/anovarmnp.h.R index 49040d82..9adb9989 100644 --- a/R/anovarmnp.h.R +++ b/R/anovarmnp.h.R @@ -105,7 +105,6 @@ anovaRMNPResults <- if (requireNamespace("jmvcore", quietly=TRUE)) R6::R6Class( options=options, name="comp", title="Pairwise Comparisons (Durbin-Conover)", - refs="PMCMR", visible="(pairs)", clearWith=list( "measures"), diff --git a/jamovi/00refs.yaml b/jamovi/00refs.yaml index 7e143e1c..188c97cf 100644 --- a/jamovi/00refs.yaml +++ b/jamovi/00refs.yaml @@ -73,14 +73,6 @@ refs: publisher: '[R package]. Retrieved from https://cran.r-project.org/package=mvnormtest' url: https://cran.r-project.org/package=mvnormtest - PMCMR: - type: 'software' - author: Pohlert, T. - year: 2021 - title: 'PMCMR: Calculate Pairwise Multiple Comparisons of Mean Rank Sums' - publisher: '[R package]. Retrieved from https://cran.r-project.org/package=PMCMR' - url: https://cran.r-project.org/package=PMCMR - ROCR: type: 'software' author: Sing, T., Sander, O., Beerenwinkel, N., & Lengauer, T., Unterthiner, T., & Ernst, F. G. M. diff --git a/jamovi/anovarmnp.r.yaml b/jamovi/anovarmnp.r.yaml index 420ffe16..55be57a9 100644 --- a/jamovi/anovarmnp.r.yaml +++ b/jamovi/anovarmnp.r.yaml @@ -30,7 +30,6 @@ items: title: Pairwise Comparisons (Durbin-Conover) type: Table description: a table of the pairwise comparisons - refs: PMCMR visible: (pairs) clearWith: - measures From 5a04cc1dc26a4c99174f67487461a3697db67b60 Mon Sep 17 00:00:00 2001 From: Sebastian Jentschke Date: Wed, 25 Feb 2026 11:59:37 +0900 Subject: [PATCH 3/5] Moved calculation of Durbin-Conover-Test into a separate function (calcDurbin), added option for p-adjustment, and added further unit tests --- R/anovarmnp.b.R | 63 ++++++++++++++++++++-------------- tests/testthat/testanovarmnp.R | 55 ++++++++++++++++++++++++++--- 2 files changed, 87 insertions(+), 31 deletions(-) diff --git a/R/anovarmnp.b.R b/R/anovarmnp.b.R index 30d9b67b..c3ae0c4b 100644 --- a/R/anovarmnp.b.R +++ b/R/anovarmnp.b.R @@ -32,34 +32,15 @@ anovaRMNPClass <- R6::R6Class( table <- self$results$get('comp') - # adapted from https://www.nist.gov/ -> search "Friedman Test" and - # https://en.wikipedia.org/wiki/Durbin_test - # NB: assumes that rows with NAs are omitted - # b denotes blocks (participant / unit of observation) - b <- nrow(mat) - # k denotes treatments (conditions / measures) - k <- ncol(mat) - - # assign ranks per block and calculate the rank sum per treatment - Rij <- apply(mat, 1, rank) - Rj <- apply(Rij, 1, sum) - - # use the formulas on the NIST web page to calculate the t-value - # for the rank sum, which is then used to determine the p-value - A1 <- sum(Rij ^ 2) - C1 <- (b * k * (k + 1) ^ 2) / 4 - T1 <- (k - 1) / (A1 - C1) * (sum(Rj ^ 2) - b * C1) - DF <- (b - 1) * (k - 1) - D <- sqrt(((2 * (A1 - C1) * b) / DF) * (1 - T1 / (b * (k - 1)))) - - rowNo <- 1 - for (i in 1:(k - 1)) { - for (j in (i + 1):k) { - T <- abs(Rj[i] - Rj[j]) / D - table$setRow(rowNo = rowNo, list(stat = T, p = 2 * (1 - pt(T, DF)))) - rowNo <- rowNo + 1 + resDurbin <- calcDurbin(mat) + + for (i in seq_len(nrow(resDurbin))) { + if (table$getCell("i1", rowNo = i)$value == resDurbin[i, "i1"] && + table$getCell("i2", rowNo = i)$value == resDurbin[i, "i2"]) { + table$setRow(rowNo = i, list(stat = resDurbin[i, "t"], p = resDurbin[i, "p"])) } } + } private$.preparePlot(data) @@ -131,3 +112,33 @@ anovaRMNPClass <- R6::R6Class( return(p) }) ) + +calcDurbin=function(mat, adjustP = "none") { + adjustP <- match.arg(adjustP, p.adjust.methods) + + # adapted from https://www.nist.gov/ -> search "Friedman Test" and + # https://en.wikipedia.org/wiki/Durbin_test + # NB: assumes that rows with NAs are omitted + # b denotes blocks (participant / unit of observation) + b <- nrow(mat) + # k denotes treatments (conditions / measures) + k <- ncol(mat) + + # assign ranks per block and calculate the rank sum per treatment + Rij <- apply(mat, 1, rank) + Rj <- apply(Rij, 1, sum) + + # use the formulas on the NIST web page to calculate the t-value + # for the rank sum, which is then used to determine the p-value + A1 <- sum(Rij ^ 2) + C1 <- (b * k * (k + 1) ^ 2) / 4 + T1 <- (k - 1) / (A1 - C1) * (sum(Rj ^ 2) - b * C1) + DF <- (b - 1) * (k - 1) + DN <- sqrt(((2 * (A1 - C1) * b) / DF) * (1 - T1 / (b * (k - 1)))) + + combPairs <- t(combn(colnames(mat), 2)) + valT <- apply(combPairs, 1, function(c) abs(diff(Rj[c])) / DN) + + data.frame(i1 = combPairs[, 1], i2 = combPairs[, 2], + t = valT, p = stats::p.adjust(2 * (1 - pt(valT, DF)), adjustP)) +} diff --git a/tests/testthat/testanovarmnp.R b/tests/testthat/testanovarmnp.R index 74fcfcf8..97670ea7 100644 --- a/tests/testthat/testanovarmnp.R +++ b/tests/testthat/testanovarmnp.R @@ -24,6 +24,11 @@ testthat::test_that('All options in the anovaRMNP work (sunny)', { testthat::expect_equal(3, mainTable[['df']]) testthat::expect_equal(0.532, mainTable[['p']], tolerance = 1e-3) + # Test descriptives table + descTable <-r$desc$asDF + testthat::expect_equal(c('x1', 'x2', 'x3', 'x4'), descTable[['level']]) + testthat::expect_equal(as.vector(sapply(data, mean)), descTable[['mean']], tolerance = 1e-3) + testthat::expect_equal(as.vector(sapply(data, median)), descTable[['median']], tolerance = 1e-3) # Test pairwise comparisons pairsTable <- r$comp$asDF @@ -36,9 +41,49 @@ testthat::test_that('All options in the anovaRMNP work (sunny)', { c(0.382, 0.95, 0.574, 0.349, 0.154, 0.617), pairsTable[['p']], tolerance = 1e-3 ) - # Test descriptives table - descTable <-r$desc$asDF - testthat::expect_equal(c('x1', 'x2', 'x3', 'x4'), descTable[['level']]) - testthat::expect_equal(as.vector(sapply(data, mean)), descTable[['mean']], tolerance = 1e-3) - testthat::expect_equal(as.vector(sapply(data, median)), descTable[['median']], tolerance = 1e-3) + # further tests of calcDurbin -------------------------------------------------------------------------------------- + # example for a complete block design from + # https://www.rdocumentation.org/packages/PMCMRplus/versions/4.0/topics/posthoc.durbin.test + y <- matrix(c(3.88, 5.64, 5.76, 4.25, 5.91, 4.33, 30.58, 30.14, 16.92, 23.19, + 26.74, 10.91, 25.24, 33.52, 25.45, 18.85, 20.45, 26.67, 4.44, 7.94, 4.04, + 4.4, 4.23, 4.36, 29.41, 30.72, 32.92, 28.23, 23.35, 12, 38.87, 33.12, 39.15, + 28.06, 38.23, 26.65), nrow = 6, ncol = 6, dimnames=list(1:6, LETTERS[1:6])) + r <- calcDurbin(y) + expect_equal(r[, 1], c("A", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "D", "D", "E")) + expect_equal(r[, 2], c("B", "C", "D", "E", "F", "C", "D", "E", "F", "D", "E", "F", "E", "F", "F")) + expect_equal(r[, 3], c(4.4821073, 5.0797216, 0.5976143, 5.6773359, 7.470179, 0.5976143, 3.8844930, 1.1952286, + 2.9880720, 4.4821073, 0.5976143, 2.390457, 5.0797216, 6.8725650, 1.792843), tolerance = 1e-6) + expect_equal(r[, 4], c(0.0001426546, 0.000030327980, 0.55547245630, 0.0000065360, 0.00000008006288, 0.55547245626, + 0.0006661548, 0.243213110008, 0.00621361995, 0.0001426546, 0.555472456255, 0.02468018576662, + 0.000030327983, 0.00000033329893, 0.08510737977997)) + # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff < 1e-15 + # using the same dataset with Holm correction for multiple comparisons + r <- calcDurbin(y, "holm") + expect_equal(r[, 4], c(0.001426545736, 0.000363935798, 1, 0.000084968305, 0.000001200943, 1, 0.005329238412, + 0.972852440031, 0.043495339635, 0.001426545736, 1, 0.148081114600, 0.000363935798, + 0.000004666185, 0.425536898900)) + # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="holm"), diff < 1e-14 + + # example from https://real-statistics.com/one-way-analysis-of-variance-anova/kruskal-wallis-test/conover-test-after-kw/ + y <- matrix(c(18, 49, 33, 19, 24, 17, 48, 22, 35, 30, 81, 32, 42, 62, 37, 44, 38, 47, 49, 31, 48, 31, 25, + 22, 30, 30, 32, 15, 64, 50), ncol = 3, dimnames = list(seq(10), c("Ctl", "New", "Old"))) + r <- calcDurbin(y) + expect_equal(r[, 1], c("Ctl", "Ctl", "New")) + expect_equal(r[, 2], c("New", "Old", "Old")) + expect_equal(r[, 3], c(2.5, 0.5, 2.0)) + expect_equal(r[, 4], c(0.02230802, 0.62313246, 0.06082147)) + # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff < 1e-15 + + # example from https://medium.com/@serurays/introduction-to-statistical-testing-in-r-part-3-non-parametric-tests-2cde06ae2893 + set.seed(9) + y <- matrix(c(sample(1:10, 30, replace = TRUE, prob = c(0.20, 0.20, 0.20, 0.15, 0.15, 0.10, 0.05, 0.05, 0, 0)), + sample(1:10, 30, replace = TRUE, prob = c(0.10, 0.10, 0.15, 0.20, 0.20, 0.10, 0.10, 0.10, 0, 0)), + sample(1:10, 30, replace = TRUE, prob = c(0.05, 0.05, 0.10, 0.15, 0.25, 0.20, 0.10, 0.05, 0, 0))), + ncol = 3, dimnames = list(seq(30), c("A", "B", "C"))) + r <- calcDurbin(y) + expect_equal(r[, 1], c("A", "A", "B")) + expect_equal(r[, 2], c("B", "C", "C")) + expect_equal(r[, 3], c(2.74883920, 2.67650132, 0.07233787)) + expect_equal(r[, 4], c(0.007959595, 0.009656408, 0.942581904)) + # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff <- 1e-16 }) From f84579a1139d9551cec6ebd7a40c79bb718e790c Mon Sep 17 00:00:00 2001 From: Sebastian Jentschke Date: Wed, 25 Feb 2026 12:23:22 +0900 Subject: [PATCH 4/5] Reverted earlier commit: 'Adapt to changes in the handling of the table function (changes in R 4.6)' --- R/descriptives.b.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/R/descriptives.b.R b/R/descriptives.b.R index 45f9c2cb..01df0be1 100644 --- a/R/descriptives.b.R +++ b/R/descriptives.b.R @@ -990,11 +990,7 @@ descriptivesClass <- R6::R6Class( next() table <- tables$get(var) - if (hasName(freqs, var)) { - freq <- freqs[[var]] - } else { - freq <- freqs - } + freq <- freqs[[var]] tableVars <- c(var, splitBy) allLevels <- lapply(jmvcore::select(self$data, tableVars), levels) From 49ffa05b6f68a9e5211a0113ffd9a6901c424c8c Mon Sep 17 00:00:00 2001 From: Sebastian Jentschke Date: Thu, 26 Feb 2026 12:06:44 +0900 Subject: [PATCH 5/5] Updated unit tests for anovarmnp (separated tests for calcDurbin) --- tests/testthat/testanovarmnp.R | 171 ++++++++++++++++++++++++++------- 1 file changed, 138 insertions(+), 33 deletions(-) diff --git a/tests/testthat/testanovarmnp.R b/tests/testthat/testanovarmnp.R index 97670ea7..3ba4e157 100644 --- a/tests/testthat/testanovarmnp.R +++ b/tests/testthat/testanovarmnp.R @@ -40,50 +40,155 @@ testthat::test_that('All options in the anovaRMNP work (sunny)', { testthat::expect_equal( c(0.382, 0.95, 0.574, 0.349, 0.154, 0.617), pairsTable[['p']], tolerance = 1e-3 ) +}) + - # further tests of calcDurbin -------------------------------------------------------------------------------------- +testthat::test_that('calcDurbin works for a complete block design without p-value adjustment', { # example for a complete block design from # https://www.rdocumentation.org/packages/PMCMRplus/versions/4.0/topics/posthoc.durbin.test - y <- matrix(c(3.88, 5.64, 5.76, 4.25, 5.91, 4.33, 30.58, 30.14, 16.92, 23.19, - 26.74, 10.91, 25.24, 33.52, 25.45, 18.85, 20.45, 26.67, 4.44, 7.94, 4.04, - 4.4, 4.23, 4.36, 29.41, 30.72, 32.92, 28.23, 23.35, 12, 38.87, 33.12, 39.15, - 28.06, 38.23, 26.65), nrow = 6, ncol = 6, dimnames=list(1:6, LETTERS[1:6])) - r <- calcDurbin(y) - expect_equal(r[, 1], c("A", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "D", "D", "E")) - expect_equal(r[, 2], c("B", "C", "D", "E", "F", "C", "D", "E", "F", "D", "E", "F", "E", "F", "F")) - expect_equal(r[, 3], c(4.4821073, 5.0797216, 0.5976143, 5.6773359, 7.470179, 0.5976143, 3.8844930, 1.1952286, - 2.9880720, 4.4821073, 0.5976143, 2.390457, 5.0797216, 6.8725650, 1.792843), tolerance = 1e-6) - expect_equal(r[, 4], c(0.0001426546, 0.000030327980, 0.55547245630, 0.0000065360, 0.00000008006288, 0.55547245626, - 0.0006661548, 0.243213110008, 0.00621361995, 0.0001426546, 0.555472456255, 0.02468018576662, - 0.000030327983, 0.00000033329893, 0.08510737977997)) # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff < 1e-15 + + # GIVEN a matrix `y` representing a complete block design + y <- matrix( + c( + 3.88, 5.64, 5.76, 4.25, 5.91, 4.33, 30.58, 30.14, 16.92, 23.19, 26.74, 10.91, + 25.24, 33.52, 25.45, 18.85, 20.45, 26.67, 4.44, 7.94, 4.04, 4.4, 4.23, 4.36, 29.41, + 30.72, 32.92, 28.23, 23.35, 12, 38.87, 33.12, 39.15, 28.06, 38.23, 26.65 + ), nrow = 6, ncol = 6, dimnames=list(1:6, LETTERS[1:6]) + ) + + # WHEN calling `calcDurbin` with default arguments (no p-value adjustment) + r <- calcDurbin(y) + + # THEN the returned pairwise comparisons are as expected + expect_equal( + r[, 1], + c("A", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "D", "D", "E") + ) + expect_equal( + r[, 2], + c("B", "C", "D", "E", "F", "C", "D", "E", "F", "D", "E", "F", "E", "F", "F") + ) + # AND the test statistics are as expected + expect_equal( + r[, 3], + c( + 4.4821073, 5.0797216, 0.5976143, 5.6773359, 7.470179, 0.5976143, 3.8844930, 1.1952286, + 2.9880720, 4.4821073, 0.5976143, 2.390457, 5.0797216, 6.8725650, 1.792843 + ), tolerance = 1e-6 + ) + # AND the p-values are as expected + expect_equal( + r[, 4], + c( + 0.0001426546, 0.000030327980, 0.55547245630, 0.0000065360, 0.00000008006288, + 0.55547245626, 0.0006661548, 0.243213110008, 0.00621361995, 0.0001426546, + 0.555472456255, 0.02468018576662, 0.000030327983, 0.00000033329893, 0.08510737977997 + ) + ) +}) + +testthat::test_that('calcDurbin works for a complete block design with Holm correction', { # using the same dataset with Holm correction for multiple comparisons - r <- calcDurbin(y, "holm") - expect_equal(r[, 4], c(0.001426545736, 0.000363935798, 1, 0.000084968305, 0.000001200943, 1, 0.005329238412, - 0.972852440031, 0.043495339635, 0.001426545736, 1, 0.148081114600, 0.000363935798, - 0.000004666185, 0.425536898900)) + # example for a complete block design from + # https://www.rdocumentation.org/packages/PMCMRplus/versions/4.0/topics/posthoc.durbin.test # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="holm"), diff < 1e-14 + + # GIVEN a matrix `y` representing a complete block design + y <- matrix( + c( + 3.88, 5.64, 5.76, 4.25, 5.91, 4.33, 30.58, 30.14, 16.92, 23.19, 26.74, 10.91, + 25.24, 33.52, 25.45, 18.85, 20.45, 26.67, 4.44, 7.94, 4.04, 4.4, 4.23, 4.36, 29.41, + 30.72, 32.92, 28.23, 23.35, 12, 38.87, 33.12, 39.15, 28.06, 38.23, 26.65 + ), nrow = 6, ncol = 6, dimnames=list(1:6, LETTERS[1:6]) + ) + # WHEN calling `calcDurbin` with "holm" p-value adjustment + r <- calcDurbin(y, "holm") + + # THEN the returned p-values are as expected + expect_equal( + r[, 4], + c( + 0.001426545736, 0.000363935798, 1, 0.000084968305, 0.000001200943, 1, 0.005329238412, + 0.972852440031, 0.043495339635, 0.001426545736, 1, 0.148081114600, 0.000363935798, + 0.000004666185, 0.425536898900 + ) + ) +}) + +testthat::test_that('calcDurbin works with data from real-statistics.com (Conover test after KW)', { # example from https://real-statistics.com/one-way-analysis-of-variance-anova/kruskal-wallis-test/conover-test-after-kw/ - y <- matrix(c(18, 49, 33, 19, 24, 17, 48, 22, 35, 30, 81, 32, 42, 62, 37, 44, 38, 47, 49, 31, 48, 31, 25, - 22, 30, 30, 32, 15, 64, 50), ncol = 3, dimnames = list(seq(10), c("Ctl", "New", "Old"))) - r <- calcDurbin(y) - expect_equal(r[, 1], c("Ctl", "Ctl", "New")) - expect_equal(r[, 2], c("New", "Old", "Old")) - expect_equal(r[, 3], c(2.5, 0.5, 2.0)) - expect_equal(r[, 4], c(0.02230802, 0.62313246, 0.06082147)) # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff < 1e-15 + + # GIVEN a smaller matrix `y` representing data across three groups + y <- matrix( + c( + 18, 49, 33, 19, 24, 17, 48, 22, 35, 30, + 81, 32, 42, 62, 37, 44, 38, 47, 49, 31, + 48, 31, 25, 22, 30, 30, 32, 15, 64, 50 + ), ncol = 3, dimnames = list(seq(10), c("Ctl", "New", "Old")) + ) + # WHEN calling `calcDurbin` on this data + r <- calcDurbin(y) + + # THEN the returned pairwise comparisons are as expected + expect_equal( + r[, 1], + c("Ctl", "Ctl", "New") + ) + expect_equal( + r[, 2], + c("New", "Old", "Old") + ) + # AND the test statistics are as expected + expect_equal( + r[, 3], + c(2.5, 0.5, 2.0) + ) + # AND the p-values are as expected + expect_equal( + r[, 4], + c(0.02230802, 0.62313246, 0.06082147) + ) +}) + +testthat::test_that('calcDurbin works with randomly generated data', { # example from https://medium.com/@serurays/introduction-to-statistical-testing-in-r-part-3-non-parametric-tests-2cde06ae2893 + # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff <- 1e-16 + + # GIVEN a matrix `y` generated with randomized probabilities using a fixed seed set.seed(9) - y <- matrix(c(sample(1:10, 30, replace = TRUE, prob = c(0.20, 0.20, 0.20, 0.15, 0.15, 0.10, 0.05, 0.05, 0, 0)), - sample(1:10, 30, replace = TRUE, prob = c(0.10, 0.10, 0.15, 0.20, 0.20, 0.10, 0.10, 0.10, 0, 0)), - sample(1:10, 30, replace = TRUE, prob = c(0.05, 0.05, 0.10, 0.15, 0.25, 0.20, 0.10, 0.05, 0, 0))), - ncol = 3, dimnames = list(seq(30), c("A", "B", "C"))) + y <- matrix( + c( + sample(1:10, 30, replace = TRUE, prob = c(0.20, 0.20, 0.20, 0.15, 0.15, 0.10, 0.05, 0.05, 0, 0)), + sample(1:10, 30, replace = TRUE, prob = c(0.10, 0.10, 0.15, 0.20, 0.20, 0.10, 0.10, 0.10, 0, 0)), + sample(1:10, 30, replace = TRUE, prob = c(0.05, 0.05, 0.10, 0.15, 0.25, 0.20, 0.10, 0.05, 0, 0)) + ), + ncol = 3, dimnames = list(seq(30), c("A", "B", "C")) + ) + + # WHEN calling `calcDurbin` on the randomly generated matrix r <- calcDurbin(y) - expect_equal(r[, 1], c("A", "A", "B")) - expect_equal(r[, 2], c("B", "C", "C")) - expect_equal(r[, 3], c(2.74883920, 2.67650132, 0.07233787)) - expect_equal(r[, 4], c(0.007959595, 0.009656408, 0.942581904)) - # results are equal to those from PMCMRplus::durbinAllPairsTest(y, p.adj="none"), diff <- 1e-16 + + # THEN the returned pairwise comparisons are as expected + expect_equal( + r[, 1], + c("A", "A", "B") + ) + expect_equal( + r[, 2], + c("B", "C", "C") + ) + # AND the test statistics are as expected + expect_equal( + r[, 3], + c(2.74883920, 2.67650132, 0.07233787) + ) + # AND the p-values are as expected + expect_equal( + r[, 4], + c(0.007959595, 0.009656408, 0.942581904) + ) })