# To report test coverage, run # devtools::test_coverage() makeData <- function(N = 40) { set.seed(1) data <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 100, sd = 50), rnorm(N, mean = 120, sd = 25), rnorm(N, mean = 80, sd = 50), rnorm(N, mean = 100, sd = 12), rnorm(N, mean = 100, sd = 50)), Group = c(rep("ZControl1", N), rep("Control2", N), rep("Group1", N), rep("Group2", N), rep("Group3", N), rep("Group4", N)), Gender = rep(c(rep('Male', N/2), rep('Female', N/2)), 6), ID = rep(1:N, 6) ) # Shuffle data[sample(nrow(data)), ] } makeES1 <- function() { set.seed(1) data <- makeData() DurgaDiff(data[data$Group %in% c("ZControl1", "Group1"),], data.col = "Measurement", group.col = "Group", R = 1000) } makePairedData <- function(addSomeNAs = FALSE, reverseGroups = FALSE) { N <- 40 set.seed(1) data <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 100, sd = 50), rnorm(N, mean = 120, sd = 25), rnorm(N, mean = 80, sd = 50), rnorm(N, mean = 100, sd = 12), rnorm(N, mean = 100, sd = 50)), Group = c(rep("ZControl1", N), rep("Control2", N), rep("Group1", N), rep("Group2", N), rep("Group3", N), rep("Group4", N)), Sex = rep(c(rep('Male', N/2), rep('Female', N/2)), 6), ID = rep(1:N, 6) ) if (addSomeNAs) data[c(1, 4, 10, 65), 1] <- NA # Shuffle data <- data[sample(nrow(data)), ] if (reverseGroups) { DurgaDiff(data[data$Group %in% c("ZControl1", "Group1"),], groups = c("ZControl1", "Group1"), na.rm = addSomeNAs, data.col = "Measurement", group.col = "Group", R = 1000, id.col = "ID") } else { DurgaDiff(data[data$Group %in% c("ZControl1", "Group1"),], na.rm = addSomeNAs, data.col = "Measurement", group.col = "Group", R = 1000, id.col = "ID") } } compareDiffs <- function(d1, d2, tolerance = 0.1) { expect_equal(d1$groups, d2$groups) expect_equal(d1$group.names, d2$group.names) expect_equal(d1$effect.type, d2$effect.type) expect_equal(length(d1$group.differences), length(d2$group.differences)) for (i in seq_len(length(d1$group.differences))) { pd1 <- d1$group.differences[[i]] pd2 <- d2$group.differences[[i]] expect_equal(pd1$groups, pd2$groups) expect_equal(pd1$t0, pd2$t0) expect_equal(pd1$R, pd2$R) expect_equal(pd1$bca[, 1], pd2$bca[, 1]) # Columns 3 and 4 can vary randomly expect_equal(pd1$bca[, c(4, 5)], pd2$bca[, c(4, 5)], tolerance = tolerance) } } plotWithBrackets <- function(d) { ylim <- extendrange(d$data[[d$data.col]]) ylim[2] <- extendrange(ylim, f = nrow(d$group.statistics) * 0.05)[2] labelFn <- function(diff) sprintf("%.2g %g%% CI [%.2g, %.2g]",diff$t0, diff$bca[1] * 100, diff$bca[4], diff$bca[5]) DurgaPlot(d, ef.size = FALSE, bty = "n", ylim = ylim) |> DurgaBrackets(labels = labelFn) } ########################################################################## # Tests start here #### test_that("multiple groups + contrasts", { n <- 30 df <- data.frame(group = c(rep(c("G1", "G2", "G3"), each = n)), group2 = c(rep(c("A", "B"), each = 3 * n / 2)), val = c(rnorm(n, 1), rnorm(n, 2), rnorm(n, 3))) d <- expect_error(DurgaDiff(df, "val", c("group", "group2")), NA) expect_equal(length(d$group.differences), 6) # Group names should be compoiste groups after combining the two group columns plotFn <- plotWithBrackets plotFn <- identity d <- expect_error(DurgaDiff(df, "val", c("group", "group2"), contrasts = c("G2, A - G1, A", "G3, B - G2, B", "G2, B - G2, A")), NA) expect_equal(length(d$group.differences), 3) # Attempting to use the comma-separated string version should fail due to ambiguity expect_error(DurgaDiff(df, "val", c("group", "group2"), contrasts = c("G2, A - G1, A, G3, B - G2, B"))) # Matrix form should work contrasts <- matrix(c("G2, A", "G1, A", "G3, B", "G2, B", "G2, B", "G2, A"), nrow = 2) d <- DurgaDiff(df, "val", c("group", "group2"), contrasts = contrasts) expect_equal(length(d$group.differences), 3) }) test_that("wide to long", { n <- 10 dfw <- data.frame(control = rnorm(n), treatment = rnorm(n, 0.2), id = seq_len(n)) dfl <- wideToLong(dfw, c("control", "treatment")) expect_equal(nrow(dfl), 2 * n) # The generated id values should be the same as ours expect_equal(dfl$id.1, dfl$id) # Group values should all be one of "control", "treatment" expect_true(all(unique(dfl$group) %in% c("control", "treatment"))) # Specifying a NULL id col should be the same as not specifying one dfln <- wideToLong(dfw, c("control", "treatment"), NULL) expect_equal(dfln, dfl) # Specifying an id column means a new one should not be generated dfli <- wideToLong(dfw, c("control", "treatment"), "id") expect_equal(names(dfli), c("id", "group", "value")) expect_equal(unique(dfli$group), c("control", "treatment")) # Unpaired in wide format n1 <- 20 n2 <- 14 dfw <- data.frame(control = rnorm(n1), treatment = c(rnorm(n2, 0.2), rep(NA, n1 - n2))) dfl <- wideToLong(dfw, c("control", "treatment"), NULL) expect_equal(unique(dfl$group), c("control", "treatment")) # Paired dp <- DurgaDiff(dfw, groups = c("control", "treatment"), na.rm = TRUE) expect_true(dp$paired.data) # When paired, any ids without matched values a removed from the data set expect_equal(sum(dp$data$group == "control"), min(n1, n2)) expect_equal(sum(dp$data$group == "treatment"), min(n1, n2)) # Unpaired du <- DurgaDiff(dfw, groups = c("control", "treatment"), id.col = NULL, na.rm = TRUE) expect_true(!du$paired.data) # When unpaired, unmatched values remain, only NAs are removed expect_equal(sum(du$data$group == "control"), n1) expect_equal(sum(du$data$group == "treatment"), n2) # 3 groups n <- 20 df3 <- data.frame(control = rnorm(n), g1 = rnorm(n, 1), g2 = rnorm(n, 2)) d3 <- DurgaDiff(df3, groups = c("control", "g1", "g2")) #DurgaPlot(d3, bty = "n") |> DurgaBrackets() expect_equal(d3$group.names, c("control", "g1", "g2")) expect_equal(unique(d3$data$group), c("control", "g1", "g2")) expect_equal(length(d3$group.differences), 3) du3 <- DurgaDiff(df3, groups = c("control", "g1", "g2"), id.col = NULL) #DurgaPlot(du3, bty = "n") |> DurgaBrackets() expect_equal(du3$group.names, c("control", "g1", "g2")) expect_equal(unique(du3$data$group), c("control", "g1", "g2")) expect_equal(length(du3$group.differences), 3) # # How can I use wide format but include another column in the groups? # dfw <- data.frame(control = rnorm(n), treatment = rnorm(n, 1.2), id = seq_len(n), # species = sample(c("Sp 1", "Sp 2", "Sp 3"), n, replace = TRUE)) # # ANSWER: There is no builtin solution. First convert wide to long, then pass # # to DurgaDiff and specify the 2nd group variable # dfl <- wideToLong(dfw, c("control", "treatment")) # DurgaDiff(dfl, "value", c("group", "species")) |> DurgaPlot() }) test_that("Combine vars", { df <- data.frame(v1 = c("A", "A", "B", "B"), v2 = c("C", "D", "C", "D")) expect_equal(combineVariables(df, c("v1", "v2")), c("A & C", "A & D", "B & C", "B & D")) df <- data.frame(a = rep(1:3, each = 4), b = rep(1:4, 3)) expect_equal(combineVariables(df, 1:2), c("1 & 1", "1 & 2", "1 & 3", "1 & 4", "2 & 1", "2 & 2", "2 & 3", "2 & 4", "3 & 1", "3 & 2", "3 & 3", "3 & 4")) # Pathological case. Ideally this test would pass, but it has not been implemented to work, so don't run the test # df <- data.frame(t = c("a", "a", "a & b", "b & c", "a"), # s = c("b", "b & c", "c", "c", "b & c")) # expect_true(all(table(combineVariables(df, 1:2)) == 1)) }) test_that("Group interactions", { set.seed(1) n <- 40 df2 <- data.frame(group = rep(c("A", "B"), each = n), sex = sample(c("Female", "Male"), 2 * n, replace = TRUE)) df2$val <- ifelse(df2$group == "A", 1, 2) + ifelse(df2$sex == "Female", 1, 0) + rnorm(2 * n) par(mfrow = c(2, 2)) # Define interaction by specifying 2 group variables da <- DurgaDiff(df2, "val", c("group", "sex")) expect_equal(length(da$groups), 4) expect_equal(length(da$group.differences), 6) DurgaPlot(da, bty = "n", ylim = c(-1.5, 6), ef.size = FALSE, main = "Auto interaction") |> DurgaBrackets() # Manual interaction df2$group_sex <- combineVariables(df2, c("group", "sex"), ", ") dm <- DurgaDiff(df2, "val", "group_sex") expect_equal(length(dm$groups), 4) expect_equal(length(dm$group.differences), 6) DurgaPlot(dm, bty = "n", ylim = c(-1.5, 6), ef.size = FALSE, main = "Manual interaction") |> DurgaBrackets() compareDiffs(da, dm) # 3 groups set.seed(1) n <- 100 df3 <- data.frame(group = rep(c("A", "B"), each = n), sex = sample(c("Female", "Male"), 2 * n, replace = TRUE), filter = sample(c("Filter", "No filter"), 2 * n, replace = TRUE)) df3$val <- ifelse(df3$group == "A", 5, 7) + ifelse(df3$sex == "Female", 1, 0) + ifelse(df3$filter == "Yes", 0.6, 0) + rnorm(2 * n) da <- DurgaDiff(df3, "val", c("group", "sex", "filter")) expect_equal(length(da$groups), 8) expect_equal(length(da$group.differences), 28) DurgaPlot(da, bty = "n", main = "Auto interaction") # Manual interaction df3$group_sex <- combineVariables(df3, c("group", "sex", "filter"), ", ") dm <- DurgaDiff(df3, "val", "group_sex") expect_equal(length(dm$groups), 8) expect_equal(length(dm$group.differences), 28) DurgaPlot(dm, bty = "n", main = "Manual interaction") compareDiffs(da, dm) # Formula interface set.seed(1) par(mfrow = c(1, 2)) # Define interaction by specifying 2 group variables f2 <- expect_error(DurgaDiff(val ~ group*sex, df2), NA) expect_equal(length(f2$groups), 4) expect_equal(length(f2$group.differences), 6) DurgaPlot(f2, bty = "n", main = "Formula interface") # Define interaction by specifying 3 group variables f3 <- expect_error(DurgaDiff(val ~ group + sex + filter, df3), NA) DurgaPlot(f3, bty = "n", main = "Formula interface") }) test_that("Paired vs unpaired", { set.seed(1) n <- 10 effSz <- 1 ind <- rnorm(n, 10, 10) df <- data.frame(value = c(ind, ind + effSz + rnorm(n, 0, 1)), group = rep(c("A", "B"), each = n), id = rep(1:n, 2)) dpd <- DurgaDiff(df, "value", "group", "id", effect.type = "cohens d") dud <- DurgaDiff(df, "value", "group", effect.type = "cohens d") # Cohen's d should be identical for paired and unpaired expect_equal(dpd$group.differences[[1]]$t0, dud$group.differences[[1]]$t0) # CI width should be less for paired than unpaired expect_lt(diff(dpd$group.differences[[1]]$bca[4:5]), diff(dud$group.differences[[1]]$bca[4:5])) # Expect bias-correction to change a bit due to different degrees of freedom dpg <- DurgaDiff(df, "value", "group", "id", effect.type = "hedges g") dug <- DurgaDiff(df, "value", "group", effect.type = "hedges g") expect_false(isTRUE(all.equal(dpg$group.differences[[1]]$t0, dug$group.differences[[1]]$t0))) # CI width should be less for paired than unpaired expect_lt(diff(dpg$group.differences[[1]]$bca[4:5]), diff(dug$group.differences[[1]]$bca[4:5])) # Bias correction should have almost no effect for large sample sizes n <- 500 ind <- rnorm(n, 10, 10) df <- data.frame(value = c(ind, ind + effSz + rnorm(n, 0, 1)), group = rep(c("A", "B"), each = n), id = rep(1:n, 2)) dpg <- DurgaDiff(df, "value", "group", "id", effect.type = "hedges g") dug <- DurgaDiff(df, "value", "group", effect.type = "hedges g") # Don't expect identical but pretty close expect_true(isTRUE(all.equal(dpg$group.differences[[1]]$t0, dug$group.differences[[1]]$t0, tolerance = 1e-4, scale = 1))) }) test_that("Effect size visible", { es <- 20 set.seed(1) vA <- rnorm(10) df <- data.frame(group = c(rep("A", 10), rep("B", 10)), val = c(vA, vA + es + rnorm(10, 0.5)), id = rep(1:10, 2)) # This requires a human to check that the plot is valid expect_error(DurgaPlot(DurgaDiff(val ~ group, df, effect.type = "cohens d"), main = "Is effect size entirely visible?"), NA) expect_error(DurgaPlot(DurgaDiff(val ~ group, df, id.col = "id", effect.type = "cohens d"), main = "Is effect size entirely visible?"), NA) #plot(cohens_d(dabestr::dabest(df, group, val, idx = c("A", "B")))) #plot(cohens_d(dabestr::dabest(df, group, val, idx = c("A", "B"), paired = TRUE, id.column = id))) }) test_that("paired with groups", { data <- structure(list(Spider = c("A1", "A11", "A12", "A13", "A14", "A2", "A3", "A4", "A5", "A6", "A7", "A9", "B2", "B3", "C1", "C2", "D3", "E1", "E10", "E11", "E2", "E3", "E4", "E5", "E6", "E7", "E8", "E9", "F1", "F2", "F3", "G01", "G02", "G03 ", "G04", "G05", "G06", "G07", "G08 ", "G09", "G10", "G11", "G12", "G13", "G14", "G15", "G16", "G17", "G18", "G19", "G20", "H01", "H02", "H03", "H04", "H05", "H06 ", "H07", "H08", "H09", "H10", "H11", "H12", "H13", "H14", "H16", "I01", "I02", "I03", "I04", "I05", "I06", "I07", "I08", "I09", "I10", "I11", "A1", "A11", "A12", "A13", "A14", "A2", "A3", "A4", "A5", "A6", "A7", "A9", "B2", "B3", "C2", "D3", "E1", "E11", "E2", "E3", "E5", "E6", "E7", "E8", "E9", "F1", "F2", "F3", "G01", "G02", "G03 ", "G04", "G05", "G06", "G07", "G08 ", "G09", "G10", "G11", "G12", "G15", "G16", "G17", "G18", "G19", "H01", "H02", "H03", "H04", "H05", "H06 ", "H07", "H08", "H09", "H10", "H11", "H13", "H14", "H15", "H16", "I01", "I02", "I03", "I04", "I05", "I06", "I07", "I09", "I10", "I11", "A1", "A11", "A12", "A13", "A2", "A3", "A4", "A5", "A6", "A7", "A9", "B2", "B3", "C2", "D3", "E1", "E3", "E5", "E6", "E7", "E8", "E9", "F1", "F2", "F3", "G03 ", "G04", "G06", "G07", "G08 ", "G09", "G10", "G12", "G13", "G14", "G15", "G16", "G17", "G18", "G19", "G20", "H01", "H02", "H03", "H04", "H05", "H06 ", "H09", "H10", "H11", "H12", "H13", "H14", "H15", "I01", "I02", "I03", "I04", "I05", "I06", "I07", "I08", "I09", "I10", "I11", "A1", "A11", "A12", "A13", "A14", "A3", "A4", "A5", "A6", "A7", "A9", "B2", "B3", "C2", "E5", "E6", "E7", "E8", "E9", "F1", "F3", "G01", "G04", "G06", "G07", "G08 ", "G09", "G10", "H01", "H02", "H03", "H04", "H06 ", "H09", "H10", "H11", "H12", "H13", "H14", "H15"), Wrap = c(185.75, 81.25, 99.07, 301.11, 111.7, 141.6, 156.4, 276.14, 49.51, 278.94, 279.22, 170.69, 93.82, 663.66, 83.94, 155.05, 69.26, 81.62, 92.87, 36.92, 143.82, 59.84, 165.04, 56.52, 156.91, 69.26, 70.28, 99.11, 523.59, 111.91, 168.2, 631.32, 341.04, 110.86, 197.9, 289.96, 372.66, 423.38, 459.65, 120.76, 82.51, 151.3, 190.98, 68.71, 20.3, 54.18, 75.85, 80.64, 74.43, 36.4, 66.74, 108.52, 62.86, 89.06, 631.77, 219.78, 144.55, 62.01, 62.08, 371.95, 233.48, 135.79, 151.9, 210.16, 239.6, 347.78, 144.8, 177.46, 145.5, 124.84, 95.25, 119.38, 277.49, 77, 129.55, 164.2, 225.38, 131.59, 376.31, 147.63, 246.17, 148.48, 158.91, 149.49, 123.72, 131.69, 121.09, 121.72, 108.95, 224.66, 74.31, 143.31, 300.31, 318.36, 43.14, 92.9, 67.48, 131.74, 42.59, 76.07, 301.75, 38.75, 316.33, 95.66, 126.45, 171.76, 304.77, 105.25, 208.41, 257.78, 32.41, 193.61, 118.46, 90.07, 420.58, 151.14, 340.76, 77.67, 95.17, 124.56, 71.13, 130.84, 160.56, 350.49, 34.07, 468.27, 257.55, 117.44, 113.99, 95.45, 181.4, 145.99, 197.04, 139.07, 380.36, 131.72, 248.03, 297.9, 77.07, 72.37, 100.41, 641.88, 70.42, 91.82, 200.94, 83.34, 257.72, 194.93, 160.14, 104.9, 157.38, 146.97, 107.91, 197.39, 125.3, 509.2, 247.17, 343.11, 293.38, 385.48, 237.19, 256.04, 263.29, 129.2, 175.92, 86.63, 97.96, 97.59, 102.42, 204.43, 398.39, 111.3, 62.04, 64.15, 73.39, 7.04, 91.56, 190.14, 86.57, 57.65, 85.45, 98.11, 52.84, 43.77, 36.91, 134.45, 59.39, 365.75, 80.97, 156.96, 80.52, 142.83, 105.88, 83.21, 154.51, 151.05, 101.14, 124.24, 71.53, 183.64, 54.16, 100.55, 58.28, 59.27, 175.66, 195.94, 203.02, 112.22, 221.76, 176.28, 201.72, 64.77, 46.92, 87.95, 41.39, 67.25, 54.22, 45.99, 51.69, 38.44, 66.77, 68.12, 70.05, 197.7, 124.19, 475.43, 153.2, 143.39, 215.45, 152.27, 393.36, 108.23, 190.35, 123.31, 142.69, 111.12, 53.67, 105.46, 123.39, 119.8, 81.22, 60.74, 83.96, 95.81, 84.77, 170.54, 163.63, 156.01, 173.24, 161.43, 127.645, 164.02), Prey_difficulty = c("difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "easy", "easy", "easy", "easy", "easy", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult", "difficult")), class = "data.frame", row.names = c(NA, -252L)) ts14 <- c("A1", "A11", "A12", "A13", "A3", "A4", "A5", "A6", "A7", "A9", "B2", "B3", "C2", "E5", "E6", "E7", "E8", "E9", "F1", "F3", "G04", "G06", "G07", "G08 ", "G09", "G10", "H01", "H02", "H03", "H04", "H06 ", "H09", "H10", "H11", "H13", "H14") # Order of groups should not affect errors expect_error(DurgaDiff(data[data$Spider %in% ts14, ], "Wrap", "Prey_difficulty", id = "Spider"), "5 Spiders") expect_error(DurgaDiff(data[data$Spider %in% ts14, ], "Wrap", "Prey_difficulty", id = "Spider", groups = c("difficult", "easy")), "5 Spiders") expect_error(DurgaDiff(data[data$Spider %in% ts14, ], "Wrap", "Prey_difficulty", id = "Spider", groups = c("easy", "difficult")), "5 Spiders") }) # TODO # test_that("No CI", { # df <- data.frame(val = rnorm(10000), grp = rep(c("G1", "G2"), exch = 5000)) # d <- DurgaDiff(df, 1, 2, R = 5000) # }) test_that("transparency", { expect_equal(DurgaTransparent("red", 0.8), "#FF000033") expect_equal(DurgaTransparent("#ff0000", 0.8), "#FF000033") expect_equal(DurgaTransparent("red", 0.8, TRUE), "#FF000033") expect_equal(DurgaTransparent("red", 0.8, FALSE), "#FF000033") expect_equal(DurgaTransparent("#ff0000cc", 0.8), "#FF000033") expect_equal(DurgaTransparent("#ff0000cc", 0.8, TRUE), "#FF000029") expect_equal(DurgaTransparent("#ff0000cc", 0.8, FALSE), "#FF000033") }) test_that("expand contrasts", { .makeX <- function(g) matrix(g, nrow = 2) expect_equal(expandContrasts("g1 - g2, g3 - g4", c("g1", "g2", "g3", "g4")), .makeX(c("g1", "g2", "g3", "g4")), ignore_attr = TRUE) expect_error(expandContrasts("g1 - g2, g3 - g5", c("g1", "g2", "g3", "g4"))) expect_error(expandContrasts("g1 - g1", c("g1", "g2", "g3", "g4"))) expect_equal(expandContrasts(" g1-g2 ", c("g1", "g2")), .makeX(c("g1", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts(" g1 -g2 ", c("g2", "g1")), .makeX(c("g1", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts(" g2 - g1 ", c("g2", "g1")), .makeX(c("g2", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts("g1 - g2", c("g4", "g1", "g2")), .makeX(c("g1", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts("g1 - g2", c("g4", "g2", "g1")), .makeX(c("g1", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts(". - g1", c("g1", "g2", "g3")), .makeX(c("g2", "g1", "g3", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts(". - g2", c("g1", "g2", "g3")), .makeX(c("g1", "g2", "g3", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts("g3-.", c("g1", "g2", "g3")), .makeX(c("g3", "g1", "g3", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts("*", c("g1", "g2")), .makeX(c("g2", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts("*", c("g2", "g1")), .makeX(c("g1", "g2")), ignore_attr = TRUE) expect_equal(expandContrasts("*", c("g1")), NULL, ignore_attr = TRUE) expect_equal(expandContrasts(" g1 - g11 ", c("g1", "g11")), .makeX(c("g1", "g11")), ignore_attr = TRUE) expect_equal(expandContrasts(" g1 - g11 ", c("g11", "g1")), .makeX(c("g1", "g11")), ignore_attr = TRUE) expect_equal(expandContrasts(" g11 - g1 ", c("g1", "g11")), .makeX(c("g11", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts(" . - g1", c("g1", "g11")), .makeX(c("g11", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts(" . - g1", c("g1", "g11", "g111")), .makeX(c("g11", "g1", "g111", "g1")), ignore_attr = TRUE) expect_equal(expandContrasts(" g1 - .", c("g1", "g11")), .makeX(c("g1", "g11")), ignore_attr = TRUE) expect_equal(expandContrasts(" g1 - . ", c("g1", "g11", "g111")), .makeX(c("g1", "g11", "g1", "g111")), ignore_attr = TRUE) expect_error(expandContrasts(" . - . ", c("g1", "g11"))) expect_equal(expandContrasts("g1 - g2", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g1", "g2"))) expect_equal(expandContrasts("G1 - G2", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g1", "g2"))) expect_equal(expandContrasts("g2 - g1", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts("G2 - G1", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts("G2 - g1", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts("g3 - g2", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g3", "g2"))) expect_equal(expandContrasts("G3 - G2", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g3", "g2"))) expect_equal(expandContrasts("g4 - g1", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g4", "g1"))) expect_equal(expandContrasts("G4 - G1, g3-G2", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g4", "g1", "g3", "g2"))) expect_equal(expandContrasts(c("G4 - G1", "g3-G2"), c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g4", "g1", "g3", "g2"))) expect_equal(expandContrasts("g2 - g1", c("g1", "g2", "g3", "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts("g2 - g1", c(G1 = "g1", G2 = "g2", G3 = "g3", G4 = "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts("G2 - G1", c(G1 = "g1", G2 = "g2", G3 = "g3", G4 = "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1"))) expect_equal(expandContrasts(c("G2 - G1", "G3 - G1"), c(G1 = "g1", G2 = "g2", G3 = "g3", G4 = "g4"), c("G1", "G2", "G3", "G4")), .makeX(c("g2", "g1", "g3", "g1"))) }) test_that("group names and contrasts", { data <- makeData() groups <- c(Ctrl = "ZControl1", `G 1` = "Group1", `G 2` = "Group2") d <- DurgaDiff(data, "Measurement", "Group", groups = groups) DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = "*") # Contrasts with data values DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = c("Group2 - Group1")) # Contrasts with group labels DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = c("G 2 - G 1")) # Contrasts with mixture of labels and data values DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = c("G 2 - Group1")) # Each pair of plots in a row should be identical. One uses group data value, the other uses group label # Seems to be a bug? in R, restoring mfrow after DurgaPlot changes the margins! mar <- par("mar") op <- par(mfrow = c(2, 2)) on.exit(par(op)) expect_error(DurgaPlot(d, contrasts = "Group1 - ZControl1", main = "1 contrast"), NA) expect_error(DurgaPlot(d, contrasts = "G 1 - Ctrl", main = "1 contrast"), NA) expect_error(DurgaPlot(d, contrasts = "Group1 - ZControl1, Group2 - ZControl1", main = "2 contrasts"), NA) expect_error(DurgaPlot(d, contrasts = "G 1 - Ctrl, G 2 - Ctrl", main = "2 contrasts"), NA) par(mar = mar) }) test_that("matrix data", { n <- 20 m <- matrix(c(val = rnorm(n), sample(1:3, n, replace = TRUE)), nrow = n) d <- DurgaDiff(m, 1, 2) expect_error(DurgaPlot(d), NA) }) test_that("violin symbology", { op <- par(mfrow = c(1, 2)) data <- makeData() d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2")) expect_error(DurgaPlot(d, main = "Default symbology"), NA) expect_error(DurgaPlot(d, main = "Custom violins", violin.params = list(lwd = 2, ljoin = 2, density = 10)), NA) par(op) }) test_that("violin no fill", { data <- makeData() d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group3", "Group4")) expect_error(DurgaPlot(d, main = "No violin fill", violin.width = 2, violin.fill = F, violin.params = list(lwd = 2), ef.size.violin.fill = F, ef.size.violin.shape = "full"), NA) }) test_that("contrast plots", { data <- makeData() groups <- c("ZControl1", "Group1", "Group2", "Group3") # Default contrasts d <- DurgaDiff(data, "Measurement", "Group", groups = groups) ng <- length(d$groups) expect_equal(length(d$group.differences), ng * (ng - 1) / 2) DurgaPlot(d, main = "Default contrasts") # 2 contrasts contrasts <- "Group1 - ZControl1, Group2 - ZControl1, Group3 - ZControl1" d <- DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = contrasts) expect_equal(length(d$group.differences), 3) expect_equal(d$group.differences[[1]]$groups[1], "Group1") expect_equal(d$group.differences[[1]]$groups[2], "ZControl1") expect_equal(d$group.differences[[2]]$groups[1], "Group2") expect_equal(d$group.differences[[2]]$groups[2], "ZControl1") expect_equal(d$group.differences[[3]]$groups[1], "Group3") expect_equal(d$group.differences[[3]]$groups[2], "ZControl1") DurgaPlot(d, violin = F, central.tendency.width = 0.1, central.tendency.symbol = "segment", central.tendency.params = , main = "Explicit contrasts") # Shorthand for same as above d <- DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = ". - ZControl1") expect_equal(length(d$group.differences), 3) expect_equal(d$group.differences[[1]]$groups[1], "Group1") expect_equal(d$group.differences[[1]]$groups[2], "ZControl1") expect_equal(d$group.differences[[2]]$groups[1], "Group2") expect_equal(d$group.differences[[2]]$groups[2], "ZControl1") expect_equal(d$group.differences[[3]]$groups[1], "Group3") expect_equal(d$group.differences[[3]]$groups[2], "ZControl1") DurgaPlot(d, points.method = "jitter", main = "Explicit contrast shorthand") # Invalid contrasts expect_error(DurgaDiff(data, "Measurement", "Group", groups = groups, contrasts = 1)) expect_error(DurgaPlot(d, contrasts = 1)) # Just 1 contrast data <- makeData() d <- DurgaDiff(data, "Measurement", "Group", groups = groups) DurgaPlot(d, contrasts = c(`Diff` = "Group3 - Group2"), main = "Plot 1 labelled diff") DurgaPlot(d, contrasts = c(`Diff` = "Group2 - Group3"), mai = "Plot negative diff") d <- DurgaDiff(data, "Measurement", "Group", groups = groups, effect.type = "cohens d") DurgaPlot(d, contrasts = c(`Diff` = "Group3 - Group2"), main = "Plot 1 Cohen's") DurgaPlot(d, contrasts = c(`Diff` = "Group2 - Group3"), mai = "Plot negative Cohen's") # Single contrast in diff d <- DurgaDiff(data, "Measurement", "Group", groups = c("Group3", "Group2")) DurgaPlot(d, main = "Restricted groups in diff") d <- DurgaDiff(data, "Measurement", "Group", groups = c("Group3", "Group2"), effect.type = "cohens d") DurgaPlot(d, main = "Restricted groups in diff Cohen's") }) test_that("group stats", { checkGroup <- function(gs, vals, popMean) { gs <- unname(gs) expect_equal(gs[1], mean(vals)) expect_equal(gs[2], median(vals)) expect_equal(gs[3], sd(vals)) expect_equal(gs[4], sd(vals) / sqrt(length(vals))) expect_lt(gs[5], popMean) expect_gt(gs[6], popMean) } # Test single group at a time set.seed(3) mean <- 50 sd <- 17 ns <- c(10, 20, 80, 200) for (n in ns) { v <- rnorm(n, mean, sd) df <- data.frame(v, rep("g", length(v))) d <- DurgaDiff(df, 1, 2) checkGroup(d$group.statistics, v, mean) } # Now combine all groups into one data set set.seed(1) df <- do.call(rbind, lapply(ns, function(n) data.frame(val = rnorm(n, mean, sd), group = rep(sprintf("g%03d", n), n)))) d <- DurgaDiff(df, 1, 2) for (i in seq_along(ns)) { gn <- sprintf(sprintf("g%03d", ns[i])) checkGroup(d$group.statistics[i, ], df$val[df$group == gn], mean) } # For debugging if (FALSE) { DurgaPlot(d, violin = F, points.method = "jitter", ef.size = F, central.tendency.dx = 0.15, error.bars.type = "CI") abline(h = mean, col = "lightgrey") DurgaPlot(d, add = T, violin = F, points = F, ef.size = F, central.tendency.dx = 0.25, error.bars.type = "SD") DurgaPlot(d, add = T, violin = F, points = F, ef.size = F, central.tendency.dx = 0.35, error.bars.type = "SE") } }) test_that("print", { d <- makeES1() # This is a regex expected <- paste0("Bootstrapped effect size\\n", " Measurement ~ Group\\n", "Groups:\\n", " mean median sd se CI\\.lower CI\\.upper n\\n", "Group1 122\\.9210 118\\.8367 21\\.31850 3\\.370750 116\\.84544 130\\.4658 40\\n", "ZControl1 102\\.3007 103\\.2276 22\\.16668 3\\.504859 95\\.36523 108\\.6278 40\\n", "Unpaired Difference of means \\(R = 1000, bootstrap CI method = bca\\):\\n", " ZControl1 - Group1: -20\\.6203, 95% CI \\[-29\\.[0-9]+, -11\\.[0-9]+\\]") expect_output(print(d), expected) # Print mismatched lines if (FALSE) { got <- capture.output(print(d)) expLines <- gsub("\\\\", "", strsplit(expected, "\\\\n")[[1]]) matched <- got == expLines print(got[!matched]) print(expLines[!matched]) } checkSummaryMatches <- function(d1, d2) { d1s <- capture.output(print(d1)) d2s <- capture.output(print(d2)) # Crude check - just remove all numbers d1s <- sub("[0-9].*", "", d1s) d2s <- sub("[0-9].*", "", d2s) expect_equal(d1s, d2s) } # Check that formula output matches non-formula set.seed(1) # Avoid white-space differences caused by different printed precisions d1 <- DurgaDiff(petunia, 1, 2) d2 <- DurgaDiff(height ~ group, petunia) checkSummaryMatches(d1, d2) # Paired set.seed(1) # Avoid white-space differences caused by different printed precisions d1 <- DurgaDiff(insulin, 1, 2, 3) d2 <- DurgaDiff(sugar ~ treatment, insulin, id.col = "id") checkSummaryMatches(d1, d2) d <- DurgaDiff(log(sugar) ~ treatment, insulin, id.col = "id") s <- capture.output(print(d)) expect_equal(s[2], " log(sugar) ~ treatment") # Spaces in names `Scapus length` <- rnorm(40, 10) `Group name` <- sample(c("Group one", "Group two", "Group three"), 40, replace = TRUE) # Just check it doesn't crash expect_error(capture.output(print(DurgaDiff(`Scapus length` ~ `Group name`))), NA) expect_error(capture.output(print(DurgaDiff(log(`Scapus length`) ~ `Group name`))), NA) }) test_that("contrasts", { data <- makeData() # Default contrasts d <- DurgaDiff(data, "Measurement", "Group") ng <- length(d$groups) expect_equal(length(d$group.differences), ng * (ng - 1) / 2) # All in 1 string contrasts = "Group1 - ZControl1, Group2 - ZControl1, Group3 - ZControl1, Group4 - ZControl1" d <- DurgaDiff(data, "Measurement", "Group", contrasts = contrasts) expect_equal(length(d$group.differences), 4) expect_equal(d$group.differences[[1]]$groups[1], "Group1") expect_equal(d$group.differences[[1]]$groups[2], "ZControl1") expect_equal(d$group.differences[[2]]$groups[1], "Group2") expect_equal(d$group.differences[[2]]$groups[2], "ZControl1") expect_equal(d$group.differences[[3]]$groups[1], "Group3") expect_equal(d$group.differences[[3]]$groups[2], "ZControl1") expect_equal(d$group.differences[[4]]$groups[1], "Group4") expect_equal(d$group.differences[[4]]$groups[2], "ZControl1") contrasts = c("Group1 - ZControl1", "Group2 - ZControl1", "Group3 - ZControl1", "Group4 - ZControl1") d2 <- DurgaDiff(data, "Measurement", "Group", contrasts = contrasts) expect_equal(length(d2$group.differences), 4) expect_equal(d2$group.differences[[1]]$groups[1], "Group1") expect_equal(d2$group.differences[[1]]$groups[2], "ZControl1") compareDiffs(d2, d) contrasts <- matrix(c("Group1", "ZControl1", "Group2", "ZControl1", "Group3", "ZControl1", "Group4", "ZControl1"), nrow = 2) d2 <- DurgaDiff(data, "Measurement", "Group", contrasts = contrasts) expect_equal(length(d2$group.differences), 4) expect_equal(d2$group.differences[[1]]$groups[1], "Group1") expect_equal(d2$group.differences[[1]]$groups[2], "ZControl1") compareDiffs(d2, d, tolerance = 1) expect_error(DurgaDiff(data, "Measurement", "Group", contrasts = "Group2:ZControl")) expect_error(DurgaDiff(data, "Measurement", "Group", contrasts = "ZControl")) expect_error(DurgaDiff(data, "Measurement", "Group", contrasts = "Group 2 - Group 1")) # Allow "" as meaning no contrasts d <- DurgaDiff(data, "Measurement", "Group", contrasts = "") # Wildcard contrasts d <- DurgaDiff(data, "Measurement", "Group", effect.type = "cohens d", contrasts = "*") ng <- length(unique(data$Group)) expect_equal(length(d$group.differences), ng * (ng - 1) / 2) }) test_that("difference effect types", { n <- 100 set.seed(1) realDiff <- 1 df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + realDiff)), group = c(rep("Control", n), rep("Group", n)), id = c(1:n, 1:n)) # Check all effect types expect_error(DurgaDiff(df, effect.type = "wrong", data.col = 1, group.col = 2)) # Should throw an error expect_error(DurgaDiff(df, effect.type = "mean", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "cohens d", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "hedges g", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "mean", id.col = "id", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "cohens d", id.col = "id", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "hedges g", id.col = "id", data.col = 1, group.col = 2), NA) # Check unstandardised diff (diff of means) d <- DurgaDiff(df, data.col = 1, group.col = 2, effect.type = "mean") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") expect_lt(pwd$t0, realDiff) expect_lt(pwd$bca[4], realDiff) expect_gt(pwd$bca[5], realDiff) # Check Cohen's D d <- DurgaDiff(df, data.col = 1, group.col = 2, effect.type = "cohens d") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") expect_equal(pwd$t0, 0.918991, tolerance = 0.0001) # Cohen's d expect_lt(pwd$bca[4], 0.918991) # Should be positive but small expect_gt(pwd$bca[5], 0.918991) # Save Cohen's d for later cohensD <- pwd$t0 # Swap groups d <- DurgaDiff(df, groups = c("Group", "Control"), data.col = 1, group.col = 2, effect.type = "cohens d") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Control") expect_equal(pwd$groups[2], "Group") expect_equal(pwd$t0, -0.918991, tolerance = 0.0001) # Cohen's d expect_lt(pwd$bca[4], -0.918991) # Should be negative but small expect_gt(pwd$bca[5], -0.918991) # Check Hedges' g d <- DurgaDiff(df, data.col = 1, group.col = 2, effect.type = "hedges g") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") # Estimate Hedges' g by correcting Cohen's d (Lakens 2013, p. 3 eqn 4) n1 <- d$group.statistics[2,"n"] n2 <- d$group.statistics[1,"n"] hedgesG <- cohensD * (1 - 3 / (4 * (n1 + n2) - 9)) # Approximate method expect_equal(pwd$t0, hedgesG, tolerance = 0.0001) expect_lt(pwd$bca[4], 0.918991) # Should be positive but small expect_gt(pwd$bca[5], 0.918991) # Swap groups d <- DurgaDiff(df, groups = c("Group", "Control"), data.col = 1, group.col = 2, effect.type = "hedges g") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Control") expect_equal(pwd$groups[2], "Group") expect_equal(pwd$t0, -hedgesG, tolerance = 0.0001) expect_lt(pwd$bca[4], -hedgesG) # Should be negative but small expect_gt(pwd$bca[5], -hedgesG) ### Paired effect sizes ### # Check unstandardised diff d <- DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3, effect.type = "mean") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") expect_lt(pwd$t0, realDiff) expect_lt(pwd$bca[4], realDiff) expect_gt(pwd$bca[5], realDiff) # Check Cohen's D d <- DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3, effect.type = "cohens d_z") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") diffs <- df$val[df$group == "Group"] - df$val[df$group == "Control"] cohensD <- mean(diffs) / sd(diffs) expect_equal(pwd$t0, cohensD) expect_lt(pwd$bca[4], cohensD) # Should be positive but small expect_gt(pwd$bca[5], cohensD) d <- DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3, effect.type = "hedges g_z") pwd <- d$group.difference[[1]] expect_equal(pwd$groups[1], "Group") expect_equal(pwd$groups[2], "Control") # Estimate Hedges' g by using Hedges' approximate method to correct Cohen's d (Cummings 2012, eqn 11.13) hedgesG <- cohensD * (1 - (3 / (4 * (n - 1) - 1))) expect_equal(pwd$t0, hedgesG, tolerance = 1.5 * 10e-6) expect_lt(pwd$bca[4], hedgesG) # Should be positive but small expect_gt(pwd$bca[5], hedgesG) }) test_that("group factors", { n <- 100 df <- data.frame(val = c(rnorm(n), rnorm(n, mean = 1)), group = factor(c(rep("Control", n), rep("Treatment", n))), id = c(1:n, 1:n)) d <- DurgaDiff(df, effect.type = "mean", data.col = 1, group.col = 2) pwd <- d$group.difference[[1]] expect_equal(d$group.names, c("Control", "Treatment")) expect_equal(pwd$groups[1], "Treatment") expect_equal(pwd$groups[2], "Control") # Check all effect types expect_error(DurgaDiff(df, effect.type = "mean", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "cohens d", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, effect.type = "hedges g", data.col = 1, group.col = 2), NA) expect_error(DurgaDiff(df, id.col = "id", data.col = 1, group.col = 2), NA) }) test_that("difference handles NA", { n <- 100 df <- data.frame(val = c(rnorm(n), rnorm(n, mean = 1)), group = c(rep("Control", n), rep("Group", n))) df[c(1, 4, 10, 65), 1] <- NA # This should throw an error (something about na.rm) expect_error(DurgaDiff(df, na.rm = FALSE, data.col = 1, group.col = 2), "na.rm") # This should NOT throw an error expect_error(DurgaDiff(df, na.rm = TRUE, data.col = 1, group.col = 2), NA) # Check column name checks expect_error(DurgaDiff(df, id.col = "id", na.rm = TRUE, data.col = 1, group.col = 2), "column name") expect_error(DurgaDiff(df, data.col = "XXX", group.col = 2), "column name") expect_error(DurgaDiff(df, data.col = 1, group.col = "XXX"), "column name") }) test_that("two groups", { N <- 40 df <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 120, sd = 25)), Group = c(rep("Control", N), rep("Treatment", N)), Gender = rep(c(rep('Male', N/2), rep('Female', N/2)), 2), ID = rep(1:N, 2) ) d2 <- DurgaDiff(df, data.col = 1, group.col = 2) # This should NOT throw an error op <- par(mar = c(5, 4, 4, 4) + 0.1) expect_error(DurgaPlot(d2, ef.size = TRUE, main = "Two groups, effect size default"), NA) expect_error(DurgaPlot(d2, ef.size.pos = "right", main = "Two groups, effect size right"), NA) par(op) expect_error(DurgaPlot(d2, ef.size = FALSE, main = "Two groups, no effect size"), NA) expect_error(DurgaPlot(d2, ef.size.position = "below", main = "Two groups, effect size below"), NA) }) test_that("three groups", { N <- 40 df <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 120, sd = 25), rnorm(N, mean = 80, sd = 50)), Group = c(rep("ZControl1", N), rep("Group1", N), rep("Group2", N)), Gender = rep(c(rep('Male', N/2), rep('Female', N/2)), 3), ID = rep(1:N, 3) ) d3 <- DurgaDiff(df, data.col = 1, group.col = 2) # This should NOT throw an error expect_error(DurgaPlot(d3, bar = FALSE, box = FALSE, main = "Three groups"), NA) expect_error(DurgaPlot(d3, violin.trunc = FALSE, main = "No violin truncation"), NA) expect_error(DurgaPlot(d3, violin.trunc = 0.05, main = "0.05 violin truncation"), NA) expect_error(DurgaPlot(d3, ef.size.violin = FALSE, main = "No effect size violin"), NA) }) test_that("three groups with factor", { N <- 40 # Use factor to control group order df <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 120, sd = 25), rnorm(N, mean = 80, sd = 50)), Group = factor(c(rep("ZControl1", N), rep("Group1", N), rep("Group2", N)), levels = c("ZControl1", "Group1", "Group2")), Gender = rep(c(rep('Male', N/2), rep('Female', N/2)), 3), ID = rep(1:N, 3) ) d3 <- DurgaDiff(df, data.col = 1, group.col = 2) # This should NOT throw an error expect_error(DurgaPlot(d3, bar = FALSE, box = FALSE, main = "Group factor"), NA) }) test_that("many groups", { n <- 12 groupMean <- round(rnorm(n, mean = 20, sd = 8)) val <- c(sapply(groupMean, function(m) rnorm(n, m, 4))) groups <- sapply(seq_len(n), function(i) paste0("G", i, "-", groupMean[i])) trt <- c(sapply(seq_along(groupMean), function(i) rep(groups[i], n))) df <- data.frame(Height = val, Treatment = trt) d <- DurgaDiff(df, "Height", "Treatment", groups = groups) op <- par(cex = 0.8) expect_error(DurgaPlot(d, main = "1/3) Many groups"), NA) expect_error(DurgaPlot(d, main = "2/3) Many groups, control-.", contrasts = paste(df$Treatment[1], "-.")), NA) expect_error(DurgaPlot(d, main = "3/3) Many groups, .-control", contrasts = paste(" . - ", df$Treatment[1])), NA) par(op) }) test_that("plots work", { #### This fails with "figure margins too large" if (FALSE) { N <- 40 set.seed(1) data <- data.frame(Measurement = c(rnorm(N, mean = 100, sd = 25), rnorm(N, mean = 100, sd = 50), rnorm(N, mean = 120, sd = 25), rnorm(N, mean = 80, sd = 50), rnorm(N, mean = 100, sd = 12), rnorm(N, mean = 100, sd = 50)), Group = c(rep("ZControl1", N), rep("Control2", N), rep("Group1", N), rep("Group2", N), rep("Group3", N), rep("Group4", N)), Gender = rep(c(rep('Male', N/2), rep('Female', N/2)), 6), ID = rep(1:N, 6) ) # Shuffle data <- data[sample(nrow(data)), ] es <- DurgaDiff(data[data$Group %in% c("ZControl1", "Group1"),], data.col = "Measurement", group.col = "Group", R = 1000) es2 <- DurgaDiff(data[data$Group %in% c("ZControl1", "Group1"),], data.col = "Measurement", group.col = "Group", R = 1000, id.col = "ID") # Generate some plots l <- layout(matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20), nrow = 4, ncol = 5, byrow = TRUE)) #layout.show(l) #par(mfrow = c(2, 4)) #a) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = FALSE, box.fill = FALSE, central.tendency = "median", error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #b) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "white", box.fill = "white", central.tendency = "median", error.bars = "SD", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5),) #c) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "white", box.fill = "white", central.tendency = "median", error.bars = "SE", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5),) #d) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "white", box.fill = "white", central.tendency.type = "mean", error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5),) #e) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "white", box.fill = "white", central.tendency.type = "mean", error.bars = "SD", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5),) #f) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "white", box.fill = "white", central.tendency.type = "mean", error.bars = "SE", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5),) #g) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = DurgaTransparent(c("red", "blue"), .5), box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) #h) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = DurgaTransparent(c("red", "blue"), .5), box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #i) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) #j) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #k) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "left-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = "white", box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #l) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "left-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) ##m) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "left-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #n) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "right-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = "white", box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #o) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "right-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) #p) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "right-half", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #q) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "full", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = DurgaTransparent(c("red", "blue"), .5), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) #r) DurgaPlot(es2, bar = FALSE, bar.fill = FALSE, violin.shape = "full", violin = DurgaTransparent(c("red", "blue"), .4), violin.fill = DurgaTransparent(c("red", "blue"), .8), box = DurgaTransparent(c("grey10"), .1), box.fill = DurgaTransparent(c("red", "blue"), .7), error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = TRUE, points = DurgaTransparent(c("red", "blue"), .5), paired = TRUE) #s) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "full", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = "white", box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5)) #t) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin.shape = "full", violin = DurgaTransparent(c("red", "blue"), .6), violin.fill = DurgaTransparent(c("red", "blue"), .6), box = "white", box.fill = "white", error.bars.type = "CI", central.tendency.type = "mean", central.tendency = FALSE, ef.size = FALSE, points = FALSE) } expect_equal(1, 1) }) test_that("box FALSE works", { es <- makeES1() DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = FALSE, box.fill = FALSE, central.tendency.type = "median", error.bars.type = "CI", error.bars.cross.width = 0.05, ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5), main = "Violin FALSE, median, no effect size") DurgaPlot(es, violin = FALSE, central.tendency = FALSE, error.bars = FALSE, ef.size = FALSE, main = "No central tendency, error bar, effect size") expect_equal(1, 1) }) test_that("central tendency FALSE works", { es <- makeES1() expect_false(es$paired.data) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "red", box.fill = "blue", central.tendency = FALSE, error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5), main = "Central tendency FALSE") expect_equal(1, 1) }) test_that("paired works", { es <- makePairedData() expect_true(es$paired.data) DurgaPlot(es, paired = DurgaTransparent("green", 0.5), paired.lty = 2, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "red", box.fill = "blue", central.tendency = FALSE, error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5), main = "Paired") DurgaPlot(es, violin = FALSE, ef.size = FALSE, paired.lwd = 3, main = "Paired, no violin, no effect size") op <- par(mar = c(5, 4, 4, 4) + 0.1) on.exit(par(op)) DurgaPlot(es, violin = FALSE, ef.size = TRUE, points = FALSE, main = "Paired, no violin, effect size, no points") DurgaPlot(es, violin.shape = c("left", "right"), violin.width = 0.2, main = "Custom") # Craft simple paired data to ensure that pair lines are correct n <- 10 control <- sort(rnorm(n, 0)) before <- sort(rnorm(n, 3)) after <- sort(rnorm(n, 2.5)) treatments <- c("Control", "Before", "After") df <- data.frame(Id = rep(seq_len(n), 3), Treatment = rep(treatments, each = n), Value = c(control, before, after)) # Shuffle set.seed(1) sortedD <- DurgaDiff(df, "Value", "Treatment", "Id", groups = treatments) df <- df[sample(nrow(df)), ] set.seed(1) d <- DurgaDiff(df, "Value", "Treatment", "Id", groups = treatments) # Shuffling the rows should not affect the results compareDiffs(sortedD, d, tolerance = 0.1) DurgaPlot(d, contrasts = "Before-Control, After - Before", ef.size = FALSE, violin = FALSE, points.dx = c(0.2, 0, -0.2), main = "No intersecting lines?") expect_equal(1, 1) }) test_that("paired (reversed groups) works", { es <- makePairedData(reverseGroups = TRUE) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "red", box.fill = "blue", central.tendency = FALSE, error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5), main = "Paired with reversed groups") expect_equal(1, 1) }) test_that("paired with NAs works", { es <- makePairedData(addSomeNAs = TRUE) DurgaPlot(es, bar = FALSE, bar.fill = FALSE, violin = FALSE, box = "red", box.fill = "blue", central.tendency = FALSE, error.bars.type = "CI", ef.size = FALSE, points = DurgaTransparent(c("red", "blue"), .5), main = "Paired with some NAs") expect_equal(1, 1) }) test_that("bar charts work", { es <- makeES1() DurgaPlot(es, bar = TRUE, violin = FALSE, box = FALSE, box.fill = "blue", error.bars = FALSE, error.bars.type = "CI", ef.size = FALSE, points = FALSE, main = "Bar chart, no error bars") DurgaPlot(es, bar = TRUE, violin = FALSE, box = FALSE, box.fill = "blue", error.bars = TRUE, error.bars.type = "CI", ef.size = FALSE, points = FALSE, main = "Bar chart, error bars") expect_equal(1, 1) }) test_that("custom effect axis", { n <- 100 df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1)), group = rep(c("Control", "Group"), each = 2), id = rep(1:n, each = 2)) ef.size.ticks = c("Large negative effect" = -0.8, "Medium negative effect" = -0.5, "Small negative effect" = -0.2, "No effect" = 0, "Small positive effect" = 0.2, "Medium positive effect" = 0.5, "Large positive effect" = 0.8) d <- DurgaDiff(df, effect.type = "cohens d", data.col = 1, group.col = 2) op <- par(mar = c(5, 4, 4, 10) + 0.1) on.exit(par(op)) expect_error(DurgaPlot(d, ef.size.ticks = ef.size.ticks, ef.size.params = list(las = 1), ef.size.label = "", main = "Cohen's with custom labels"), NA) }) test_that("Axis las", { n <- 100 df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1), rnorm(n, mean = 10 + 1.4)), group = rep(c("Group1", "Group2", "Group3"), each = n)) d2 <- DurgaDiff(df, groups = c("Group1", "Group2"), data.col = 1, group.col = 2) op <- par(mar = c(5, 4, 4, 4)) expect_error(DurgaPlot(d2, las = 1, ef.size.params = list(las = 1), main = "las horizontal"), NA) par(op) d3 <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d3, las = 1, ef.size.params = list(las = 1), main = "las horizontal"), NA) }) test_that("Other data frame classes", { n <- 40 df <- tibble::tibble(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1), rnorm(n, mean = 10 + 1.4)), group = rep(c("Group1", "Group2", "Group3"), each = n)) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, main = "Tibble"), NA) df <- data.table::data.table(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1), rnorm(n, mean = 10 + 1.4)), group = rep(c("Group1", "Group2", "Group3"), each = n)) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, main = "Data.table"), NA) }) test_that("point colours", { n <- 40 df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1)), group = rep(c("Group1", "Group2"), each = n), sex = sample(factor(c("Male", "Female")), n, replace = TRUE)) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, points = 1:2, main = "Group colours"), NA) expect_error(DurgaPlot(d, points = 1:5, main = "Group colours (truncated palette)"), NA) expect_error(DurgaPlot(d, points = as.numeric(df$sex) + 1, main = "Sex colours"), NA) }) test_that("detect missing paired data", { n <- 40 # Test whether missing paired data is detected correctly set.seed(1) df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1)), group = rep(c("Group1", "Group2"), each = n), id = rep(1:n, 2)) # If we take a random sample of rows, some will be missing their paired data df <- df[sample(seq_len(nrow(df)), round(n / 2)), ] expect_error(DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3), "paired data") }) test_that("plot contrasts", { n <- 40 set.seed(1) df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 1)), group = rep(c("Group1", "Group2"), each = n), sex = sample(factor(c("Male", "Female")), n, replace = TRUE)) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, points = 1:2, main = "1/3) Default contrast"), NA) d <- DurgaDiff(df, data.col = 1, group.col = 2, contrasts = ". - Group1") expect_error(DurgaPlot(d, points = 1:2, main = "2/3) DurgaDiff contrast"), NA) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, points = 1:2, main = "3/3) DurgaPlot contrast", contrasts = ". - Group1"), NA) }) test_that("effect size position", { n <- 20 set.seed(1) df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + 0.8), rnorm(n, mean = 10 + 0.4)), group = rep(c("Group1", "Group2", "Group3"), each = n), sex = sample(factor(c("Male", "Female", "Juvenile")), n, replace = TRUE)) d <- DurgaDiff(df, data.col = 1, group.col = 2) expect_error(DurgaPlot(d, contrasts = Filter(function(d) d$bca[4] > 0 || d$bca[5] < 0, d$group.differences), main = "1/3) ef.size position default, filtered contrasts"), NA) expect_error(DurgaPlot(d, ef.size.position = "below", contrasts = "Group3 - Group1", main = "2/3) ef.size below, 1 contrast"), NA) expect_error(DurgaPlot(d, main = "3/3) ef.size default position, shorthand contrasts", contrasts = ". - Group1"), NA) # Check missing contrast d <- DurgaDiff(df, data.col = 1, group.col = 2, contrasts = "Group2 - Group1") expect_error(DurgaPlot(d, contrasts = "Group3 - Group1")) }) test_that("Grouped plot", { # Attempt to reproduce Figure 1 in # Ostap-Chec, M., Opalek, M., Stec, D., & Miler, K. (2021). Discontinued alcohol consumption elicits withdrawal symptoms in honeybees. Biology Letters, 17(6), 20210182. doi:10.1098/rsbl.2021.0182 op <- par(mar = c(4, 5, 1, 10.5) + 0.1, mfrow = c(2, 1)) on.exit(par(op)) cols <- RColorBrewer::brewer.pal(4, "Dark2") set.seed(1) for (i in 1:2) { n <- 6 meansA <- c("water" = 5, "0.125%" = 5, "1.25%" = 5, "2%" = 12, "sugar" = 100) groups <- c("no exposure", "short exposure", "exposure withheld", "constant exposure") df <- data.frame(group = rep(groups, each = n), prob = sapply(names(meansA), function(nm) pmin(100, rnorm(n * length(groups), meansA[nm], meansA[nm] / 2))), check.names = FALSE) # Convert wide to long format rows <- lapply(names(meansA), function(nm) cbind.data.frame(group = df$group, treatment = rep(nm, nrow(df)), prob = df[[paste0("prob.", nm)]])) long <- as.data.frame(do.call(rbind, rows)) long$combined <- paste(long$group, long$treatment) cg <- sapply(names(meansA), function(t) paste(groups, t)) d <- DurgaDiff(long, "prob", "combined", groups = cg, contrasts = NULL) ps <- expect_error(DurgaPlot(d, x.axis = FALSE, left.las = 1, left.ylab = "probability of\nresponding (%)", frame.plot = FALSE, ef.size = FALSE, violin = FALSE, points = FALSE, group.dx = c(0.6, 0.2, -0.2, -0.6), central.tendency = cols, error.bars = cols), NA) # Force axis lines to cover the plot extents axis(1, at = par("usr")[1:2], labels = c("", ""), lwd.ticks = 0) axis(2, at = par("usr")[3:4], labels = c("", ""), lwd.ticks = 0) # Calculate centre of each group at <- colMeans(matrix(ps$extents[,1], nrow = 4)) axis(1, at = at, labels = names(meansA)) title(xlab = "solution") legend(par("usr")[2] + 0.2, mean(par("usr")[3:4]), yjust = 0.5, xpd = NA, c("1, no exposure", "2, short exposure", "3, exposure withheld", "4, constant exposure"), col = cols, lwd = 3, pch = 19, bty = "n") } }) test_that("group labels etc", { # Add in some fake sex data to the insulin data set data <- cbind(insulin, Sex = sample(c("Male", "Female"), nrow(insulin), replace = TRUE)) # Thin it the data so individual symbols are visible data <- data[data$id %in% 1:15, ] d <- DurgaDiff(data, "sugar", "treatment", "id", groups = c("Before insulin" = "before", "After insulin" = "after"), na.rm = TRUE) expect_error(DurgaPlot(d, left.ylab = "Blood sugar level", violin.shape = c("left", "right"), violin.dx = c(-0.055, 0.055), violin.width = 0.3, points = "black", points.params = list(bg = ifelse(data$Sex == "Female", "red", "blue"), pch = ifelse(data$Sex == "Female", 21, 24)), ef.size.pch = 19, ef.size.violin = "blue", ef.size.violin.shape = "full", central.tendency = FALSE, main = "Customised plot"), NA) d <- DurgaDiff(iris, data.col = "Sepal.Length", group.col = "Species") expect_error(DurgaPlot(d, bar = TRUE, error.bars.type = "SD", points = FALSE, main = "Bar chart with std. deviation"), NA) expect_error(DurgaPlot(d, box = TRUE, error.bars = TRUE, central.tendency.type = "median", error.bars.type = "CI", points = FALSE, main = "Box plot with 95% CI"), NA) expect_error(DurgaPlot(d, bar = TRUE, central.tendency.symbol = "segment", error.bars.type = "SE", points = FALSE, main = "Bar chart with SE"), NA) }) test_that("plot miscellanea", { d <- DurgaDiff(damselfly, "length", "maturity", groups = c("Immature" = "juvenile", "Mature" = "adult")) # Axis text is smaller when there are multiple columns op <- par(mfrow = c(1, 3)) on.exit(par(op)) DurgaPlot(d, ef.size.position = "below", main = "Text size consistent") par(mar = c(5, 4, 4, 6) + 0.1) DurgaPlot(d, bar = T) DurgaPlot(d, box = T, xlim = c(0, 5), ylim = c(28, 40), main = "Explicit limits") expect_equal(1, 1) expect_error(DurgaDiff(petunia, 1, 2, groups = c("self-fertilised" = "self_fertilised", "intercrossed" = "inter_cross", "Westerham-crossed" = "westerham_cross"), contrasts = c("Westerham-crossed - self-fertilised", "Westerham-crossed - intercrossed", "intercrossed - self-fertilised")), NA) }) test_that("CI", { x <- rnorm(200) CI90 <- meanCI(x, .9) CI95 <- meanCI(x, .95) expect_lt(CI95[1], mean(x)) expect_lt(CI95[1], CI90[1]) expect_gt(CI95[2], mean(x)) expect_gt(CI95[2], CI90[2]) set.seed(1) d90 <- DurgaDiff(petunia, "height", "group", ci.conf = .9) d95 <- DurgaDiff(petunia, "height", "group", ci.conf = .95) d99 <- DurgaDiff(petunia, "height", "group", ci.conf = .99) # CI should be smaller for lower levels for (g in seq_along(d90$groups)) { # Ensure we are comparing the same things expect_equal(d90$group.differences[[g]]$groupLabels, d95$group.differences[[g]]$groupLabels) expect_equal(d90$group.differences[[g]]$groupLabels, d99$group.differences[[g]]$groupLabels) # Check CI of mean differences # Lower limits expect_gt(d90$group.differences[[g]]$bca[4], d95$group.differences[[g]]$bca[4]) expect_gt(d95$group.differences[[g]]$bca[4], d99$group.differences[[g]]$bca[4]) # Upper limits expect_lt(d90$group.differences[[g]]$bca[5], d95$group.differences[[g]]$bca[5]) expect_lt(d95$group.differences[[g]]$bca[5], d99$group.differences[[g]]$bca[5]) # Check CI of group means expect_gt(d90$group.statistics[g, "CI.lower"], d95$group.statistics[g, "CI.lower"]) expect_gt(d90$group.statistics[g, "CI.lower"], d95$group.statistics[g, "CI.lower"]) expect_lt(d95$group.statistics[g, "CI.upper"], d99$group.statistics[g, "CI.upper"]) expect_lt(d95$group.statistics[g, "CI.upper"], d99$group.statistics[g, "CI.upper"]) } }) test_that("Diff error detection", { n <- 40 set.seed(1) realDiff <- 1 df <- data.frame(val = c(rnorm(n, mean = 10), rnorm(n, mean = 10 + realDiff)), group = c(rep("Control", n), rep("Group", n)), id = c(1:n, 1:n)) # Shuffle df <- df[sample(nrow(df)), ] # Missing group expect_error(DurgaDiff(df, data.col = 1, group.col = 2, groups = c("Control", "Group", "Nonexistent")), "Nonexistent") # Non factor ID d1 <- DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3) expect_true(d1$paired.data) # Factor ID df$id <- factor(df$id) d2 <- DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3) expect_true(d2$paired.data) compareDiffs(d2, d1) # Missing IDs - remove some of the Group values removeIdxs <- sample(which(df$group == "Group"), 10) df <- df[-removeIdxs, ] expect_error(DurgaDiff(df, data.col = 1, group.col = 2, id.col = 3), "id") # Non-numeric value expect_error(DurgaDiff(df, data.col = 3, group.col = 2, id.col = 3), "numeric") # Empty data frame expect_error(DurgaDiff(df[numeric(0), ], data.col = 1, group.col = 2, id.col = 3), "No data") }) test_that("single group", { op <- par(mfrow = c(1, 2)) on.exit(par(op)) n <- 40 df <- data.frame(val = rnorm(n, 10), group = rep("G1", each = n)) d <- DurgaDiff(df, data.col = 1, group.col = 2) p <- expect_error(DurgaPlot(d, main = "1 group in data"), NA) df <- data.frame(val = c(rnorm(n, 10), rnorm(n, 11)), group = rep(c("G1", "G2"), each = n)) d <- DurgaDiff(df, data.col = 1, group.col = 2, groups = "G1") expect_error(DurgaPlot(d, main = "1 group in diff"), NA) }) test_that("group colours", { op <- par(mfrow = c(2, 2)) on.exit(par(op)) d <- DurgaDiff(petunia, 1, 2) DurgaPlot(d, ef.size = FALSE, group.colour = "Set1", main = "Group colours Set1") DurgaPlot(d, ef.size = FALSE, group.colour = "blue", main = "group.colours blue") DurgaPlot(d, ef.size = FALSE, group.colour = c("red", "Green", "blue"), main = "group.colours RGB, points fill coloured", points = "black", points.params = list(pch = 21)) DurgaPlot(d, ef.size = FALSE, group.colour = c("red", "Green", "#0000ff"), main = "group.colours RGB, points red fill", points.params = list(pch = 21, bg = "red")) # Invalid colour/palette expect_error(DurgaPlot(d, group.colour = "Fred")) expect_error(DurgaPlot(d, group.colour = c("Set1", "Dark1"))) }) test_that("minor formatting", { d <- DurgaDiff(damselfly, 1, 3) expect_error(DurgaPlot(d, main = "Offset lines, styled error bars, ef size lines", error.bars = "red", error.bars.lwd = 5, error.bars.lty = 4, ef.size.mean.line.dx = 0.2, ef.size.line.col = "blue", ef.size.line.lwd = 3, ef.size.line.lty = 3, ef.size.violin = "blue", ef.size.violin.fill = "pink"), NA) expect_error(DurgaPlot(d, main = "Offset lines, styled error bars, ef size line", ef.size.position = "below", error.bars = "red", error.bars.lwd = 5, error.bars.lty = 4, ef.size.mean.line.dx = 0.2, ef.size.line.col = "blue", ef.size.line.lwd = 3, ef.size.line.lty = 1, ef.size.violin = "blue", ef.size.violin.fill = "pink"), NA) }) test_that("Spaces in column names", { n <- 40 df <- data.frame(`Genital length` = c(rnorm(n, mean = 10)), `Cannibalism y/n` = sample(c("Yes", "No"), n, replace = TRUE), check.names = FALSE) d <- DurgaDiff(df, "Genital length", "Cannibalism y/n") expect_error(DurgaPlot(d, xlab = "Cannibalism?", main = "Spaces in names, xlab"), NA) }) test_that("Formula", { n <- 40 set.seed(1) valVar <- c(rnorm(n, mean = 10), rnorm(n, mean = 11)) df <- data.frame(`Scapus length` = valVar, val = valVar, group = c(rep("Control", n), rep("Group", n)), id = c(1:n, 1:n), check.names = FALSE) `Group cat` <- df$group op <- par(mfrow = c(1, 2)) on.exit(par(op)) DurgaPlot(DurgaDiff(valVar ~ group, df, id.col = "id", groups = c("Control", "Group")), main = "Formula interface, paired") DurgaPlot(DurgaDiff(df, "val", "group", "id", groups = c("Control", "Group")), main = "Standard interface, paired") # More complicated formulae DurgaPlot(DurgaDiff(log(valVar) ~ group, df), main = "Formula interface") DurgaPlot(DurgaDiff(`Scapus length` ~ group, data = df), main = "Formula interface") DurgaPlot(DurgaDiff(log(`Scapus length`) ~ group, data = df), main = "Formula interface") DurgaPlot(DurgaDiff(log(`Scapus length`) ~ `Group cat`, data = df), main = "Formula interface") expect_error(DurgaDiff(`Scapus lengh` ~ group, data = df)) }) test_that("custom stat", { # Example from http://www.estimationstats.com/#/analyze/shared-control ds <- "A1 A2 B1 B2 C1 C2 D1 D2 8.885 10.135 8 -35 3.375 6.625 0.54 -0.54 14.38 11.94 7 -30 -0.3 2.3 1.98 0.02 8.015 6.025 17 -25 10.025 11.975 1.1 0.9 5.835 3.045 15 -20 2.35 3.65 3.42 0.58 5.47 1.87 12 -15 7.675 8.325 2.54 1.46 12.06 12.64 5 -10 9 9 1.655 2.345 11.72 9.66 6 -5 7.325 6.675 4.865 1.135 10.315 9.265 19 0 6.65 5.35 3.98 2.02 5.065 6.155 16 5 4.975 3.025 3.1 2.9 8.235 10.785 11 10 3.3 0.7 2.215 3.785 15.08 12.36 18 15 11.625 8.375 6.305 1.695 13.485 10.175 9 20 17.765 8.235 5.42 2.58 11.3 12.38 14 25 17.09 6.91 4.54 3.46 9.82 9.66 13 30 19.41 8.59 3.655 4.345 9.565 6.955 10 35 20.735 9.265 2.775 5.225" dw <- read.delim(text = ds) # Convert to long format l <- lapply(colnames(dw), function(c) data.frame(value = as.numeric(dw[[c]]), group = c, id = seq_along(c))) df <- do.call(rbind, l) # Single difference dd1 <- DurgaDiff(df, 1, 2, groups = c("A1", "A2")) dc1 <- DurgaDiff(df, 1, 2, groups = c("A1", "A2"), effect.type = function(x1, x2) { median(x2) - median(x1) }) # Many groups set.seed(1) dd <- DurgaDiff(df, 1, 2, contrasts = ". - A1") set.seed(1) dc <- DurgaDiff(df, 1, 2, contrasts = ". - A1", effect.type = function(x1, x2) { median(x2) - median(x1) }) # Check that group details are the same, pairwise differences are different od <- capture.output(print(dd)) oc <- capture.output(print(dc)) expect_equal(which(od == oc), 1:12) # Check it can be plotted # One contrast op <- par(mfrow = c(1, 2)) expect_error(DurgaPlot(dd1, main = "Defaults"), NA) expect_error(DurgaPlot(dc1, main = "Custom median differences"), NA) par(op) # Many contrasts op <- par(mfrow = c(2, 1)) expect_error(DurgaPlot(dd, main = "Defaults"), NA) expect_error(DurgaPlot(dc, main = "Custom median differences"), NA) par(op) }) test_that("custom labels", { data <- makeData() d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2")) # Override effect size label d$group.differences[[2]]$label.plot <- expression(italic("Sp. 1") ~ "-" ~ italic("Sp. 2")) # Override effect name label expect_error(DurgaPlot(d, ef.size.label = expression(bold("Bold name")), x.axis = F), NA) # Override print label d$group.differences[[2]]$label.print <- "Sp. 1 : Sp. 2" s <- capture.output(print(d)) expect_match(s[10], "^ *Sp. 1 : Sp. 2:") }) test_that("ef range plot bug", { n <- c(20, 30, 50, 30) seasons <- c("Spring", "Summer", "Autumn", "Winter") df <- data.frame(proportion = c( rnorm(n[1], mean = 0.1, sd = 0.01), rnorm(n[2], mean = 0.15, sd = 0.015), rnorm(n[3], mean = 0.25, sd = 0.02), rnorm(n[4], mean = 0.6, sd = 0.03) ), season = rep(seasons, n) ) d <- DurgaDiff(df, 1, 2, contrasts = ". - Winter", groups = seasons) expect_error(DurgaPlot(d, main = "Effect size ylim correct"), NA) }) test_that("ef below layout", { data <- makeData() d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2")) op <- par(mfrow = c(1, 2), cex = 0.7) on.exit(par(op)) expect_error(DurgaPlot(d, main = "Default EF layout"), NA) expect_error(DurgaPlot(d, ef.size.top.pad = 1.5, ef.size.height = 0.6, main = "Smaller gap, larger height"), NA) }) test_that("ef size ticks", { data <- makeData() op <- par(mfrow = c(2, 2), mar = c(5, 4, 4, 1) + 0.1) on.exit(par(op)) d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2")) DurgaPlot(d, main = "Custom ef ticks", ef.size.ticks = c(30, 0, -50)) DurgaPlot(d, contrasts = "Group2 - ZControl1", main = "Custom ef ticks", ef.size.ticks = c(30, 0, -50)) d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2"), effect.type = "cohens d") DurgaPlot(d, main = "Custom ef ticks", ef.size.ticks = c(-2, 0, 1)) expect_error(DurgaPlot(d, contrasts = "Group2 - ZControl1", main = "Custom ef ticks", ef.size.ticks = c(-2, 0, 1)), NA) }) test_that("ef size labels", { data <- makeData() op <- par(mfrow = c(2, 2)) on.exit(par(op)) d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2")) DurgaPlot(d, main = "Custom ef labels", ef.size.ticks = c("Big" = 30, "None" = 0, "Huge" = -50), ef.size.params = list(las = 1)) DurgaPlot(d, contrasts = "Group2 - ZControl1", main = "Custom ef labels", ef.size.ticks = c("Big" = 30, "None" = 0, "Huge" = -50), ef.size.params = list(las = 1)) d <- DurgaDiff(data, 1, 2, groups = c("ZControl1", "Group1", "Group2"), effect.type = "cohens d") DurgaPlot(d, main = "Custom ef labels", ef.size.ticks = c("Huge" = -2, "None" = 0, "Big" = 1), ef.size.params = list(las = 1)) expect_error(DurgaPlot(d, contrasts = "Group2 - ZControl1", main = "Custom ef labels", ef.size.ticks = c("Huge" = -2, "None" = 0, "Big" = 1), ef.size.params = list(las = 1)), NA) }) test_that("ef size symbology", { data <- makeData() op <- par(mfrow = c(1, 2)) on.exit(par(op)) d <- DurgaDiff(data, 1, 2, groups = c(Control = "ZControl1", "Group1", "Group2")) expect_error(DurgaPlot(d, main = "Custom ef symbology", ef.size = 2:3, ef.size.lty = 3:2, ef.size.lwd = c(4, 2)), NA) d <- DurgaDiff(data, 1, 2, groups = c(Control = "ZControl1", "Group1"), effect.type = "cohens d") expect_error(DurgaPlot(d, main = "Custom ef symbology", ef.size.lty = 2, ef.size.lwd = 4), NA) }) test_that("pathological data", { # There is a weird case when the density does not extend past the data n <- 10000 df <- data.frame(Value = c(rnorm(n), 500, rnorm(n), -500), group = rep(c("G1", "G2"), each = n + 1)) d <- DurgaDiff(df, 1, 2, contrasts = NULL, R = NA) # Don't bootstrap anything expect_error(DurgaPlot(d, main = "Pathological case - don't crash!"), NA) }) test_that("points layout", { # Make some bimodal data n1 <- 300 n2 <- 100 n3 <- 50 n <- n1 + n2 + n3 data <- data.frame(Measurement = c(rnorm(n1, mean = 1), rnorm(n2, mean = 4), rnorm(n3, mean = 7, sd = 0.2), rnorm(n1, mean = 4), rnorm(n2, mean = 0.5), rnorm(n3, mean = -1, sd = 0.5)), Group = rep(c("Group 1", "Group 2"), each = n)) op <- par(mfrow = c(2, 2)) on.exit(par(op)) d <- DurgaDiff(data, 1, 2, contrasts = "") expect_error(DurgaPlot(d, main = "Default point layout", ef.size = FALSE), NA) expect_error(DurgaPlot(d, main = "Adjust = 0.7", ef.size = FALSE, points.adjust = 0.7, violin.adj = 0.7), NA) expect_error(DurgaPlot(d, main = "Method = tukey", ef.size = FALSE, points.method = "tukey"), NA) expect_error(DurgaPlot(d, main = "Method = overplot", ef.size = FALSE, points.method = "overplot"), NA) }) test_that("missing right axis", { # Bug: Can fail with very small group sizes, because bootstrap wasn't stratified stats <- structure(list(cat = c("RelearningwalksOld", "RelearningwalksOld", "RelearningwalksOld", "RelearningwalksOld", "RelearningwalksOld", "RelearningwalksOld", "RelearningNew", "RelearningNew", "RelearningNew", "RelearningNew", "RelearningNew"), `Convex hull area` = c(5749.50993182344, 12052.2167954142, 11692.8373439029, 5327.12611957551, 434.590489399867, 19118.6301846389, 7934.00774651262, 28646.8325343417, 25339.7006681887, 11880.8673679235, 26535.4117772832), `Max displacement` = c(8.41844123108767, 6.34767201343218, 25.188225965698, 6.20186762835405, 2.10714057756224, 5.98568895263556, 11.3602307617808, 15.453027281305, 11.3599509921636, 18.5884235814962, 21.898981986908), Duration = c(25.02, 43.14, 25.6, 44.76, 13.1, 55.48, 8.28, 9.52, 22, 9.28, 11.96)), row.names = c(NA, -11L), class = "data.frame") learningSheets <- c(`Old nest relearning walks` = "RelearningwalksOld", `New nest relearning walks` = "RelearningNew") set.seed(1) d <- DurgaDiff(stats, "Convex hull area", "cat", groups = learningSheets) expect_error(DurgaPlot(d, main = "Effect size axis present?"), NA) data <- data.frame(Val = rnorm(10), Group = c(1, 1, rep(2, 8))) expect_error(DurgaDiff(data, 1, 2), NA) }) test_that("Wide format data", { n <- 10 wide <- data.frame(Control = rnorm(n, 10), Treatment = rnorm(n, 11), Mass = rnorm(n, 2, .4)) l <- stats::reshape(wide, direction = "long", varying = list(c("Control", "Treatment")), idvar = "Specimen", v.names = "Measurement", # Name of the value column timevar = "Condition", # Name of the group column times = c("Control", "Treatment")) d1 <- DurgaDiff(l, "Measurement", "Condition") d2 <- DurgaDiff(wide, id.col = "Specimen", groups = c("Control", "Treatment")) par(mfrow = c(2, 2)) DurgaPlot(d1, main = "Manual wide to long") expect_error(DurgaPlot(d2, main = "Auto wide to long"), NA) # 3 groups n <- 20 df3 <- data.frame(Control = rnorm(n, 10), Treatment1 = rnorm(n, 11), Treatment2 = rnorm(n, 11.1, .4), Treatment3 = rnorm(n, 9, 2)) d <- DurgaDiff(df3, id.col = "Id", groups = c("Control", "Treatment1", "Treatment2", "Treatment3")) DurgaPlot(d, main = "Paired - multiple groups", left.ylab = "Mass (g)") DurgaPlot(d, main = "Paired - multiple groups, contrasts", left.ylab = "Mass (g)", contrasts = "Treatment1 - Control, Treatment2 - Treatment1, Treatment3 - Treatment2") # Id column that isn't unique df <- data.frame(id = rep(1:5, each = 2), Control = rnorm(10, 10), Treatment = rnorm(10, 11)) expect_error(DurgaDiff(df, id.col = "id", groups = c("Control", "Treatment"))) # Id column not specified df <- data.frame(Control = rnorm(10, 10), Treatment = rnorm(10, 11)) d <- DurgaDiff(df, groups = c("Control", "Treatment")) expect_true(d$paired.data) # Existing id column not specified df <- data.frame(id = sample(10), Control = rnorm(10, 10), Treatment = rnorm(10, 11)) d <- DurgaDiff(df, groups = c("Control", "Treatment")) expect_true(d$paired.data) }) test_that("mapY fn", { my <- buildMapYFn(c(0, 1), c(0, 1)) expect_equal(my(0), 0) expect_equal(my(1), 1) expect_equal(my(2), 2) expect_equal(my(-2), -2) my <- buildMapYFn(c(1, 2), c(2, 1)) expect_equal(my(1), 2) expect_equal(my(2), 1) expect_equal(my(3), 0) my <- buildMapYFn(c(1, 0), c(1, 0)) expect_equal(my(0), 0) expect_equal(my(1), 1) expect_equal(my(2), 2) expect_equal(my(-2), -2) }) test_that("group means", { # Should a 2-sample group generate an error because the group mean CI can't be estimated? # No, because sometimes people have 2 sample groups and can't get more data df <- data.frame(val = c(rnorm(2), rnorm(4, 1)), group = c("A", "A", "B", "B", "B", "B")) d <- expect_error(DurgaDiff(df, "val", "group"), NA) expect_true(is.na(d$group.statistics[1, "CI.lower"])) expect_true(is.na(d$group.statistics[1, "CI.upper"])) expect_true(!is.na(d$group.statistics[2, "CI.lower"])) expect_true(!is.na(d$group.statistics[2, "CI.upper"])) })