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..c3ae0c4b 100644 --- a/R/anovarmnp.b.R +++ b/R/anovarmnp.b.R @@ -31,19 +31,16 @@ 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 - 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] - )) - 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) @@ -115,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/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/R/descriptives.b.R b/R/descriptives.b.R index b5efbfb2..01df0be1 100644 --- a/R/descriptives.b.R +++ b/R/descriptives.b.R @@ -1000,7 +1000,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 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 diff --git a/tests/testthat/testanovarmnp.R b/tests/testthat/testanovarmnp.R index 74fcfcf8..3ba4e157 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 @@ -35,10 +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 ) +}) - # 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) + +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 + # 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 + # 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/ + # 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")) + ) + + # WHEN calling `calcDurbin` on the randomly generated matrix + r <- calcDurbin(y) + + # 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) + ) })