---
title: "RMarkdown version of lab4-1"
author: "Constantin T Yiannoutsos"
date: "October 9, 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = F)
```

```{r, packages, echo=FALSE, results="hide", message=FALSE}
  options(width=72)
```

# Introduction

First, we import the library ```gsDesign``` into R.
```{r, enrollment, message=FALSE}
################################################
### Sequential monitoring of clinical trials ###
################################################
library(gsDesign)
```

# The R package ```gsDesign```

The function ```gsDesign``` can carry out lots of analyses but here we just focus on the calculation of boundaries. We tell $\verb|R|$ that we want 5 analyses $k=5$ (including interim and final). By default, the function ```gsDesign``` assumes equally spaced analyses. Also, we specify that we want two sided test (i.e. both lower and upper boundaries) and type-I error $\alpha=5$\%. Be careful here as the function always used one-sided type-I error. Thus, if you want a 5\% alpha, you need to specify ```alpha=0.025```!

## Pocock boundaries (equally-spaced analyses)

The following are the bounds proposed by Pocock for $k=5$ equally spaced interim analyses.
```{r, Pocockboundaries, message=FALSE}
# Pocock boundaries
x1 = gsDesign(k = 5,test.type = 2, alpha = 0.025, sfu = "Pocock")
gsBoundSummary(x1)
```

The boundaries are constant across the analyses $Z=\pm2.41$ and the p-value we have to use in a two-sided test is $2\times 0.0079=0.0158$. You can ignore the remaining output at this stage. We can also use the $\verb|plot|$ to produce a graph.


```{r, figure1, fig.width=4.5, fig.height=4.5, message=FALSE, fig.align='center', fig.cap='<b>Figure 1. </b>Pocock boundaries assuming $k=5$ looks and &alpha=0.05$.'}
par(mar=c(4,4,1,1))
ss1<-plot(x1, main = "")
ss2<-ss1+theme(legend.position = "none")+
  theme(axis.text=element_text(size=10))+
  theme(axis.title=element_text(size=15,face="bold"))
ss2
```

The following are the $Z$ statistics and the associated $p$ value boundaries
```{r, Pocockboundarylist, message=FALSE}
# Z's and p-values
cbind(Z = gsBoundSummary(x1)$Eff[gsBoundSummary(x1)$Val=="Z"],
      P = 2*gsBoundSummary(x1)$Eff[gsBoundSummary(x1)$Val=="p (1-sided)"])
```
Notice how far the p-value at the last (fifth) analysis is from the nominal 5\% alpha level!

Another issue that you need to realize is that the introduction of the Pocock interim review impacted the sample size. 

Recall the following line from the output above:

```{r, hiddenPocock, echo=FALSE, tidy=TRUE}
tmp<-capture.output(gsBoundSummary(x1))
cat(tmp[23],"\n")
```

In fact, the expected sample size is inflated by about 21\% compared to the fixed design (more on that later).

## O'Brien-Fleming boundaries (equally-spaced analyses)

To specify the O'Brien-Fleming method, just use ```sfu = "OF"```. Note that according to the O'Brien-Fleming method it's much more difficult to reject the null hypothesis in the beginning of the study (i.e. you require a higher level of evidence to do so). The advantage is that the p-value is close to the nominal level at the final analysis.

```{r, OFboundaries, message=FALSE}
# O'Brien-Fleming boundaries
x2 = gsDesign(k = 5,test.type = 2, alpha = 0.025, sfu = "OF")
gsBoundSummary(x2)
```

The following are the $Z$ statistics and the associated $p$ value boundaries

```{r, OFboundarylist, message=FALSE}
# Z's and p-values
cbind(Z = gsBoundSummary(x2)$Eff[gsBoundSummary(x2)$Val=="Z"],
      P = 2*gsBoundSummary(x2)$Eff[gsBoundSummary(x2)$Val=="p (1-sided)"])
```
Note that these are twice those listed in the previous output.  Also note that the p value for the final analysis is not too far from the nominal 5\% alpha level despite four interim analyses. We will see why in a minute.

Another issue has to do with the samle inflation caused by the O'Brien-Fleming procedure.  Recall the output from above:

```{r, hiddenOF, echo=FALSE, tidy=TRUE}
tmp<-capture.output(gsBoundSummary(x2))
cat(tmp[23],"\n")
```

We see that the sample size is inflated less than 3\% after four interim analysis (!) another advantage of the O'Brien-Fleming procedure.

```{r, figure2, fig.width=4.5, fig.height=4.5, message=FALSE, fig.align='center', fig.cap="<b>Figure 2.</b> O'Brien-Fleming boundaries assuming $k=5$ looks and &alpha=0.05$."}
par(mar=c(4,4,1,1))
ss1<-plot(x2,main = "")
ss2<-ss1+theme(legend.position = "none")+
  theme(axis.text=element_text(size=10))+
  theme(axis.title=element_text(size=15,face="bold"))
ss2
```

# Alpha spending functions
Because we want to be able to carry out an analysis at any point, including adding ad hoc analyses at any point during the execution of the study, we adopt the \textit{spending function} approach.  This approach literally governs how the alpha (Type-I error) will be spent during the study. 

## Example: Spending equal amount of error at each analysis
Consider the following situation where we expend the alpha level over $k=5$ analyses, so that $\alpha_1=0.01, \cdots, \alpha_5=0.01$.  

The situation is handled by the following \texttt{R} code:
```{r, linearalphaspend, message=FALSE}
t<-c(0.2,0.4,0.6,0.8,1.0)
x <- gsDesign(k=5, test.type=2, beta=0.2, 
              sfu=sfLinear, sfupar=c(t,t))
x
```

This is shown pictorially as follows:

```{r, linalphaspendfig, fig.width=4.5, fig.height=4.5, message=FALSE, warning=FALSE, fig.align='center', fig.cap="<b>Figure 3.</b> Bounds for equally spending &alpha=0.05$ over $k=5$ analyses."}
par(mar=c(4,4,1,1))
ss1<- plot(x, main="")
ss2<-ss1+theme(legend.position = "none")+
  theme(axis.text=element_text(size=10))+
  theme(axis.title=element_text(size=15,face="bold"))
ss2
```

Note here that, despite the fact that an equal amount of alpha is spent at each analysis, the bounds are not constant. This is because, while, for the first analysis

$$
P(|Z_1|>c_1)=\alpha_1=0.01
$$

so that $c_1=z_{1-a_1/2}$, for the second interim analysis, $c_2$ must satisfy the #conditional probability#

$$
P(|Z_2|>c_2|-c_1<Z_1<c_1)=\alpha_2=0.01
$$

Thus, the boundary points $c_1\neq c_2$; and so on.

## Continuous spending functions

In their 1983 paper,~\cite{KimDeMets1983} introduced the concept of a \textit{continuous} (alpha) spending function. This is a function that is defined in all points $t \in [0,1]$ (i.e., in all trial fractions from the beginning to the end of the study.)  There are several families of such functions:

### Approximate Pocock and O'Brien-Fleming spending functions

These are 
$$
\alpha(t)=\alpha\ln\left\{1+(e-1)t\right\}
$$
and
$$
\alpha(t)=2-2\Phi(z_{\alpha/4}/\sqrt{t})
$$
for the case of the Pocock and O'Brien-Fleming-like functions respectively.

### The _Power_ spending function family

These spending functions have the general form
$$
\alpha(t)=\min\left (\alpha t^\rho, \alpha\right )
$$

where $\rho>0$. Spending functions with $\rho=0.8$ approximate the Pocock spending function. (Note that this is close to the linear function employed in the example of Section~\ref{sec:linear} above), while $\rho=3$ approximates the O'Brien-Fleming function .  Power spending functions for several choices of $\rho$ are shown in the following Figure.

```{r, powerspendingfunfigno, fig.keep="none", message=FALSE, fig.align='center'}
# Plot of power family error spending functions
t<-seq(0,1,.05)

p08=0.05*t^0.8 # Pocock
p1=0.05*t # Linear
p3=0.05*t^3.0 # O-F

plot(t,p08, lty=2, type="l", ylab="Alpha spent", xlab="Trial fraction", main="")
lines(t, p1, lty=1, type="l", lwd=2)
lines(t, p3, lty=3, type="l", lwd=2)

legend(0,.05, legend = c(expression(paste(rho, " = ", 0.8, " (Pocock)")),
                         expression(paste(rho, " = ", 1, " (linear)")),
                         expression(paste(rho, " = ", 3, " (O-F)"))), 
       lty=c(2,1,3), lwd=c(2,2,2), title="Power spending functions")


```

```{r, powerspendingfunfigyes, echo=FALSE, fig.width=4.5, fig.height=4.5, message=FALSE, fig.align='center', fig.cap="<b> Figure 4. </b>Power family of spending functions for various choices of $\rho$."}
par(mar=c(4,4,1,1))
# Plot of power family error spending functions
t<-seq(0,1,.05)

p08=0.05*t^0.8 # Pocock
p1=0.05*t # Linear
p3=0.05*t^3.0 # O-F

par(mar=c(4,4,1,1))
plot(t,p08, lty=2, type="l", ylab="Alpha spent", xlab="Trial fraction", main="")
lines(t, p1, lty=1, type="l", lwd=2)
lines(t, p3, lty=3, type="l", lwd=2)

legend(0,.05, legend = c(expression(paste(rho, " = ", 0.8, " (Pocock)")),
                         expression(paste(rho, " = ", 1, " (linear)")),
                         expression(paste(rho, " = ", 3, " (O-F)"))), 
       lty=c(2,1,3), lwd=c(2,2,2), title="Power spending functions")


```

### The Hwang, Shi \& DeCani family of spending functions

Hwang, Shih \& DeCani (1990) proposed the following family of spending functions:
\begin{eqnarray*}
\alpha(t)=\left \{
\begin{array}{cc}
  \alpha(1-e^{-\gamma t})/(1-e^\gamma) & {\rm if} \gamma\neq 0 \\
  \alpha t & {\rm if} \gamma=0
\end{array}
\right .
\end{eqnarray*}
Alpha is spent more quickly the larger the value of $\gamma$.  This is shown in Figure~5.

```{r, HSDspendingfunfigno, fig.keep="none", message=FALSE}
# Plot of HSD error spending functions
t<-seq(0,1,.05)

g0=0.05*t # Pocock
g4=0.05*(1-exp(4*t))/(1-exp(4)) # O-F
g2=0.05*(1-exp(-2*t))/(1-exp(-2))

plot(t,g2, lty=2, lwd=2, type="l", 
     ylab="Alpha spent", xlab="Trial fraction", main="")
lines(t,g0, lty=1, lwd=2)
lines(t,g4, lty=3, lwd=2)

legend(0,.05, legend = c(expression(paste(gamma, " =   ", 2)),
                         expression(paste(gamma, " =   ", 0, " (Pocock)")),
                         expression(paste(gamma, " = ", -4, " (O-F)"))), 
       lty=c(2,1,3), lwd=c(2,2,2), title="HSD spending functions")
```

```{r, HSDspendingfunfigyes, echo=FALSE, fig.width=4.5, fig.height=4.5, message=TRUE, fig.align='center', fig.cap='<b>Figure 5.</b> HSD family of spending functions for various choices of &gamma.'}
# Plot of HSD error spending functions
t<-seq(0,1,.05)

g0=0.05*t # Pocock
g4=0.05*(1-exp(4*t))/(1-exp(4)) # O-F
g2=0.05*(1-exp(-2*t))/(1-exp(-2))

par(mar=c(4,4,1,1))
plot(t,g2, lty=2, lwd=2, type="l", 
     ylab="Alpha spent", xlab="Trial fraction", main="")
lines(t,g0, lty=1, lwd=2)
lines(t,g4, lty=3, lwd=2)

legend(0,.05, legend = c(expression(paste(gamma, " =   ", 2)),
                         expression(paste(gamma, " =   ", 0, " (Pocock)")),
                         expression(paste(gamma, " = ", -4, " (O-F)"))), 
       lty=c(2,1,3), lwd=c(2,2,2), title="HSD spending functions")
```

## Applications of spending functions

Spending functions were a ``game changer" in the effort to monitor a clinical trial.  The following simple example underlines this fact.

### Example: Inserting an interim analysis after the start of a study

Suppose that you have planned a study with $k=3$ total analyses (two interim and one final), to be carried out at the $\tau_1=20\%$, $\tau_2=50\%$ and $\tau_3=100\%$ trial fractions.  The analyses are going to be carried out at the $\alpha=0.025$ (one-sided) with power 90\% (i.e., $\beta=0.1$) and we will use a power spending functions with $\rho=1$ (i.e., a linear spending of the Type-1 error). The study design, as it stands at the start is as follows:

```{r, insertanalysis1, message=FALSE}
# Inserting an additional interim analysis
t<-c(0.2,0.5,1.0)
gsDesign(k=3, test.type=1, sfu=sfLinear, sfupar=c(t,t), timing=t)
```

We can see that the alpha level will be spent in the three analyses as

\begin{center}
\begin{tabular}{ccc}
\hline
$k$   & $z$ boundary & $\alpha$ spent\\
\hline
1 & 2.58 & 0.0050\\
2 & 2.38 & 0.0125\\
3 & 2.14 & 0.0250\\
\hline
\end{tabular}
\end{center}

and the final analysis will be carried out at a p value $p=0.0161$.

Now suppose that, after the second interim analysis at the $\tau_2=0.5$ trial fraction, we (or, more likely, the DSMB reviewing this study) decide to add one more interim analysis at the point $\tau^*_3=0.75$ (making now the final analysis at the time point $\tau^*_4=1.00$).  Spending functions allow us to do this almost effortlessly!  Consider the following output:

```{r, insertanalysis2, message=FALSE}
# adding an interim look at t=0.75
t<-c(0.2,0.5, 0.75, 1.0)
gsDesign(k=4, test.type=1, sfu=sfLinear, sfupar=c(t,t), timing=t)
```

There are a couple of remarkable things in the previous output that must be remarked on.  First of all, the addition of the third analysis only depended on the alpha level spent cumulatively up to the second analysis (had a revision been necessary, this would be have been an unworkable situation). On the other hand, adding the third interim analysis, and the additional alpha spent prior to the final analysis, required the final analysis to be carried out at an alpha level of 0.062 (or a p value of 0.0124) rather than 0.0125 like in the previous situation (where the final analysis was the third analysis). Thus, expending more alpha level has a downside: It requires a higher threshold of evidence at the final analysis if you do not interrupt the study during any of the interim analyses.

# Effect of group sequential monitoring on sample size

The effect of introducing a number of interim analysis on top of a fixed design invariably results in increasing the required sample size to preserve the power.

Consider the three following cases shown in Figure 6
```{r, revisedss, fig.width=4.5, fig.height=4.5, message=FALSE, fig.align='center', fig.cap='<b>Figure 6.</b> Sample size inflation (compared to the fixed design) according to &gamma in the HSD family of spending functions'}
x<-array(dim=11)
for(i in 1:11) {
x[i]<-gsDesign(k=5, test.type=2, alpha=0.025, sfu=sfHSD, sfupar=-i+1)$en[1]
}

par(mar=c(4,4,1,1))
plot(seq(0,-10,-1), x, type="b", ylab="Sample size inflation", 
     xlab=expression(paste(gamma, " parameter in the HSD spending function")))
```

To figure out why there is this significant inflation of the sample size, we can consider the rates of spending (recall Figure 5 above). Long story short, the earlier the alpha is spent, the higher the sample size inflation compared to the fixed design.  This is the reason why the Pocock prodedure resulted in such a large sample size inflation compared to the O'Brien-Fleming