### test-BuyseTest-previousBug.R --- ##---------------------------------------------------------------------- ## Author: Brice Ozenne ## Created: apr 17 2018 (16:46) ## Version: ## Last-Updated: sep 23 2024 (18:12) ## By: Brice Ozenne ## Update #: 244 ##---------------------------------------------------------------------- ## ### Commentary: ## ### Change Log: ##---------------------------------------------------------------------- ## ### Code: if(FALSE){ library(testthat) library(BuyseTest) library(data.table) library(survival) library(prodlim) } context("Check that bugs that have been reported are fixed \n") ## * settings BuyseTest.options(check = TRUE, keep.pairScore = TRUE, method.inference = "none", pool.strata = "Buyse", add.1.presample = FALSE, trace = 0) ## * Joris: jeudi 5 avril 2018 à 14:57 dt.sim <- data.table( ttt=c(rep(0,3),rep(1,3)), timeOS = c(10,20,30,15,20,35), eventOS = c(1,1,0,0,1,1), Mgrade.tox = -c(1,2,3,2,4,2) ) test_that("number of pairs - argument neutral.as.uninf", { for(iCorrection in c(FALSE,TRUE)){ ## iCorrection <- TRUE ; iCorrection <- FALSE BT.T <- BuyseTest(ttt~TTE(timeOS,threshold=0,status=eventOS) + cont(Mgrade.tox,threshold=0), data = dt.sim, neutral.as.uninf = TRUE, scoring.rule = "Gehan", correction.uninf = iCorrection) BTS.T <- as.data.table(model.tables(BT.T, percentage = FALSE)) BT.F <- BuyseTest(ttt~TTE(timeOS,threshold=0,status=eventOS) + cont(Mgrade.tox,threshold=0), data = dt.sim, neutral.as.uninf = FALSE, scoring.rule = "Gehan", correction.uninf = iCorrection) BTS.F <- as.data.table(model.tables(BT.F, percentage = FALSE)) ## neutral.as.uninf does not impact the results for first endpoint expect_equal(BTS.T[1,c("favorable","unfavorable","neutral","uninf","delta","Delta")], BTS.F[1,c("favorable","unfavorable","neutral","uninf","delta","Delta")]) ## check consistency of the number of pairs ## neutral.as.uninf = TRUE ## summary(BT.T) expect_equal(BTS.T[endpoint == "Mgrade.tox", favorable+unfavorable+neutral+uninf], BTS.T[endpoint == "Mgrade.tox", total]) expect_equal(BTS.T[endpoint == "timeOS", neutral+uninf], BTS.T[endpoint == "Mgrade.tox", total]) expect_equal(BTS.T[endpoint == "Mgrade.tox", total], BTS.T[endpoint == "Mgrade.tox", favorable+unfavorable+neutral+uninf]) ## neutral.as.uninf = FALSE expect_equal(BTS.F[endpoint == "Mgrade.tox", favorable+unfavorable+neutral+uninf], BTS.F[endpoint == "Mgrade.tox", total]) expect_equal(BTS.F[endpoint == "timeOS", uninf], BTS.F[endpoint == "Mgrade.tox", total]) expect_equal(BTS.F[endpoint == "Mgrade.tox", total], BTS.F[endpoint == "Mgrade.tox", favorable+unfavorable+neutral+uninf]) ## compared to known value if(iCorrection == FALSE){ keep.col <- c("endpoint","total","favorable","unfavorable","neutral","uninf","delta","Delta") test <- as.data.table(model.tables(BT.T, percentage = TRUE)[,keep.col]) GS <- data.table("endpoint" = c("timeOS", "Mgrade.tox"), "total" = c(100.00000, 44.44444), "favorable" = c(44.44444, 22.22222), "unfavorable" = c(11.11111, 11.11111), "neutral" = c(11.11111, 11.11111), "uninf" = c(33.33333, 0.00000), "delta" = c(0.3333333, 0.1111111), "Delta" = c(0.3333333, 0.4444444)) ## butils::object2script(test) attr(test,"index") <- NULL expect_equal(test, GS, tol = 1e-6) ## class(BTS.T[["n.resampling"]]) ## class(GS[["n.resampling"]]) } } }) ## * Emeline T: samedi 26 mai 2018 à 14:39 (Version 1.0) ## ERROR: Error in xy.coords(x, y, setLab = FALSE) : 'x' and 'y' lengths differ ## butils:::object2script(data[175:325,], digits = 8) data <- data.frame("X" = c(175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325), "trt" = c(1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0), "time" = c(0.34972271, 0.15919528, 0.27021802, 0.95994001, 0.46312769, 0.31997409, 0.23637537, 0.89707932, 0.01708799, 0.39366592, 0.50017673, 0.446804, 0.50040844, 0.32638758, 0.6262995, 0.14755089, 0.12050248, 0.07458989, 0.04339593, 0.59982912, 1.41101136, 0.2675838, 0.52521586, 1.14631933, 0.74191272, 1.2196955, 0.02732667, 1.3869635, 0.75430971, 0.3780356, 0.42434206, 1.28254783, 0.65964535, 0.80568326, 1.09058069, 0.14099648, 1.30095204, 0.69223441, 1.33892841, 0.73062582, 0.28980283, 1.74724314, 0.85952631, 0.40828457, 1.26493484, 0.96396552, 0.75849828, 0.70308743, 1.71091642, 1.01266995, 0.29350899, 0.79999462, 0.90685983, 0.2697463, 0.92647206, 0.00936012, 0.69425291, 0.82894713, 0.28051478, 1.40047767, 0.83924557, 0.61605441, 0.56216195, 0.68796769, 1.83362936, 0.45955409, 1.381266, 1.34455702, 0.30326241, 2.42955884, 0.53467431, 1.00931952, 1.11490004, 0.72048666, 0.07125682, 0.34582823, 0.33357166, 0.47453535, 0.27259304, 0.60673207, 0.95520791, 0.05198433, 0.82662585, 1.15532297, 0.87506277, 1.37889663, 0.12846039, 0.68540728, 0.77377909, 0.81177511, 0.29095231, 2.02666276, 0.21531326, 0.45024274, 1.43151175, 0.46492612, 0.14985886, 0.22205914, 1.59582145, 0.76701798, 1.23825982, 0.33712561, 1.07747869, 0.06973708, 1.27342747, 0.42610371, 1.0686674, 2.03964558, 0.5787245, 1.05125486, 0.24393524, 1.02678662, 0.2725943, 0.59435986, 0.32627314, 0.39337226, 0.71167895, 0.58597973, 0.3605633, 1.24886565, 0.43183396, 0.75826836, 0.22063575, 0.28832416, 0.16407274, 0.91388552, 0.62053192, 2.46164696, 0.28193246, 0.33575549, 0.51327929, 0.90610562, 0.43071919, 1.392834, 0.69855789, 0.81717857, 0.46312768, 0.11466708, 0.42909682, 0.29334352, 0.76480274, 0.80197241, 0.40497033, 0.68113025, 0.98833506, 0.58629864, 0.00627822, 0.35254414, 0.52416901, 0.67108879, 0.49179438), "event" = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), stringsAsFactors = FALSE) BT_tau0 <- BuyseTest(data=data, treatment="trt", endpoint="time", type="timeToEvent", threshold=as.numeric(0), status="event", scoring.rule="Peron", method.inference = "none", cpus=1, trace = 0) ## * Brice: 09/06/18 6:51 (Tied event with tte endpoint) ## when computing the integral for peron with double censoring ## the ordering of the data modified the ouput ## this has been correct with version 1.4 data(cancer, package = "survival") test_that("ordering of tied event does not affect BuyseTest", { ## veteran2[veteran2$time==100,] BT.all <- BuyseTest(trt ~ tte(time, threshold = 0, status = "status"), data = veteran, scoring.rule = "Peron", method.inference = "none", correction.uninf = FALSE) veteran1 <- veteran[order(veteran$time,veteran$status),c("time","status","trt")] ## veteran1[veteran2$time==100,] BT1.all <- BuyseTest(trt ~ tte(time, threshold = 0, status = "status"), data = veteran1, scoring.rule = "Peron", method.inference = "none", correction.uninf = FALSE) veteran2 <- veteran[order(veteran$time,-veteran$status),c("time","status","trt")] ## ## veteran2[veteran2$time==100,] BT2.all <- BuyseTest(trt ~ tte(time, threshold = 0, status = "status"), data = veteran2, scoring.rule = "Peron", method.inference = "none", correction.uninf = FALSE) ## effect of the ordering expect_equal(coef(BT.all, statistic = "winRatio"), coef(BT1.all, statistic = "winRatio")) expect_equal(coef(BT.all, statistic = "winRatio"), coef(BT2.all, statistic = "winRatio")) ## number of pairs expect_equal(as.double(BT.all@n.pairs), prod(table(veteran$trt)), tol = 1e-5) expect_equal(as.double(BT1.all@n.pairs), prod(table(veteran$trt)), tol = 1e-5) expect_equal(as.double(BT2.all@n.pairs), prod(table(veteran$trt)), tol = 1e-5) ## values of the pairs expect_true(all(getPairScore(BT.all, endpoint = 1)[["favorable"]]>=0)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["favorable"]]<=1)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["unfavorable"]]>=0)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["unfavorable"]]<=1)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["neutral"]]>=0)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["neutral"]]<=1)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["uninformative"]]>=0)) expect_true(all(getPairScore(BT.all, endpoint = 1)[["uninformative"]]<=1)) expect_true(all(getPairScore(BT.all, endpoint = 1)[,favorable + unfavorable]<=1+1e-12)) ## tolerance ## survival ## getSurvival(BT.all, endpoint = 1, strata = 1)$lastSurv: only 0 so no uninformative paris expect_equal(as.double(coef(BT.all, statistic = "count.uninf", cumulative = FALSE)), 0.0, tol = 1e-12) ## result expect_equal(as.double(coef(BT.all, statistic = "winRatio")), 0.8384569, tol = 1e-5) }) ## * Brice: 26/09/18 x:xx (Multiple thresholds in Julien's simulations) HR1 <- 0.65 TpsFin <- 60 #values for Taux.Censure HazC <- 0.1 set.seed(10) HazT <- 0.1*(HR1) n.Treatment <- 100 n.Control <- 100 n <- n.Treatment+n.Control group <- c(rep(1, n.Treatment),rep(0, n.Control)) TimeEvent.Ctr <- rexp(n.Control,HazC) TimeEvent.Tr <- rexp(n.Control,HazT) TimeEvent<-c(TimeEvent.Tr,TimeEvent.Ctr) Time.Cens<-runif(n,0,TpsFin) Time<-pmin(Time.Cens,TimeEvent) Event<-Time==TimeEvent Event<-as.numeric(Event) tab <- data.frame(group,Time,Event, stringsAsFactors = FALSE) test_that("Multiple thresholds",{ BuyseresPer <- BuyseTest(data=tab, endpoint=c("Time","Time","Time","Time","Time","Time","Time","Time","Time","Time","Time","Time","Time","Time","Time"), treatment="group", type=c("TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE","TTE"), status=c("Event","Event","Event","Event","Event","Event","Event","Event","Event","Event","Event","Event","Event","Event","Event"), threshold=c(42,39,36,33,30,27,24,21,18,15,12,9,6,3,0), n.resampling=500, trace=0, scoring.rule="Peron", correction.uninf=F, method.inference="none") resS <- as.data.table(model.tables(BuyseresPer)) ## pairs are correctly transfered from one endpoint to another expect_equal(resS[threshold > tail(threshold,1), neutral + uninf], resS[threshold < threshold[1], total], tol = 1e-2) ## butils::object2script(as.double(BuyseresPer@count.favorable), digit = 2) GS <- c(260.64, 35.93, 37.33, 147.32, 272.14, 263.6, 235.7, 213.21, 390.29, 408.73, 514.7, 514.34, 744.78, 865.21, 1095.26) expect_equal(as.double(coef(BuyseresPer, statistic = "count.favorable", cumulative = FALSE)), GS, tol = 1e-5) ## butils::object2script(as.double(BuyseresPer@count.unfavorable), digit = 2) GS <- c(0, 0, 6.97, 25.66, 43.89, 34.8, 46.38, 105.42, 199.85, 338.55, 407.72, 521.83, 548.02, 782.94, 938.8) expect_equal(as.double(coef(BuyseresPer, statistic = "count.unfavorable", cumulative = FALSE)), GS, tol = 1e-5) ## butils::object2script(as.double(BuyseresPer@count.neutral), digit = 2) GS <- c(9580.24173298, 9580.24173298, 9573.26856493, 9433.45578036, 9140.84153639, 8860.00546388, 8577.9282735, 8259.29654759, 7675.01664374, 6933.58435144, 6011.16786378, 4975.00117239, 3682.20846492, 2034.06256111, 0) expect_equal(as.double(coef(BuyseresPer, statistic = "count.neutral", cumulative = FALSE)), GS, tol = 1e-5) ## butils::object2script(as.double(BuyseresPer@count.uninf), digit = 2) GS <- c(159.12095011, 123.19041299, 85.85998481, 52.68680886, 29.27044937, 11.70817975, 11.70817975, 11.70817975, 5.85408987, 0, 0, 0, 0, 0, 0) expect_equal(as.double(coef(BuyseresPer, statistic = "count.uninf", cumulative = FALSE)), GS, tol = 1e-1) ## butils::object2script(as.double(BuyseresPer@delta.netBenefit), digit = 5) GS <- c(0.02606, 0.00359, 0.00304, 0.01217, 0.02282, 0.02288, 0.01893, 0.01078, 0.01904, 0.00702, 0.0107, -0.00075, 0.01968, 0.00823, 0.01565) expect_equal(as.double(coef(BuyseresPer, statistic = "netBenefit", cumulative = FALSE)), GS, tol = 1e-3) ## butils::object2script(as.double(BuyseresPer@delta.winRatio), digit = 5) GS <- c(Inf, Inf, 5.35344, 5.74093, 6.19986, 7.57457, 5.08161, 2.02241, 1.95291, 1.2073, 1.2624, 0.98564, 1.35904, 1.10508, 1.16666) expect_equal(as.double(coef(BuyseresPer, statistic = "winRatio", cumulative = FALSE)), GS, tol = 1e-3) }) ## * Brice: 30/10/18 4:36 Neutral pairs with 0 threshold df <- data.frame("survie" = c(2.1, 4.1, 6.1, 8.1, 4, 6, 8, 10), "event" = c(1, 1, 1, 0, 1, 0, 0, 1), "group" = c(0, 0, 0, 0, 1, 1, 1, 1), "score" = 1, stringsAsFactors = FALSE) test_that("1 TTE endpoint - Gehan (no correction)", { Peron <- BuyseTest(group ~ tte(survie, status = event, threshold = 0), data = df, scoring.rule = "Peron", correction.uninf = FALSE) expect_equal(as.double(coef(Peron, statistic = "count.neutral", cumulative = FALSE)),0) ## should not be any neutral pair with a threshold of 0 }) ## * Hickey, Graeme: 8 mars 2019 14:54 p-value permutation ## I have one question, which I hope you can help with. ## If using method.inference = “permutation”, the P-values are slightly different for the net benefit and win ratio summary methods. ## However, if you use using method.inference = “bootstrap”, the P-values are identical, as I would expect. ## Can you explain why they differ with the permutation test? set.seed(1) dt <- simBuyseTest(50) test_that("same p.value (permutation test) for winRatio and net Benefit", { e.perm <- BuyseTest(treatment ~ bin(toxicity), data = dt, method.inference = "permutation", n.resampling = 100, trace = 0) netBenefit.perm <- suppressWarnings(confint(e.perm, statistic = "netBenefit")) winRatio.perm <- suppressWarnings(confint(e.perm, statistic = "winRatio")) Delta.netBenefit <- coef(e.perm, statistic = "netBenefit") Delta.winRatio <- coef(e.perm, statistic = "winRatio") DeltaResampling.netBenefit <- e.perm@DeltaResampling[,1,"netBenefit"] DeltaResampling.winRatio <- e.perm@DeltaResampling[,1,"winRatio"] manual <- c(netBenefit = mean(abs(DeltaResampling.netBenefit) >= abs(Delta.netBenefit)), netBenefit.atanh = mean(abs(atanh(DeltaResampling.netBenefit)) >= abs(atanh(Delta.netBenefit))), winRatio = mean(abs(DeltaResampling.winRatio-1) >= abs(Delta.winRatio-1)), winRatio.log = mean(abs(log(DeltaResampling.winRatio)) >= abs(log(Delta.winRatio))) ) expect_equal(netBenefit.perm[,"p.value"], winRatio.perm[,"p.value"]) expect_equal(unname(manual["netBenefit"]), netBenefit.perm[,"p.value"]) expect_equal(unname(manual["winRatio.log"]), winRatio.perm[,"p.value"]) ## note CI are not agreeing with p-values suppressWarnings(confint(e.perm, statistic = "netBenefit", conf.level = 1-0.48)) suppressWarnings(confint(e.perm, statistic = "winRatio", conf.level = 1-0.48)) }) ## * Alice, Brouquet-Laglair: 3 avril 2019 p-value bootstrap df <- rbind(data.frame(score = rep(1,5), tox = 0, group = 1, stringsAsFactors = FALSE), data.frame(score = rep(0,5), tox = 0, group = 0, stringsAsFactors = FALSE) ) test_that("BuyseTest without variability", { e.BT_ustat <- BuyseTest(group ~ bin(tox) + cont(score), data = df, method.inference = "u-statistic", trace = 0) e.BT_boot <- BuyseTest(group ~ bin(tox) + cont(score), data = df, method.inference = "studentized bootstrap", n.resampling = 10, trace = 0) confintTempo <- confint(e.BT_ustat) expect_equal(unname(confintTempo[,"p.value"]),1:0) confintTempo <- suppressMessages(confint(e.BT_boot, transformation = FALSE)) expect_equal(unname(confintTempo[,"p.value"]),1:0) }) ## * graemeleehickey (issue #2 on Github): 8 september 2019 p-value bootstrap test_that("Boostrap - issue in the summary", { BT.keep <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status") + cont(karno), data = veteran, keep.pairScore = TRUE, scoring.rule = "Gehan", trace = 0, method.inference = "bootstrap", n.resampling = 20, seed = 10) expect_error(capture.output(summary(BT.keep, statistic = "winRatio")), regexp = NA) ## no error }) ## * graemeleehickey (issue #3 on Github): 22 september 2019 BuysePower test_that("BuysePower - error in print", { simFCT <- function(n.C, n.T){ out <- data.table(Y=rnorm(n.C+n.T), T=c(rep(1,n.C),rep(0,n.T)) ) return(out) } ## the error was when setting trace to 4 tempo <- capture.output({ xx <- powerBuyseTest(sim = simFCT, sample.size = 100, n.rep = 2, formula = T ~ cont(Y), method.inference = "u-statistic", trace = 4, seed = 10) }) yy <- powerBuyseTest(sim = function(n.C, n.T){ out <- data.table(Y=rnorm(n.C+n.T), T=c(rep(1,n.C),rep(0,n.T)) ) return(out) }, sample.size = 100, n.rep = 2, formula = T ~ cont(Y), method.inference = "u-statistic", trace = 0, seed = 10) expect_equal(xx,yy) ## xx <- powerBuyseTest(sim = simFCT, ## sample.sizeC = c(100), ## sample.sizeT = c(100), ## n.rep = 10, ## cpus = 3, ## formula = T ~ cont(Y), ## method.inference = "u-statistic", ## trace = 4) }) ## * brice ozenne: 11/13/19 4:11 hierachical in BuyseTest test_that("BuyseTest - hierarchical", { BuyseTest.options(order.Hprojection = 1) BT.nH1 <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status", weight = 1) + cont(karno, threshold = 0, weight = 1), hierarchical = FALSE, data = veteran, method.inference = "u-statistic", trace = 0) BT.nH.5 <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status") + cont(karno, threshold = 0), hierarchical = FALSE, data = veteran, method.inference = "u-statistic", trace = 0) expect_equal(unname(coef(BT.nH.5)), c(-0.04382918, -0.05949414), tol = 1e-6) expect_equal(coef(BT.nH.5)*2, coef(BT.nH1), tol = 1e-6) expect_equal(BT.nH1@covariance[,"netBenefit"],4*BT.nH.5@covariance[,"netBenefit"], tol = 1e-6) expect_equal(as.double(confint(BT.nH.5)$se), c(0.04880450,0.08700807), tol = 1e-6) expect_equal(confint(BT.nH.5)$se*2, confint(BT.nH1)$se, tol = 1e-6) }) ## * graemeleehickey (issue #4 on Github): 6 october 2019 simBuyseTest test_that("simBuyseTest - rate vs. scale", { scale <- 2 args <- list(scale.T = scale, scale.censoring.T = scale+1, scale.C = scale, scale.censoring.C = scale+1, scale.CR = scale) set.seed(10) test <- simBuyseTest(1e4, argsBin = NULL, argsCont = NULL, argsTTE = args, latent = TRUE) expect_equal(scale,mean(test[treatment == "C", mean(eventtimeUncensored)]), tol = 1e-2) expect_equal(scale,mean(test[treatment == "T", mean(eventtimeUncensored)]), tol = 1e-2) expect_equal(scale+1,mean(test[treatment == "C", mean(eventtimeCensoring)]), tol = 1e-1) expect_equal(scale+1,mean(test[treatment == "T", mean(eventtimeCensoring)]), tol = 1e-1) }) ## * graemeleehickey (issue #6 on Github): 15 march 2020 powerBuyseTest args <- list(scale.T = c((3:5) / 10), scale.censoring.T = rep(1, 3)) simFCT <- function(n.C, n.T) { simBuyseTest(100, argsBin = NULL, argsCont = NULL, argsTTE = args) } test_that("powerBuyseTest - status vs. censoring", { valid <- powerBuyseTest(sim = simFCT, sample.size = 100, n.rep = 2, formula = treatment ~ tte(eventtime1, status = status1), method.inference = "u-statistic", scoring.rule = "Gehan", trace = 0) expect_error(powerBuyseTest(sim = simFCT, sample.size = 100, n.rep = 2, formula = treatment ~ tte(eventtime1, censoring = status1), method.inference = "u-statistic", scoring.rule = "Gehan")) valid <- capture.output(powerBuyseTest(sim = simFCT, sample.size = 100, n.rep = 2, formula = treatment ~ tte(eventtime1, status = status1), method.inference = "u-statistic", scoring.rule = "Gehan", trace = 4)) }) ## * brice ozenne : 04/26/20 2:36 uninformative pairs Peron dt.prodlim <- rbind(data.table(treat=0, time = c(1:8,rep(9,12)), status = c(rep(1,8),rep(0,12)) ), data.table(treat=1, time = c(1:8,rep(9,12)), status = c(0,rep(1,7),rep(0,12)) )) e.prodlim <- prodlim(Hist(time, status) ~ treat, data = dt.prodlim) ## plot(e.prodlim) dt.sim <- data.table(treat = c(0:1), time = 8, status = 0) test_that("uniformative pair after last observation",{ ## warning because only uninformative expect_warning(e.BP <- BuyseTest(treat ~ tte(time, status, threshold=2), model.tte = e.prodlim, data = dt.sim, method.inference = "none")) expect_equal(as.double(e.BP@count.neutral), 0) expect_equal(as.double(e.BP@count.uninf), 1) }) ## * brice ozenne : 10/08/20 3:26 last time tie (event/censor) ## butils::object2script(mydata, digit = 3) dt <- data.table("bras" = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), "OS" = c(0.427, 1.708, 2.004, 2.792, 3.088, 3.384, 3.417, 3.647, 3.778, 3.844, 5.092, 5.355, 5.453, 6.012, 6.209, 6.209, 6.307, 6.702, 7.786, 8.049, 8.739, 9.461, 11.367, 11.728, 11.925, 11.991, 12.648, 12.746, 13.042, 13.338, 13.436, 13.666, 13.798, 16.097, 16.097, 0.854, 1.84, 3.055, 3.515, 4.172, 5.059, 5.158, 5.223, 5.519, 5.585, 6.307, 6.34, 6.373, 6.767, 6.899, 6.965, 7.129, 7.589, 7.589, 7.589, 7.753, 8.18, 9.133, 9.198, 9.855, 10.315, 11.498, 13.141, 13.239, 13.305, 13.568, 13.929, 15.21, 16.459, 19.087, 20.237, 20.532, 21.846, 22.273, 26.445, 27.989), "etat" = c(0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1)) test_that("last time is a tie with both event and censor",{ test <- BuyseTest(bras ~ tte(OS, status = etat), data = dt, method.inference = "u-statistic", scoring.rule = "Peron", trace = 0) expect_equal(as.double(c(coef(test, statistic = "count.favorable"),coef(test, statistic = "count.unfavorable"),coef(test, statistic = "count.neutral"))), c(892.6111, 520.3092, 0.0000 ), tol = 1e-3) ## dt[c(1,36)] }) ## * brice ozenne : 10/12/20 9:46 only censored event in one group test_that("one group with only censoring, one group with no censoring",{ dt <- data.table("treatment" = c(rep("C",10),rep("T",10)), "time" = c(1:10,1:10), "status" = c(rep(1,10),rep(0,10))) e.Peron <- BuyseTest(treatment ~ tte(time, status = status, threshold = 0), data = dt, scoring.rule = "Peron") expect_equal(as.double(coef(e.Peron, statistic = "netBenefit")),1) dt2 <- data.table("treatment" = c(rep("C",10),rep("T",10)), "time" = c(1:10,1:10), "status" = c(c(rep(1,9),0),rep(0,10))) e2.Peron <- BuyseTest(treatment ~ tte(time, status = status, threshold = 0), data = dt2, scoring.rule = "Peron") expect_equal(as.double(coef(e2.Peron, statistic = "netBenefit")),0.9) dt3 <- data.table("treatment" = c("C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"), "time" = c(0.302, 0.307, 0.336, 0.347, 0.348, 0.459, 0.494, 0.525, 0.587, 0.588, 0.098, 0.116, 0.180, 0.229, 0.306, 0.318, 0.452, 0.485, 1.025, 1.339), "status" = c(0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) e3.Peron <- BuyseTest(treatment ~ tte(time, status = status, threshold = 0), data = dt3, scoring.rule = "Peron") expect_equal(as.double(coef(e3.Peron, statistic = "netBenefit")),0.733333333) }) ## * brice ozenne : 02/18/21 12:00 subset factor strata test_that("subset factor strata",{ dt <- data.table("treatment" = c(rep("C",100),rep("T",100)), "time" = rnorm(200, mean = 100), "status" = 1, "strata" = factor(1:5)) dtR <- dt[strata %in% 1:3] dtR[, strata := droplevels(strata)] test <- BuyseTest(treatment ~ tte(time, status = status, threshold = 0) + strata, data = dt[strata %in% 1:3], trace = FALSE) GS <- BuyseTest(treatment ~ tte(time, status = status, threshold = 0) + strata, data = dtR, trace = FALSE) expect_equal(confint(test),confint(GS), tol = 1e-6) }) ## * Mickaël De Backer : Oct 11, 2021 11:04:31 transformation and variance strata test_that("backtransformation after permutation",{ gpc_ex1 <- BuyseTest(trt~cont(time), data=veteran, seed = 10, n.resampling = 100, method.inference = "permutation") test <- suppressWarnings(confint(gpc_ex1, statistic='winratio', conf.level = 0.95) ) expect_true(test$estimate < test$upper.ci) expect_true(test$estimate > test$lower.ci) }) test_that("U-stat with stratification",{ for(iOrder in 1:2){ ## iOrder <- 2 BuyseTest.options(order.Hprojection = iOrder) GPC.stratified <- BuyseTest(trt~cont(time) + celltype, data=veteran, method.inference = "u-statistic") ls.GPC <- lapply(split(veteran, veteran$celltype), function(iData){ BuyseTest(trt~cont(time), data=iData, method.inference = "u-statistic") }) weights.strata <- GPC.stratified@n.pairs/sum(GPC.stratified@n.pairs) expect_equivalent(sum(weights.strata*coef(GPC.stratified,statistic="netBenefit",strata = TRUE)), coef(GPC.stratified,statistic="netBenefit"), tol = 1e-6) expect_equivalent(sum(weights.strata^2*sapply(ls.GPC,function(iGPC){iGPC@covariance[,"netBenefit"]})), GPC.stratified@covariance[,"netBenefit"], tol = 1e-6) } }) ## * SamSalvaggio (issue #9 on Github): 21 december 2021 permutation test_that("p-value with permutation",{ set.seed(1) dt <- simBuyseTest(50, argsCont = list(mu.T = 100, mu.C = 1)) ## extremely large difference so always in favor of treatment GPC.perm <- BuyseTest(treatment~cont(score), data=dt, trace = FALSE, method.inference = "permutation") BuyseTest.options(add.1.presample = FALSE) expect_equal(suppressWarnings(confint(GPC.perm)$p.value), 0) BuyseTest.options(add.1.presample = TRUE) expect_equal(suppressWarnings(confint(GPC.perm)$p.value), 1/1001) }) ## * SamSalvaggio (issue #10 on Github): 6 june 2023 restricted test_that("restriction via formula interface",{ set.seed(1) dt <- simBuyseTest(n.T = 50, n.C = 50, names.strata = "strat_column", n.strata = 6, argsTTE = list(name = c("tte1","tte2","tte3", "tte4", "tte5", "tte6", "tte7", "tte8", "tte9"), name.censoring = c("cnsr1","cnsr2","cnsr3", "cnsr4", "cnsr5", "cnsr6", "cnsr7", "cnsr8", "cnsr9"), scale.T = c(200,100,250,300,150,200,350,400,450), scale.censoring.T = c(10^5,10^5,10^5,10^5,10^5,10^5,10^5,10^5,10^5), scale.C = c(200,100,250,300,150,200,350,400,450), scale.censoring.C = c(10^5,10^5,10^5,10^5,10^5,10^5,10^5,10^5,10^5)), argsBin = list(name = "bin_var")) formula1 <- treatment ~ strata(strat_column) + tte(tte1, status = cnsr1, threshold = 10, restriction = 365) + tte(tte2, status = cnsr2, threshold = 10, restriction = 365) + tte(tte3, status = cnsr3, threshold = 10, restriction = 365) + tte(tte4, status = cnsr4, threshold = 10, restriction = 365) + bin(bin_var, operator = "<0") + tte(tte5, status = cnsr5, threshold = 10, restriction = 365) + tte(tte6, status = cnsr6, threshold = 10, restriction = 365) + tte(tte7, status = cnsr7, threshold = 10, restriction = 365) + tte(tte8, status = cnsr8, threshold = 10, restriction = 365) + tte(tte9, status = cnsr9, threshold = 10, restriction = 365) GPC.v1 <- BuyseTest(formula1, data = dt, trace = FALSE) GPC.v2 <- BuyseTest(treatment = "treatment", strata = "strat_column", endpoint = c("tte1","tte2","tte3","tte4","bin_var","tte5","tte6","tte7","tte8","tte9"), status = c("cnsr1","cnsr2","cnsr3","cnsr4","cnsr5","cnsr6","cnsr7","cnsr8","cnsr9"), type = c("tte","tte","tte","tte","bin","tte","tte","tte","tte","tte"), operator = c(">0",">0",">0",">0","<0",">0",">0",">0",">0",">0"), threshold = c(10,10,10,10,NA,10,10,10,10,10), restriction = c(365,365,365,365,NA,365,365,365,365,365), data = dt, trace = FALSE) }) ## * Brice: mandag 23-07-03 at 15:54 add.halfNeutral in summary test_that("number of pairs with add.halfNeutral",{ set.seed(1) dt <- simBuyseTest(n.T = 50, n.C = 50) BT <- BuyseTest(treatment ~ bin(toxicity), add.halfNeutral = TRUE, data = dt, trace = FALSE) expect_equal(100,model.tables(BT)$total, tol = 0.0001) }) ## * SamSalvaggio: (issue #2 on Github): 27 september 2023 restriction with multiple endpoint test_that("restriction with multiple endpoints",{ set.seed(1) dt <- simBuyseTest(n.T = 50, n.C = 50, names.strata = "strat_column", n.strata = 2, argsTTE = list(name = c("tte1"), name.censoring = c("cnsr1"), scale.T = c(200), scale.censoring.T = c(10^5), scale.C = c(200), scale.censoring.C = c(10^5)), argsBin = list(p.T = list(c(0.4,0.6))), argsCont = list(mu.T = c(0.3), sigma.T = 1) ) formula1 <- treatment ~ strat_column + bin(toxicity, operator = "<0") + cont(score, threshold = 0.5) + tte(tte1, status = cnsr1, threshold = 10, restriction = 365) GPC <- BuyseTest(formula1, data = dt, trace = FALSE, method.inference = "u-statistic") test <- confint(GPC) rownames(test) <- NULL GS <- data.frame("estimate" = c(-0.04265403, -0.00552923, -0.00592417), "se" = c(0.09983531, 0.1159374, 0.11753888), "lower.ci" = c(-0.2342771, -0.22865458, -0.23200199), "upper.ci" = c(0.15215946, 0.21814805, 0.22076088), "null" = c(0, 0, 0), "p.value" = c(0.66957926, 0.9619629, 0.95980315)) expect_equivalent(test,GS,tol = 1e-6) }) ## * graemeleehickey: (issue #13 on Github): fredag 24-02-09 at 12:20 Handling ties with Gehan's scoring rule test_that("Handling ties with Gehan's scoring rule",{ dat <- data.frame(time = c(10, 10), event = c(0, 1), treat = c(0, 1)) test <- BuyseTest(treat ~ TTE(time, status = event), data = dat, trace = FALSE, method.inference = "none", scoring.rule = "Gehan") expect_equal(c(0,1,0,0), c(coef(test, statistic = "count.favorable"), coef(test, statistic = "count.unfavorable"), coef(test, statistic = "count.neutral"), coef(test, statistic = "count.uninf")) ) testR <- BuyseTest(treat ~ TTE(time, status = event, restriction = 10), data = dat, trace = FALSE, method.inference = "none", scoring.rule = "Gehan") expect_equal(c(0,0,1,0), c(coef(testR, statistic = "count.favorable"), coef(testR, statistic = "count.unfavorable"), coef(testR, statistic = "count.neutral"), coef(testR, statistic = "count.uninf")) ) dat <- data.frame(time = c(10, 10), event = c(0, 1), treat = c(1, 0)) test <- BuyseTest(treat ~ TTE(time, status = event), data = dat, trace = FALSE, method.inference = "none", scoring.rule = "Gehan") expect_equal(c(1,0,0,0), c(coef(test, statistic = "count.favorable"), coef(test, statistic = "count.unfavorable"), coef(test, statistic = "count.neutral"), coef(test, statistic = "count.uninf")) ) testR <- BuyseTest(treat ~ TTE(time, status = event, restriction = 10), data = dat, trace = FALSE, method.inference = "none", scoring.rule = "Gehan") expect_equal(c(0,0,1,0), c(coef(testR, statistic = "count.favorable"), coef(testR, statistic = "count.unfavorable"), coef(testR, statistic = "count.neutral"), coef(testR, statistic = "count.uninf")) ) }) ## * aghaynes: (issue #14 on Github): Mi-september 2024 Bug with win odds test_that("Handling neutral pairs with win odds",{ set.seed(1) ddd <- simBuyseTest(1e2) ddd$sss <- ddd$id %% 2 ## ** estimate eFALSE.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 0.1), data = ddd, add.halfNeutral = FALSE, trace = FALSE, method.inference = "u-statistic") expect_equivalent(coef(eFALSE.BT, cumulative = TRUE, statistic = "winRatio"), cumsum(eFALSE.BT@count.favorable)/cumsum(eFALSE.BT@count.unfavorable), tol = 1e-6) expect_equivalent(coef(eFALSE.BT, cumulative = FALSE, statistic = "winRatio"), eFALSE.BT@count.favorable/eFALSE.BT@count.unfavorable, tol = 1e-6) eTRUE.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 0.1), data = ddd, add.halfNeutral = TRUE, trace = FALSE, method.inference = "u-statistic") expect_equivalent(coef(eTRUE.BT, cumulative = TRUE, statistic = "winRatio"), (cumsum(eTRUE.BT@count.favorable)+0.5*eTRUE.BT@count.neutral)/(cumsum(eTRUE.BT@count.unfavorable)+0.5*eTRUE.BT@count.neutral), tol = 1e-6) expect_equivalent(coef(eTRUE.BT, cumulative = FALSE, statistic = "winRatio"), (eTRUE.BT@count.favorable+0.5*eTRUE.BT@count.neutral)/(eTRUE.BT@count.unfavorable+0.5*eTRUE.BT@count.neutral), tol = 1e-6) eSFALSE.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 0.1) + sss, data = ddd, add.halfNeutral = FALSE, trace = FALSE, method.inference = "u-statistic") expect_equivalent(coef(eSFALSE.BT, cumulative = TRUE, statistic = "winRatio"), cumsum(colSums(eSFALSE.BT@count.favorable))/cumsum(colSums(eSFALSE.BT@count.unfavorable)), tol = 1e-6) expect_equivalent(coef(eSFALSE.BT, cumulative = FALSE, statistic = "winRatio"), colSums(eSFALSE.BT@count.favorable)/colSums(eSFALSE.BT@count.unfavorable), tol = 1e-6) expect_equivalent(coef(eSFALSE.BT, cumulative = TRUE, strata = TRUE, statistic = "winRatio"), riskRegression::rowCumSum(eSFALSE.BT@count.favorable)/riskRegression::rowCumSum(eSFALSE.BT@count.unfavorable), tol = 1e-6) expect_equivalent(coef(eSFALSE.BT, cumulative = FALSE, strata = TRUE, statistic = "winRatio"), eSFALSE.BT@count.favorable/eSFALSE.BT@count.unfavorable, tol = 1e-6) eSTRUE.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 0.1) + sss, data = ddd, add.halfNeutral = TRUE, trace = FALSE, method.inference = "u-statistic") expect_equivalent(coef(eSTRUE.BT, cumulative = TRUE, statistic = "winRatio"), (cumsum(colSums(eSTRUE.BT@count.favorable))+0.5*colSums(eSTRUE.BT@count.neutral))/(cumsum(colSums(eSTRUE.BT@count.unfavorable))+0.5*colSums(eSTRUE.BT@count.neutral)), tol = 1e-6) expect_equivalent(coef(eSTRUE.BT, cumulative = FALSE, statistic = "winRatio"), (colSums(eSTRUE.BT@count.favorable)+0.5*colSums(eSTRUE.BT@count.neutral))/(colSums(eSTRUE.BT@count.unfavorable)+0.5*colSums(eSTRUE.BT@count.neutral)), tol = 1e-6) expect_equivalent(coef(eSTRUE.BT, cumulative = TRUE, strata = TRUE, statistic = "winRatio"), (riskRegression::rowCumSum(eSTRUE.BT@count.favorable)+0.5*eSTRUE.BT@count.neutral)/(riskRegression::rowCumSum(eSTRUE.BT@count.unfavorable)+0.5*eSTRUE.BT@count.neutral), tol = 1e-6) expect_equivalent(coef(eSTRUE.BT, cumulative = FALSE, strata = TRUE, statistic = "winRatio"), (eSTRUE.BT@count.favorable+0.5*eSTRUE.BT@count.neutral)/(eSTRUE.BT@count.unfavorable+0.5*eSTRUE.BT@count.neutral), tol = 1e-6) ## ** iid (outcome specific) expect_equivalent(eTRUE.BT@iidAverage$favorable, eFALSE.BT@iidAverage$favorable + 0.5 * eTRUE.BT@iidAverage$neutral, tol = 1e-6) expect_equivalent(eTRUE.BT@iidAverage$unfavorable, eFALSE.BT@iidAverage$unfavorable + 0.5 * eTRUE.BT@iidAverage$neutral, tol = 1e-6) expect_equivalent(eTRUE.BT@iidAverage$neutral, eFALSE.BT@iidAverage$neutral, tol = 1e-6) expect_equivalent(eSTRUE.BT@iidAverage$favorable, eSFALSE.BT@iidAverage$favorable + 0.5 * eSTRUE.BT@iidAverage$neutral, tol = 1e-6) expect_equivalent(eSTRUE.BT@iidAverage$unfavorable, eSFALSE.BT@iidAverage$unfavorable + 0.5 * eSTRUE.BT@iidAverage$neutral, tol = 1e-6) expect_equivalent(eSTRUE.BT@iidAverage$neutral, eSFALSE.BT@iidAverage$neutral, tol = 1e-6) ## ** se (cumulated over outcomes) if(attr(eSTRUE.BT@method.inference,"hprojection")==1){ expect_equivalent(colSums(getIid(eFALSE.BT, statistic = "favorable")^2), eFALSE.BT@covariance[,"favorable"], tol = 1e-6) expect_equivalent(colSums(getIid(eTRUE.BT, statistic = "favorable")^2), eTRUE.BT@covariance[,"favorable"], tol = 1e-6) expect_equivalent(colSums(getIid(eSFALSE.BT, statistic = "favorable")^2), eSFALSE.BT@covariance[,"favorable"], tol = 1e-6) expect_equivalent(colSums(getIid(eSTRUE.BT, statistic = "favorable")^2), eSTRUE.BT@covariance[,"favorable"], tol = 1e-6) } xxx <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.1) + cont(score, threshold = 3) + sss, data = ddd, add.halfNeutral = TRUE, method.inference = "u-statistic") summary(xxx, statistic = "winRatio") summary(xxx, statistic = "netBenefit") })