---
title: "Test of the function gs_power_npe"
output:
  rmarkdown::html_document:
    toc: true
    toc_float: true
    toc_depth: 2
    number_sections: true
    highlight: "textmate"
    css: "custom.css"
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Test of the function gs_power_npe}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(gt)
library(dplyr)
library(tibble)
library(testthat)
library(gsDesign)
devtools::load_all()

setwd("./..")
my_path <- "./old_function/"
source_files <- list.files(my_path, "*.R$")
sapply(paste0(my_path, source_files), source)
```


# Test 1: Examples from spec

## Default

The default of `gs_power_npe` is a single analysis with type I error controlled.

```{r}
x1 <- gs_power_npe(theta = 0) %>% filter(bound == "upper")
x2 <- gs_power_npe_(theta = 0) %>% filter(Bound == "Upper")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  mutate(`Computated from` = c("new version", "old version")) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Fixed bound

```{r}
x1 <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_b,
  upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b,
  lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_b,
  upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b,
  lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Same fixed efficacy bounds, no futility bound (i.e., non-binding bound), null hypothesis

```{r}
x1 <- gs_power_npe(
  theta = rep(0, 3),
  info = (1:3) * 40,
  upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lpar = rep(-Inf, 3)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = rep(0, 3),
  info = (1:3) * 40,
  upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lpar = rep(-Inf, 3)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Fixed bound with futility only at analysis 1; efficacy only at analyses 2, 3

```{r}
x1 <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_b,
  upar = c(Inf, 3, 2),
  lower = gs_b,
  lpar = c(qnorm(.1), -Inf, -Inf)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_b,
  upar = c(Inf, 3, 2),
  lower = gs_b,
  lpar = c(qnorm(.1), -Inf, -Inf)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Spending function bounds

```{r}
# Lower spending based on non-zero effect
x1 <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Same bounds, but power under different theta

```{r}
x1 <- gs_power_npe(
  theta = c(.15, .25, .35),
  info = (1:3) * 40,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = c(.15, .25, .35),
  info = (1:3) * 40,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -1, timing = NULL)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Two-sided symmetric spend, O'Brien-Fleming spending

Typically, 2-sided bounds are binding

```{r}
x1 <- gs_power_npe(
  theta = rep(0, 3),
  info = (1:3) * 40,
  binding = TRUE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
) %>% mutate(`Computated from` = "new version")

x2 <- gs_power_npe_(
  theta = rep(0, 3),
  info = (1:3) * 40,
  binding = TRUE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
x1 %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

## Re-use these bounds under alternate hypothesis

Always use binding = TRUE for power calculations

```{r}
xx1 <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  binding = TRUE,
  upar = (x1 %>% filter(bound == "upper"))$z,
  lpar = -(x1 %>% filter(bound == "upper"))$z
) %>% mutate(`Computated from` = "new version")

xx2 <- gs_power_npe_(
  theta = c(.1, .2, .3),
  info = (1:3) * 40,
  binding = TRUE,
  upar = (x1 %>% filter(bound == "upper"))$z,
  lpar = -(x1 %>% filter(bound == "upper"))$z
) %>% mutate(`Computated from` = "old version")
```

```{r, echo=FALSE}
xx1 %>%
  union_all(xx2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

# Test 2: Fixed Design

The power at the following analysis is expected at 0.975.

```{r}
gs_power_npe(
  theta = 0,
  upper = gs_b, upar = qnorm(0.025),
  lower = gs_b, lpar = -Inf
)
```

# Test 3: `info` != `info0` != `info1`

If one inputs `info` in `upar`

```{r}
x1_a <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 80, info0 = (1:3) * 90 + 10, info1 = (1:3) * 70 - 5, info_scale = 0,
  upper = gs_b, upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b, lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "new version", `Info scale` = 0)
x1_b <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 80, info0 = (1:3) * 90 + 10, info1 = (1:3) * 70 - 5, info_scale = 1,
  upper = gs_b, upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b, lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "new version", `Info scale` = 1)
x1_c <- gs_power_npe(
  theta = c(.1, .2, .3),
  info = (1:3) * 80, info0 = (1:3) * 90 + 10, info1 = (1:3) * 70 - 5, info_scale = 2,
  upper = gs_b, upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b, lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "new version", `Info scale` = 2)
x2 <- gs_power_npe_(
  theta = c(.1, .2, .3),
  info = (1:3) * 80, info0 = (1:3) * 90 + 10, info1 = (1:3) * 70 - 5,
  upper = gs_b, upar = gsDesign::gsDesign(k = 3, sfu = gsDesign::sfLDOF)$upper$bound,
  lower = gs_b, lpar = c(-1, 0, 0)
) %>% mutate(`Computated from` = "old version")
```

```{r}
x1_a %>%
  union_all(x1_b) %>%
  union_all(x1_c) %>%
  union_all(x2 %>% rename(bound = Bound, probability = Probability, z = Z, analysis = Analysis)) %>%
  arrange(analysis) %>%
  group_by(analysis, bound) %>%
  gt() %>%
  tab_style(
    style = list(cell_fill(color = "#d3edeb")),
    locations = cells_body(rows = `Computated from` == "old version")
  )
```

# Test 3: Developer Tests

## 1-sided test

```{r}
r <- 80
x <- gs_power_npe(
  theta = 0,
  info = (1:3) * 400,
  binding = FALSE, r = r,
  upper = gs_b, # gs_spending_bound,
  upar = gsDesign(k = 3, test.type = 1, sfu = sfLDOF)$upper$bound,
  lower = gs_b,
  lpar = rep(-Inf, 3)
) %>% filter(bound == "upper")

y <- gs_power_npe_(
  theta = 0,
  info = (1:3) * 400,
  binding = FALSE, r = r,
  upper = gs_b, # gs_spending_bound,
  upar = gsDesign(k = 3, test.type = 1, sfu = sfLDOF)$upper$bound,
  # list(par = list(sf = gsDesign::sfLDOF, param = NULL, total_spend = 0.025)),
  lower = gs_b,
  lpar = rep(-Inf, 3)
) %>% filter(Bound == "Upper")

z <- gsProbability(
  k = 3,
  theta = 0,
  n.I = (1:3) * 400,
  b = gsDesign(k = 3, test.type = 1, sfu = sfLDOF)$upper$bound, a = rep(-20, 3), r = r
)
```

```{r, echo=FALSE}
tibble(
  `Calculated from` = rep(c("new version", "old version", "gsDesign"), each = 3),
  Analysis = rep(1:3, 3),
  `upper bound` = c(
    x %>% filter(bound == "upper") %>% select(z) %>% unlist() %>% as.numeric(),
    y %>% filter(Bound == "Upper") %>% select(Z) %>% unlist() %>% as.numeric(),
    z$upper$bound
  ),
  `upper prob` = c(
    x %>% filter(bound == "upper") %>% select(probability) %>% unlist() %>% as.numeric(),
    y %>% filter(Bound == "Upper") %>% select(Probability) %>% unlist() %>% as.numeric(),
    cumsum(z$upper$prob)
  )
) %>%
  arrange(Analysis) %>%
  group_by(Analysis) %>%
  gt()
```

# Test 4: Independent Tests

## Expect equal with mvtnorm for efficacy and futility bounds

```{r}
info <- c(40, 100)
r <- info[1] / info[2]

test <- gs_power_npe(
  theta = 0,
  info = info,
  info0 = NULL,
  binding = FALSE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, param = NULL, total_spend = 0.025),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfLDOF, param = NULL, total_spend = 0.02)
)

test1 <- test %>% filter(bound == "upper")
test2 <- test %>% filter(bound == "lower")

alpha.t <- 0.025
b.ia <- gsDesign::sfLDOF(alpha = alpha.t, t = r)
alpha.ia <- b.ia$spend

Pb <- function(alpha.t, alpha.ia, r, b) {
  temp <- mvtnorm::pmvnorm(
    lower = c(-Inf, b),
    upper = c(qnorm(1 - alpha.ia), Inf),
    corr = rbind(c(1, sqrt(r)), c(sqrt(r), 1))
  )
  return(alpha.t - alpha.ia - temp)
}

b <- uniroot(Pb, c(1.96, 4), alpha.t = alpha.t, alpha.ia = alpha.ia, r = r)

pb <- 1 - pnorm(b$root)

expect_equal(
  object = test1$z,
  expected = c(qnorm(1 - alpha.ia), b$root),
  tolerance = 0.001
)
expect_equal(
  object = test1$probability,
  expected = cumsum(c(b.ia$spend, pb)),
  tolerance = 0.001
)

beta.t <- 0.02
a.ia <- gsDesign::sfLDOF(alpha = beta.t, t = r)
beta.ia <- a.ia$spend

Pa <- function(beta.t, beta.ia, r, a) {
  temp <- mvtnorm::pmvnorm(
    lower = c(-Inf, qnorm(beta.ia)),
    upper = c(a, Inf),
    corr = rbind(c(1, sqrt(r)), c(sqrt(r), 1))
  )
  return(beta.t - beta.ia - temp)
}

a <- uniroot(Pa, c(-4, 1.96), beta.t = beta.t, beta.ia = beta.ia, r = r)

pa <- pnorm(a$root)

expect_equal(
  object = test2$z,
  expected = c(qnorm(beta.ia), a$root),
  tolerance = 0.001
)
expect_equal(
  object = test2$probability,
  expected = cumsum(c(a.ia$spend, pa)),
  tolerance = 0.001
)
```

## Expect equal with `gsDesign::gsProbability` outcome for efficacy bounds

```{r}
info <- c(40, 150, 200)

x <- gs_power_npe(
  theta = .1,
  info = info, binding = FALSE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, param = NULL, total_spend = 0.025),
  lower = gs_b,
  lpar = rep(-Inf, 3)
) %>% filter(bound == "upper")

y <- gs_power_npe(
  theta = .1,
  info = info, binding = FALSE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, param = NULL, total_spend = 0.025),
  lower = gs_b,
  lpar = rep(-Inf, 3)
) %>% filter(bound == "upper")

z <- gsDesign::gsProbability(
  k = 3, theta = .1,
  n.I = info,
  a = rep(-20, 3),
  b = gsDesign(k = 3, test.type = 1, sfu = sfLDOF, n.I = info)$upper$bound
)
```

```{r, echo=FALSE}
tibble(
  `Calculated from` = rep(c("new version", "old version", "gsDesign"), each = 3),
  Analysis = rep(1:3, 3),
  `upper bound` = c(
    x %>% filter(bound == "upper") %>% select(z) %>% unlist() %>% as.numeric(),
    y %>% filter(bound == "upper") %>% select(z) %>% unlist() %>% as.numeric(),
    z$upper$bound
  ),
  `upper prob` = c(
    x %>%
      filter(bound == "upper") %>%
      select(probability) %>%
      unlist() %>%
      as.numeric(),
    y %>%
      filter(bound == "upper") %>%
      select(probability) %>%
      unlist() %>%
      as.numeric(),
    cumsum(z$upper$prob)
  )
) %>%
  arrange(Analysis) %>%
  group_by(Analysis) %>%
  gt()
```

# Test 5: Compare with `gsDesign` under information-based design

Information-based design is useful when testing for a natural parameter $\delta$ (e.g., treatment difference on a relevant scale such as risk difference) where the variance of the estimate of $\delta$ is unknown. 
The basic *canonical form* of represents information-based design, so it is a particularly simple way to check corresponding basic calculations for sample size, bounds and power in `gs_power_npe()` and `gs_design_npe()`.

## Step 1: set the design assumptions

```{r}
k <- 2 # Number of analyses
test.type <- 4
alpha <- 0.025 # 1-sided Type I error
beta <- 0.15 # Type 2 error (1 - targeted power)
astar <- .1
timing <- 0.4 # Timing (information fraction) at interim analyses
sfu <- sfHSD # Efficacy bound spending function
sfupar <- -1 # Upper bound spending function parameters, if any
sfl <- sfLDPocock # Lower bound spending function, if used (test.type > 2)
sflpar <- 0 # Lower bound spending function parameters, if any
delta <- 0.1 # Natural parameter difference (assumed value - H0 value)
delta1 <- 0.1 # Natural parameter assumed value
delta0 <- 0 # Natural parameter difference under H0
endpoint <- "info"
n.fix <- 0
```

# References
