### test-BuyseTest-iid.R --- ##---------------------------------------------------------------------- ## Author: Brice Ozenne ## Created: jan 8 2019 (11:54) ## Version: ## Last-Updated: sep 23 2024 (12:08) ## By: Brice Ozenne ## Update #: 220 ##---------------------------------------------------------------------- ## ### Commentary: ## ### Change Log: ##---------------------------------------------------------------------- ## ### Code: if(FALSE){ library(testthat) library(BuyseTest) library(data.table) } ## options("stringsAsFactors" = TRUE) context("Check correct computation of the variance \n") var2 <- function(x){var(x)*(length(x)-1)/length(x)} cov2 <- function(x,y){cov(x,y)*(length(x)-1)/length(x)} coef2 <- function(x){ do.call(cbind,setNames(lapply(c("favorable","unfavorable","netBenefit"), function(iStat){coef(x, statistic = iStat)}),c("favorable","unfavorable","netBenefit"))) } ## * Settings BuyseTest.options(check = TRUE, keep.pairScore = FALSE, method.inference = "u-statistic", pool.strata = "Buyse", trace = 0) ## * iid average ## ** 1 binary variable ## *** no strata ## equal number in each group d <- data.table(id = 1:4, group = c("C","C","T","T"), toxicity = c(1,0,1,0)) test_that("iid: binary and no strata (balanced groups)", { ## first order BuyseTest.options(order.Hprojection = 1) e.BT <- BuyseTest(group ~ bin(toxicity), data = d, method.inference = "u-statistic") e2.BT <- BuyseTest(group ~ bin(toxicity), data = d, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e2.BT@covariance) expect_equal(as.double(e.BT@covariance), c(1/16,1/16,-1/16, 1/4, 4) ) expect_equal(as.double(getIid(e.BT, statistic = "favorable")), c(-1/8,1/8,1/8,-1/8)) expect_equal(as.double(getIid(e.BT, statistic = "unfavorable")), c(1/8,-1/8,-1/8,1/8)) ## second order BuyseTest.options(order.Hprojection = 2) e.BT <- BuyseTest(group ~ bin(toxicity), data = d, keep.pairScore = FALSE, method.inference = "u-statistic") e1.BT <- BuyseTest(group ~ bin(toxicity), data = d, keep.pairScore = TRUE, method.inference = "u-statistic") e2.BT <- BuyseTest(group ~ bin(toxicity), data = d, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e1.BT@covariance) ## assumes vs. does not assumes binary score when computing second order terms expect_equal(e.BT@covariance, e2.BT@covariance) expect_equal(as.double(e.BT@covariance), c(5/64, 5/64, -3/64, 1/4, 4) ) }) ## unequal number in each group d.bis <- data.table(id = 1:4, group = c("C","T","T","T"), toxicity = c(1,1,1,0)) test_that("iid: binary and no strata (unbalanced groups)", { ## first order BuyseTest.options(order.Hprojection = 1) suppressWarnings(e.BT <- BuyseTest(group ~ bin(toxicity), data = d.bis, method.inference = "u-statistic")) suppressWarnings(e2.BT <- BuyseTest(group ~ bin(toxicity), data = d.bis, keep.pairScore = TRUE, method.inference = "u-statistic-bebu")) expect_equal(e.BT@covariance, e2.BT@covariance) expect_equal(as.double(e.BT@covariance), c(0, 2/27, 0, 2/27, 0) ) expect_equal(as.double(getIid(e.BT, statistic = "favorable")), c(0,0,0,0)) expect_equal(as.double(getIid(e.BT, statistic = "unfavorable")), c(0,-1/9,-1/9,2/9)) ## second order BuyseTest.options(order.Hprojection = 2) suppressWarnings(e.BT <- BuyseTest(group ~ bin(toxicity), data = d.bis, keep.pairScore = FALSE, method.inference = "u-statistic")) suppressWarnings(e1.BT <- BuyseTest(group ~ bin(toxicity), data = d.bis, keep.pairScore = TRUE, method.inference = "u-statistic")) suppressWarnings(e2.BT <- BuyseTest(group ~ bin(toxicity), data = d.bis, keep.pairScore = TRUE, method.inference = "u-statistic-bebu")) expect_equal(e.BT@covariance, e1.BT@covariance) ## assumes vs. does not assumes binary score when computing second order terms expect_equal(e.BT@covariance, e2.BT@covariance) expect_equal(as.double(e.BT@covariance), c(0, 2/27, 0, 2/27, 0) ) }) ## *** strata d2 <- rbind(cbind(d, strata = 1), cbind(d, strata = 2), cbind(d, strata = 3)) test_that("iid: binary with strata (balanced groups)", { for(iOrder in 1:2){ ## iOrder <- 2 BuyseTest.options(order.Hprojection = iOrder) e.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2, keep.pairScore = FALSE, method.inference = "u-statistic") e2.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2, keep.pairScore = TRUE, method.inference = "u-statistic") ebebu.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance,e2.BT@covariance) expect_equal(e.BT@covariance,ebebu.BT@covariance) ls.BT <- lapply(split(d2,d2$strata), function(iData){ BuyseTest(group ~ bin(toxicity), data = iData, method.inference = "u-statistic") }) M.estimate <- do.call(rbind,lapply(ls.BT, function(iBT){coef2(iBT)})) M.covariance <- do.call(rbind,lapply(ls.BT, function(iBT){iBT@covariance[,c(1:2,4)]})) n.pair <- sapply(ls.BT, function(iBT){iBT@n.pairs}) ntot.pair <- sum(n.pair) weight <- n.pair/ntot.pair expect_equal(colSums(riskRegression::colMultiply_cpp(M.estimate,weight)),as.double(coef2(e.BT))) expect_equal(colSums(riskRegression::colMultiply_cpp(M.covariance,weight^2)),as.double(e.BT@covariance[,c(1:2,4)])) expect_equal( rowSums(getIid(e.BT, statistic = "netBenefit", scale = FALSE, center = FALSE, strata = TRUE)), as.double(unlist(lapply(ls.BT, getIid, statistic = "netBenefit", scale = FALSE, center = FALSE))) ) ## getIid(e.BT, statistic = "netBenefit", scale = FALSE, center = FALSE, strata = TRUE) ## getIid(ls.BT[[1]], statistic = "netBenefit", scale = FALSE, center = FALSE) ## getIid(ls.BT[[2]], statistic = "netBenefit", scale = FALSE, center = FALSE) if(iOrder==1){ expect_equal(as.double(e.BT@covariance), c(1/16, 1/16, -1/16, 1/4, 4)/3 ) }else{ expect_equal(as.double(e.BT@covariance), c(5/64, 5/64, -3/64, 1/4, 4)/3 ) } expect_equal(as.double(getIid(e.BT, statistic = "favorable")), rep(c(-1/8,1/8,1/8,-1/8),3)/3) expect_equal(as.double(getIid(e.BT, statistic = "unfavorable")), rep(c(1/8,-1/8,-1/8,1/8),3)/3) } }) d2.bis <- rbind(cbind(d.bis, strata = 1), cbind(d, strata = 2), cbind(d.bis, strata = 3)) test_that("iid: binary and strata (unbalanced groups)", { for(iOrder in 1:2){ ## iOrder <- 2 BuyseTest.options(order.Hprojection = iOrder) e.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2.bis, keep.pairScore = FALSE, method.inference = "u-statistic") e2.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2.bis, keep.pairScore = TRUE, method.inference = "u-statistic") ebebu.BT <- BuyseTest(group ~ bin(toxicity) + strata, data = d2.bis, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") ## 0 0.1875 0 expect_equal(e.BT@covariance, e2.BT@covariance) expect_equal(e.BT@covariance, ebebu.BT@covariance) ls.BT <- lapply(split(d2.bis,d2.bis$strata), function(iData){ suppressWarnings(BuyseTest(group ~ bin(toxicity), data = iData, method.inference = "u-statistic")) }) M.estimate <- do.call(rbind,lapply(ls.BT, function(iBT){coef2(iBT)})) M.covariance <- do.call(rbind,lapply(ls.BT, function(iBT){iBT@covariance[,c(1:2,4)]})) n.pair <- sapply(ls.BT, function(iBT){iBT@n.pairs}) ntot.pair <- sum(n.pair) weight <- n.pair/ntot.pair expect_equal(colSums(riskRegression::colMultiply_cpp(M.estimate,weight)),as.double(coef2(e.BT))) expect_equal(colSums(riskRegression::colMultiply_cpp(M.covariance,weight^2)),as.double(e.BT@covariance[,c(1:2,4)])) expect_equal( rowSums(getIid(e.BT, statistic = "netBenefit", scale = FALSE, center = FALSE, strata = TRUE)), as.double(unlist(lapply(ls.BT, getIid, statistic = "netBenefit", scale = FALSE, center = FALSE))) ) if(iOrder==1){ expect_equal(as.double(e.BT@covariance), c(0.01, 0.02333333, -0.01, 0.05333333, 0.21399177), tol = 1e-7 ) }else{ expect_equal(as.double(e.BT@covariance), c(0.0125, 0.02583333, -0.0075, 0.05333333, 0.22633745), tol = 1e-7 ) } expect_equal(as.double(getIid(e.BT, statistic = "favorable")), c(0, 0, 0, 0, -0.05, 0.05, 0.05, -0.05, 0, 0, 0, 0), tol = 1e-7) expect_equal(as.double(getIid(e.BT, statistic = "unfavorable")), c(0, -0.03333333, -0.03333333, 0.06666667, 0.05, -0.05, -0.05, 0.05, 0, -0.03333333, -0.03333333, 0.06666667), tol = 1e-7) } }) ## ** 1 continuous variable n <- 5 set.seed(10) dt <- simBuyseTest(n, name.cluster = NULL) dtS <- rbind(cbind(S = 1, dt), cbind(S = 2, dt), cbind(S = 3, dt)) dtS$score1 <- dtS$score dtS$score2 <- dtS$score dtS[S == 1, score1 := 1] dtS[S > 1, score2 := 1] test_that("Manual calculation of second order H projection (no strata)",{ BuyseTest.options(order.Hprojection = 1) e.BT_c1 <- BuyseTest(treatment ~ cont(score), data = dt, trace = 0, method.inference = "u-statistic") BuyseTest.options(order.Hprojection = 2) e.BT_c2 <- BuyseTest(treatment ~ cont(score), data = dt, trace = 0, method.inference = "u-statistic") e.BT_c3 <- BuyseTest(treatment ~ cont(score), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") ## manual calculation dt.pair <- getPairScore(e.BT_c3)[,.(index.C,index.T,favorable,unfavorable)] dt.pair[, H1C.favorable := mean(favorable), by = index.C] dt.pair[, H1T.favorable := mean(favorable), by = index.T] dt.pair[, H1C.favorable := H1C.favorable - mean(favorable)] dt.pair[, H1T.favorable := H1T.favorable - mean(favorable)] dt.pair[, H1C.unfavorable := mean(unfavorable), by = index.C] dt.pair[, H1T.unfavorable := mean(unfavorable), by = index.T] dt.pair[, H1C.unfavorable := H1C.unfavorable - mean(unfavorable)] dt.pair[, H1T.unfavorable := H1T.unfavorable - mean(unfavorable)] dt.pair[, H2.favorable := favorable - H1C.favorable - H1T.favorable] dt.pair[, H2.unfavorable := unfavorable - H1C.unfavorable - H1T.unfavorable] ## check H1 expect_true(all(abs(dt.pair[!duplicated(index.C),H1C.favorable/.N]-getIid(e.BT_c3, statistic = "favorable")[1:n])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.C),H1C.unfavorable/.N]-getIid(e.BT_c3, statistic = "unfavorable")[1:n])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.T),H1T.favorable/.N]-getIid(e.BT_c3, statistic = "favorable")[n+(1:n)])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.T),H1T.unfavorable/.N]-getIid(e.BT_c3, statistic = "unfavorable")[n+(1:n)])<1e-6)) ## check H2 manual <- dt.pair[,.(favorable = var2(H2.favorable)/.N, unfavorable = var2(H2.unfavorable)/.N, covariance = cov2(H2.favorable,H2.unfavorable)/.N)] expect_equal(as.double(manual), as.double((e.BT_c2@covariance - e.BT_c1@covariance)[1,c("favorable","unfavorable","covariance")])) expect_equal(as.double(manual), c(0.00256, 0.00256, -0.00256), tol = 1e-5) }) test_that("Manual calculation of second order H projection (strata)",{ BuyseTest.options(order.Hprojection = 1) e.BT_c1 <- BuyseTest(treatment ~ cont(score) + S, data = dtS, trace = 0, method.inference = "u-statistic") BuyseTest.options(order.Hprojection = 2) e.BT_c2 <- BuyseTest(treatment ~ cont(score) + S, data = dtS, trace = 0, method.inference = "u-statistic") e.BT_c3 <- BuyseTest(treatment ~ cont(score) + S, data = dtS, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") ## manual calculation dt.pair <- getPairScore(e.BT_c3)[,.(strata,index.C,index.T,favorable,unfavorable)] dt.pair[, H1C.favorable := mean(favorable), by = index.C] dt.pair[, H1T.favorable := mean(favorable), by = index.T] dt.pair[, H1C.favorable := H1C.favorable - mean(favorable), by = "strata"] dt.pair[, H1T.favorable := H1T.favorable - mean(favorable), by = "strata"] dt.pair[, H1C.unfavorable := mean(unfavorable), by = index.C] dt.pair[, H1T.unfavorable := mean(unfavorable), by = index.T] dt.pair[, H1C.unfavorable := H1C.unfavorable - mean(unfavorable), by = "strata"] dt.pair[, H1T.unfavorable := H1T.unfavorable - mean(unfavorable), by = "strata"] dt.pair[, H2.favorable := favorable - H1C.favorable - H1T.favorable] dt.pair[, H2.unfavorable := unfavorable - H1C.unfavorable - H1T.unfavorable] ## check H1 expect_true(all(abs(dt.pair[!duplicated(index.C),H1C.favorable/.N]-getIid(e.BT_c3, statistic = "favorable")[which(dtS$treatment=="C")])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.C),H1C.unfavorable/.N]-getIid(e.BT_c3, statistic = "unfavorable")[which(dtS$treatment=="C")])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.T),H1T.favorable/.N]-getIid(e.BT_c3, statistic = "favorable")[which(dtS$treatment=="T")])<1e-6)) expect_true(all(abs(dt.pair[!duplicated(index.T),H1T.unfavorable/.N]-getIid(e.BT_c3, statistic = "unfavorable")[which(dtS$treatment=="T")])<1e-6)) ## check H2 manual <- dt.pair[,.(favorable = var2(H2.favorable)/.N, unfavorable = var2(H2.unfavorable)/.N, covariance = cov2(H2.favorable,H2.unfavorable)/.N)] expect_equal(as.double(manual), as.double((e.BT_c2@covariance - e.BT_c1@covariance)[1,c("favorable","unfavorable","covariance")])) expect_equal(as.double(manual), c(0.00085333, 0.00085333, -0.00085333), tol = 1e-5) }) ## ** 1 TTE variable test_that("iid: TTE and no strata",{ BuyseTest.options(order.Hprojection = 1) e.BT_tte1 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, scoring.rule = "Gehan", method.inference = "u-statistic") e.BT_tte2 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1e5) + tte(eventtime, status, threshold = 1-1e-5), data = dt, trace = 0, keep.pairScore = TRUE, scoring.rule = "Gehan", method.inference = "u-statistic") e.BT_tte3 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1) + tte(eventtime, status, threshold = 1-1e-5), data = dt, trace = 0, keep.pairScore = TRUE, scoring.rule = "Gehan", method.inference = "u-statistic") expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte2)[2,], tol = 1e-6) expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte3)[1,], tol = 1e-6) expect_equivalent(confint(e.BT_tte3)[1,], confint(e.BT_tte3)[2,], tol = 1e-6) }) ## ** Two endpoints ## *** no strata test_that("iid: two endpoints (no strata - first order)", { ## different endpoints BuyseTest.options(order.Hprojection = 1) e.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 1), data = dtS, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e2.BT@covariance) GS <- matrix(c(0.01066667, 0.01258667, 0, 0.00085333, 0, -0.00138667, 0.01066667, 0.01621333, NaN, 136.66666667), nrow = 2, ncol = 5, dimnames = list(c("toxicity", "score_t1"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) ## same endpoint e.BT <- BuyseTest(treatment ~ cont(score, threshold = 2) + cont(score, threshold = 1), data = dtS, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ cont(score, threshold = 2) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e2.BT@covariance) GS <- matrix(c(0.01194667, 0.01408, 0, 0.00085333, 0, -0.00149333, 0.01194667, 0.01792, NaN, 108), nrow = 2, ncol = 5, dimnames = list(c("score_t2", "score_t1"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) ## same endpoint tte e.BT <- BuyseTest(treatment ~ tte(eventtime, threshold = 1, status = status) + tte(eventtime, threshold = 0, status = status), data = dtS, scoring.rule = "Gehan", method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ tte(eventtime, threshold = 1, status = status) + tte(eventtime, threshold = 0, status = status), data = dtS, keep.pairScore = TRUE, scoring.rule = "Gehan", method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e2.BT@covariance) GS <- matrix(c(0, 0.00234667, 0.00234667, 0.00832, 0, 0.00064, 0.00234667, 0.00938667, 0, 0.04938272), nrow = 2, ncol = 5, dimnames = list(c("eventtime_t1", "eventtime"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) ## cluster argument expect_equal(unname(getIid(e.BT, cluster = 1:NROW(dtS))), unname(getIid(e.BT)), tol = 1e-6) expect_equivalent(confint(e.BT, cluster = 1:NROW(dtS)), confint(e.BT), tol = 1e-6) }) test_that("iid: two endpoints (no strata - second order)", { BuyseTest.options(order.Hprojection = 2) ## different endpoints e.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 1), data = dtS, keep.pairScore = FALSE, method.inference = "u-statistic") e1.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ bin(toxicity) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") remaining.pairs <- e1.BT@tablePairScore[[2]]$index.pair expect_equal(e1.BT@tablePairScore[[1]][index.pair %in% remaining.pairs,.(index.C,index.T)], e1.BT@tablePairScore[[2]][,.(index.C,index.T)]) expect_equal(e.BT@covariance, e1.BT@covariance) expect_equal(e.BT@covariance, e2.BT@covariance) GS <- matrix(c(0.01066667, 0.01284267, 0, 0.00096711, 0, -0.00139378, 0.01066667, 0.01659733, NaN, 150.88888889), nrow = 2, ncol = 5, dimnames = list(c("toxicity", "score_t1"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) ## same endpoint e.BT <- BuyseTest(treatment ~ cont(score, threshold = 2) + cont(score, threshold = 1), data = dtS, keep.pairScore = FALSE, method.inference = "u-statistic") e1.BT <- BuyseTest(treatment ~ cont(score, threshold = 2) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ cont(score, threshold = 2) + cont(score, threshold = 1), data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e1.BT@covariance) expect_equal(e.BT@covariance, e2.BT@covariance) remaining.pairs <- e1.BT@tablePairScore[[2]]$index.pair expect_equal(e1.BT@tablePairScore[[1]][index.pair %in% remaining.pairs,.(index.C,index.T)], e1.BT@tablePairScore[[2]][,.(index.C,index.T)]) GS <- matrix(c(0.01211733, 0.01425067, 0, 0.00096711, 0, -0.00147911, 0.01211733, 0.018176, NaN, 118.13333333), nrow = 2, ncol = 5, dimnames = list(c("score_t2", "score_t1"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) ## same endpoint tte e.BT <- BuyseTest(treatment ~ tte(eventtime, threshold = 1, status = status) + tte(eventtime, threshold = 0, status = status), data = dtS, scoring.rule = "Gehan", keep.pairScore = FALSE, method.inference = "u-statistic") e1.BT <- BuyseTest(treatment ~ tte(eventtime, threshold = 1, status = status) + tte(eventtime, threshold = 0, status = status), data = dtS, scoring.rule = "Gehan", keep.pairScore = TRUE, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ tte(eventtime, threshold = 1, status = status) + tte(eventtime, threshold = 0, status = status), data = dtS, keep.pairScore = TRUE, scoring.rule = "Gehan", method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance, e1.BT@covariance) expect_equal(e.BT@covariance, e2.BT@covariance) remaining.pairs <- e1.BT@tablePairScore[[2]]$index.pair expect_equal(e1.BT@tablePairScore[[1]][index.pair %in% remaining.pairs,.(index.C,index.T)], e1.BT@tablePairScore[[2]][,.(index.C,index.T)]) GS <- matrix(c(0, 0.00251733, 0.00251733, 0.008576, 0, 0.000512, 0.00251733, 0.01006933, 0, 0.05432099), nrow = 2, ncol = 5, dimnames = list(c("eventtime_t1", "eventtime"),c("favorable", "unfavorable", "covariance", "netBenefit", "winRatio")) ) expect_equal(e.BT@covariance, GS, tol = 1e-6 ) }) ## *** strata test_that("iid: two endpoints (strata)", { for(iOrder in 1:2){ ## iOrder <- 1 BuyseTest.options(order.Hprojection = iOrder) e.BT <- BuyseTest(treatment ~ cont(score1, threshold = 1) + cont(score2, threshold = 1) + S, data = dtS, method.inference = "u-statistic") e2.BT <- BuyseTest(treatment ~ cont(score1, threshold = 1) + cont(score2, threshold = 1) + S, data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic") ebebu.BT <- BuyseTest(treatment ~ cont(score1, threshold = 1) + cont(score2, threshold = 1) + S, data = dtS, keep.pairScore = TRUE, method.inference = "u-statistic-bebu") expect_equal(e.BT@covariance,e2.BT@covariance) expect_equal(e.BT@covariance,ebebu.BT@covariance) ls.BT <- lapply(split(dtS,dtS$S), function(iData){ BuyseTest(treatment ~ cont(score1, threshold = 1) + cont(score2, threshold = 1), data = iData, method.inference = "u-statistic") }) ls.estimate <- lapply(ls.BT, function(iBT){coef2(iBT)}) ls.covariance <- lapply(ls.BT, function(iBT){iBT@covariance[,c(1:2,4)]}) n.pair <- sapply(ls.BT, function(iBT){iBT@n.pairs}) ntot.pair <- sum(n.pair) weight <- n.pair/ntot.pair expect_equal(Reduce("+",lapply(1:length(weight), function(iS){ls.estimate[[iS]]*weight[iS]})),coef2(e.BT), tol = 1e-7) expect_equal(Reduce("+",lapply(1:length(weight), function(iS){ls.covariance[[iS]]*weight[iS]^2})),e.BT@covariance[,c(1:2,4)]) expect_equal( Reduce("+",getIid(e.BT, statistic = "netBenefit", scale = FALSE, center = FALSE, strata = TRUE)), do.call(rbind,lapply(ls.BT, getIid, statistic = "netBenefit", scale = FALSE, center = FALSE)) ) if(iOrder==1){ expect_equal(as.double(e.BT@covariance[2,]), c(0.01408, 0.00085333, -0.00149333, 0.01792, 108), tol = 1e-6 ) }else{ expect_equal(as.double(e.BT@covariance[2,]), c(0.014592, 0.00119467, -0.00145067, 0.018688, 138.4), tol = 1e-6 ) } } }) ## *** cluster dtS <- rbind(cbind(S = 1, dt, id = 1:NROW(dt)), cbind(S = 2, dt, id = 1:NROW(dt)), cbind(S = 3, dt, id = 1:NROW(dt))) test_that("cluster option", { BuyseTest.options(order.Hprojection = 1) e.BT <- BuyseTest(treatment ~ cont(score, threshold = 1) + bin(toxicity), data = dt, method.inference = "u-statistic") e3.BT <- BuyseTest(treatment ~ cont(score, threshold = 1) + bin(toxicity) + S, data = dtS, method.inference = "u-statistic") expect_equal(unname(getIid(e.BT)),unname(getIid(e3.BT, cluster = dtS$id)), tol = 1e-6) expect_equivalent(as.data.frame(confint(e.BT)), as.data.frame(confint(e3.BT, cluster = dtS$id)), tol = 1e-6) }) ## * iid Peron ## ** 1 TTE variable n <- 5 set.seed(10) dt <- simBuyseTest(n) dt$X0 <- 0 dt$treatment2 <- as.numeric(dt$treatment=="C") ## data.table("treatment" = c("C", "C", "C", "C", "C", "T", "T", "T", "T", "T"), ## "eventtime" = c(0.60539304, 0.31328027, 0.03946623, 0.32147489, 1.57044952, 0.29069131, 0.19522131, 0.04640668, 0.05277335, 0.43062009), ## "status" = c(0, 1, 0, 1, 0, 0, 0, 0, 1, 1), ## "toxicity" = c("yes", "yes", "yes", "yes", "yes", "no", "yes", "yes", "yes", "yes"), ## "score" = c(-1.85374045, -0.07794607, 0.96856634, 0.18492596, -1.37994358, 1.10177950, 0.75578151, -0.23823356, 0.98744470, 0.74139013), ## "X0" = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), ## "treatment2" = c(1, 1, 1, 1, 1, 0, 0, 0, 0, 0)) test_that("iid with nuisance parameters: 1 TTE",{ BuyseTest.options(order.Hprojection = 1) e.BT_tte1 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte1.bis <- BuyseTest(treatment2 ~ tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte2 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1e5) + tte(eventtime, status, threshold = 1-1e-5), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte3 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1) + tte(eventtime, status, threshold = 1-1e-5), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte4 <- BuyseTest(treatment ~ bin(X0) + tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte5 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 2) + bin(X0) + tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") e.BT_tte5.bis <- BuyseTest(treatment2 ~ tte(eventtime, status, threshold = 2) + bin(X0) + tte(eventtime, status, threshold = 1), data = dt, trace = 0, keep.pairScore = TRUE, method.inference = "u-statistic") ## unchanged when switching around the groups expect_equal(e.BT_tte1@count.unfavorable,e.BT_tte1.bis@count.favorable) expect_equal(e.BT_tte1@covariance["unfavorable"],e.BT_tte1.bis@covariance["favorable"]) expect_equal(e.BT_tte5@count.unfavorable,e.BT_tte5.bis@count.favorable) expect_equal(e.BT_tte5@covariance["unfavorable"],e.BT_tte5.bis@covariance["favorable"]) ## results does not depend on previously used thresholds expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte2)[2,], tol = 1e-6) expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte3)[1,], tol = 1e-6) expect_equivalent(confint(e.BT_tte2)[2,], confint(e.BT_tte3)[2,], tol = 1e-6) expect_equivalent(confint(e.BT_tte3)[1,], confint(e.BT_tte3)[2,], tol = 1e-6) expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte4)[2,], tol = 1e-6) expect_equivalent(confint(e.BT_tte1)[1,], confint(e.BT_tte5)[3,], tol = 1e-6) }) ## ** 1 TTE variable and 1 binary n <- 5 set.seed(10) dt <- simBuyseTest(n, argsTTE = list(scale.T = 2, scale.censoring.T = 1)) ## data.table("treatment" = c("C", "C", "C", "C", "C", "T", "T", "T", "T", "T"), ## "eventtime" = c(0.60539304, 0.98681219, 0.03946623, 1.28589957, 1.57044952, 0.29069131, 0.19522131, 0.04640668, 0.21109340, 0.69214123), ## "status" = c(0, 0, 0, 1, 0, 0, 0, 0, 1, 0), ## "toxicity" = c("yes", "yes", "yes", "yes", "yes", "no", "yes", "yes", "yes", "yes"), ## "score" = c(-1.85374045, -0.07794607, 0.96856634, 0.18492596, -1.37994358, 1.10177950, 0.75578151, -0.23823356, 0.98744470, 0.74139013)) test_that("iid with nuisance parameters: 1 TTE + 1 binary",{ BuyseTest.options(order.Hprojection = 1) e.BT_ttebin <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1) + bin(toxicity), data = dt, keep.pairScore = TRUE, method.inference = "u-statistic") test <- confint(e.BT_ttebin) GS <- matrix(c(-0.33333333, -0.13333333, 0.24130536, 0.36004622, -0.70573842, -0.69241599, 0.18339631, 0.52579681, 0, 0, 0.20172157, 0.7144262), nrow = 2, ncol = 6, dimnames = list(c("eventtime_t1", "toxicity"),c("estimate", "se", "lower.ci", "upper.ci","null", "p.value")) ) expect_equivalent(test, as.data.frame(GS), tol = 1e-6) ## exponential approximation of the survival when computing the influence function e.TTEM <- BuyseTTEM(Hist(eventtime,status)~treatment, data = dt, iid=TRUE, iid.surv="prodlim") attr(e.TTEM, "iidNuisance") <- TRUE e.BT_ttebin <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1) + bin(toxicity), data = dt, keep.pairScore = TRUE, model.tte = e.TTEM, method.inference = "u-statistic") test <- confint(e.BT_ttebin) GS <- matrix(c(-0.33333333, -0.13333333, 0.23587679, 0.35518499, -0.69967949, -0.68733243, 0.17180423, 0.51874253, 0, 0, 0.19153767, 0.71069249), nrow = 2, ncol = 6, dimnames = list(c("eventtime_t1", "toxicity"),c("estimate", "se", "lower.ci", "upper.ci", "null","p.value")) ) expect_equivalent(test, as.data.frame(GS), tol = 1e-6) ## GS <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 1) + bin(toxicity), ## data = dt, ## keep.pairScore = TRUE, ## method.inference = "bootstrap") }) ## ** 2 TTE variables n.patients <- c(20,20) set.seed(10) dt.sim <- simBuyseTest(n.T = n.patients[1], n.C = n.patients[2], argsBin = list(p.T = list(c(0.5,0.5),c(0.25,0.75))), argsCont = list(mu.T = 1:3, sigma.T = rep(1,3)), argsTTE = list(scale.T = 1:3, scale.censoring.T = rep(1,3))) setkeyv(dt.sim,c("treatment","eventtime1")) ## dt.sim[,status1.bis := c(status1[1:(.N-1)],1),by="treatment"] ## make sure last observation is a case test_that("iid with nuisance parameters: 2 TTE",{ ## plot(prodlim(Hist(eventtime1,status1.bis) ~ treatment, data = dt.sim)) e.BT_tte1 <- BuyseTest(treatment ~ tte(eventtime2, status2, threshold = 1), data = dt.sim, method.inference = "u-statistic") ## e.BT_tte2 <- BuyseTest(treatment ~ tte(eventtime1, status1.bis, threshold = 1e5) + tte(eventtime2, status2, threshold = 1), ## data = dt.sim, ## method.inference = "u-statistic") e.BT_tte3 <- BuyseTest(treatment ~ tte(eventtime1, status1, threshold = 1e5) + tte(eventtime2, status2, threshold = 1), data = dt.sim, method.inference = "u-statistic") expect_equal(e.BT_tte1@count.favorable[1], e.BT_tte3@count.favorable[2]) expect_equal(e.BT_tte1@count.unfavorable[1], e.BT_tte3@count.unfavorable[2]) expect_equal(e.BT_tte1@count.neutral[1], e.BT_tte3@count.neutral[2]) expect_equal(e.BT_tte1@count.uninf[1], e.BT_tte3@count.uninf[2]) ## expect_equal(e.BT_tte1@covariance[1,],e.BT_tte3@covariance[2,], tol = 1e-6) e.BT_tte <- BuyseTest(treatment ~ tte(eventtime1, status1, threshold = 1) + tte(eventtime2, status2, threshold = 1) + bin(toxicity1), data = dt.sim, method.inference = "u-statistic") test <- confint(e.BT_tte) GS <- matrix(c(0.26401345, 0.179801, 0.00853608, 0.19937703, 0.2148585, 0.2408939, -0.14852611, -0.24811831, -0.43304747, 0.59828277, 0.54900838, 0.4468153, 0, 0, 0, 0.20703017, 0.41296871, 0.97173422), nrow = 3, ncol = 6, dimnames = list(c("eventtime1_t1", "eventtime2_t1", "toxicity1"),c("estimate", "se", "lower.ci", "upper.ci", "null", "p.value")) ) expect_equivalent(test, as.data.frame(GS), tol = 1e-3) ## exponential approximation of the survival when computing the influence function e.TTEM <- list(BuyseTTEM(Hist(eventtime1,status1)~treatment, data = dt.sim, iid=TRUE, iid.surv="prodlim"), BuyseTTEM(Hist(eventtime2,status2)~treatment, data = dt.sim, iid=TRUE, iid.surv="prodlim")) attr(e.TTEM, "iidNuisance") <- TRUE e.BT_tte <- BuyseTest(treatment ~ tte(eventtime1, status1, threshold = 1) + tte(eventtime2, status2, threshold = 1) + bin(toxicity1), data = dt.sim, model.tte = e.TTEM, method.inference = "u-statistic") test <- confint(e.BT_tte) GS <- matrix(c(0.26401345, 0.179801, 0.00853608, 0.17433724, 0.19140473, 0.21888023, -0.09657678, -0.20304112, -0.39734512, 0.5633411, 0.51495999, 0.41162398, 0, 0, 0, 0.14902033, 0.35809689, 0.96889279), nrow = 3, ncol = 6, dimnames = list(c("eventtime1_t1", "eventtime2_t1", "toxicity1"),c("estimate", "se", "lower.ci", "upper.ci", "null", "p.value")) ) expect_equivalent(test, as.data.frame(GS), tol = 1e-3) }) ## ** 1 TTE with strata test_that("iid with nuisance parameters: 1 TTE + strata",{ BuyseTest.options(order.Hprojection = 1) e.BT_tteS <- BuyseTest(treatment ~ tte(eventtime1, status1, threshold = 1) + toxicity1, data = dt.sim, method.inference = "u-statistic") ## e.BT_tteS@covariance ls.BT_tteS <- lapply(split(dt.sim,dt.sim$toxicity1), function(iData){ BuyseTest(treatment ~ tte(eventtime1, status1, threshold = 1), data = iData, method.inference = "u-statistic") }) M.estimate <- do.call(rbind,lapply(ls.BT_tteS, function(iBT){coef2(iBT)})) M.covariance <- do.call(rbind,lapply(ls.BT_tteS, function(iBT){iBT@covariance[,c(1:2,4)]})) n.pair <- sapply(ls.BT_tteS, function(iBT){iBT@n.pairs}) ntot.pair <- sum(n.pair) weight <- n.pair/ntot.pair expect_equal(colSums(riskRegression::colMultiply_cpp(M.estimate,weight)),as.double(coef2(e.BT_tteS))) expect_equal(colSums(riskRegression::colMultiply_cpp(M.covariance,weight^2)),as.double(e.BT_tteS@covariance[,c(1:2,4)])) test <- rowSums(getIid(e.BT_tteS, statistic = "netBenefit", scale = FALSE, center = FALSE, strata = TRUE)) GS <- unlist(lapply(ls.BT_tteS, getIid, statistic = "netBenefit", scale = FALSE, center = FALSE)) expect_equal(test[unlist(split(1:NROW(dt.sim),dt.sim$toxicity1))], unname(GS), tol = 1e-7) }) ## * normalization iid n <- 200 set.seed(10) dt <- simBuyseTest(n, n.strata = 3) test_that("iid - remove normalization", { e.all <- BuyseTest(treatment~bin(toxicity)+cont(score)+strata, method.inference = "u-statistic", data = dt, trace = 0) ## summary(e.all) e.strata <- BuyseTest(treatment~bin(toxicity)+cont(score), method.inference = "u-statistic", data = dt[strata=="a"], trace = 0) ## summary(e.strata) iid.all <- getIid(e.all, scale = FALSE, center = FALSE, strata = "a")[which(dt$strata=="a"),"score"] iid.strata <- as.double(getIid(e.strata, scale = FALSE, center = FALSE)[,"score"]) expect_equal(unname(iid.all),unname(iid.strata), tol = 1e-9) GS <- c(0.37313, -0.70149, 0.40299, -0.13433, -0.76119, 0.43284, 0.28358, -0.79104, -0.73134, 0.28358, -0.43284, -0.13433, -0.16418, -0.22388, -0.31343, -0.67164, -0.97015, 0.43284, 0.28358, -0.64179, -0.9403, 0.19403, 0.43284, 0.01493, -0.91045, -0.13433, 0.04478, 0.52239, 0.91045, 0.10448, 0.07463, 0.40299, 0.28358, -0.43284, 0.07463, -0.79104, -0.73134, 0.07463, -0.64179, -0.67164, 0.64179, -0.22388, -0.76119, 1, -0.34328, 0.70149, -0.31343, -0.58209, -0.49254, -0.58209, -0.31343, 0.49254, 0.43284, -0.73134, 0.22388, 0.91045, 0.73134, 0.28358, 0.49254, 0.70149, -0.07463, 0.43284, -0.61194, -0.8806, 0.91045, -0.2, -0.01538, 0.32308, -0.96923, -0.84615, -0.41538, -0.96923, 0.66154, -0.66154, 0.13846, 0.38462, -0.72308, -0.75385, 0.41538, -0.75385, 0.2, 0.81538, -0.87692, -0.26154, 0.38462, 0.32308, 1, -0.23077, -0.41538, -0.78462, -0.07692, 0.50769, 0.41538, -0.96923, -0.87692, 0.90769, 0.93846, -0.16923, -0.26154, 0.75385, 0.13846, -0.01538, 0.10769, -0.01538, -0.87692, 0.41538, -0.66154, -0.75385, 0.01538, 0.63077, 0.32308, 0.29231, 0.2, -0.78462, 0.87692, 0.87692, 0.47692, -0.41538, 0.56923, -0.2, 0.96923, 0.87692, -0.44615, 0.2, -0.2, -0.50769, -0.87692, 0.01538, -0.87692, -0.87692, -0.75385, -0.04615) expect_equal(as.double(iid.all),GS, tol = 1e-4) }) ###################################################################### ### test-BuyseTest-iid.R ends here