---
title: "Constraints, Interim Estimates and Conditional Power Functions for the Optimal Conditional Error Function"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using Constraints and Interim Estimates for the Optimal Conditional Error Function}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.dim = c(6,4)
)
```

```{r setup}
library(optconerrf)
```

# Summary

This vignette provides examples on some slightly more advanced methods which can be applied in the framework of the optimal conditional error function.
Firstly, the use of constraints for the optimal conditional error function is illustrated and afterwards, the use of an interim estimate for the effect at which the conditional power should be achieved as well as the implementation of a data-dependent conditional power function are discussed.

# Constraints

In the framework of two-stage adaptive designs, it may not be feasible to allow the second-stage information to fall below a certain minimum $I_{2,\min}$, e.g., if a certain amount of data on safety or secondary endpoints are required.
Conversely, budgeting considerations may impose an upper limit $I_{2, \max}$ on the information that can be collected in the second-stage.
Instead of ignoring these constraints in the planning phase and applying them "post-hoc", they can be accounted for in the minimisation procedure of the optimal conditional error function.
For illustration purposes, we use a trial design with overall type I error rate $\alpha=0.025$, early decision boundaries $\alpha_1 = 0.000158$ and $\alpha_0 = 0.5$ with conditional power $CP=0.9$ at the fixed alternative $\Delta_1 = 0.25$.
After collecting a first-stage information of $I_1=85$, the information of the second-stage of the trial should lie between $I_{2, \min} = 85/2$ and $I_{2, \max} = 85\cdot 2$.
To implement this design using the optimal conditional error function with the maximum likelihood ratio approach, the following code is required:

```{r}
trialDesignConstraints <- getDesignOptimalConditionalErrorFunction(
  alpha = 0.025,
  alpha1 = 0.000158,
  alpha0 = 0.5,
  conditionalPower = 0.9,
  delta1 = 0.25,
  useInterimEstimate = FALSE,
  likelihoodRatioDistribution = "maxlr",
  firstStageInformation = 85,
  minimumSecondStageInformation = 85/2,
  maximumSecondStageInformation = 85*2
)
```

The constraints are supplied by the arguments `minimumSecondStageInformation` and `maximumSecondStageInformation`.
Plotting the second-stage information of this design highlights the application of the constraints (represented by the gray dashed lines):

```{r, fig.dim = c(6,4)}
# type = 2 for second-stage information
plot(trialDesignConstraints, type = 2) +
  ggplot2::geom_hline(yintercept = c(85/2, 85*2), linetype = "dashed", colour = "gray")
```

Alternatively, the constraints may be directly specified on the conditional error function as $\alpha_{2,\min}$ and $\alpha_{2,\max}$ due to

$$
\alpha_{2,\min} = \Phi(\Phi^{-1}(CP)- \Delta_1\sqrt{I_{2,\max}})
$$
and
$$
\alpha_{2,\max} = \Phi(\Phi^{-1}(CP)- \Delta_1\sqrt{I_{2,\min}}).
$$
Note that constraining the *maximum* second-stage information corresponds to imposing a constraint on the *minimum* conditional error and that constraining the *minimum* second-stage information corresponds to imposing a constraint on the *maximum* conditional error.
Plugging in $I_{2, \max} = 85\cdot 2$ and $I_{2,\min} = 85/2$, we get $\alpha_{2,\min} \approx 0.02396$ and $\alpha_{2,\max} \approx 0.36382$ and can create the same design, now using `minimumConditionalError` and `maximumConditionalError` for the constraints:

```{r}
trialDesignConstraintsConditionalError <- getDesignOptimalConditionalErrorFunction(
  alpha = 0.025,
  alpha1 = 0.000158,
  alpha0 = 0.5,
  conditionalPower = 0.9,
  delta1 = 0.25,
  useInterimEstimate = FALSE,
  likelihoodRatioDistribution = "maxlr",
  firstStageInformation = 85,
  minimumConditionalError =  0.02396,
  maximumConditionalError = 0.36382
)
```

Comparing the expected second-stage information of the designs for $\Delta=0.25$ yields essentially the same values, with a minor difference resulting from rounding of $\alpha_{2,\min}$ and $\alpha_{2,\max}$. 

```{r}
getExpectedSecondStageInformation(
  trialDesignConstraints, 
  likelihoodRatioDistribution = "fixed", 
  deltaLR = 0.25
)
getExpectedSecondStageInformation(
  trialDesignConstraintsConditionalError, 
  likelihoodRatioDistribution = "fixed",
  deltaLR = 0.25
)
```

It is also possible to specify constraints on *both* scales simultaneously, i.e., on the second-stage information and on the conditional error function at the same time.
If this is done, the more restrictive constraint is applied. 
However, when using a fixed effect size $\Delta_1$ and constant conditional power $CP$, the scales have a one-to-one correspondence (as seen above in the formulas for $\alpha_{2,\min}$ and $\alpha_{2,\max}$) and supplying one of them is sufficient. 

Note that the specification of inappropriate constraints can make it impossible to fully exhaust (and not exceed) the $\alpha$ level, i.e., find a level constant for the design.
(As an example, consider a design with $\alpha=0.025$, $\alpha_1=0$ and $\alpha_0=1$. Since the full $\alpha$ level is to be spent in the second stage, it follows from the level condition (see `?getLevelConstant`) that $0.025 = \int_0^1 \alpha_2(p_1)dp_1$. The exemplary constraint $\alpha_{2,\max}=0.01$ is inappropriate since when applying it, $\int_0^1 \min(0.01, \alpha_2(p_1))dp_1 \leq \int_0^1 0.01 dp_1 = 0.01 \neq 0.025$ and the level condition cannot be fulfilled.)

# Interim Estimate for Treatment Effect

In certain applications, using a fixed effect $\Delta_1$ at which the conditional power (be it constant or a function, more details below) should be achieved may not be appropriate.
For example, we may wish to estimate this effect from the first-stage data if only little or unreliable evidence regarding the treatment effect is available during initial trial planning.
A simple estimate for the interim effect is given by:

$$
\hat\Delta_1(p_1) = \Phi^{-1}(1-p_1)/ \sqrt{I_1}.
$$
An issue with this estimate is the fact that it may become negative for $p_1 > 0.5$.
However, negative effects are part of the null hypothesis $H_0: \Delta \leq 0$ and cannot reasonably be powered for.
Furthermore, a clinical trial typically operates with a minimally clinically relevant effect size $\Delta_{1, \min} > 0$, which marks the lower boundary of effects considered worthwhile detecting.
For this reason, it is generally advisable to cut off interim effects that fall below $\Delta_{1,\min}$ (and instead use $\Delta_{1, \min}$).
An upper limit $\Delta_{1, \max}$ may also be put in place to deselect overly optimistic interim effects (and instead use $\Delta_{1, \max}$).
Implementing both cutoffs leads to the truncated interim estimate:

$$
\tilde\Delta_1(p_1) = \min(\max(\hat\Delta_1(p_1), \, \Delta_{1,\min}),\, \Delta_{1,\max}).
$$
When creating a design object with `getDesignOptimalConditionalErrorFunction()`, the use of an interim estimate for the treatment effect at which the conditional power should be achieved is specified by setting `useInterimEstimate = TRUE`.
When an interim estimate is used, the lower limit $\Delta_{1,\min}$ must be specified by the argument `delta1Min`.
An upper limit is not required, but can optionally be specified via `delta1Max`.
Both arguments may also be specified on the non-centrality parameter scale with $\vartheta_{1,\min} = \Delta_{1, \min}\cdot\sqrt{I_1}$ as `ncp1Min` and $\vartheta_{1,\max} = \Delta_{1, \max}\cdot\sqrt{I_1}$ as `ncp1Max`. 

As an example, consider the same trial setting as above with $\alpha=0.025$, $\alpha_1=0.000158$, $\alpha_0=0.5$, $CP=0.9$, $I_1 = 85$ and the maximum likelihood ratio.
Now, instead of using a fixed effect $\Delta_1$, we wish to estimate this effect from the data with a minimally clinically relevant effect size of $\Delta_{1, \min}=0.15$.
Additionally, we wish to deselect interim estimate effects above $\Delta_{1,\max} = 0.35$.

The design object is then created as:

```{r}
trialDesignInterimEstimate <- 
  getDesignOptimalConditionalErrorFunction(
    alpha = 0.025,
    alpha1 = 0.000158,
    alpha0 = 0.5,
    conditionalPower = 0.9,
    useInterimEstimate = TRUE,
    delta1Min = 0.15,
    delta1Max = 0.35,
    firstStageInformation = 85,
    likelihoodRatioDistribution = "maxlr"
  )
```

Plotting the optimal conditional error function and second-stage information from this design yields:

```{r}
plot(trialDesignInterimEstimate, type = 1)
plot(trialDesignInterimEstimate, type = 2)
```

Note that for a design with an early efficacy boundary $\alpha_1>0$, an upper limit for the effect is always implicitly given by $\Delta_{1,\max} = \Phi^{-1}(1-\alpha_1)/\sqrt{I_1}$, which is $\Phi^{-1}(1-0.000158)/\sqrt{85}=0.390672$ in our example.
While this always applies for the upper limit, a futility boundary $\alpha_0 < 1$ may not always correspond to a lower limit $\Delta_{1, \min}>0$ (more precisely, the implicit boundary $\Delta_{1,\min} \leq 0$ would occur if $0.5 \leq \alpha_0 < 1$).
For this reason, it is mandatory that $\Delta_{1,\min}>0$ (or, alternatively, $\vartheta_{1,\min}>0$) are explicitly specified, regardless of the futility boundary $\alpha_0$.

It is also possible to specify constraints on either or both the optimal conditional error scale or the second-stage information scale for designs implementing an interim estimate.
Rescaling the constraints, as opposed to the above case using a fixed effect size $\Delta_1$, must now take into account the interim estimate $\tilde{\Delta}_1(p_1)$ and requires the use of the constraint functions:


$$
\alpha_{2,\min}(p_1) = \Phi(\Phi^{-1}(CP)- \tilde{\Delta}_1(p_1)\sqrt{I_{2,\max}})
$$
and

$$
\alpha_{2,\max}(p_1) = \Phi(\Phi^{-1}(CP)- \tilde{\Delta}_1(p_1)\sqrt{I_{2,\min}}).
$$

Since $I_{2, \max} \geq I_{2,\min}$, it follows that $\alpha_{2,\min}(p_1) \leq \alpha_{2,\max}(p_1)$.
Note that both constraint functions are *non-decreasing* in $p_1$, since $\tilde{\Delta}_1(p_1)$ is non-increasing in $p_1$.
To ensure that the resulting conditional error function is nevertheless non-increasing, we may exploit a property of $\tilde{\Delta}_1(p_1)$:
The smallest effect size of the interim estimate is always given by $\Delta_{1,\min}$ (or may occur for $p_1 = \alpha_0$).
As the smallest effect size corresponds to the largest second-stage information, we can plug in the smallest effect in the formula for $\alpha_{2,\min}(p_1)$ and achieve a constant constraint.
If an upper limit for the effect size is provided, either via $\Delta_{1,\max}$ or at $p_1 = \alpha_1$, the same shortcut is applicable for $\alpha_{2,\max}(p_1)$.

# Using a Conditional Power Function

The conditional power describes the probability of a successful trial outcome at the final analysis (i.e., rejection of $H_0$) given the interim data and an effect assumption, which can either be pre-defined or derived from the interim data (see the section above).
Aiming to achieve the same value for the conditional power, e.g., 80%, regardless of the interim results of the trial may not be intuitive.
Instead, one may aim for a higher conditional power if the interim results are promising and settle for a lower conditional power for poor interim results.
Such an approach can be realised through the use of a conditional power function which depends on the interim data.
Generally, a small first-stage p-value may warrant a large conditional power, whereas for weaker first-stage evidence, a lower conditional power may be more appropriate.
Thus, the conditional power function should be decreasing (or at least non-increasing) in the first-stage p-value.

Importantly, the use of a (decreasing) conditional power function may lead to a non-monotone second-stage information function.
For a fixed conditional power, the second-stage information required to achieve that conditional power will always be increasing in $p_1$ (provided the conditional error function is non-increasing in $p_1$).
A decreasing conditional power function may impact the second-stage information function in a way that it increases for small first-stage p-values, reaches a maximum and decreases for larger first-stage p-values.
In `optconerrf`, the user may freely specify a conditional power function to `getDesignOptimalConditionalErrorFunction()` via the argument `conditionalPowerFunction`.
In the current implementation, the first-stage p-value is passed as the first and only argument to the conditional power function.
For illustration purposes, we consider the exemplary conditional power function $CP(p_1) = \Phi(1-p_1)$, which is decreasing in the first-stage p-value $p_1$:

```{r}
myConditionalPowerFunction <- function(firstStagePValue) {
  return(pnorm(1 - firstStagePValue))
}
```

A simple plot shows the properties of this function:

```{r}
p1 <- seq(0, 1, 0.01)
plot(p1, myConditionalPowerFunction(p1), type = "l")
abline(h = myConditionalPowerFunction(c(0, 1)), col = "blue")
```

For small first-stage p-values, the conditional power function leads to values just above 0.8 and approaches 0.5 as the first-stage p-value becomes larger.
As above, we will use $\alpha = 0.025$, $\alpha_1 = 0.000158$ and $\alpha_0 = 0.5$.
Implementing this in a design object is straightforward:

```{r}
trialDesignConditionalPowerFunction <- 
  getDesignOptimalConditionalErrorFunction(
    alpha = 0.025,
    alpha1 = 0.000158,
    alpha0 = 0.5,
    conditionalPowerFunction = myConditionalPowerFunction,
    delta1 = 0.25,
    useInterimEstimate = FALSE,
    firstStageInformation = 100,
    likelihoodRatioDistribution = "maxlr"
  )
```

Plotting the optimal conditional error function resulting from this design yields a seemingly normal function:

```{r}
plot(trialDesignConditionalPowerFunction, type = 1)
```

The merits of using a conditional power function become more apparent when plotting the second-stage information resulting from the design:

```{r}
plot(trialDesignConditionalPowerFunction, type = 2)
```

This second plot shows that the second-stage information required to achieve the conditional power (function) first increases in the first-stage p-value, reaches a maximum and then decreases as the first-stage evidence becomes weaker.
This is different to a design with a constant conditional power, for which the second-stage information function is always non-decreasing (and typically, increasing) in the first-stage p-value.
Using a decreasing conditional power function thus allows for potential savings in first-stage information in case of weak first-stage evidence.

The application of *either* constraints *or* an interim estimate for the treatment effect is possible in combination with a conditional power function, however the combination of all three methods may lead to complications further elaborated in the next section.

# Limitations

The `optconerrf` package provides support for each of the three methods detailed above and for each combination of up to two methods.
However, using all three of them simultaneously may lead to issues.

To illustrate this, presume we wish to implement a design with:

 * $\alpha = 0.025$, $\alpha_1 = 0.000158$, $\alpha_0 = 0.5$
 * maximum likelihood ratio
 * first-stage information $I_1 = 100$
 * second-stage information constraint $I_{2,\min} = 50$
 * conditional power function $CP(p_1)=\Phi(1-p_1)$
 * an interim estimate with $\Delta_{1,\min}=0.25$
 
 
```{r}
combinedDesign <- 
  getDesignOptimalConditionalErrorFunction(
    alpha = 0.025,
    alpha1 = 0.000158,
    alpha0 = 0.5,
    conditionalPowerFunction = myConditionalPowerFunction,
    delta1Min = 0.25,
    useInterimEstimate = TRUE,
    firstStageInformation = 100,
    likelihoodRatioDistribution = "maxlr",
    minimumSecondStageInformation = 50
  )
```

Plotting the optimal conditional error function resulting from this design shows a non-monotone function:

```{r}
plot(combinedDesign)
```

This comes from the fact that now, the constraint functions are more difficult to handle than before.

We see that:

$$
\alpha_{2,\min}(p_1) = \Phi(\Phi^{-1}(CP(p_1)) - \tilde{\Delta}_1(p_1)\sqrt{I_{2,\max}})
$$
and

$$
\alpha_{2,\max}(p_1) = \Phi(\Phi^{-1}(CP(p_1)) - \tilde{\Delta}_1(p_1)\sqrt{I_{2,\min}})
$$

may not be monotone in $p_1$, because both $CP(p_1)$ and $\tilde{\Delta}_1(p_1)$ are non-increasing in $p_1$.
We can no longer use the workaround of plugging in $\Delta_{1,\min}$ (respectively, $\Delta_{1, \max}$) for the interim estimate, because it is no longer guaranteed that the maximum (respectively, minimum) second-stage information is achieved for this effect size (due to the conditional power function).
