%\VignetteIndexEntry{BuyseTest: overview}
%\VignetteEngine{R.rsp::asis}
%\VignetteKeyword{PDF}
%\VignetteKeyword{vignette}
%\VignetteKeyword{package}

## 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.halfboot <- 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.halfboot, statistic = "netBenefit"),
               winRatio = confint(BT.halfboot, statistic = "winRatio"),
               favorable = confint(BT.halfboot, statistic = "favorable"))
Mstat

## ** Stratified GPC

## chunk 21
ffstrata <- trt ~ tte(time, threshold = 20, status = "status") + strata(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)

## ** Standardization

## chunk 30
BTstd <- BuyseTest(ffstrata, data = veteran, trace = 0, pool.strata = "standardization")
model.tables(BTstd)[,c("strata","total","delta","Delta","lower.ci","upper.ci","p.value")]

## chunk 31
BuyseTest(trt ~ tte(time, threshold = 20, status = "status"),
          data = rbind(veteran[veteran$celltype == "smallcell" & veteran$trt == "Pl",],
                       veteran[veteran$celltype == "squamous" & veteran$trt == "Exp",]),
          trace = 0)

## chunk 32
BTstd.boot <- BuyseTest(ffstrata, data = veteran, trace = 0, pool.strata = "standardization",
                        method.inference = "bootstrap", n.resampling = 100, seed = 10)
confint(BTstd.boot, strata = TRUE)[1:6,]

## ** Using multiple endpoints (hierarchical)

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

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


## ** Using multiple endpoints (non-hierarchical)

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

## chunk 37
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 38
confint(BuyseTest(trt ~ cont(karno, threshold = 0), data = veteran, trace = 0))

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

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

## ** Statistical inference

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

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

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

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

## ** What if smaller is better?

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

## ** What about a multiplicative instead of additive threshold?

## chunk 46
ffstar <- trt ~ tte(time, status = "status", threshold = 20, operator = "*")
BTstar <- BuyseTest(ffstar, data = veteran, trace = 0)
print(BTstar)

## ** Stopping comparison for neutral pairs

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

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

## chunk 50
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 51
BTse <- sensitivity(BT.ustat, threshold = seq(0,500, length.out=10),
                          trace = FALSE)
BTse

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

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

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

## chunk 55
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 56
library(ggplot2)
autoplot(BTse.ustat)

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

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

## chunk 60
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 61
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 (long format)

## chunk 63
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 64
getPairScore(BT.keep, endpoint = 1)

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

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

## chunk 67
BT.keep

## ** Extracting the contribution of each pair to the statistic (matrix format)

## chunk 68
Mscore <- getMatrixScore(BT.keep, endpoint = 1)
dim(Mscore)
Mscore[1:4,1:4]

## chunk 69
MscoreAll <- getMatrixScore(BT.keep, endpoint = 1, within = TRUE)
dim(MscoreAll)

## chunk 70
rbind(cbind(MscoreAll[1:4,1:4],NA,MscoreAll[70:73,1:4]),
      NA,
      cbind(MscoreAll[1:4,70:73],NA,MscoreAll[70:73,70:73]))

## ** Extracting the survival probabilities

## chunk 71
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 72
outSurv <- getSurvival(BT.keep2, endpoint = 1, strata = 1)
str(outSurv)

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

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

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

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

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

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

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

## chunk 79
Sc132/Sc97

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

## chunk 80
head(outSurv$survJumpT)

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

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

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

## chunk 84
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 85
outSurv$lastSurv

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

## chunk 86
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 87
100*dt[,mean(status==0)]

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

## chunk 89
names(survival::survreg.distributions)

## ** Gehan's scoring rule

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

## ** Latta's scoring rule

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

## chunk 92
dt[,.SD[which.max(eventtime)]]

## chunk 93
e.prodlimLatta <- prodlim(Hist(eventtime, status) ~ 1, data = dt)
attr(e.prodlimLatta,"iidNuisance") <- TRUE

## chunk 94
e.LattaManual <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                           model.tte = e.prodlimLatta,
                           data = dt, scoring.rule = "Peron", trace = 0)
confint(e.LattaManual)

## ** Peron's scoring rule

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

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

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

## chunk 98
e.PeronManual <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                           model.tte = e.prodlimPeron,
                           data = dt, scoring.rule = "Peron", trace = 0)
confint(e.PeronManual)

## chunk 99
dt2 <- dt[order(dt$eventtime)]
e.PeronManual2 <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                            model.tte = prodlim(Hist(eventtime, status) ~ treatment, data = dt2),
                            data = dt, scoring.rule = "Peron", trace = 0)
confint(e.PeronManual2)

## ** Efron's scoring rule

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

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

## chunk 102
dtEfron <- copy(dt)
dtEfron[treatment=="T" & id == 154, c("eventtime","status") := .(eventtime+1e-10,1)]
## dt2[treatment=="T" & eventtime > 1.58865, status := 1]

## chunk 103
e.prodlimEfron <- prodlim(Hist(eventtime, status) ~ treatment, data = dtEfron)
attr(e.prodlimEfron,"iidNuisance") <- TRUE

## chunk 104
e.Efron <- BuyseTest(treatment ~ tte(eventtime, status, threshold = 0.5),
                     model.tte = e.prodlimEfron,
                     data = dt, scoring.rule = "Peron", trace = 0)
confint(e.Efron)

## ** Scoring rule based on a parametric survival model

## chunk 105
names(survreg.distributions)

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

## chunk 108
library(survival)
e.survregWeibull <- survreg(Surv(eventtime, status) ~ treatment, data = dt, 
                            dist = "weibull")

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

## ** Correction via re-weighting citep:peron2021correcting
## *** At the trial level

## chunk 110
BT.trialW <- 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.trialW)

## chunk 111
iScore <- getPairScore(BT.trialW, endpoint = 1)

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

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

## *** At the pair level

## chunk 114
BT.pairI <- 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.pairI)

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

## chunk 116
c(favorable = unname(coef(BT.pairI, statistic = "favorable")),
  unfavorable = unname(coef(BT.pairI, statistic = "unfavorable")),
  neutral = unname(coef(BT.pairI, statistic = "neutral")))

## *** Note on the use of the corrections

## chunk 117
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 119
e.ref <- BuyseTest(group ~ tte(time,status), data = df, method.inference = "none",
                   trace = FALSE)
s.ref <- model.tables(e.ref, column = c("favorable","unfavorable","neutral","uninf","Delta"))
s.ref

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

## chunk 121
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 122
set.seed(10)
simBuyseTest(n.T = 5, n.C = 5)

## chunk 123
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 124
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 125
null <- c("netBenefit" = 0)

## chunk 126
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 = 10, trace = 0)

## chunk 127
summary(powerW)

## chunk 128
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 129
summary(nW)

## * Modifying default options

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

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

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

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

## * References