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

## * Paired design

## chunk 2
data(diabetic, package = "survival")
head(diabetic)

## chunk 3
diabetic$juvenile <- diabetic$age <= 19
library(LMMstar)
summarize(age ~ juvenile, data = diabetic[!duplicated(diabetic$id),])

## chunk 4
diabeticJ <- diabetic[diabetic$juvenile,]

## ** Wald methods (Gehan scoring rule)

## chunk 5
library(BuyseTest)
e.BTjuv <- BuyseTest(trt ~ tte(time,status) + strata(id, match = TRUE), 
                     data = diabeticJ, trace = FALSE,
                     scoring.rule =  "Gehan")
model.tables(e.BTjuv, percentage = FALSE)

## chunk 6
mean(coef(e.BTjuv, strata = TRUE))

## chunk 7
N <- nobs(e.BTjuv)["pairs"]
pw <- coef(e.BTjuv, statistic = "favorable")
pl <- coef(e.BTjuv, statistic = "unfavorable")
sqrt((pw + pl - (pw - pl)^2)/N)

## chunk 8
confint(e.BTjuv)

## chunk 9
confint(e.BTjuv, transform = FALSE)

## chunk 10
sqrt(var(coef(e.BTjuv, strata = TRUE))/N)

## chunk 11
sqrt(var(coef(e.BTjuv, strata = TRUE))/N) * sqrt((N-1)/N)

## ** MOVER method (Gehan scoring rule)

## chunk 12
BuyseTest:::mover(e.BTjuv)

## ** Wald methods (Peron scoring rule)

## chunk 13
library(prodlim)
e.BTjuv2 <- BuyseTest(trt ~ tte(time,status) + strata(id, match = TRUE), 
                      data = diabeticJ, trace = FALSE,
                      model.tte = prodlim(Hist(time,status)~ trt, data = diabeticJ))
model.tables(e.BTjuv2, percentage = FALSE)

## chunk 14
c(sqrt(var(coef(e.BTjuv2, strata = TRUE))/N),
  sqrt(var(coef(e.BTjuv2, strata = TRUE))/N) * sqrt((N-1)/N)
  )

## chunk 15
model.tte <- prodlim(Hist(time,status)~ trt, data = diabeticJ)
attr(model.tte, "iidNuisance") <- FALSE
confint(BuyseTest(trt ~ tte(time,status) + strata(id, match = TRUE), 
                  data = diabeticJ, trace = FALSE,
                  model.tte = model.tte))

## chunk 16
pw.peron <- coef(e.BTjuv2, statistic = "favorable")
pl.peron <- coef(e.BTjuv2, statistic = "unfavorable")
sqrt((pw.peron + pl.peron - (pw.peron - pl.peron)^2)/N)

## chunk 17
confint(e.BTjuv2)

## chunk 18
pw.peronS <- coef(e.BTjuv2, statistic = "favorable", strata = TRUE)
pl.peronS <- coef(e.BTjuv2, statistic = "unfavorable", strata = TRUE)
Hterm1 <- (pw.peronS - pl.peronS) - (pw.peron - pl.peron)

## chunk 19
Hterm2.obs <- e.BTjuv2@iidNuisance$favorable - e.BTjuv2@iidNuisance$unfavorable
Hterm2.pair <- Hterm2.obs[diabeticJ$trt==0]
table(Hterm2.obs[diabeticJ$trt==1])

## chunk 20
c(average = sqrt(crossprod((Hterm1/N))),
  nuisance = sqrt(crossprod((Hterm2.pair/N))),
  all = sqrt(crossprod((Hterm1/N + Hterm2.pair/N))))

## * Multiple cross-over design

## chunk 21
data(profil, package = "BuyseTest")
profil <- profil[order(profil$id),]
profil[profil$id==1 & profil$period==1,]

## chunk 22
profil.L <- reshape(profil, direction = "long", idvar = c("id","time"),
                    varying = c("rcs","number","duration"), v.names = "value",
                    timevar = "outcome", times = c("rcs","number","duration"))

## chunk 23
ggRCS <- ggplot(profil.L[profil.L$id %in% 1:5,],
                aes(x = time, group = id, color = treatment, y = value))
ggRCS <- ggRCS + geom_point() + geom_line() 
ggRCS <- ggRCS + facet_grid(outcome~id, labeller = label_both, scales = "free_y")
ggRCS <- ggRCS + labs(y="")
ggRCS

## ** With 2 treatment groups

## chunk 25
lowProfil <- profil[profil$treatment %in% c("placebo","lowDose"),]
lowProfil$treatment <- droplevels(lowProfil$treatment)

## chunk 26
fff <- treatment ~ cont(rcs, threshold = 1.45, operator = "<0") + cont(number, threshold = 0.35, operator = "<0") + cont(duration, threshold = 3, operator = "<0")

## chunk 27
ls.lowGPC <- by(lowProfil, lowProfil$id, BuyseTest, formula = fff, trace = FALSE)
df.lowGPC <- data.frame(do.call(rbind,lapply(ls.lowGPC, coef)),
                        do.call(rbind,lapply(ls.lowGPC, nobs)),
                        Buyse = as.numeric(NA), CMH = as.numeric(NA))
df.lowGPC$Buyse <- with(df.lowGPC, pairs/sum(pairs))
df.lowGPC$CMH <- with(df.lowGPC, (pairs/(placebo+lowDose))/sum(pairs/(placebo+lowDose)))
head(df.lowGPC)

## chunk 28
rbind(average = colMeans(df.lowGPC[,1:3]),
      Buyse = apply(df.lowGPC[,1:3], 2, weighted.mean, w = df.lowGPC$Buyse),
      CMH = apply(df.lowGPC[,1:3], 2, weighted.mean, w = df.lowGPC$CMH))

## chunk 29
lowGPC <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                    data = lowProfil, trace = FALSE)
confint(lowGPC)

## chunk 30
lowGPC_Buyse <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                          data = lowProfil, pool.strata = "Buyse", trace = FALSE)
confint(lowGPC_Buyse)

## chunk 31
lowGPC_CMH <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                        data = lowProfil, pool.strata = "CMH", trace = FALSE)
confint(lowGPC_CMH)

## chunk 32
sqrt(apply(df.lowGPC[,1:3],2,var)/NROW(df.lowGPC))

## chunk 33
BuyseTest.options(order.Hprojection = 2)
lowGPC2 <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                     data = lowProfil, trace = FALSE)
confint(lowGPC2)

## chunk 34
df.lowGPC_center <- sweep(df.lowGPC[,1:3], MARGIN = 2, FUN = "-", STATS = coef(lowGPC_Buyse))
df.lowGPC_Wcenter <- sweep(df.lowGPC_center, MARGIN = 1, FUN = "*", STATS = df.lowGPC$Buyse)
sqrt(colSums(df.lowGPC_Wcenter^2))

## ** Accounting for time (WORK IN PROGRESS!)

## chunk 35
ls.lowGPC_period <- by(lowProfil, lowProfil$id, BuyseTest,
                       formula = update(fff,.~.+strata(period)), trace = FALSE)

## chunk 36
rbind(noStrata = nobs(ls.lowGPC[[1]]),
      PeriodStrata = nobs(ls.lowGPC_period[[1]]))

## chunk 37
nobs(ls.lowGPC_period[[1]],strata=TRUE)

## chunk 38
df.lowGPC_period <- data.frame(do.call(rbind,lapply(ls.lowGPC_period, coef)),
                               do.call(rbind,lapply(ls.lowGPC_period, nobs)),
                               Buyse = as.numeric(NA), CMH = as.numeric(NA))
df.lowGPC_period$Buyse <- with(df.lowGPC_period, pairs/sum(pairs))
df.lowGPC_period$CMH <- with(df.lowGPC_period, (pairs/(placebo+lowDose))/sum(pairs/(placebo+lowDose)))

## chunk 39
rbind(average = colMeans(df.lowGPC_period[,1:3]),
      Buyse = apply(df.lowGPC_period[,1:3], 2, weighted.mean, w = df.lowGPC_period$Buyse),
      CMH = apply(df.lowGPC_period[,1:3], 2, weighted.mean, w = df.lowGPC_period$CMH))

## ** With 3 or more treatment groups (WORK IN PROGRESS!) 

## chunk 40
CasinoTest(update(fff, . ~ . + strata(id)), data = profil,
           method.inference = "none", type = "unweighted")
## do not trust inference since it would require accounting for matching

## * References