library(data.table)
library(prodlim)
## run BuyseTest
library(survival) ## import veteran
BT.keep <- BuyseTest(trt ~ tte(time, threshold = 20, status = "status") + cont(karno),
data = veteran, keep.pairScore = TRUE,
trace = 0, method.inference = "none")
## Extract scores
pScore <- getPairScore(BT.keep, endpoint = 1)
## look at one pair
indexPair <- intersect(which(pScore$index.1 == 22),
which(pScore$index.2 == 71))
pScore[indexPair]
## retrive pair in the original dataset
pVeteran <- veteran[pScore[indexPair,c(index.1,index.2)],]
pVeteran
## the observation from the control group is censored at 97
## the observation from the treatment group has an event at 112
## since the threshold is 20, and (112-20)<97
## we know that the pair is not in favor of the treatment
## the formula for probability in favor of the control is
## Sc(97)/Sc(112+20)
## where Sc(t) is the survival at time t in the control arm.
## we first estimate the survival in each arm
e.KM <- prodlim(Hist(time,status)~trt, data = veteran)
## and compute the survival
iSurv <- predict(e.KM, times = c(97,112+20),
newdata = data.frame(trt = 1, stringsAsFactors = FALSE))[[1]]
## the probability in favor of the control is then
pUF <- iSurv[2]/iSurv[1]
pUF
## and the complement to one of that is the probability of being neutral
pN <- 1 - pUF
pN
if(require(testthat)){
testthat::expect_equal(pUF, pScore[indexPair, unfavorable])
testthat::expect_equal(pN, pScore[indexPair, neutral])
}
Run the code above in your browser using DataLab