#+TITLE: Performing GPC with repeated measurements
#+Author: Brice Ozenne
#+BEGIN_SRC R :exports none :results quiet :session *R* :cache no
options(width = 95)
if(system("whoami",intern=TRUE)=="bozenne"){  
  setwd("~/Documents/GitHub/BuyseTest/inst/doc-paired/")
}else if(system("whoami",intern=TRUE)=="sund\\hpl802"){  
  setwd("c:/Users/hpl802/Documents/Github/BuyseTest/inst/doc-paired/")
}
library(survival) ## avoid messages when loading the package later on
library(ggplot2) ## avoid messages when loading the package later on
library(prodlim) ## avoid messages when loading the package later on
#+END_SRC

#+RESULTS:

This vignette describes how to use Generalized Pairwise comparisons
(GPC) in a design where each subject undergo several
measurements. Traditional GPC assumes independent observations which
is likely not the case in such design since observations from the same
subject are typically correlated. This is an area of development for
the BuyseTest package so the content of this vignette may change in
newer versions.

* Paired design

The Diabetic Retinopathy Study (DRS), which data is available from
 \Rlogo package survival, included 197 patients who had one of their
 eye randomized to laser treatment while the other did not receive any
 treatment:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
data(diabetic, package = "survival")
head(diabetic)
#+END_SRC

#+RESULTS:
:   id laser age   eye trt risk  time status
: 1  5 argon  28  left   0    9 46.23      0
: 2  5 argon  28 right   1    9 46.23      0
: 3 14 xenon  12  left   1    8 42.50      0
: 4 14 xenon  12 right   0    6 31.30      1
: 5 16 xenon   9  left   1   11 42.27      0
: 6 16 xenon   9 right   0   11 42.27      0

The outcome was time to blindness (visual acuity drop below a certain
threshold). In the real study =status= equal to 0 mixes death and
censoring (due to drop-out or end of study) but this complication will
be neglected here for simplicity.


\bigskip

We will replicate some of the analyzes presented in
cite:matsouaka2022robust. In this paper they split the dataset into
juvenile and adult patients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
diabetic$juvenile <- diabetic$age <= 19
library(LMMstar)
summarize(age ~ juvenile, data = diabetic[!duplicated(diabetic$id),])
#+END_SRC

#+RESULTS:
:   juvenile observed missing     mean        sd min q1 median    q3 max
: 1    FALSE       83       0 35.30120 11.242054  20 25     34 45.00  58
: 2     TRUE      114       0 10.21053  4.713892   1  7     10 13.75  19

and we will focus on the juvenile patients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
diabeticJ <- diabetic[diabetic$juvenile,]
#+END_SRC

#+RESULTS:

** Wald methods (Gehan scoring rule)

To mimic the methodology underlying the results presented in Table 1
of cite:matsouaka2022robust, we perform GPC stratified by patient
using the Gehan scoring rule:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)
#+END_SRC

#+RESULTS:
:   endpoint total favorable unfavorable neutral uninf     Delta   lower.ci  upper.ci    p.value
: 1     time   114        39          21       3    51 0.1578947 0.02591623 0.2844633 0.01922741

Indeed this scoring rule does not involve any extra-modeling, only
evaluating the patient specific net benefit and averaging them:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
mean(coef(e.BTjuv, strata = TRUE))
#+END_SRC

#+RESULTS:
: [1] 0.1578947

cite:matsouaka2022robust propose to estimate the standard error as:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
N <- nobs(e.BTjuv)["pairs"]
pw <- coef(e.BTjuv, statistic = "favorable")
pl <- coef(e.BTjuv, statistic = "unfavorable")
sqrt((pw + pl - (pw - pl)^2)/N)
#+END_SRC

#+RESULTS:
:       time 
: 0.06631828

which matches what =BuyseTest= output:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
confint(e.BTjuv)
#+END_SRC

#+RESULTS:
:       estimate         se   lower.ci  upper.ci null    p.value
: time 0.1578947 0.06631828 0.02591623 0.2844633    0 0.01922741

By default =confint= uses a hyperbolic tangent to compute confidence
intervals (CIs), which will then differ from the 'Wald' discussed in
cite:matsouaka2022robust. These 'untransformed Wald' CIs can be
retrieved by setting the argument =transform= to =FALSE=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
confint(e.BTjuv, transform = FALSE)
#+END_SRC

#+RESULTS:
:       estimate         se   lower.ci  upper.ci null    p.value
: time 0.1578947 0.06631828 0.02791329 0.2878762    0 0.01727214

\clearpage

_Note:_ naively one may think to estimate the standard error as:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sqrt(var(coef(e.BTjuv, strata = TRUE))/N)
#+END_SRC

#+RESULTS:
:      pairs 
: 0.06661108

This is equivalent (in large samples to the previous formula). Indeed:
#+BEGIN_EXPORT latex
\begin{align*}
&\Prob[X>Y] + \Prob[Y>X] - (\Prob[X>Y] - \Prob[Y>X])^2 \\
=& \Prob[X>Y] + \Prob[Y>X] - \Prob[X>Y]^ - \Prob[Y>X]^2 + 2 \Prob[X>Y] \Prob[Y>X] \\
=& \Prob[X>Y](1-\Prob[X>Y]) + \Prob[Y>X](1-\Prob[Y>X]) + 2 \Prob[X>Y] \Prob[Y>X] \\
=& \Prob[X>Y](1-\Prob[X>Y]) + \Prob[Y>X](1-\Prob[Y>X]) \\
 & - 2 (0 - \Prob[X>Y] \Prob[Y>X] - \Prob[X>Y] \Prob[Y>X] + \Prob[X>Y] \Prob[Y>X] \\
=& \Var\left[\Ind[X>Y]\right] + \Var\left[\Ind[X<Y]\right] - 2 \Cov\left(\Ind[X>Y],\Ind[X<Y]\right) \\
=& \Var\left[\Ind[X>Y]-\Ind[X<Y]\right] \\
\end{align*}
#+END_EXPORT

There is only a factor =N/(N-1)= difference between the two:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sqrt(var(coef(e.BTjuv, strata = TRUE))/N) * sqrt((N-1)/N)
#+END_SRC

#+RESULTS:
:      pairs 
: 0.06631828

** MOVER method (Gehan scoring rule)

The method recommended by cite:matsouaka2022robust is the MOVER
approach, which has been developped for a binary scoring rule
(e.g. Gehan). An experimental function with the same name has been
implemented in the BuyseTest package:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
BuyseTest:::mover(e.BTjuv)
#+END_SRC
#+RESULTS:
:   estimate      lower      upper     pvalue 
: 0.15789474 0.02540421 0.28317729 0.01967878

leading to the same results as those of the table 1 in the original article, up to rounding.

\clearpage

** Wald methods (Peron scoring rule)

To better account for censoring one could use the Peron scoring rule
where the survival is estimated across all subjects within a treatment
group. One has to specify the survival model, otherwise, the BuyseTest
function will estimate a treatmnet and strata specific survival curve
which is impossible to perform here. The =model.tte= argument can be
used to specify such survival model:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)
#+END_SRC

#+RESULTS:
:   endpoint total favorable unfavorable neutral    uninf    Delta   lower.ci  upper.ci
: 1     time   114  47.36525    24.29552       3 39.33923 0.202366 0.05045454 0.3451254
:       p.value
: 1 0.009329589

Ignoring the uncertainty of the survival model, the standard would be:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
c(sqrt(var(coef(e.BTjuv2, strata = TRUE))/N),
  sqrt(var(coef(e.BTjuv2, strata = TRUE))/N) * sqrt((N-1)/N)
  )
#+END_SRC

#+RESULTS:
:      pairs      pairs 
: 0.06595510 0.06566518

depending on whether a small sample correction is used or not. This
matches the output of =BuyseTest= when ignoring the uncertainty of the
survival model:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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))
#+END_SRC

#+RESULTS:
:      estimate         se   lower.ci  upper.ci null     p.value
: time 0.202366 0.06566518 0.07088227 0.3269375    0 0.002726979

\Warning Because the pairwise win and loss score are no more binary, the
previous formula no more simplifies into the formula presented in
cite:matsouaka2022robust:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)
#+END_SRC

#+RESULTS:
:       time 
: 0.07179718

\clearpage 

To account for the uncertainty of the survival model, =BuyseTest=
outputs a slightly higher standard error:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
confint(e.BTjuv2)
#+END_SRC

#+RESULTS:
:      estimate         se   lower.ci  upper.ci null     p.value
: time 0.202366 0.07569815 0.05045454 0.3451254    0 0.009329589

This is achieved by considering two sources of uncertainty:
- average of a finite number of pairs:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)
#+END_SRC

#+RESULTS:

- propage the uncertainty of the survival model to the net
  benefit. Because each pair appear twice (control and treatment) the
  impact of removing a pair on the net benefit is stored in the
  control and the treated is set to 0:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
Hterm2.obs <- e.BTjuv2@iidNuisance$favorable - e.BTjuv2@iidNuisance$unfavorable
Hterm2.pair <- Hterm2.obs[diabeticJ$trt==0]
table(Hterm2.obs[diabeticJ$trt==1])
#+END_SRC  

#+RESULTS:
: 
:   0 
: 114

After rescaling the terms by a factor N (number of pairs, to account
for the pooling) we retrieve the uncertainty output by =BuyseTest=:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
c(average = sqrt(crossprod((Hterm1/N))),
  nuisance = sqrt(crossprod((Hterm2.pair/N))),
  all = sqrt(crossprod((Hterm1/N + Hterm2.pair/N))))
#+END_SRC

#+RESULTS:
:    average   nuisance        all 
: 0.06566518 0.02084622 0.07569815

\clearpage

* Multiple cross-over design
# Can be compared to https://onlinelibrary.wiley.com/doi/10.1002/sim.9648

A more complex design is a cross-over where each patient may
repeteadly experience each treatment. As an example, we will consider
the PROFIL trial whose dataset is available in the BuyseTest package:
#+BEGIN_SRC R :exports both :results none :session *R* :cache no
data(profil, package = "BuyseTest")
profil <- profil[order(profil$id),]
profil[profil$id==1 & profil$period==1,]
#+END_SRC

#+RESULTS:
#+begin_example
   id age male period time treatment rcs number duration
1   1  23    0      1    1  highDose 1.8      2       16
2   1  23    0      1    2  highDose 2.3      2       26
3   1  23    0      1    3  highDose 0.0      0        0
4   1  23    0      1    4  highDose 0.0      0        0
5   1  23    0      1    5  highDose 2.7      1       16
6   1  23    0      1    6  highDose 1.9      1       10
7   1  23    0      1    7  highDose 0.0      0        0
8   1  23    0      1    8   placebo 1.5      1       11
9   1  23    0      1    9   placebo 2.9      2       39
10  1  23    0      1   10   placebo 2.3      1       22
11  1  23    0      1   11   placebo 0.0      0        0
12  1  23    0      1   12   placebo 0.0      0        0
13  1  23    0      1   13   placebo 2.4      1       20
14  1  23    0      1   14   placebo 1.9      1        8
15  1  23    0      1   15   lowDose 3.1      1       13
16  1  23    0      1   16   lowDose 3.0      2      161
17  1  23    0      1   17   lowDose 0.0      0        0
18  1  23    0      1   18   lowDose 0.0      0        0
19  1  23    0      1   19   lowDose 2.4      1       35
20  1  23    0      1   20   lowDose 0.0      0        0
21  1  23    0      1   21   lowDose 0.0      0        0
#+end_example

The software output displays the information of the first patient
relative to the first period (out of 4) during which the patient is
sequentially assigned one of three treatments and her outcomes (=rcs=,
=number=, and =duration=) are monitored daily. To obtain a graphical
display of the outcomes over time we first reshape the data:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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"))
#+END_SRC

#+RESULTS:
and make a spaghetti plot for the first 5 patients:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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
#+END_SRC

#+RESULTS:

#+BEGIN_SRC R :exports none :results output :session *R* :cache no
ggsave(ggRCS + theme(text = element_text(size=15), 
                     axis.line = element_line(linewidth = 1.25),
                     axis.ticks = element_line(linewidth = 2),
                     axis.ticks.length=unit(.25, "cm"),
                     legend.key.size = unit(2,"line"),
                     legend.position="bottom",
                     legend.direction = "horizontal"),
       filename = file.path("inst","doc-paired","figures","spaghetti-Nof1.pdf"),
       width = 12, height = 7)
#+END_SRC

#+RESULTS:

#+ATTR_LaTeX: :width 1\textwidth :options trim={0 0 0 0} :placement [!h]
[[./figures/spaghetti-Nof1.pdf]]

** With 2 treatment groups

\noindent Since =BuyseTest= can only handle two treatment group, we restrict the
dataset to placebo and low dose:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
lowProfil <- profil[profil$treatment %in% c("placebo","lowDose"),]
lowProfil$treatment <- droplevels(lowProfil$treatment)
#+END_SRC

#+RESULTS:

# From Joris paper on patient preferences and multiple treatments into consideration for individualized Net Benefit 
# Thresholds were taken from either latest meta-analysis or large trial in RP on the daily number of attacks, their duration and the RCS VAS[23,24].
# They were set to 0.35 attacks per day, 3 minutes, and 1.45 points respectively.
We will use the following hierarchy and threshold of clinical
relevance:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
fff <- treatment ~ cont(rcs, threshold = 1.45, operator = "<0") + cont(number, threshold = 0.35, operator = "<0") + cont(duration, threshold = 3, operator = "<0")
#+END_SRC

#+RESULTS:



\noindent One could either run a separate GPC for each patient:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)
#+END_SRC

#+RESULTS:
:      rcs_t1.45 number_t0.35  duration_t3 placebo lowDose pairs      Buyse       CMH
: 1 -0.016581633 -0.015306122 -0.021683673      28      28   784 0.04694049 0.0364368
: 2  0.000000000  0.153061224  0.183673469      28      28   784 0.04694049 0.0364368
: 3 -0.001632653  0.003265306 -0.008979592      35      35  1225 0.07334451 0.0455460
: 4  0.117346939  0.225765306  0.154336735      28      28   784 0.04694049 0.0364368
: 5 -0.043083900 -0.047619048 -0.029478458      21      21   441 0.02640402 0.0273276
: 6  0.102040816  0.092970522  0.077097506      21      21   441 0.02640402 0.0273276

\noindent and pool the patient-specific Net Treatment
Benefits. Different weighting scheme are possible, e.g. same weight
for all patients, weight dependent on the number of blocks experienced
by the patient:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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))
#+END_SRC

#+RESULTS:
:          rcs_t1.45 number_t0.35 duration_t3
: average 0.02741742   0.02755903  0.03397497
: Buyse   0.01628547   0.03730092  0.04215064
: CMH     0.02145018   0.03266743  0.03869945

This can be replicated using a single call to =BuyseTest= specifying a
strata in the formula:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
lowGPC <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                    data = lowProfil, trace = FALSE)
confint(lowGPC)
#+END_SRC

#+RESULTS:
:                estimate         se    lower.ci   upper.ci null   p.value
: rcs_t1.45    0.02741742 0.02047690 -0.01273920 0.06748574    0 0.1808076
: number_t0.35 0.02755903 0.03139999 -0.03401048 0.08892015    0 0.3803606
: duration_t3  0.03397497 0.02801978 -0.02099009 0.08873527    0 0.2256649

By default, equal weights are given to each patient but other
weighting schemes can be used by specifying the =pool.strata= argument:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
lowGPC_Buyse <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                          data = lowProfil, pool.strata = "Buyse", trace = FALSE)
confint(lowGPC_Buyse)
#+END_SRC

#+RESULTS:
:                estimate         se     lower.ci   upper.ci null    p.value
: rcs_t1.45    0.01628547 0.01771680 -0.018444486 0.05097618    0 0.35807018
: number_t0.35 0.03730092 0.02668638 -0.015057839 0.08945568    0 0.16257765
: duration_t3  0.04215064 0.02400508 -0.004957166 0.08907178    0 0.07946061

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
lowGPC_CMH <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                        data = lowProfil, pool.strata = "CMH", trace = FALSE)
confint(lowGPC_CMH)
#+END_SRC

#+RESULTS:
:                estimate         se    lower.ci   upper.ci null   p.value
: rcs_t1.45    0.02145018 0.01855984 -0.01493878 0.05778241    0 0.2479363
: number_t0.35 0.03266743 0.02784903 -0.02195881 0.08709921    0 0.2411232
: duration_t3  0.03869945 0.02491698 -0.01019050 0.08740482    0 0.1207617

\clearpage

In term of uncertainty quantification, it is important to specify
=match=TRUE= when using a single call so =BuyseTest= does not treat
each line of the dataset as an independent replicate. An intuitive way
to evaluate the standard error of the pooled estimator is to use the
across subject variability:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
sqrt(apply(df.lowGPC[,1:3],2,var)/NROW(df.lowGPC))
#+END_SRC

#+RESULTS:
:    rcs_t1.45 number_t0.35  duration_t3 
:   0.02075177   0.03182148   0.02839590

\noindent This is, up to a factor =N/(N-1)= exactly what the single call
approach returns. Actually we can retrieve this value by modifying the default options:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
BuyseTest.options(order.Hprojection = 2)
lowGPC2 <- BuyseTest(update(fff, . ~ . + strata(id, match = TRUE)),
                     data = lowProfil, trace = FALSE)
confint(lowGPC2)
#+END_SRC

#+RESULTS:
: estimate         se    lower.ci   upper.ci null   p.value
: rcs_t1.45    0.02741742 0.02075177 -0.01327825 0.06802241    0 0.1866527
: number_t0.35 0.02755903 0.03182148 -0.03483625 0.08974030    0 0.3867027
: duration_t3  0.03397497 0.02839590 -0.02172779 0.08946745    0 0.2318709

\noindent  Similarly when using other weighting scheme. For instance we can
retrieve the results of the Buyse pooling scheme doing:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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))
#+END_SRC

#+RESULTS:
: rcs_t1.45 number_t0.35  duration_t3 
:   0.01771680   0.02668638   0.02400508

\clearpage

** Accounting for time (WORK IN PROGRESS!)

In previous approaches, all pairwise comparisons are performed within
each patient, not only within-block comparisons. As of version
3.3.3. it is not yet possible to match at the patient-level and only
perform comparisons within the same period. One would have to perform
a separate GPC for each patient, with period as a strata variable:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
ls.lowGPC_period <- by(lowProfil, lowProfil$id, BuyseTest,
                       formula = update(fff,.~.+strata(period)), trace = FALSE)
#+END_SRC

#+RESULTS:

We can check that fewer pairs were made thanks to stratification:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
rbind(noStrata = nobs(ls.lowGPC[[1]]),
      PeriodStrata = nobs(ls.lowGPC_period[[1]]))
#+END_SRC

#+RESULTS:
:              placebo lowDose pairs
: noStrata          28      28   784
: PeriodStrata      28      28   196

In particular 49 pairs where made in each period:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
nobs(ls.lowGPC_period[[1]],strata=TRUE)
#+END_SRC

#+RESULTS:
:   placebo lowDose pairs
: 1       7       7    49
: 2       7       7    49
: 3       7       7    49
: 4       7       7    49

One would then combine the results in a single data.frame:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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)))
#+END_SRC

#+RESULTS:

and decide on a weighting scheme to pool the results:
#+BEGIN_SRC R :exports both :results output :session *R* :cache no
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))
#+END_SRC

#+RESULTS:
:          rcs_t1.45 number_t0.35 duration_t3
: average 0.02338177   0.03351172  0.04072124
: Buyse   0.01731834   0.04085114  0.04858818
: CMH     0.02333549   0.03361112  0.04077570

Resampling methods like a non-parametric bootstrap can then be used to
quantify uncertainty and obtain confidence intervals and p-values.

\clearpage

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

Handling more than two treatment groups is still an area of
development for the BuyseTest package. Several approaches have been
proposed in the litterature and here we focus on one that aim at
handling the non-transitivity issues that comes with Wilcoxon-like
tests citep:lumley2016characterising. This approach compare the
treatment-specific distribution to a pooled distribution over all
treatment groups citep:thangavelu2007wilcoxon:

#+BEGIN_SRC R :exports both :results output :session *R* :cache no
CasinoTest(update(fff, . ~ . + strata(id)), data = profil,
           method.inference = "none", type = "unweighted")
## do not trust inference since it would require accounting for matching
#+END_SRC

#+RESULTS:
#+begin_example
                   endpoint treatment     estimate se lower.ci upper.ci null p.value
placebo: rcs            rcs   placebo -0.012193158 NA       NA       NA   NA      NA
lowDose: rcs            rcs   lowDose  0.008015351 NA       NA       NA   NA      NA
highDose: rcs           rcs  highDose  0.004177807 NA       NA       NA   NA      NA
placebo: number      number   placebo -0.020434169 NA       NA       NA   NA      NA
lowDose: number      number   lowDose  0.011566310 NA       NA       NA   NA      NA
highDose: number     number  highDose  0.008867860 NA       NA       NA   NA      NA
placebo: duration  duration   placebo -0.024148505 NA       NA       NA   NA      NA
lowDose: duration  duration   lowDose  0.013113939 NA       NA       NA   NA      NA
highDose: duration duration  highDose  0.011034566 NA       NA       NA   NA      NA
#+end_example

Its main drawback is that the assessment of say =placebo=
vs. =lowDose= is now influenced by =highDose=.


* References
:PROPERTIES:
:UNNUMBERED: t
:END:

#+BEGIN_EXPORT latex
\begingroup
\renewcommand{\section}[2]{}
#+END_EXPORT

bibliographystyle:apalike
[[bibliography:bibliography.bib]]

#+BEGIN_EXPORT latex
\endgroup
#+END_EXPORT

* CONFIG                                                           :noexport:
#+LANGUAGE:  en
#+LaTeX_CLASS: org-article
#+LaTeX_CLASS_OPTIONS: [12pt]
#+OPTIONS:   title:t author:t toc:nil todo:nil
#+OPTIONS:   H:3 num:t 
#+OPTIONS:   TeX:t LaTeX:t
** Display of the document
# ## space between lines
#+LATEX_HEADER: \RequirePackage{setspace} % to modify the space between lines - incompatible with footnote in beamer
#+LaTeX_HEADER:\renewcommand{\baselinestretch}{1.1}
# ## margins
#+LaTeX_HEADER: \geometry{a4paper, left=10mm, right=10mm, top=10mm}
# ## personalize the prefix in the name of the sections
#+LaTeX_HEADER: \usepackage{titlesec}
# ## fix bug in titlesec version
# ##  https://tex.stackexchange.com/questions/299969/titlesec-loss-of-section-numbering-with-the-new-update-2016-03-15
#+LaTeX_HEADER: \usepackage{etoolbox}
#+LaTeX_HEADER: 
#+LaTeX_HEADER: \makeatletter
#+LaTeX_HEADER: \patchcmd{\ttlh@hang}{\parindent\z@}{\parindent\z@\leavevmode}{}{}
#+LaTeX_HEADER: \patchcmd{\ttlh@hang}{\noindent}{}{}{}
#+LaTeX_HEADER: \makeatother
** Color
# ## define new colors
#+LATEX_HEADER: \RequirePackage{colortbl} % arrayrulecolor to mix colors
#+LaTeX_HEADER: \definecolor{myorange}{rgb}{1,0.2,0}
#+LaTeX_HEADER: \definecolor{mypurple}{rgb}{0.7,0,8}
#+LaTeX_HEADER: \definecolor{mycyan}{rgb}{0,0.6,0.6}
#+LaTeX_HEADER: \newcommand{\lightblue}{blue!50!white}
#+LaTeX_HEADER: \newcommand{\darkblue}{blue!80!black}
#+LaTeX_HEADER: \newcommand{\darkgreen}{green!50!black}
#+LaTeX_HEADER: \newcommand{\darkred}{red!50!black}
#+LaTeX_HEADER: \definecolor{gray}{gray}{0.5}
# ## change the color of the links
#+LaTeX_HEADER: \hypersetup{
#+LaTeX_HEADER:  citecolor=[rgb]{0,0.5,0},
#+LaTeX_HEADER:  urlcolor=[rgb]{0,0,0.5},
#+LaTeX_HEADER:  linkcolor=[rgb]{0,0,0.5},
#+LaTeX_HEADER: }
** Font
# https://tex.stackexchange.com/questions/25249/how-do-i-use-a-particular-font-for-a-small-section-of-text-in-my-document
#+LaTeX_HEADER: \newenvironment{note}{\small \color{gray}\fontfamily{lmtt}\selectfont}{\par}
#+LaTeX_HEADER: \newenvironment{activity}{\color{orange}\fontfamily{qzc}\selectfont}{\par}
** Symbols
# ## valid and cross symbols
#+LaTeX_HEADER: \RequirePackage{pifont}
#+LaTeX_HEADER: \RequirePackage{relsize}
#+LaTeX_HEADER: \newcommand{\Cross}{{\raisebox{-0.5ex}%
#+LaTeX_HEADER:		{\relsize{1.5}\ding{56}}}\hspace{1pt} }
#+LaTeX_HEADER: \newcommand{\Valid}{{\raisebox{-0.5ex}%
#+LaTeX_HEADER:		{\relsize{1.5}\ding{52}}}\hspace{1pt} }
#+LaTeX_HEADER: \newcommand{\CrossR}{ \textcolor{red}{\Cross} }
#+LaTeX_HEADER: \newcommand{\ValidV}{ \textcolor{green}{\Valid} }
# ## warning symbol
#+LaTeX_HEADER: \usepackage{stackengine}
#+LaTeX_HEADER: \usepackage{scalerel}
#+LaTeX_HEADER: \newcommand\Warning[1][3ex]{%
#+LaTeX_HEADER:   \renewcommand\stacktype{L}%
#+LaTeX_HEADER:   \scaleto{\stackon[1.3pt]{\color{red}$\triangle$}{\tiny\bfseries !}}{#1}%
#+LaTeX_HEADER:   \xspace
#+LaTeX_HEADER: }
# # R Software
# ## R logo
#+LATEX_HEADER:\definecolor{grayR}{HTML}{8A8990}
#+LATEX_HEADER:\definecolor{grayL}{HTML}{C4C7C9}
#+LATEX_HEADER:\definecolor{blueM}{HTML}{1F63B5}
#+LATEX_HEADER: \newcommand{\Rlogo}[1][0.07]{
#+LATEX_HEADER: \begin{tikzpicture}[scale=#1]
#+LATEX_HEADER: \shade [right color=grayR,left color=grayL,shading angle=60] 
#+LATEX_HEADER: (-3.55,0.3) .. controls (-3.55,1.75) 
#+LATEX_HEADER: and (-1.9,2.7) .. (0,2.7) .. controls (2.05,2.7)  
#+LATEX_HEADER: and (3.5,1.6) .. (3.5,0.3) .. controls (3.5,-1.2) 
#+LATEX_HEADER: and (1.55,-2) .. (0,-2) .. controls (-2.3,-2) 
#+LATEX_HEADER: and (-3.55,-0.75) .. cycle;
#+LATEX_HEADER: 
#+LATEX_HEADER: \fill[white] 
#+LATEX_HEADER: (-2.15,0.2) .. controls (-2.15,1.2) 
#+LATEX_HEADER: and (-0.7,1.8) .. (0.5,1.8) .. controls (2.2,1.8) 
#+LATEX_HEADER: and (3.1,1.2) .. (3.1,0.2) .. controls (3.1,-0.75) 
#+LATEX_HEADER: and (2.4,-1.45) .. (0.5,-1.45) .. controls (-1.1,-1.45) 
#+LATEX_HEADER: and (-2.15,-0.7) .. cycle;
#+LATEX_HEADER: 
#+LATEX_HEADER: \fill[blueM] 
#+LATEX_HEADER: (1.75,1.25) -- (-0.65,1.25) -- (-0.65,-2.75) -- (0.55,-2.75) -- (0.55,-1.15) -- 
#+LATEX_HEADER: (0.95,-1.15)  .. controls (1.15,-1.15) 
#+LATEX_HEADER: and (1.5,-1.9) .. (1.9,-2.75) -- (3.25,-2.75)  .. controls (2.2,-1) 
#+LATEX_HEADER: and (2.5,-1.2) .. (1.8,-0.95) .. controls (2.6,-0.9) 
#+LATEX_HEADER: and (2.85,-0.35) .. (2.85,0.2) .. controls (2.85,0.7) 
#+LATEX_HEADER: and (2.5,1.2) .. cycle;
#+LATEX_HEADER: 
#+LATEX_HEADER: \fill[white]  (1.4,0.4) -- (0.55,0.4) -- (0.55,-0.3) -- (1.4,-0.3).. controls (1.75,-0.3) 
#+LATEX_HEADER: and (1.75,0.4) .. cycle;
#+LATEX_HEADER: 
#+LATEX_HEADER: \end{tikzpicture}
#+LATEX_HEADER: }

** Code
:PROPERTIES:
:ID:       2ec77c4b-f83d-4612-9a89-a96ba1b7bf70
:END:
# Documentation at https://org-babel.readthedocs.io/en/latest/header-args/#results
# :tangle (yes/no/filename) extract source code with org-babel-tangle-file, see http://orgmode.org/manual/Extracting-source-code.html 
# :cache (yes/no)
# :eval (yes/no/never)
# :results (value/output/silent/graphics/raw/latex)
# :export (code/results/none/both)
#+PROPERTY: header-args :session *R* :tangle yes :cache no ## extra argument need to be on the same line as :session *R*
# Code display:
#+LATEX_HEADER: \RequirePackage{fancyvrb}
#+LATEX_HEADER: \DefineVerbatimEnvironment{verbatim}{Verbatim}{fontsize=\small,formatcom = {\color[rgb]{0.5,0,0}}}
# ## change font size input (global change)
# ## doc: https://ctan.math.illinois.edu/macros/latex/contrib/listings/listings.pdf
# #+LATEX_HEADER: \newskip kipamount    kipamount =6pt plus 0pt minus 6pt
# #+LATEX_HEADER: \lstdefinestyle{code-tiny}{basicstyle=\ttfamily\tiny, aboveskip =  kipamount, belowskip =  kipamount}
# #+LATEX_HEADER: \lstset{style=code-tiny}
# ## change font size input (local change, put just before BEGIN_SRC)
# ## #+ATTR_LATEX: :options basicstyle=\ttfamily\scriptsize
# ## change font size output (global change)
# ## \RecustomVerbatimEnvironment{verbatim}{Verbatim}{fontsize=\tiny,formatcom = {\color[rgb]{0.5,0,0}}}
** Lists
#+LATEX_HEADER: \RequirePackage{enumitem} % better than enumerate
** Image and graphs
#+LATEX_HEADER: \RequirePackage{epstopdf} % to be able to convert .eps to .pdf image files
#+LATEX_HEADER: \RequirePackage{capt-of} % 
#+LATEX_HEADER: \RequirePackage{caption} % newlines in graphics
#+LaTeX_HEADER: \RequirePackage{tikz-cd} % graph
# ## https://tools.ietf.org/doc/texlive-doc/latex/tikz-cd/tikz-cd-doc.pdf
** Table
#+LATEX_HEADER: \RequirePackage{booktabs} % for nice lines in table (e.g. toprule, bottomrule, midrule, cmidrule)
** Inline latex
# @@latex:any arbitrary LaTeX code@@
** Algorithm
#+LATEX_HEADER: \RequirePackage{amsmath}
#+LATEX_HEADER: \RequirePackage{algorithm}
#+LATEX_HEADER: \RequirePackage[noend]{algpseudocode}
** Math
#+LATEX_HEADER: \RequirePackage{dsfont}
#+LATEX_HEADER: \RequirePackage{amsmath,stmaryrd,graphicx}
#+LATEX_HEADER: \RequirePackage{prodint} % product integral symbol (\PRODI)
# ## lemma
# #+LaTeX_HEADER: \RequirePackage{amsthm}
# #+LaTeX_HEADER: \newtheorem{theorem}{Theorem}
# #+LaTeX_HEADER: \newtheorem{lemma}[theorem]{Lemma}
*** Template for shortcut
#+LATEX_HEADER: \usepackage{ifthen}
#+LATEX_HEADER: \usepackage{xifthen}
#+LATEX_HEADER: \usepackage{xargs}
#+LATEX_HEADER: \usepackage{xspace}
#+LATEX_HEADER: \newcommand\defOperator[7]{%
#+LATEX_HEADER:	\ifthenelse{\isempty{#2}}{
#+LATEX_HEADER:		\ifthenelse{\isempty{#1}}{#7{#3}#4}{#7{#3}#4 \left#5 #1 \right#6}
#+LATEX_HEADER:	}{
#+LATEX_HEADER:	\ifthenelse{\isempty{#1}}{#7{#3}#4_{#2}}{#7{#3}#4_{#1}\left#5 #2 \right#6}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommand\defUOperator[5]{%
#+LATEX_HEADER: \ifthenelse{\isempty{#1}}{
#+LATEX_HEADER:		#5\left#3 #2 \right#4
#+LATEX_HEADER: }{
#+LATEX_HEADER:	\ifthenelse{\isempty{#2}}{\underset{#1}{\operatornamewithlimits{#5}}}{
#+LATEX_HEADER:		\underset{#1}{\operatornamewithlimits{#5}}\left#3 #2 \right#4}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommand{\defBoldVar}[2]{	
#+LATEX_HEADER:	\ifthenelse{\equal{#2}{T}}{\boldsymbol{#1}}{\mathbf{#1}}
#+LATEX_HEADER: }
**** Probability
#+LATEX_HEADER: \newcommandx\Esp[2][1=,2=]{\defOperator{#1}{#2}{E}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Prob[2][1=,2=]{\defOperator{#1}{#2}{P}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Qrob[2][1=,2=]{\defOperator{#1}{#2}{Q}{}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Var[2][1=,2=]{\defOperator{#1}{#2}{V}{ar}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Cov[2][1=,2=]{\defOperator{#1}{#2}{C}{ov}{\lbrack}{\rbrack}{\mathbb}}
#+LATEX_HEADER: \newcommandx\Binom[2][1=,2=]{\defOperator{#1}{#2}{B}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Gaus[2][1=,2=]{\defOperator{#1}{#2}{N}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Wishart[2][1=,2=]{\defOperator{#1}{#2}{W}{ishart}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Likelihood[2][1=,2=]{\defOperator{#1}{#2}{L}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\logLikelihood[2][1=,2=]{\defOperator{#1}{#2}{\ell}{}{(}{)}{}}
#+LATEX_HEADER: \newcommandx\Information[2][1=,2=]{\defOperator{#1}{#2}{I}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Score[2][1=,2=]{\defOperator{#1}{#2}{S}{}{(}{)}{\mathcal}}
**** Operators
#+LATEX_HEADER: \newcommandx\Vois[2][1=,2=]{\defOperator{#1}{#2}{V}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\IF[2][1=,2=]{\defOperator{#1}{#2}{IF}{}{(}{)}{\mathcal}}
#+LATEX_HEADER: \newcommandx\Ind[1][1=]{\defOperator{}{#1}{1}{}{(}{)}{\mathds}}
#+LATEX_HEADER: \newcommandx\Max[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{min}}
#+LATEX_HEADER: \newcommandx\Min[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{max}}
#+LATEX_HEADER: \newcommandx\argMax[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{argmax}}
#+LATEX_HEADER: \newcommandx\argMin[2][1=,2=]{\defUOperator{#1}{#2}{(}{)}{argmin}}
#+LATEX_HEADER: \newcommandx\cvD[2][1=D,2=n \rightarrow \infty]{\xrightarrow[#2]{#1}}
#+LATEX_HEADER: \newcommandx\Hypothesis[2][1=,2=]{
#+LATEX_HEADER:         \ifthenelse{\isempty{#1}}{
#+LATEX_HEADER:         \mathcal{H}
#+LATEX_HEADER:         }{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#2}}{
#+LATEX_HEADER: 		\mathcal{H}_{#1}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\mathcal{H}^{(#2)}_{#1}
#+LATEX_HEADER:         }
#+LATEX_HEADER:         }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommandx\dpartial[4][1=,2=,3=,4=\partial]{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#3}}{
#+LATEX_HEADER: 		\frac{#4 #1}{#4 #2}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\left.\frac{#4 #1}{#4 #2}\right\rvert_{#3}
#+LATEX_HEADER: }
#+LATEX_HEADER: }
#+LATEX_HEADER: \newcommandx\dTpartial[3][1=,2=,3=]{\dpartial[#1][#2][#3][d]}
#+LATEX_HEADER: \newcommandx\ddpartial[3][1=,2=,3=]{
#+LATEX_HEADER: 	\ifthenelse{\isempty{#3}}{
#+LATEX_HEADER: 		\frac{\partial^{2} #1}{\partial #2^2}
#+LATEX_HEADER: 	}{
#+LATEX_HEADER: 	\frac{\partial^2 #1}{\partial #2\partial #3}
#+LATEX_HEADER: }
#+LATEX_HEADER: } 
**** General math
#+LATEX_HEADER: \newcommand\Real{\mathbb{R}}
#+LATEX_HEADER: \newcommand\Rational{\mathbb{Q}}
#+LATEX_HEADER: \newcommand\Natural{\mathbb{N}}
#+LATEX_HEADER: \newcommand\trans[1]{{#1}^\intercal}%\newcommand\trans[1]{{\vphantom{#1}}^\top{#1}}
#+LATEX_HEADER: \newcommand{\independent}{\mathrel{\text{\scalebox{1.5}{$\perp\mkern-10mu\perp$}}}}
#+LaTeX_HEADER: \newcommand\half{\frac{1}{2}}
#+LaTeX_HEADER: \newcommand\normMax[1]{\left|\left|#1\right|\right|_{max}}
#+LaTeX_HEADER: \newcommand\normTwo[1]{\left|\left|#1\right|\right|_{2}}
#+LATEX_HEADER: \newcommand\Veta{\boldsymbol{\eta}}
#+LATEX_HEADER: \newcommand\VX{\mathbf{X}}
** Notations



