---
title: "harmBounds"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{harmBounds}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

\code{\link[harmBounds:plot.harmbound]{plot()}}

```{r, echo = FALSE}
has_ggplot <- requireNamespace("ggplot2", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = has_ggplot
)

library(harmBounds)
library(ggplot2)
```

The harmBounds package calculates stopping probabilities, defines stopping boundaries and generates plots for safety monitoring using an event based approach.

The idea is that under the null hypothesis of no safety concern, safety events are expected to occur at a frequency proportional to the randomization ratio. Assuming a single event per patient, we can do one sample binomial exact tests on the proportion of events in the treatment arm. In case of several events, we can use a beta-binomial framework, assuming a certain intra-class correlation coefficient. Evidence that the proportion of events in the treatment arm is higher than what is expected from randomization indicates a safety problem.

The safety monitoring can be done continuously after every event or at a pre-specified total number of events.
The nominal test-wise alpha would be calibrated to obtain desired properties, such as a certain proportion of stopping under alternative scenarios of some degree of safety problems (the power). An overall type I error control can also be implemented but we would usually not recommend that for safety testing, as the consequence of a type II error (not stopping for safety if the treatment is not safe) is arguably worse than of a type I error (stopping if the treatment is safe).

The procedure is computationally efficient and fully reproducible, as all calculations are exact and derived in closed form from the binomial or the beta-binomial distribution.


## Calculation of boundaries

Stopping boundaries can be calculated with function `getHarmBound`.
Let's assume a trial with 10 safety interim analyses after every 10 events (combined over both arms) up to a total of 100 events and a nominal test-wise alpha of 0.025.

```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, pH0 = 0.5)
hb
```


The data frame *bounds* specifies the boundaries at each interim analysis, *events_treatment* is the minimal number of events in the treatment arm that would lead to a stopping of the trial. The list *stopprob* contains stopping probability at each analysis for the null and optionally alternative hypotheses (see below). Data frame *opchar* shows the cumulative stopping probability and the expected number of events for each hypothesis.

In the example above, the overall type I error of the safety testing is `r round(hb$opchar$cum_stop_prob*100,1)`%. The expected number of events (i.e. the average if the trial would be repeated) is `r round(hb$opchar$expected_events,0)`.
Typically, the trial might not be stopped at the last interim analysis, which can be specified using the *maxevents* option. The option only influences the expected events.

```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, pH0 = 0.5,
  maxevents = 150)
hb$opchar
```

## Alternative specification of rejection criteria

Instead of using an alpha for each test, we can also control the family-wise type I error rate 
  at a specific level (even though we do not necessarily recommend that).
  Refer to the section on [Choosing the test-wise alpha](#sec-alpha) for more details.
  
  
```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_total = 0.05, pH0 = 0.5,
  maxevents = 150)
hb$opchar
```  

The overall stopping probability under H0 (the type I error) is `r round(hb$opchar$cum_stop_prob*100,1)`%, i.e. as close to 5% as possible (given the discrete nature of the test).If a strict control at 5% would required a *alpha_total* of 0.049 would have to be chosen.

Another more sensible option would be to target the type II error (or the power) for a specific alternative. This option requires the specification of an alternative, e.g. assuming that 60% of the events occuring in the experimental arm would indicate a safety problem (see the section on [Operating characteristics](#sec-opchar) for more details about the alternatives).

```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), power = 0.8, 
  pH0 = 0.5, pH1 = 0.6, maxevents = 150)
hb$opchar
```  

Specifying *alpha_total* or *power* is much slower as it relies on an iterative procedure ro determine the alpha used at each step. It could make sense to first determine the *alpha_test* and then use it directly in `getHarmBound`.

This can be done separately using function `getAlphaPerTest`:

```{r}
alphaPerTest<-getAlphaPerTest(nevents = seq(10, 100, by = 10), 
  alpha_total = 0.05, pH0 = 0.5)

alphaPerTest

alphaPerTest<-getAlphaPerTest(nevents = seq(10, 100, by = 10),
  power = 0.8,  pH0 = 0.5, pH1 = 0.6)

alphaPerTest

```


## Plotting of the boundaries

The boundaries can be plotted using `harmboundPlot` (or the harmbound.plot-method):  

```{r message=FALSE, warning=FALSE}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, pH0 = 0.5,
  maxevents = 150)
plot(hb)
```

Where the bars indicate the time points of the interim analysis with the rejection region in red,
  and the dashed line represents the expectation (can be removed with 'H0line = FALSE')

Continuous monitoring could also be implemented:

```{r message=FALSE, warning=FALSE}
hb<-getHarmBound(nevents = 0:100, alpha_test = 0.025, pH0 = 0.5, 
  maxevents = 150)
plot(hb) 
```

Observed data can be added as vector with 0 and 1, indicated the sequence of the arms in which events occurred (0 being the control and 1 the treatment arm). In this example, the boundary is not breached.

```{r message=FALSE, warning=FALSE}
set.seed(123)
eventgroups<-rbinom(n = 100, size = 1, prob = 0.5)
plot(hb,observed=eventgroups) 
```

## Operating characteristics  {#sec-opchar}

Stopping probabilities and expected number of events can be obtained for alternative scenarios with the `getHarmBound` function.
The alternative hypothesis can be specified as 

- pH1: the proportion of the events in the treatment arm, 
  with 0.5 being the null scenario for a 1:1 rando, and numbers from 0.5 to 1 indicating a 
  safety problem.
  
- rrH1: the risk ratio (treatment / control), with 1 being the null scenario and numbers >1 indicating a safety problem.

- rdH1: the risk difference (treatment minus control), with 0 being the null scenario and numbers >0 indicating a safety problem. Here the control proportion (r0) has to be specified.

- orH1: the odds ratio (treatment / control), with 1 being the null scenario and number >1 indicating a safety problem. Here the control proportion (r0) has to be specified.

To stay in the event-based setting, we recommend using the first two options (pH1 or rrH1).

```{r}
#with proportion of events in the treatment arm:
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, pH0 = 0.5,
  pH1 = 0.6, maxevents = 150)
hb
```
Note that *stopprob* is now of length two with a data frame for each condition.

We would stop in `r round(100*hb$opchar$cum_stop_prob[hb$opchar$hyp=="H0"],1)`% of the trials under the null (no safety problem) and in `r round(100*hb$opchar$cum_stop_prob[hb$opchar$hyp=="H1"],1)`% under the alternative (safety problem with 60% of the events in the treatment arm).
Because of the higher probability of early stopping, the expected number of events is reduced from `r round(hb$opchar$expected_events[hb$opchar$hyp=="H0"],0)` under H0 to `r round(hb$opchar$expected_events[hb$opchar$hyp=="H1"],0)` under H1.

The (cumulative) stopping probabilities at each step can be plotted using *absstopPlot* or *cumstopPlot* (or the harmbound.plot-methods). We now using more than one alternative to make it more interesting:

```{r, fig.height=6, fig.width=8}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, 
  pH0 = 0.5, pH1 = c(0.6, 0.7, 0.8) , maxevents = 150)
plot(hb, which = "abs_stopping")
```

```{r}
plot(hb, which = "cum_stopping")
```

We can also specify a vector of alternatives and plot the cumulative stopping probabilities or the expected number of events for each alternative.

```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025, pH0 = 0.5, 
  pH1 = seq(0,1,l=100), maxevents = 150)
plot(hb, which = "opchar_stop")
plot(hb, which = "opchar_n")
```

## Choosing the test-wise alpha {#sec-alpha}

To choose the test-wise alpha, we recommend to control the power rather than the type I error.
  I.e. to check that a sufficient proportion of trials are stopped under 
  an appropriate alternative hypothesis (reflecting a safety problem).

A target power can be specified directly in `getHarmbound` (see above).

To get a more complete picture, we would check the stopping probabilities for a grid of alphas:

```{r, warning=FALSE}
alphaPerTest <- getAlphaPerTest(nevents = seq(10, 100, by = 10), pH0 = 0.5,
  alpha_total = 0.05)
alist<-c(0.001,0.01,alphaPerTest,0.025,0.05)

hbl<-lapply(alist,function(x)
  getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = x,
  pH0 = 0.5, pH1 = seq(0,1,l=100), maxevents = 150)$opchar)

hbd<-data.frame(do.call(rbind,hbl),alpha=rep(alist,each=nrow(hbl[[1]])))
hbd$alpha<-as.factor(round(hbd$alpha,4))

ggplot(hbd, aes(x = p, y = cum_stop_prob, colour=alpha)) + 
  	geom_line() +
    ylab("Stopping probability") +
    xlab("Proportion of events in treatment arm") +
  	scale_x_continuous(breaks=seq(0,1,by=0.1),limits=c(0.2,0.8)) +
  	scale_y_continuous(breaks=seq(0,1,by=0.2))

ggplot(hbd, aes(x = p, y = expected_events, colour=alpha)) + 
  	geom_line() +
    ylab("Expected number of events") +
    xlab("Proportion of events in treatment arm") +
  	scale_x_continuous(breaks=seq(0,1,by=0.1),limits=c(0.2,0.8)) +
  	scale_y_continuous(limits=c(0, 150))
```

or on the transformed risk ratio scale:

```{r, warning=FALSE}
alphaPerTest <- getAlphaPerTest(nevents = seq(10, 100, by = 10), pH0 = 0.5,
  alpha_total = 0.05)
alist<-c(0.001,0.01,alphaPerTest,0.025,0.05)

hbl<-lapply(alist,function(x)
  getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = x, pH0 = 0.5,
  rrH1 = seq(0,10,l=100), maxevents = 150)$opchar)

hbd<-data.frame(do.call(rbind,hbl),alpha=rep(alist,each=nrow(hbl[[1]])))
hbd$alpha<-as.factor(round(hbd$alpha,4))
hbd$lrr<-log(hbd$rr)

ggplot(hbd, aes(x = rr, y = cum_stop_prob, colour=alpha)) + 
    geom_line() +
    ylab("Stopping probability") +
    xlab("Risk ratio (treatment/control)") +
    scale_x_continuous(trans='log', breaks=c(0.5,0.75,1,1.5,2,3,4),
      minor_breaks=NULL,limits=c(0.5,5)) +
    scale_y_continuous(breaks=seq(0,1,by=0.2))
```


The test-wise alpha that leads to a reasonable stopping probability under an appropriate alternative can be selected.
However, the steepness of the curve is not influenced by the alpha. We can only choose between making more type I or II errors but not reduce the overall error rate.

```{r, echo=FALSE, include=FALSE}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.05, pH0 = 0.5,
  rrH1 = 1.5)$opchar

hb2<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = alphaPerTest, pH0 = 0.5,
  rrH1 = 1.5)$opchar
```

We could e.g. specify that at least 70% of the trials should be stopped if the risk of having a safety event increase by half (i.e. a risk ratio of 1.5), i.e. a test-wise alpha of close to 0.05. However, this would also lead to a stopping of `r round(hb$cum_stop_prob[hb$p==0.5]*100,1)`% of the trials under the null.

On the other hand, if we control the type I error at 5% using a test-wise alpha of `r round(alphaPerTest,4)`, only `r round(hb2$cum_stop_prob[hb2$p==0.6]*100,1)`% of the trials are stopped if there is a safety problem with a risk ratio of 1.5.


## Further paramaters

The number of interim analyses has some impact on steepness and location of the curve, but it is rather limited, as long as the last interim analysis is at the same time point.

```{r, warning=FALSE}
nelist<-list(c(100),c(50,100),c(10,20,50,100),
             seq(10, 100, by = 10),seq(10, 100, by = 1))

hbl<-lapply(nelist,function(x)
  getHarmBound(nevents = x, alpha_test = 0.025, pH0 = 0.5,
    pH1 = seq(0,1,l=100), maxevents = 150)$opchar)

nis<-unlist(lapply(nelist,length))
net<-unlist(lapply(nelist,function(x) paste(x,collapse = "/")))
net[[length(net)-1]]<-"every 10"
net[[length(net)]]<-"continous"
nism<-paste0(nis," (",net,")")
hbd<-data.frame(do.call(rbind,hbl),ni=rep(nism,each=nrow(hbl[[1]])))
hbd$ni<-factor(hbd$ni,levels=nism)

ggplot(hbd, aes(x = p, y = cum_stop_prob, colour=ni)) + 
    geom_line() +
    ylab("Stopping probability") +
    xlab("Proportion of events in treatment arm") +
    scale_x_continuous(breaks=seq(0,1,by=0.1),limits=c(0.2,0.8)) +
    scale_y_continuous(breaks=seq(0,1,by=0.2)) +
    labs(colour="Number of IAs")


```

To substantially decrease the total error rates (i.e. get steeper curves), we would need a higher total number of events: 

```{r, warning=FALSE}

nelist<-list(seq(10, 20, by = 10),seq(10, 50, by = 10), seq(10, 100, by = 10), 
  seq(10, 200, by = 10))

hbl<-lapply(nelist,function(x)
  getHarmBound(nevents = x, alpha_test = 0.025, pH0 = 0.5,
  pH1 = seq(0,1,l=100))$opchar)

nis<-unlist(lapply(nelist,max))
hbd<-data.frame(do.call(rbind,hbl),ni=rep(nis,each=nrow(hbl[[1]])))
hbd$ni<-as.factor(hbd$ni)

ggplot(hbd, aes(x = p, y = cum_stop_prob, colour=ni)) + 
  	geom_line() +
    ylab("Stopping probability") +
     xlab("Proportion of events in treatment arm") +
    scale_x_continuous(breaks=seq(0,1,by=0.1),limits=c(0.2,0.8)) +
    scale_y_continuous(breaks=seq(0,1,by=0.2)) +
    labs(colour="Total number of events")

```


## Correlated data 

If patients can have more than one event, the binomial boundaries are no longer 
valid because the actual variance is increased.
However, we can use a beta-binomial framework. That needs the specification
  of a intraclass correlation coefficient (*icc*).

The overdispersion factor (or design effect) by which the beta-binomial 
  variance exceeds the regular binomial variance corresponds to 
  1 + (n - 1) * icc and is printed by the function.
Note that the effective sample size is reduced by that factor which 
  results in a substantial loss of power, even for small ICCs.

```{r, echo=FALSE, results='hide', warning=FALSE}
capture.output(hb1<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025,
  pH0 = 0.5, pH1 = 0.6, maxevents = 150))
capture.output(hb2<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025,
  pH0 = 0.5, pH1 = 0.6, maxevents = 150, icc=0.05))

ps1<-hb1$opchar$cum_stop_prob[hb1$opchar$hyp=="H1"]
ps2<-hb2$opchar$cum_stop_prob[hb1$opchar$hyp=="H1"]

```  

For example, with an ICC of 0.05, we would only stop `r round(100*ps2)`% of the trials with 60% events in the treatment arm, compared to `r round(100*ps1)`% without correlation. 
```{r}
hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025,
  pH0 = 0.5, pH1 = 0.6, maxevents = 150)
hb$opchar

hb<-getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025,
  pH0 = 0.5, pH1 = 0.6, maxevents = 150, icc = 0.05)
hb$opchar
```  

The effective sample size can be calculated using the overdispersion factor,
  which captures the increase in variance due to the correlation (the design effect).
With an ICC of 0.05, the effective sample size is much smaller than
  the number of events: 

```{r}
nevents<-seq(10, 100, by = 10)
#overdispersion
ods<-1 + (nevents - 1) * 0.05
ods 
#effective sample size:
round(nevents/ods)
```
We can plot this loss of power for different ICCs:

```{r, warning=FALSE, results="hide"}
iccs<-c(0,0.01,0.02,0.05,0.1)

hbl<-lapply(iccs,function(x)
  getHarmBound(nevents = seq(10, 100, by = 10), alpha_test = 0.025,
  pH0 = 0.5, pH1 = seq(0,1,l=100), maxevents = 150, icc = x)$opchar)

hbd<-data.frame(do.call(rbind,hbl),icc=rep(iccs,each=nrow(hbl[[1]])))
hbd$icc<-as.factor(round(hbd$icc,2))

ggplot(hbd, aes(x = p, y = cum_stop_prob, colour=icc)) + 
  	geom_line() +
    ylab("Stopping probability") +
    xlab("Proportion of events in treatment arm") +
  	scale_x_continuous(breaks=seq(0,1,by=0.1),limits=c(0.2,0.8)) +
  	scale_y_continuous(breaks=seq(0,1,by=0.2))

```