BuyseTest: Generalized Pairwise Comparisons

## chunk 2
library(BuyseTest)
library(data.table)

## chunk 3
data(cancer, package = "survival")
veteran <- cbind(id = 1:NROW(veteran), veteran)
veteran$trt <- factor(veteran$trt,1:2,c("Pl","Exp"))
head(veteran)

## chunk 4
utils::packageVersion("BuyseTest")

## chunk 5
sessionInfo()

## * Performing generalized pairwise comparisons (GPC)

## chunk 6
BT <- BuyseTest(data = veteran, 
                endpoint = "time", 
                type = "timeToEvent", 
                treatment = "trt", 
                status = "status", 
                threshold = 20)

## chunk 7
BT.f <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                  data = veteran, trace = FALSE)

## chunk 8
BT.f@call <- list(); BT@call <- list();
testthat::expect_equal(BT.f,BT)

## ** Displaying the results

## chunk 9
summary(BT)

## chunk 10
print(BT, percentage = FALSE)

## chunk 11
model.tables(BT, percentage = FALSE)

## chunk 12
confint(BT)

## chunk 13
coef(BT)

## ** What about other summary statistics?

## chunk 14
summary(BT, statistic = "winRatio")

## chunk 15
confint(BT, statistic = "favorable")

## chunk 16
BT.perm <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                     data = veteran, trace = FALSE,
                     method.inference = "permutation", seed = 10)
confint(BT.perm, statistic = "favorable")

## chunk 17
rbind(confint(BT, statistic = "favorable", null = 0.42),
      confint(BT, statistic = "favorable", null = 0.5))

## chunk 18
BT.half <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                     data = veteran, trace = FALSE, add.halfNeutral = TRUE)
confint(BT.half, statistic = "favorable")

## chunk 19
confint(BT.half, statistic = "winRatio")

## chunk 20
BT.halfperm <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                         data = veteran, trace = FALSE, add.halfNeutral = TRUE,
                         method.inference = "bootstrap", seed = 10)
Mstat <- rbind(netBenefit = confint(BT.halfperm, statistic = "netBenefit"),
               winRatio = confint(BT.halfperm, statistic = "winRatio"),
               favorable = confint(BT.halfperm, statistic = "favorable"))
Mstat

## ** Stratified GPC

## chunk 21
ffstrata <- trt ~ tte(time, threshold = 20, status = "status") + celltype
BTstrata <- BuyseTest(ffstrata, data = veteran, trace = 0)

## chunk 22
keep.colStrata <- c("endpoint","strata", "total",
                    "favorable","unfavorable","neutral","uninf","delta","Delta")
summary(BTstrata, type.display = keep.colStrata)

## chunk 23
strata.obs <- as.data.frame(nobs(BTstrata, strata = TRUE))
strata.obs

## chunk 24
dfStrata <- model.tables(BTstrata, percentage = FALSE,
                         strata = c("squamous","smallcell","adeno","large"),
                         columns = c("strata","total","favorable","unfavorable"))
dfStrata

## chunk 25
delta <- (dfStrata$favorable - dfStrata$unfavorable)/strata.obs$pairs
delta

## chunk 26
weightCMH <- strata.obs$pairs/(strata.obs$Pl + strata.obs$Exp)

list(estimate = sum(delta * weightCMH/sum(weightCMH)),
     weight = 100*weightCMH/sum(weightCMH))

## chunk 27
BTstrata2 <- BuyseTest(ffstrata, data = veteran, trace = 0, pool.strata = "buyse")
summary(BTstrata2, type.display = keep.colStrata)

## chunk 28
confint(BTstrata2)

## chunk 29
confint(BTstrata, strata = TRUE)

## ** Using multiple endpoints

## chunk 30
ff2 <- trt ~ tte(time, threshold = 20, status = "status") + cont(karno, threshold = 0)
BT.H <- BuyseTest(ff2, data = veteran, trace = 0)
summary(BT.H)

## chunk 31
plot(BT.H, type = "hist")
plot(BT.H, type = "pie")
plot(BT.H, type = "racetrack")

## chunk 33
BT.nH <- BuyseTest(ff2, hierarchical = FALSE, data = veteran, trace = 0)
summary(BT.nH)

## chunk 34
ff2w <- trt ~ tte(time, threshold = 20, status = "status", weight = 0.8)
ff2w <- update.formula(ff2w, . ~ . + cont(karno, threshold = 0, weight = 0.2))
BT.nHw <- BuyseTest(ff2w, hierarchical = FALSE, data = veteran, trace = 0)
model.tables(BT.nHw)

## chunk 35
confint(BuyseTest(trt ~ cont(karno, threshold = 0), data = veteran, trace = 0))

## chunk 36
confint(BT.nHw, cumulative = FALSE)

## chunk 37
confint(BT.nHw, transform = FALSE)

## ** Statistical inference

## chunk 38
BT.perm <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                     data = veteran, trace = 0, method.inference = "permutation",
                     seed = 10) 
summary(BT.perm)

## chunk 39
BT.boot <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                     data = veteran, trace = 0, method.inference = "bootstrap",
                     seed = 10) 
summary(BT.boot)

## chunk 40
BT.ustat <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
                      data = veteran, trace = 0, method.inference = "u-statistic") 
summary(BT.ustat)

## chunk 41
set.seed(10)
sapply(1:10, function(i){mean(rbinom(1e4, size = 1, prob = 0.05))})

## ** What if smaller is better?

## chunk 42
ffop <- trt ~ tte(time, status = "status", threshold = 20, operator = "<0")
BTinv <- BuyseTest(ffop, data = veteran, trace = 0)
summary(BTinv)

## ** Stopping comparison for neutral pairs

## chunk 43
dt.sim <- data.table(Id = 1:2,
                     treatment = c("Yes","No"),
                     tumor = c("Yes","Yes"),
                     size = c(15,20))
dt.sim

## chunk 44
BT.pair <- BuyseTest(treatment ~ bin(tumor) + cont(size, operator = "<0"), data = dt.sim,
                     trace = 0, method.inference = "none")
summary(BT.pair)

## chunk 45
BT.pair2 <- BuyseTest(treatment ~ bin(tumor) + cont(size, operator = "<0"), data = dt.sim,
                     trace = 0, method.inference = "none", neutral.as.uninf = FALSE)
summary(BT.pair2)

## ** Is multiple testing a concern with GPC?

## chunk 46
BTse <- sensitivity(BT.ustat, threshold = seq(0,500, length.out=10),
                          trace = FALSE)
BTse

## chunk 47
BuyseMultComp(BT.H, endpoint = 1:2)

## chunk 48
BuyseMultComp(BT.nH, cumulative = FALSE, endpoint = 1:2)

## chunk 49
BuyseMultComp(list(hierarchical = BT.H, Obrien = BT.nH), cluster = "id")

## chunk 50
BTse.ustat <- sensitivity(BT.ustat, threshold = seq(0,500, length.out=10),
                          band = TRUE, adj.p.value = TRUE, seed = 10, trace = FALSE)
BTse.ustat[,c("time","estimate",
              "lower.ci","upper.ci","p.value",
              "lower.band","upper.band","adj.p.value")]

## chunk 51
library(ggplot2)
autoplot(BTse.ustat)

## chunk 53
BTse.cor <- cor(lava::iid(BTse.ustat))
range(BTse.cor[lower.tri(BTse.cor)])

## chunk 54
BTse.H <- sensitivity(BT.H, trace = FALSE,
                      threshold = list(time = seq(0,500,length = 10), karno = c(0,40,80)))
head(BTse.H)

## chunk 55
grid <- expand.grid(list("time_t20" = seq(0,500,length = 10), "karno" = c(0,40,80)))
cbind(head(grid)," " = "  ...   ",tail(grid))
BTse.H2 <-sensitivity(BT.H, threshold = grid, trace = FALSE)
range(BTse.H-BTse.H2)

## chunk 56
autoplot(BTse.H, col = NA)
##  alternative display:
## autoplot(BTse.H, position  = position_dodge(width = 15))

## * Getting additional inside: looking at the pair level
## ** Extracting the contribution of each pair to the statistic

## chunk 58
form <- trt ~ tte(time, threshold = 20, status = "status") + cont(karno)
BT.keep <- BuyseTest(form,
                     data = veteran, keep.pairScore = TRUE, 
                     trace = 0, method.inference = "none")

## chunk 59
getPairScore(BT.keep, endpoint = 1)

## chunk 60
veteran[c(70,1),]

## chunk 61
getPairScore(BT.keep, endpoint = 1)[, mean(favorable) - mean(unfavorable)]

## chunk 62
BT.keep

## ** Extracting the survival probabilities

## chunk 63
BuyseTest.options(keep.survival = TRUE, precompute = FALSE)
BT.keep2 <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status") + cont(karno),
                      data = veteran, keep.pairScore = TRUE, scoring.rule = "Peron",
                      trace = 0, method.inference = "none")

## chunk 64
outSurv <- getSurvival(BT.keep2, endpoint = 1, strata = 1)
str(outSurv)

## *** Computation of the score with only one censored event

## chunk 65
getPairScore(BT.keep2, endpoint = 1, rm.withinStrata = FALSE)[91]

## chunk 66
veteran[c(22,71),]

## chunk 67
iSurv <- outSurv$survTimeC[22,] 
iSurv

## chunk 68
Sc97 <- iSurv["survivalC_0"] 
Sc97

## chunk 69
iSurv <- outSurv$survTimeT[2,] ## survival at time 112+20
iSurv

## chunk 70
Sc132 <- iSurv["survivalC+threshold"] 
Sc132

## chunk 71
Sc132/Sc97

## *** Computation of the score with two censored events

## chunk 72
head(outSurv$survJumpT)

## chunk 73
getPairScore(BT.keep2, endpoint = 1, rm.withinStrata = FALSE)[148]

## chunk 74
veteran[c(10,72),]

## chunk 75
calcInt <- function(...){ ## no need for the functionnal derivative of the score 
    BuyseTest:::.calcIntegralSurv_cpp(..., 
                                      returnDeriv = FALSE, 
                                      derivSurv = matrix(0), 
                                      derivSurvD = matrix(0))
}

## chunk 76
denom <- as.double(outSurv$survTimeT[3,"survivalT_0"] * outSurv$survTimeC[10,"survivalC_0"])
M <- cbind("favorable" = -calcInt(outSurv$survJumpC, start = 100, 
                                  lastSurv = outSurv$lastSurv[2],
                                  lastdSurv = outSurv$lastSurv[1])/denom,
           "unfavorable" = -calcInt(outSurv$survJumpT, start = 87, 
                                    lastSurv = outSurv$lastSurv[1],
                                    lastdSurv = outSurv$lastSurv[2])/denom)
rownames(M) <- c("lowerBound","upperBound")
M

## chunk 77
outSurv$lastSurv

## * Dealing with missing values or/and right censoring 

## chunk 78
set.seed(10)
dt <- simBuyseTest(1e2, latent = TRUE, argsCont = NULL,
                   argsTTE = list(scale.T = 1/2, scale.C = 1,
                                  scale.censoring.C = 1, scale.censoring.T = 1))
dt[, eventtimeCensoring := NULL]
dt[, status1 := 1]
head(dt)

## chunk 79
100*dt[,mean(status==0)]

## chunk 80
BuyseTest(treatment ~ tte(eventtimeUncensored, status1, threshold = 0.5),
          data = dt,
          scoring.rule = "Gehan", method.inference = "none", trace = 0)

## ** Gehan's scoring rule

## chunk 81
e.G <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
          data = dt, scoring.rule = "Gehan", trace = 0)
model.tables(e.G)

## ** Peron's scoring rule

## chunk 82
e.P <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
          data = dt, scoring.rule = "Peron", trace = 0)
model.tables(e.P)

## chunk 83
dt[,.SD[which.max(eventtime)],by="treatment"]

## chunk 84
library(prodlim)
e.prodlim <- prodlim(Hist(eventtime, status) ~ treatment, data = dt)

## chunk 85
e.P1 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                  model.tte = e.prodlim,
                  data = dt, scoring.rule = "Peron", trace = 0)
model.tables(e.P1)

## chunk 86
attr(e.prodlim, "iidNuisance") <- TRUE
e.P2 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                  model.tte = e.prodlim,
                  data = dt, scoring.rule = "Peron", trace = 0)
model.tables(e.P2)

## chunk 87
library(survival)
e.survreg <- survreg(Surv(eventtime, status) ~ treatment, data = dt, 
                     dist = "weibull")
attr(e.survreg, "iidNuisance") <- TRUE

## chunk 88
e.P3 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                  model.tte = e.survreg,
                  data = dt, scoring.rule = "Peron", trace = 0)
model.tables(e.P3)

## chunk 89
e.TTEM <- BuyseTTEM(e.survreg, treatment = "treatment", iid = TRUE, n.grid = 2500)
attr(e.TTEM, "iidNuisance") <- TRUE
str(e.TTEM$peron$jumpSurvHaz[[1]][[1]])

## chunk 90
e.P4 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                  model.tte = e.TTEM,
                  data = dt, scoring.rule = "Peron", trace = 0)
model.tables(e.P4)

## ** Correction via inverse probability-of-censoring weights (IPCW)

## chunk 91
BT <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                data = dt, keep.pairScore = TRUE, trace = 0,
                scoring.rule = "Gehan", method.inference = "none", correction.uninf = 2)
summary(BT)

## chunk 92
iScore <- getPairScore(BT, endpoint = 1)

## chunk 93
iScore[,.SD[1], 
       .SDcols = c("favorableC","unfavorableC","neutralC","uninfC"),
       by = c("favorable","unfavorable","neutral","uninf")]

## chunk 94
iScore[,1/mean(favorable + unfavorable + neutral)]

## ** Correction at the pair level

## chunk 95
BT <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                data = dt, keep.pairScore = TRUE, trace = 0,
                scoring.rule = "Gehan", method.inference = "none", correction.uninf = TRUE)
summary(BT)

## chunk 96
iScore <- getPairScore(BT, endpoint = 1)
iScore[,.SD[1], 
       .SDcols = c("favorableC","unfavorableC","neutralC","uninfC"),
       by = c("favorable","unfavorable","neutral","uninf")]

## chunk 97
iScore[, .(favorable = weighted.mean(favorable, w = 1-uninf), 
           unfavorable = weighted.mean(unfavorable, w = 1-uninf), 
           neutral = weighted.mean(neutral, w = 1-uninf))]

## ** Note on the use of the corrections

## chunk 98
set.seed(10);n <- 250; 
df <- rbind(data.frame(group = "T1", time = rweibull(n, shape = 1, scale = 2), status = 1),
            data.frame(group = "T2", time = rweibull(n, shape = 2, scale = 1.8), status = 1))
df$censoring <- runif(NROW(df),0,2)
df$timeC <- pmin(df$time,df$censoring)
df$statusC <- as.numeric(df$time<=df$censoring)
plot(prodlim(Hist(time,status)~group, data = df)); title("complete data");
plot(prodlim(Hist(timeC,statusC)~group, data = df)); title("right-censored data");

## chunk 100
BuyseTest.options(method.inference = "none")
e.ref <- BuyseTest(group ~ tte(time,status), data = df, trace = FALSE)
s.ref <- model.tables(e.ref, column = c("favorable","unfavorable","neutral","uninf","Delta"))
s.ref

## chunk 101
e.correction <- BuyseTest(group ~ tte(timeC,statusC), data = df, trace = FALSE, correction.uninf = TRUE)
s.correction <- model.tables(e.correction, column = c("favorable","unfavorable","neutral","uninf","Delta"))

## chunk 102
e.Peron <- BuyseTest(group ~ tte(timeC,statusC), data = df, trace = FALSE)
s.Peron <- model.tables(e.Peron, column = c("favorable","unfavorable","neutral","uninf","Delta"))
rbind("reference" = s.ref,
      "no correction" = s.Peron,
      "correction" = s.correction)

## * Simulating data using =simBuyseTest=

## chunk 103
set.seed(10)
simBuyseTest(n.T = 5, n.C = 5)

## chunk 104
set.seed(10)
argsCont <- list(mu.T = c(5,5), mu.C = c(10,10), 
                 sigma.T = c(1,1), sigma.C = c(1,1),
                 name = c("tumorSize","score"))
dt <- simBuyseTest(n.T = 5, n.C = 5,
                   argsCont = argsCont)
dt

## * Power calculation using =powerBuyseTest=

## chunk 105
simFCT <- function(n.C, n.T){
     out <- rbind(cbind(Y=stats::rt(n.C, df = 5), group=0),
                  cbind(Y=stats::rt(n.T, df = 5) + 1/2, group=1))
     return(data.table::as.data.table(out))
}
set.seed(10)
simFCT(101,101)

## chunk 106
null <- c("netBenefit" = 0)

## chunk 107
powerW <- powerBuyseTest(sim = simFCT, method.inference = "u-statistic", null = null,
                         sample.size = c(5,10,20,30,50,100),                         
                         formula = group ~ cont(Y), 
                         n.rep = 1000, seed = 10, cpus = 6, trace = 0)

## chunk 108
summary(powerW)

## chunk 109
nW <- powerBuyseTest(sim = simFCT, method.inference = "u-statistic", 
                     power = 0.8, max.sample.size = 1000,                     
                     formula = group ~ cont(Y), null = c("netBenefit" = 0),
                     n.rep = c(1000,10), seed = 10, cpus = 5, trace = 0)

## chunk 110
summary(nW)

## * Modifying default options

## chunk 111
BuyseTest.options("trace")

## chunk 112
BuyseTest.options(trace = 0)

## chunk 113
BuyseTest.options(summary.display = list(c("endpoint","threshold","delta","Delta","information(%)")))
summary(BT)

## chunk 114
BuyseTest.options(reinitialise = TRUE)

## * References