adj-seq-p.Rmd

---
title: "Adjusted sequential p-values"
author: "Yujie Zhao, Qi Liu, Linda Sun, Keaven Anderson"
output:
  rmarkdown::html_document:
    toc: true
    toc_float: true
    toc_depth: 2
    number_sections: true
    highlight: "textmate"
    css: "custom.css"
# bibliography: "example.bib"
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Adjusted sequential p-values}
---

```{r, message = FALSE, warning = FALSE}
library(wpgsd)
library(gsDesign)
library(tibble)
library(gt)
library(dplyr)
```

# Example overview

In a 2-arm controlled clinical trial example with one primary endpoint, there are 3 patient populations defined by the status of two biomarkers A and B:

- biomarker A positive,
- biomarker B positive,
- overall population.

The 3 primary elementary hypotheses are:

- $H_1$: the experimental treatment is superior to the control in the biomarker A positive population;
- \item $H_2$: the experimental treatment is superior to the control in the biomarker B positive population;
- \item $H_3$: the experimental treatment is superior to the control in the overall population.

Assume an interim analysis and a final analysis are planned for the study and the number of events are listed as

```{r}
event_tbl <- tribble(
  ~population, ~analysis, ~event,
  "A positive", 1, 100,
  "B positive", 1, 110,
  "AB positive", 1, 80,
  "overall", 1, 225,
  "A positive", 2, 200,
  "B positive", 2, 220,
  "AB positive", 2, 160,
  "overall", 2, 450,
)
```

The observed p-values are

```{r}
obs_tbl <- tribble(
  ~hypothesis, ~analysis, ~obs_p,
  "H1", 1, 0.02,
  "H2", 1, 0.01,
  "H3", 1, 0.012,
  "H1", 2, 0.015,
  "H2", 2, 0.012,
  "H3", 2, 0.010
) %>%
  mutate(obs_Z = -qnorm(obs_p))

obs_tbl %>%
  gt() %>%
  tab_header(title = "Nominal p-values")

p_obs_IA <- (obs_tbl %>% filter(analysis == 1))$obs_p
p_obs_FA <- (obs_tbl %>% filter(analysis == 2))$obs_p
```

The information fraction of $H_1, H_2, H_3$ at IA is

```{r}
IF_IA <- c(
  ((event_tbl %>% filter(analysis == 1, population == "A positive"))$event + (event_tbl %>% filter(analysis == 1, population == "overall"))$event) /
    ((event_tbl %>% filter(analysis == 2, population == "A positive"))$event + (event_tbl %>% filter(analysis == 2, population == "overall"))$event),
  ((event_tbl %>% filter(analysis == 1, population == "B positive"))$event + (event_tbl %>% filter(analysis == 1, population == "overall"))$event) /
    ((event_tbl %>% filter(analysis == 2, population == "B positive"))$event + (event_tbl %>% filter(analysis == 2, population == "overall"))$event),
  ((event_tbl %>% filter(analysis == 1, population == "AB positive"))$event + (event_tbl %>% filter(analysis == 1, population == "overall"))$event) /
    ((event_tbl %>% filter(analysis == 2, population == "AB positive"))$event + (event_tbl %>% filter(analysis == 2, population == "overall"))$event)
)

IF_IA
```

We assign the initial weights of $H_1, H_2, H_3$ as$\left(w_1(I), w_2(I), w_3(I) \right) = (0.3, 0.3, 0.4).$ And its multiplicity strategy is visualized in below. If $H_1$ is rejected, then $3/7$ local significance level $\alpha_1$ will be propagated to $H_2$, and $4/7$ will go to $H_3$. If $H_3$ is rejected, then half of $\alpha_3$ goes to $H_1$, and half goes to $H_2$.

```{r}
# Transition matrix in Figure A1
m <- matrix(c(
  0, 3 / 7, 4 / 7,
  3 / 7, 0, 4 / 7,
  0.5, 0.5, 0
), nrow = 3, byrow = TRUE)
# Initial weights
w <- c(0.3, 0.3, 0.4)
```

```{r}
name_hypotheses <- c("H1: Biomarker A positive", "H2: Biomarker B positive", "H3: Overall Population")

hplot <- gMCPLite::hGraph(
  3,
  alphaHypotheses = w, m = m,
  nameHypotheses = name_hypotheses, trhw = .2, trhh = .1,
  digits = 5, trdigits = 3, size = 5, halfWid = 1, halfHgt = 0.5,
  offset = 0.2, trprop = 0.4,
  fill = as.factor(c(2, 3, 1)),
  palette = c("#BDBDBD", "#E0E0E0", "#EEEEEE"),
  wchar = "w"
)
hplot
```

The correlation of the 6 statistisc (2 analysis $\times$ 3 hypothesis) are

```{r}
# Event count of intersection of paired hypotheses - Table 2
# H1, H2: Hypotheses intersected.
# (1, 1) represents counts for hypothesis 1
# (1, 2) for counts for the intersection of hypotheses 1 and 2
event <- tribble(
  ~H1, ~H2, ~Analysis, ~Event,
  1, 1, 1, event_tbl %>% filter(analysis == 1, population == "A positive") %>% select(event) %>% as.numeric(),
  2, 2, 1, event_tbl %>% filter(analysis == 1, population == "B positive") %>% select(event) %>% as.numeric(),
  3, 3, 1, event_tbl %>% filter(analysis == 1, population == "overall") %>% select(event) %>% as.numeric(),
  1, 2, 1, event_tbl %>% filter(analysis == 1, population == "AB positive") %>% select(event) %>% as.numeric(),
  1, 3, 1, event_tbl %>% filter(analysis == 1, population == "A positive") %>% select(event) %>% as.numeric(),
  2, 3, 1, event_tbl %>% filter(analysis == 1, population == "B positive") %>% select(event) %>% as.numeric(),
  1, 1, 2, event_tbl %>% filter(analysis == 2, population == "A positive") %>% select(event) %>% as.numeric(),
  2, 2, 2, event_tbl %>% filter(analysis == 2, population == "B positive") %>% select(event) %>% as.numeric(),
  3, 3, 2, event_tbl %>% filter(analysis == 2, population == "overall") %>% select(event) %>% as.numeric(),
  1, 2, 2, event_tbl %>% filter(analysis == 2, population == "AB positive") %>% select(event) %>% as.numeric(),
  1, 3, 2, event_tbl %>% filter(analysis == 2, population == "A positive") %>% select(event) %>% as.numeric(),
  2, 3, 2, event_tbl %>% filter(analysis == 2, population == "B positive") %>% select(event) %>% as.numeric()
)
event

# Generate correlation from events
gs_corr <- wpgsd::generate_corr(event)
gs_corr %>% round(2)
```

# Sequential p-value

## IA

```{r, wpgsdSeqIA, message=FALSE}
seq_p_IA_H123 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H1, H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H12 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H1, H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H13 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H1, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H23 <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H1 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H1",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H2 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H3 <- calc_seq_p(
  test_analysis = 1,
  test_hypothesis = "H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ),
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)
```

```{r, weightedBonferroniSeqIA, message=FALSE}
seq_p_IA_H123_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H1, H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)

seq_p_IA_H12_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H1, H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_IA_H13_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H1, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)

seq_p_IA_H23_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)

seq_p_IA_H1_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H1",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)

seq_p_IA_H2_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)

seq_p_IA_H3_B <- calc_seq_p(
  test_analysis = 1, # stage of interest
  test_hypothesis = "H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.3)
)
```

## FA

```{r, wpgsdSeqFA, message=FALSE}
seq_p_FA_H123 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.15)
)

seq_p_FA_H12 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.15)
)

seq_p_FA_H13 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.15)
)

seq_p_FA_H23 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.15)
)

seq_p_FA_H1 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H2 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H3 <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 2,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = gsDesign::sfHSD,
  spending_fun_par = -4,
  info_frac = c(min(IF_IA), 1),
  interval = c(1e-4, 0.2)
)
```

```{r, weightedBonferroniSeqFA, message=FALSE}
seq_p_FA_H123_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H12_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H13_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H23_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H2, H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H1_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H1",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H2_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H2",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)

seq_p_FA_H3_B <- calc_seq_p(
  test_analysis = 2, # stage of interest
  test_hypothesis = "H3",
  p_obs = tibble(
    analysis = 1:2,
    H1 = c(p_obs_IA[1], p_obs_FA[1]),
    H2 = c(p_obs_IA[2], p_obs_FA[2]),
    H3 = c(p_obs_IA[3], p_obs_FA[3])
  ), # observed p-value
  alpha_spending_type = 0,
  n_analysis = 2,
  initial_weight = w,
  transition_mat = m,
  z_corr = gs_corr,
  spending_fun = list(gsDesign::sfHSD, gsDesign::sfHSD, gsDesign::sfHSD),
  spending_fun_par = list(-4, -4, -4),
  info_frac = list(c(IF_IA[1], 1), c(IF_IA[2], 1), c(IF_IA[2], 1)),
  interval = c(1e-4, 0.2)
)
```

# Adjusted-Sequential p-value

## IA

```{r, wpgsdAdjSeqIA}
adj_seq_IA_H1 <- max(seq_p_IA_H123, seq_p_IA_H12, seq_p_IA_H13, seq_p_IA_H1)
adj_seq_IA_H2 <- max(seq_p_IA_H123, seq_p_IA_H12, seq_p_IA_H23, seq_p_IA_H2)
adj_seq_IA_H3 <- max(seq_p_IA_H123, seq_p_IA_H13, seq_p_IA_H23, seq_p_IA_H3)

cat("The adjusted-sequential p-value of H1, H2, H3 in IA via WPGSD is", adj_seq_IA_H1, adj_seq_IA_H2, adj_seq_IA_H3, "\n")
```

```{r, weightedBonferroniAdjSeqIA}
adj_seq_IA_H1_B <- max(seq_p_IA_H123_B, seq_p_IA_H12_B, seq_p_IA_H13_B, seq_p_IA_H1_B)
adj_seq_IA_H2_B <- max(seq_p_IA_H123_B, seq_p_IA_H12_B, seq_p_IA_H23_B, seq_p_IA_H2_B)
adj_seq_IA_H3_B <- max(seq_p_IA_H123_B, seq_p_IA_H13_B, seq_p_IA_H23_B, seq_p_IA_H3_B)

cat("The adjusted-sequential p-value of H1, H2, H3 in FA via weighted Bonferroni is", adj_seq_IA_H1_B, adj_seq_IA_H2_B, adj_seq_IA_H3_B, "\n")
```

## FA

### WPGSD
```{r, wpgsdAdjSeqFA}
adj_seq_FA_H1 <- max(seq_p_FA_H123, seq_p_FA_H12, seq_p_FA_H13, seq_p_FA_H1)
adj_seq_FA_H2 <- max(seq_p_FA_H123, seq_p_FA_H12, seq_p_FA_H23, seq_p_FA_H2)
adj_seq_FA_H3 <- max(seq_p_FA_H123, seq_p_FA_H13, seq_p_FA_H23, seq_p_FA_H3)

cat("The adjusted-sequential p-value of H1, H2, H3 in FA via WPGSD is", adj_seq_FA_H1, adj_seq_FA_H2, adj_seq_FA_H3, "\n")
```

```{r, weightedBonferroniAdjSeqFA}
adj_seq_FA_H1_B <- max(seq_p_FA_H123_B, seq_p_FA_H12_B, seq_p_FA_H13_B, seq_p_FA_H1_B)
adj_seq_FA_H2_B <- max(seq_p_FA_H123_B, seq_p_FA_H12_B, seq_p_FA_H23_B, seq_p_FA_H2_B)
adj_seq_FA_H3_B <- max(seq_p_FA_H123_B, seq_p_FA_H13_B, seq_p_FA_H23_B, seq_p_FA_H3_B)

cat("The adjusted-sequential p-value of H1, H2, H3 in FA via weighted Bonferroni is", adj_seq_FA_H1_B, adj_seq_FA_H2_B, adj_seq_FA_H3_B, "\n")
```

# Summary

```{r}
ans <- tribble(
  ~Analysis, ~Hypothesis, ~`Sequential p-values of WPGSD`, ~`Sequential p-values of Weighted Bonferroni`, ~`Adjusted-sequential p-values of WPGSD`, ~`Adjusted-sequential p-values of Weighted Bonferroni`,
  "IA", "H123", seq_p_IA_H123, seq_p_IA_H123_B, NA, NA,
  "IA", "H12", seq_p_IA_H12, seq_p_IA_H12_B, NA, NA,
  "IA", "H13", seq_p_IA_H13, seq_p_IA_H13_B, NA, NA,
  "IA", "H23", seq_p_IA_H23, seq_p_IA_H23_B, NA, NA,
  "IA", "H1", seq_p_IA_H1, seq_p_IA_H1_B, adj_seq_IA_H1, adj_seq_IA_H1_B,
  "IA", "H2", seq_p_IA_H2, seq_p_IA_H2_B, adj_seq_IA_H2, adj_seq_IA_H2_B,
  "IA", "H3", seq_p_IA_H3, seq_p_IA_H3_B, adj_seq_IA_H3, adj_seq_IA_H3_B,
  "FA", "H123", seq_p_FA_H123, seq_p_FA_H123_B, NA, NA,
  "FA", "H12", seq_p_FA_H12, seq_p_FA_H12_B, NA, NA,
  "FA", "H13", seq_p_FA_H13, seq_p_FA_H13_B, NA, NA,
  "FA", "H23", seq_p_FA_H23, seq_p_FA_H23_B, NA, NA,
  "FA", "H1", seq_p_FA_H1, seq_p_FA_H1_B, adj_seq_FA_H1, adj_seq_FA_H1_B,
  "FA", "H2", seq_p_FA_H2, seq_p_FA_H2_B, adj_seq_FA_H2, adj_seq_FA_H2_B,
  "FA", "H3", seq_p_FA_H3, seq_p_FA_H3_B, adj_seq_FA_H3, adj_seq_FA_H3_B
)

ans %>%
  select(
    Analysis, Hypothesis,
    `Sequential p-values of Weighted Bonferroni`, `Adjusted-sequential p-values of Weighted Bonferroni`,
    `Sequential p-values of WPGSD`, `Adjusted-sequential p-values of WPGSD`
  ) %>%
  gt() %>%
  tab_spanner(
    label = "Weighted Bonferroni",
    columns = c(`Sequential p-values of Weighted Bonferroni`, `Adjusted-sequential p-values of Weighted Bonferroni`)
  ) %>%
  tab_spanner(
    label = "WPGSD",
    columns = c(`Sequential p-values of WPGSD`, `Adjusted-sequential p-values of WPGSD`)
  ) %>%
  tab_style_body(
    columns = where(is.numeric),
    style = cell_fill(color = "pink"),
    fn = function(x) x <= 0.025
  ) %>%
  fmt_number(columns = 3:6, decimals = 4) %>%
  tab_header(
    title = "(Adjusted-) sequential p-values",
    subtitle = "Multiple populations"
  ) # %>% as_latex()
```