---
title: "ALL-IN meta-analysis with Safe Logrank Testing"
author: "[Judith ter Schure](https://www.cwi.nl/people/judith-ter-schure), [Alexander Ly](https://www.cwi.nl/people/alexander-ly) & [Peter Grünwald](https://www.cwi.nl/people/peter-gruenwald)"
date: "`r Sys.Date()`"
output:
  html_document
link-citations: yes
bibliography: References - Safe logrank testing.bib
---
```{r, include = F}
numSim <- 1000
numPlot <- ifelse(numSim > 100, 100, numSim)
# set to 200 (needed for a proper output, gives a 1-2% simulation error in power and type-I error)
# set to 200: it takes 17 minutes to run all R chuncks (so at least that to knit)
# set to 1000 (better): it takes 1 hour and 25 minutes, 
# set to  4 (the minimum): it takes 21 seconds to run all R chuncks (so at least that to knit)
```


Suppose we have a medical question that urgently needs a reliable answer. Suppose further that this question is considered by multiple clinical trials running at the same time or at least overlapping in time. To provide a reliable answer as fast as possible, we want to combine the results, use interim analyses and/or encourage trials that finish early to inform whether others should be stopped. In principle, this can provide more power, and in case of time-to-event analysis, more early events to support sufficient evidence at interim. Furthermore, early stopping might prevent many events that otherwise occur if all trials are continued to their prespecified completion. However, in such an interim meta-analysis scenario, standard analysis is not able to smoothly retain type-I error guarantees.

Take the example in which trial results are only evaluated in isolation and a response follows the first positive result of a single trial. Such multiple testing inflates the type-I error dramatically if individual trials are evaluated with a standard methods. Part of the solution to this problem lies in pooling the trials: always evaluate (interim) trial results in unison and regain type-I error control over the meta-analysis. But new problems arise, because ideally, we would meta-analyze as soon as an interim result looks encouraging, but would like to test again if the first meta-analysis does not yet conclude the research effort. So, multiple testing issues also arise in meta-analysis, and standard analysis of isolated trials as well as meta-analysis end up with an unreliable inference and, subsequently, poor decisions.

The solution lies in a new statistical method that allows to flexibly pool trials, increase power and efficiency by unlimited interim analyses, while always retaining type-I error guarantees. The theory underlying such an *ALL-IN meta-analysis* is that of Safe Testing  (@grunwald2019safe) with its basic tool, the e-value. These e-values are analogous to p-values, or rather their inverse: a large e-value indicates evidence against the null hypothesis. But, completely unlike p-values, e-values of several trials can be combined effortlessly -- by multiplication -- without compromising type-I error. Such combination does not have to be prespecified and can be done even if the trials are dependent, e.g if the decision to start or stop an extra trial depends on (interim) results of earlier trials.

The Safe Testing methodology is described in detail by @grunwald2019safe. Its relation to dependency problems in meta-analysis, called accumulation bias, is outlined in @ter_schure_accumulation_2019. Safe Tests are an extension of earlier work by Vovk, Shafer and others on so-called martingale tests, and Robbins and Lai on always-valid tests, e.g. @shafer_test_2011, @robbins_statistical_1970. This latter work was a precursor to both Safe Testing and group sequential testing with alpha-spending approaches.

It could, in fact, be possible to combine interim test results via an alpha-spending or a group-sequential approach. But in contexts of multiple trials with no top-down planning these have considerable disadvantages compared to Safe Testing. Most importantly, they are much less flexible. Suppose, for example, that a new study is started halfway through our analysis; a study we did not know about when designing our meta-analyses procedure. Because alpha-spending requires us to know and fix the number of trials beforehand, it would not allow us to include this additional study. Fortunately, with Safe Testing accounting for the new study is no problem at all.

Below we illustrate, via examples in R, the Safe logrank test. Like the classical logrank test, this is a test for time-to-event data with censoring. But unlike the standard approach, the Safe logrank test outputs an e-value as measure of evidence.



## The Safe logrank test

We first introduce the Safe logrank test in its use for a single trial with only one underlying stratum. One stratum means that the baseline expected number of events per time point is equal for all enrolled participants. In case participants enter the study from different locations or hospitals, this assumption might be unrealistic. We therefore secondly introduce the Safe logrank test for the use on stratified time-to-event data. We conclude our discussion by an illustration of its use on different studies, each in its own stage of analysis. Due to the simplicity of the meta-analysis procedure -- pooling by multiplication -- neither of these extensions introduce any more complexity.


### Using the Safe logrank test: one study, one stratum, many interim analyses

Please also find the complete R markdown file [here](https://projects.cwi.nl/safestats/SafeLogrank.Rmd).

```{r, warning = F, message = F, results = "hide"}

library(devtools)
devtools::install_github("AlexanderLyNL/Safestats", ref = "logrank")
library(safestats)
library(survival)
library(knitr)
```


Consider the following survival data of a first study:

```{r}
enrollment.study1 <- 1500     # 750 treatment, 750 placebo, so:
ratio.study1 <- 1             # ratio = nT/nP = 1
fup.study1 <- 180             # folow up of 180 days
nEventsC.study1 <- 150        # 150 events in 180 days among 1500
                              # in case of no treatment effect,
                              # this defines the baseline hazard
hr1.study1 <- 0.5             # hazard ratio between treatment en placebo group
                              # (hr1: alternative hypothesis is true)
                              # 113 events anticipated within 180 days 
                              # if treatment reduces risk with hr1.study1:
anticipNevents.study1 <- nEventsC.study1*(1/(1 + ratio.study1)) + 
                         nEventsC.study1*(ratio.study1/(1 + ratio.study1))*hr1.study1  

# assuming constant hazard, we simulate the following baseline hazard and data:
lambdaC.study1 <- 1 - (1 - (nEventsC.study1/enrollment.study1))^(1/fup.study1)  
data.study1 <- generateSurvData(nP = enrollment.study1*(1/(1 + ratio.study1)),
                                nT = enrollment.study1*(ratio.study1/(1 + ratio.study1)),
                                lambdaP = lambdaC.study1,
                                lambdaT = hr1.study1*lambdaC.study1,
                                endTime = fup.study1,
                                seed = 52020)  # set seed to 52020 make the result reproducible
```

```{r, echo = F}
plot(survfit(Surv(data.study1$time, data.study1$status) ~ data.study1$group),
     xlim = c(0, fup.study1), ylim = c(0.88, 1), las = 1,
     ylab = "Survival", xlab = "d (days)", lwd = 2,
     main = paste("Study 1: hazard ratio = ", hr1.study1))
  
kable(head(data.study1, n = 10), row.names = F)
```

We see events occurring in either the treatment (`T`) or the placebo control group (`P`) at day `r unique(head(data.study1$time, n = 10))`,... (events have `status = 2`, censoring has `status = 1`). These are the possible days for an interim analysis. An analysis at day 75 (`d = 75`) would look as follows:

```{r, include = F}
d <- 75
```

```{r}
alpha <- 0.025

dataSoFar <- data.study1
dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d

designObj.study1 <- designSafeLogrank(hrMin = hr1.study1,
                                      alpha = alpha,
                                      alternative = "less",  # one-sided test hr < 1
                                      ratio = ratio.study1)
dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)

test <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group, 
                        designObj = designObj.study1)
test
```

The analysis returns an e-value of ``r round(test$eValue, digits = 1)``, that has an interpretation in terms of evidence: a larger e-value means there is more evidence against the null hypothesis. Just like a p-value, we can use a threshold based on \(\alpha\) to obtain a test that has \(\alpha\) type-I error control (has smaller than \(\alpha\) probability to cross the threshold if the null hypothesis is true). The threshold to use for the e-value is \(\frac{1}{\alpha}\), meaning that for the one-sided test above with \(\alpha = `r alpha`\), our e-value needs to be bigger than \(`r 1/alpha`\).

If we track the available evidence using a Safe logrank test, our progress over all days with possible interim analyses is as follows:

```{r}
eValues.study1 <- rep(NA, times = fup.study1)
interimDays.study1 <- unique(data.study1$time[data.study1$status == 2])

eValues.study1[1:interimDays.study1[1]] <- 1  # before you observe any event, your e-value is 1

for (d in interimDays.study1[1]:fup.study1) {
  if (!(d %in% interimDays.study1)) {           # at days on which you do not observe an event, 
    eValues.study1[d] <- eValues.study1[d - 1]  # evidence stays the same as the day before
  } else {
    dataSoFar <- data.study1
    dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d
    
    dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)
    eValues.study1[d] <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group,
                                         designObj = designObj.study1
                                         )$eValue
  }
  if (eValues.study1[d] > (1/alpha)) break()
}
```

```{r, echo = F}
plot(NULL, xlim = c(0, fup.study1), ylim = c(-6, 6),
     ylab = "e-value", xlab = "d (days)", yaxt = "n",
     main = paste("Study 1: hazard ratio = ", hr1.study1))
axis <- mapply(axis, side = 2, at = c(-6:5, log(1/alpha)/log(2), 6), cex.axis = 0.65,
               labels = c(paste0("1/", 2^(6:1)), 2^(0:4), "", 1/alpha, 2^6), las = 1)
abline(h = log(1/alpha)/log(2), lty = 2)

lines(x = 1:fup.study1, y = log(eValues.study1)/log(2), 
      'l', col = "darkgreen", lwd = 2)
```

Note that the evidence in terms of e-value is doubling on the 2-log y-axis. We observe that the threshold of \(e > `r 1/alpha`\) is crossed at the interim analysis around ``r which(eValues.study1 > (1/alpha))[1]` days`. This means that we can stop the trial or choose this moment to initiate a pooled analyses. We discuss a pooled analyses in the final paragraph.



#### Power and type-I error control of the Safe logrank test

What to expect before starting a trial -- how many days we need before reaching \(e > `r 1/alpha`\) -- depends of course on how many events we observe early on and how large the difference is in hazard between our two groups. The simulation below shows the power of our analysis, if the data is generated following the assumptions of our first study above: 

```{r}
simeValues.study1 <- matrix(NA, nrow = numSim, ncol = fup.study1)
simEvents.study1 <- rep(NA, times = numSim)

for (j in 1:numSim){
    data <- generateSurvData(nP = enrollment.study1*(1/(1 + ratio.study1)),
                             nT = enrollment.study1*(ratio.study1/(1 + ratio.study1)),
                             lambdaP = lambdaC.study1,
                             lambdaT = hr1.study1*lambdaC.study1,
                             endTime = fup.study1, 
                             seed = j)

    interimDays <- unique(data$time[data$status == 2])

    for (d in interimDays) {
      # to speed up the simulation, we do not fill all days with e-values,
      # but leave those at which evidence does not change NA
      dataSoFar <- data
      dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d
      
      dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)
      simeValues.study1[j, d] <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group,
                                                 designObj = designObj.study1
                                                 )$eValue
      if (simeValues.study1[j, d] > (1/alpha)) {
        simEvents.study1[j] <- sum(dataSoFar$status == 2)
        break()
      }
    }
  }

power.optimal <- mean(sapply(1:numSim, function(i) any(simeValues.study1[i, ] > (1/alpha), na.rm = T)))*100
```

```{r, echo = F}
plot(NULL, xlim = c(-40, fup.study1), ylim = c(-6, 6), 
     ylab = "e-value", xlab = "d (days)", yaxt = "n",
     main = paste("Study 1: hazard ratio = ", hr1.study1))
axis <- mapply(axis, side = 2, at = c(-6:5, log(1/alpha)/log(2), 6), cex.axis = 0.65,
               labels = c(paste0("1/", 2^(6:1)), 2^(0:4), "", 1/alpha, 2^6), las = 1)
abline(h = log(1/alpha)/log(2), lty = 2)
legend("bottomleft", fill = c("darkgreen", "darkred"), 
       legend = c(paste0("ever e > ", 1/alpha), paste0("never e > ", 1/alpha)), cex = 0.8)

simeValues.study1.sample <- simeValues.study1[1:(numPlot/2), ]
for (i in 1:nrow(simeValues.study1.sample)) {
    x <- which(!is.na(simeValues.study1.sample[i, ]))
    y <- na.omit(simeValues.study1.sample[i, ])
    myColour <- ifelse(any(y > (1/alpha)), "darkgreen", "darkred")
    lines(x, log(y)/log(2), col = myColour, lwd = 2)
}
for (j in 1:numSim) {
  x <- which(!is.na(simeValues.study1[j, ]))
  y <- na.omit(simeValues.study1[j, ])
  if (any(y > (1/alpha))) {
    abline(v = max(x), col = "darkgreen", lty = 3)
  }
}
```

```{r, echo = F}
hr0.study1 <- 1  # hazard ratio

simeValuesNULL.study1 <- matrix(NA, nrow = numSim, ncol = fup.study1)

for (j in 1:numSim){
    data <- generateSurvData(nP = enrollment.study1*(1/(1 + ratio.study1)),
                             nT = enrollment.study1*(ratio.study1/(1 + ratio.study1)),
                             lambdaP = lambdaC.study1,
                             lambdaT = hr0.study1*lambdaC.study1,
                             endTime = fup.study1,
                             seed = j + 2021)
    interimDays.study1 <- unique(data$time[data$status == 2])
    
    # before you observe any event, your e-value is 1
    simeValuesNULL.study1[j, 1:interimDays.study1[1]] <- 1  

    for (d in interimDays.study1[1]:fup.study1) {
      if (!(d %in% interimDays.study1)) {                        
        # at days on which you do not observe an event,
        # evidence stays the same as the day before:
        simeValuesNULL.study1[j, d] <- simeValuesNULL.study1[j, d - 1]  
      } else {
        dataSoFar <- data
        dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d
        
        dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)
        simeValuesNULL.study1[j, d] <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group,
                                                       designObj = designObj.study1
                                                       )$eValue
  
        if (simeValuesNULL.study1[j, d] > (1/alpha)) {
          break()
        }
      }
    }
}

typeIerror <- mean(sapply(1:numSim, 
                          function(i) any(simeValuesNULL.study1[i, ] > (1/alpha), na.rm = T)))*100
```

The vertical lines in the above plot show all the possible days of interim analysis at which a simulated analysis reaches an e-value of \(`r 1/alpha`\). Of those ``r numSim`` simulated paths, a sample of ``r numPlot/2`` are shown. Green e-value paths indicate that the study reaches \(e > `r 1/alpha`\) at some interim analysis, red e-value paths indicate that it does not.

``r power.optimal`%` of all simulated e-value paths are green, which means that we have an approximate power of ``r power.optimal`%` within our follow up of ``r fup.study1`` days to correctly conclude that there is a difference between the two groups (we have set that difference to be a hazard ratio of \(0.5\) (`hr1 = 0.5`). If we set the hazard ratio to \(1\) (`hr0 = 1`) and run the same simulation, we observe that our type-I error, so hitting our threshold of `r 1/alpha` within 180 days if there is no difference between the groups, is only ``r typeIerror`%`.

Unlike alpha-spending, the type-I error control does not depend on any specification of a spending function. Type-I error control is preserved by comparing the e-value to \(`r 1/alpha`\), regardless of the days of interim analysis to do so. So the threshold of \(`r 1/alpha`\) allows for all possible interim analyses, without spending all alpha, as is shown by our simulated type-I error ``r typeIerror`%`. In fact, not all \(\alpha = `r alpha*100`\%\) is spent in these possible interim analyses: the Safe test always reserves some \(\alpha\) for pooled analyses with data from other studies, and can in fact keep combining data indefinitely.

```{r, echo = F}
plot(NULL, xlim = c(-40, fup.study1), ylim = c(-6, 6), 
     ylab = "e-value", xlab = "d (days)", yaxt = "n",
     main = paste("Study 1: hazard ratio = ", hr0.study1))
axis <- mapply(axis, side = 2, at = c(-6:5, log(1/alpha)/log(2), 6), cex.axis = 0.65,
               labels = c(paste0("1/", 2^(6:1)), 2^(0:4), "", 1/alpha, 2^6), las = 1)
abline(h = log(1/alpha)/log(2), lty = 2)
legend("bottomleft", fill = c("darkgreen", "darkred"), 
       legend = c(paste0("ever e > ", 1/alpha), paste0("never e > ", 1/alpha)), cex = 0.8)

simeValuesNULL.study1.sample <- simeValuesNULL.study1[1:numPlot, ] 
for (i in 1:nrow(simeValuesNULL.study1.sample)) {
    x <- which(!is.na(simeValuesNULL.study1.sample[i, ]))
    y <- na.omit(simeValuesNULL.study1.sample[i, ])
    myColour <- ifelse(any(y > (1/alpha)), "darkgreen", "darkred")
    lines(x, log(y)/log(2), col = myColour, lwd = 2)
}
```


#### Designing the Safe logrank test

Type-I error control does not depend on what e-value or threshold we use. The power and the number of events before stopping the trial, on the other hand, do depend on how the Safe logrank test is designed. In the above simulation a design object `designObj.study1` was defined based on a two-sided test and hazard ratio `hr1 = 0.5` that exactly matched the simulation. This is the best possible scenario, and results in the following simulated events before stopping at an e-value of \(`r 1/alpha`\), out of ``r numSim`` simulated runs:

```{r, echo = F}
hist(simEvents.study1, breaks = seq(0, 150, by = 5), freq = T, xlab = "Number of events before stopping",
     ylab = paste(paste("Frequency out of ", numSim), " simulations"),
     xlim = c(0, 140), ylim = c(0, numSim/10), las = 1,
     main = paste(paste("Simulated studies: hazard ratio = ", hr1.study1), 
                  paste("\n Design Safe test: hazard ratio = ", hr1.study1)))
```

However, in reality, we can supply a hazard ratio that is of minimal clinical interest, but can never truly know it in advance. If we design a test based on a difference between the groups that is more extreme than the true value, type-I error is preserved, but we obtain suboptimal power and need more events to detect the true difference. We are underpowered, because the true effect is smaller than the effect we set to be of minimal clinical relevance.

The hazard ratio used to design the Safe test should be prespecified at the same stage as would a power analysis: before seeing any data. Adjusting the `designObject` to the observed hazard ratio is like post-hoc changing your alpha-spending function to capitalize on many early or late events. It invalidates the analysis and loses type-I error control.

```{r}
hr1.study1 <- 0.5
hr1.suboptimal <- 0.3
designObj.study1.suboptimal <- designSafeLogrank(hrMin = hr1.suboptimal,
                                                 alpha = alpha,
                                                 alternative = "less",  # one-sided test hr < 1
                                                 ratio = ratio.study1)
  
designObj.study1.suboptimal
```

The consequences of a too optimistic hazard ratio are shown in the simulation below. While the data is simulated assuming a hazard ratio of ``r hr1.study1``, the test is designed based on a hazard ratio of ``r hr1.suboptimal``.

Note that if the number of events is sufficiently large, such as in pooling time-to-event data across multiple trials, a Safe test designed with a too optimistic effect size can still have sufficient power.

```{r, echo = F}
simeValues.suboptimal <- matrix(NA, nrow = numSim, ncol = fup.study1)
simEvents.study1.suboptimal <- rep(NA, times = numSim)

for (j in 1:numSim){
    data <- generateSurvData(nP = enrollment.study1*(1/(1 + ratio.study1)),
                             nT = enrollment.study1*(ratio.study1/(1 + ratio.study1)),
                             lambdaP = lambdaC.study1,
                             lambdaT = hr1.study1*lambdaC.study1,
                             endTime = fup.study1,
                             seed = j + 3006)
    interimDays <- unique(data$time[data$status == 2])

    for (d in interimDays) {
      # to speed up the process, we do not fill all days with e-values,
      # but leave those at which evidence does not change NA
      dataSoFar <- data
      dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d
     
      dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)
      simeValues.suboptimal[j, d] <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group,
                                                     designObj = designObj.study1.suboptimal
                                                     )$eValue

      if (simeValues.suboptimal[j, d] > (1/alpha)) {
        simEvents.study1.suboptimal[j] <- sum(dataSoFar$status == 2)
        break()
      }
    }
}

power.suboptimal <- mean(sapply(1:numSim, function(i) any(simeValues.suboptimal[i, ] > (1/alpha), na.rm = T)))*100
```

The power obtained now with 1500 participants over 180 days is `power.suboptimal = `r power.suboptimal`%` instead of `power.optimal = `r power.optimal`%`. And the distribution of events needed before stopping the trial looks as follows:

```{r, echo = F}
hist(simEvents.study1.suboptimal, breaks = seq(0, 150, by = 5), freq = T, 
     xlab = "Number of events before stopping", ylab = paste(paste("Frequency out of ", numSim), " simulations"),
     xlim = c(0, 140), ylim = c(0, numSim/10), las = 1,
     main = paste(paste("Simulated studies: hazard ratio = ", hr1.study1), 
                  paste("\n Design Safe test: hazard ratio = ", hr1.suboptimal)))
```


### Combining the Safe logrank test: one study, multiple strata, many interim analyses

If our trial is composed of several strata, such as hospitals with different characteristics, it is optimal to account for this structure. After all, the underlying baseline hazard, at a specific day after randomization, might be different from hospital to hospital. So we want to allow the events that are observed on that day to tell a different story per hospital.

The standard logrank test as well as the Safe logrank test assume that at each day that observes an event, there is one underlying baseline hazard of events. If we assume the null hypothesis of no difference between placebo and treatment group, then this baseline hazard underlies observations in both groups. When the observations deviate from that pattern, the e-value grows and indicates evidence against the null. Consequently, it is not optimal to calculate an e-value based on combined events from different strata. In that case the e-value might not be able to detect a difference between treatment and control group, because the differences between the strata make it more uncertain about the baseline.

So accounting for strata -- or more general: baseline characteristics -- makes the analysis more powerful. But the type-I error is always controlled, whether we correct for baseline characteristics or not.

To account for differences between strata, we calculate an e-value per stratum and multiply these e-values. As we observed in the beginning of this text, this preserves type-I error control: for every fixed $\alpha$, the probability that a product of e-values *ever* crosses $1/\alpha$, is bounded by $\alpha$.

Intuitively, if the null hypothesis of no difference between the groups is true, each e-value will only get smaller. So if we do not specify when we multiply two or more e-values, we expect to multiply something small with something small, and do not expect the product to be anything but small. Because each e-value has a probability of only \(1\) in \(`r 1/alpha`\) to ever become larger than \(`r 1/alpha`\) if we keep adding data, so does the combined e-value if we add data by multiplying e-values.

This procedure does not inflate the type-I errors if constantly multiply our e-value evidence. But also if we are not able to do this constantly due to practical limitations, the type-I error is retained for any selection of possible interim meta-analyses. All these properties rely on the assumption that before we start the data collection, all analysis times are possible, so if we control type-I error overall, we control type-I error over the analyses we actually perform.


### Meta-analyzing the Safe logrank test: many studies, multiple stages, many interim analyses

Just like tracking an individual study e-value, as shown above, we can also track e-values from multiple studies at different stages, and pool the results into one meta-analysis. The following code describes the properties of a second study:

```{r}
enrollment.study2 <- 4000     # 2000 treatment, 2000 placebo, so:
ratio.study2 <- 1             # ratio = nT/nP = 1
fup.study2  <- 140            # follow up of 140 days
nEventsC.study2 <- 300        # 300 events in 140 days among 4000
                              # in case of no treatment effect,
                              # this defines the baseline hazard
hr1.study2 <- 0.7             # hazard ratio between treatment en placebo group
                              # (hr1: alternative hypothesis is true)
                              # 255 events anticipated within 140 days 
                              # if treatment reduces risk with hr1.study2:
anticipNevents.study2 <- nEventsC.study2*(1/(1 + ratio.study2)) + 
                         nEventsC.study2*(ratio.study2/(1 + ratio.study2))*hr1.study2  

designObj.study2 <- designSafeLogrank(hrMin = hr1.study2,
                                      alpha = alpha,
                                      alternative = "less",  # one-sided test hr < 1
                                      ratio = ratio.study2)

# assuming constant hazard, we simulate the following baseline hazard and data:
lambdaC.study2 <- 1 - (1 - (nEventsC.study2/enrollment.study2))^(1/fup.study2) 
data.study2 <- generateSurvData(nP = enrollment.study2*(1/(1 + ratio.study2)),
                                nT = enrollment.study2*(ratio.study2/(1 + ratio.study2)),
                                lambdaP = lambdaC.study2,
                                lambdaT = hr1.study2*lambdaC.study2,
                                endTime = fup.study2,
                                seed = 20)  # set seed to 20 to make the result reproducible
# Assume that this new study starts 40 days later than the previous one
data.study2$time <- data.study2$time + 40
```

```{r, echo = F}
interimDays.study2 <- unique(data.study2$time[data.study2$status == 2])

eValues.study2 <- rep(NA, times = fup.study1)
eValues.study2[1:interimDays.study2[1]] <- 1  # before you observe any event, your eValue is 1

for (d in interimDays.study2[1]:fup.study1) {
  if (!(d %in% interimDays.study2)) {           # at days on which you do not observe an event, 
    eValues.study2[d] <- eValues.study2[d - 1]  # evidence stays the same as the day before
  } else {
    dataSoFar <- data.study2
    dataSoFar$status[dataSoFar$time > d] <- 1  # censor everything after day d
   
    dataSoFar$survObj <- Surv(dataSoFar$time, dataSoFar$status)
    eValues.study2[d] <- safeLogrankTest(formula = dataSoFar$survObj ~ dataSoFar$group,
                                         designObj = designObj.study2
                                         )$eValue
  }
  if (eValues.study2[d] > (1/alpha)) break()
}
```

```{r, echo = F}
survivalObj <- survfit(Surv(data.study2$time, data.study2$status) ~ data.study2$group)
plot(survivalObj, las = 1,
     xlim = c(0, fup.study2 + 40), ylim = c(0.88, 1), 
     ylab = "Survival", xlab = "d (days)", lwd = 2,
     main = paste("Study 2: hazard ratio = ", hr1.study2))
  
kable(head(data.study2, n = 10), row.names = F)
```

Because we assume that this second study starts 40 days later than the first, the days for interim analyses on the timeline of the first study are: `r head(interimDays.study2, 10)`

We can pool the evidence in the second study with the first one by multiplying their e-values:


```{r}
eValueMeta <- sapply(1:fup.study1, function(day) prod(eValues.study1[day], eValues.study2[day]))
```


```{r, echo = F}
plot(NULL, xlim = c(0, fup.study1), ylim = c(-6, 6), 
     ylab = "e-value", xlab = "d (days)", yaxt = "n",
     main = paste(paste(paste("hazard ratio = ", hr1.study1),
                              " (study 1) and "), 
                              paste(hr1.study2, " (study 2)")))
axis <- mapply(axis, side = 2, at = c(-6:5, log(1/alpha)/log(2), 6), cex.axis = 0.65,
               labels = c(paste0("1/", 2^(6:1)), 2^(0:4), "", 1/alpha, 2^6), las = 1)
abline(h = log(1/alpha)/log(2), lty = 2)
legend("bottomleft", fill = c("darkgreen", "green", "blue"), 
       legend = c(paste0("study1: ever e > ", 1/alpha), 
                  paste0("study2: ever e > ", 1/alpha),
                  paste0("meta: ever e > ", 1/alpha)), cex = 0.8)

lines(x = 1:fup.study1, y = log(eValues.study1)/log(2), 
      'l', col = "darkgreen", lwd = 2)
lines(x = 1:length(eValues.study2),
      y = log(eValues.study2)/log(2), col = "green", lwd = 2)
lines(x = 1:length(eValueMeta),
      y = log(eValueMeta)/log(2) + 0.07, col = "blue", lwd = 2)
abline(v = which(eValues.study1 > (1/alpha))[1], col = "darkgreen", lty = 2)
abline(v = which(eValueMeta > (1/alpha))[1], col = "blue", lty = 2)
abline(v = which(eValues.study2 > (1/alpha))[1], col = "green", lty = 2)
```

We observe that the second study in its first ``r which(eValueMeta > (1/alpha))[1]` days` does not have much power individually to reach \(`r 1/alpha`\), due to very low number of events. Fortunately, pooling the little evidence it does provide with that already provided by the first study enables the meta-analysis e-value to reach the threshold sooner.


#### Type-I error control over all possible analysis times

The error control of multiplied e-values is preserved over all interim analyses and all possible timings for the pooled meta-analysis. This means that if you design to pool, it can also be that the first pooled analysis contains only one study. While such an analysis would not conventionally be called a meta-analysis, our e-value meta-analysis allows for this situation. Before observing any event, the e-value of a study is 1, indicated by the flat line of study 2 in the first 40 days of the e-value plot above. Therefore, pooling by multiplication is always possible, also when e-values of some studies just provide  \(1\).

Note that error control and power over all possible pooled analyses therefore inherently assume that stopping the pooled analysis at a stage with only single study contribution is possible. Disallowing such a scenario does not necessarily compromise type-I error control (for details see @ter_schure_accumulation_2019 and its definition of error control *surviving over time*). However, only allowing for pooled analyses of more than one trial does decrease power because it excludes the possibility of reaching \(e > `r 1/alpha`\) with only one study.

We show that type-I error control is preserved over a pooled analysis, in which sometimes the first study that provides data is sufficient, and sometimes two studies combined provide the evidence.

```{r, echo = F}
hr0.study2 <- 1

simeValuesNULL.study2 <- matrix(NA, nrow = numSim, ncol = fup.study1)
simeValueMetaNULL <- matrix(NA, nrow = numSim, ncol = fup.study1)

fup.study2.pooled <- fup.study2 + 40  # study 2 starts 40 days later

for (j in 1:numSim){
  simeValueMetaNULL[j, ] <- simeValuesNULL.study1[j, ]  # available from previous simulation
  
  # assuming constant hazard
  lambdaC.study2 <- 1 - (1 - (anticipNevents.study2/enrollment.study2))^(1/fup.study2)  
  data.study2 <- generateSurvData(nP = enrollment.study2*(1/(1 + ratio.study2)),
                                  nT = enrollment.study2*(ratio.study2/(1 + ratio.study2)),
                                  lambdaP = lambdaC.study2,
                                  lambdaT = hr0.study2*lambdaC.study2,
                                  endTime = fup.study2,
                                  seed = j + 630)
  data.study2$time <- data.study2$time + 40  # study 2 starts 40 days later
  interimDays.study2 <- unique(data.study2$time[data.study2$status == 2])
  
  # before you observe any event, your e-value is 1
  simeValuesNULL.study2[j, 1:interimDays.study2[1]] <- 1 

  for (d in interimDays.study2[1]:fup.study1) {
    # at days on which you do not observe an event,
    # evidence stays the same as the day before
    if (!(d %in% interimDays.study2)) {                        
      simeValuesNULL.study2[j, d] <- simeValuesNULL.study2[j, d - 1]
      } else {
        dataSoFar.study2 <- data.study2
        dataSoFar.study2$status[dataSoFar.study2$time > d] <- 1  # censor after day d
        
        dataSoFar$survObj <- Surv(dataSoFar.study2$time, dataSoFar.study2$status)
        simeValuesNULL.study2[j, d] <- safeLogrankTest(formula = 
                                                         dataSoFar$survObj ~ dataSoFar.study2$group,
                                                       designObj = designObj.study2
                                                       )$eValue
      }
      simeValueMetaNULL[j, d] <- prod(simeValuesNULL.study1[j, d], simeValuesNULL.study2[j, d])
      if (simeValuesNULL.study2[j, d] > (1/alpha)) break()
  }
}
```

```{r, echo = F}
plot(NULL, xlim = c(-40, fup.study1), ylim = c(-11, 6), 
     ylab = "e-value", xlab = "d (days)", yaxt = "n",
     main = paste(paste(paste("hazard ratio = ", hr0.study1),
                              " (study 1) and "), 
                              paste(hr0.study2, " (study 2)")))
axis <- mapply(axis, side = 2, at = c(-11:5, log(1/alpha)/log(2), 6), cex.axis = 0.55,
               labels = c(paste0("1/", 2^(11:1)), 2^(0:4), "", 1/alpha, 
                          paste0("  ", 2^6)), las = 1)
abline(h = log(1/alpha)/log(2), lty = 2)
legend("bottomleft", fill = c("darkgreen", "green", "blue", "darkred", "red", "pink"), 
       legend = c(paste0("study1: ever e > ", 1/alpha), 
                  paste0("study2: ever e > ", 1/alpha),
                  paste0("meta: ever e > ", 1/alpha),
                  paste0("study1: never e > ", 1/alpha),
                  paste0("study2: never e > ", 1/alpha),
                  paste0("meta: never e > ", 1/alpha)), cex = 0.6)

simeValuesNULL.study1.sample <- simeValuesNULL.study1[1:numPlot, ]
simeValuesNULL.study2.sample <- simeValuesNULL.study2[1:numPlot, ]
simeValueMetaNULL.sample <- simeValueMetaNULL[1:numPlot, ]
for (i in 1:nrow(simeValueMetaNULL.sample)){
  x <- which(!is.na(simeValueMetaNULL.sample[i, ]))
  y <- na.omit(simeValueMetaNULL.sample[i, ])
  myColour.simeValueMetaNULL <- ifelse(any(y > (1/alpha)), "blue", "pink")
  lines(x, log(y)/log(2) - 0.07, col = myColour.simeValueMetaNULL, lwd = 2) 
  x <- which(!is.na(simeValuesNULL.study1.sample[i, ]))
  y <- na.omit(simeValuesNULL.study1.sample[i, ])
  myColour.study1 <- ifelse(any(y > (1/alpha)), "darkgreen", "darkred")
  lines(x, log(y)/log(2), col = myColour.study1, lwd = 2)
  x <- which(!is.na(simeValuesNULL.study2.sample[i, ]))
  y <- na.omit(simeValuesNULL.study2.sample[i, ])
  myColour.study2 <- ifelse(any(y > (1/alpha)), "green", "red")
  lines(x, log(y)/log(2), col = myColour.study2, lwd = 2)
}
```

The plot depicts the evidence of the meta-analysis e-value. The first 40 days of each path equal the e-values of the first study. After results from the second study come in, the meta-analytic e-value starts pooling the evidence. Although the rate at which either the first or the second study reaches \(e > `r 1/alpha`\) will be larger than \(\alpha\) due to multiple testing (shown in dark green en light green), no more than \(\alpha\)\% of the pooled meta-analysis e-value paths will reach the threshold (shown in blue).

We observe that adding data will only decrease the e-value under the null hypothesis. Pooling evidence means adding data by multiplying e-values. So, also pooled e-values decrease, and due to the larger amounts of data available at each time, even more than single study e-values (indicated by the many pink e-value paths in the bottom of the graph).

```{r, echo = F}
meta.typeIerror <- mean(sapply(1:numSim, function(i) any(simeValueMetaNULL[i, ] > (1/alpha), na.rm = T)))*100
```

The type-I error of the pooled e-value in the above simulation is therefore controlled at ``r meta.typeIerror`%`.


## Summary: General recommendations Safe logrank test

We recommend to pool the evidence available at all interim times at which pooling is practically possible. The analysis per trial stratum should be performed using a Safe logrank test and results -- from multiple strata and trials at different interim stages -- should be pooled by multiplying their e-values. Type-I error is controlled if the multiplied e-values are compared to a prespecified threshold of \(\frac{1}{\alpha}\). The available evidence could be assessed based on the pooled e-values and the individual trial e-values. The pooled product of e-values should be leading, but additional requirements on individual trial e-values can be posed since those will only decrease type-I error.

So to preserve type-I error control, the statistical analysis plan needs to prespecify the design of the Safe test: minimally setting \(\alpha\), whether the test is two-sided, or left- or right-sided, and what is the effect size of minimum clinical relevance in terms of hazard ratio.

Larger e-values indicate more evidence against the null hypothesis of no difference between the treated and control groups. Trends in evidence can be inspected by plotting the e-values on a logarithmic scale at the various interim analyses. An increasing trend is a good indication of conditional power to cross the threshold in the future. Visual inspection of this trend can be used to substantiate preparations needed to act as soon as the threshold is crossed. Formal conditional power analyses can also be performed, as well as futility analyses, and will soon be available as straightforward functions in the `safestats` R package.

The Safe logrank test enables unlimited interim analyses and pooling of strata and trials. Any possible dependency between the timing of the analysis, the number of trials included and the results of those trials is allowed for. This means that to achieve optimal efficiency, we encourage intermediate decisions: Use the interim analyses to decide whether new trials are needed or should be discouraged. Use the observation of a large difference between the two groups in one stratum to increase enrollment in that stratum. Use the number of events observed in one stratum to optimize the power-design of the e-value in the next stratum. Note that the \(\frac{1}{\alpha}\) threshold as well as the design of the e-value need to be specified before observing the data that is tested. But if that is the case, anything is possible.



## References

