Multi-Arm Sample Size Estimation

# Dependencies
library(MAMS) # MAMS Version 2.0.2
library(multiarm) # Multiarm 0.13.0, devtools::install_github("mjg211/multiarm")
library(tidyverse)
library(ggplot2)
library(gt)
library(ggpubr)

# Table settings
gt_theme <- function(data){
  data |>
    tab_options(table.width = pct(100),
                heading.align = "left",
                table.align = "left",
                table.font.size = px(16),
                column_labels.font.size = px(16),
                column_labels.font.weight = "bold",
                heading.title.font.size = px(16),
                heading.title.font.weight = "bold",
                table.border.top.style = "none",
                data_row.padding = px(2),
                column_labels.padding = px(2),
                heading.padding = px(5)
    ) |>
    cols_align(align = "left", columns = everything())
}

MAMS vs. Multiarm Packages

MAMS is primarily used for group-sequential, multiarm, multistage designs, but can be used to estimate fixed sample, single stage designs (current grant). As demonstrated below, when specified correctly, the two packages perform comparably for estimating the sample size for an outcome with a normal distribution. I prefer multiarm as more transparent and flexible for specifying multiple testing correction and power types. Additionally, the package offers functions for calculating Bernoulli and Poisson distributions (not described here). See: https://mjg211.github.io/multiarm/.

Multiple Testing Correction:

Bonferroni multiple testing correction, although more conservative, is typically preferred when outcomes tend to be highly correlated (as will certainly be the case among related psychometric scales).

Power Type:

The disjunctive power (or minimal power) is the probability of finding at least one true intervention effect across all of the outcomes.
The conjunctive power (or maximal power) is the probability of finding a true intervention effect on all outcomes.
The marginal (or individual) power is the probability of finding a true intervention effect on a particular outcome and is calculated separately for each outcome.

“When the clinical objective is to detect an intervention effect for at least one of the outcomes the disjunctive power and marginal power are recommended whereas the conjunctive power is recommended when the clinical objective is to detect an intervention effect on all the outcomes.” See: https://doi.org/10.1186/s12874-019-0754-4

Method	MTC	N	n_site	N_cor¹	n_site_cor¹	Alpha	Power	Delta1	Delta0	SD	PowerType
Table 1: Comparison of MAMS vs. multiarm Packages for Sample Size Estimation (Delta1=0.5 at 80% Power)
MAMS	Dunnett	276	28	386	39	0.05	0.8	0.5	0.1	1	Marginal
Multiarm	Dunnett	272	27	381	38	0.05	0.8	0.5	0.1	1	Marginal
Multiarm	Bonferroni	284	28	398	40	0.05	0.8	0.5	0.1	1	Marginal
¹ Corrected assuming 40% drop-out rate.

The more conservative Bonferroni correction only modestly increases the necessary sample size. Current study is recruiting for each arm from 10 different sites, as such the n_site calculated by the package(s) should be divided by 10 to derive the actual sample size per site (obviously, the total study sample size N remains the same).

Multiarm Power x N x Delta1 Estimates

Dunnett MTC

ma_est_dunnett <- function() {
  delta1 <- seq(0.2, 0.8, by = 0.1)
  beta <- seq(0.05, 0.95, by = 0.05)
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "dunnett", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Dunnett"
  )
  return(df)
}

t_ma_est_dunnett <- ma_est_dunnett()

t_ma_est_dunnett |> 
  filter(Power_pct %in% c(70, 80, 90)) |>
  filter(Delta1 %in% c(0.2, 0.5, 0.8)) |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 2: Multiarm Power x N x Delta1; MTC: Dunnett**")) |>   
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction
Table 2: Multiarm Power x N x Delta1; MTC: Dunnett
0.2	0.1	90	224	2240	314	3136	Marginal	Dunnett
0.2	0.2	80	169	1688	237	2364	Marginal	Dunnett
0.2	0.3	70	134	1340	188	1876	Marginal	Dunnett
0.5	0.1	90	36	360	51	504	Marginal	Dunnett
0.5	0.2	80	28	272	39	381	Marginal	Dunnett
0.5	0.3	70	22	216	31	303	Marginal	Dunnett
0.8	0.1	90	14	140	20	196	Marginal	Dunnett
0.8	0.2	80	11	108	16	152	Marginal	Dunnett
0.8	0.3	70	9	84	12	118	Marginal	Dunnett
¹ Corrected assuming 40% drop-out rate.

ggplot(t_ma_est_dunnett) +
  aes(x = Power_pct, y = N, colour = as.factor(Delta1), group = as.factor(Delta1)) +
  geom_line() +
  scale_color_hue(direction = 1) +
  scale_y_continuous(trans = "log10") +
  labs(
    x = "Power%",
    y = "log10(Estimated N_cor)",
    title = "Figure 1: Multi-Arm, Power x N_cor x Delta1",
    subtitle = "MTC: Dunnett",
    color = "Delta1") +
  theme_classic2() +
  theme(plot.title = element_text(face = "bold", size = 14L),
    plot.subtitle = element_text(size = 14L))

Bonferroni MTC

ma_est_bon <- function() {
  delta1 <- seq(0.2, 0.8, by = 0.1)
  beta <- seq(0.05, 0.95, by = 0.05)
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Bonferroni"
  )
  return(df)
}

t_ma_est_bon <- ma_est_bon()

t_ma_est_bon |> 
  filter(Power_pct %in% c(70, 80, 90)) |>
  filter(Delta1 %in% c(0.2, 0.5, 0.8)) |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 3: Multiarm Power x N x Delta1; MTC: Bonferroni**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction
Table 3: Multiarm Power x N x Delta1; MTC: Bonferroni
0.2	0.1	90	233	2328	326	3260	Marginal	Bonferroni
0.2	0.2	80	177	1764	247	2470	Marginal	Bonferroni
0.2	0.3	70	141	1408	198	1972	Marginal	Bonferroni
0.5	0.1	90	38	376	53	527	Marginal	Bonferroni
0.5	0.2	80	29	284	40	398	Marginal	Bonferroni
0.5	0.3	70	23	228	32	320	Marginal	Bonferroni
0.8	0.1	90	15	148	21	208	Marginal	Bonferroni
0.8	0.2	80	12	112	16	157	Marginal	Bonferroni
0.8	0.3	70	9	88	13	124	Marginal	Bonferroni
¹ Corrected assuming 40% drop-out rate.

ggplot(t_ma_est_bon) +
  aes(x = Power_pct, y = N_cor, colour = as.factor(Delta1), group = as.factor(Delta1)) +
  geom_line() +
  scale_color_hue(direction = 1) +
  scale_y_continuous(trans = "log10") +
  labs(
    x = "Power%",
    y = "log10(Estimated N_cor)",
    title = "Figure 2 : Multi-Arm, Power x N_cor x Delta1",
    subtitle = "MTC: Bonferroni",
    color = "Delta1"
  ) +
  theme_classic2() + 
    theme(plot.title = element_text(face = "bold", size = 14L),
    plot.subtitle = element_text(size = 14L))

Example Estimation: All Power Types - Bonferroni NTC

Supposing the range of treatment effects anticipated is 0.3, 0.4, and 0.5. Calculating the estimated sample size range needed overall and by site for 70% power.

Marginal Power Type

ma_est_bon_ex <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N/10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_ex <- ma_est_bon_ex()


t_ma_est_bon_ex |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 4: Example Estimation - Marginal Power Type**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 4: Example Estimation - Marginal Power Type
0.2	0.2	80	177	1764	247	2470	Marginal	Bonferroni	0.017
0.3	0.2	80	79	784	110	1098	Marginal	Bonferroni	0.017
0.4	0.2	80	45	444	63	622	Marginal	Bonferroni	0.017
0.5	0.2	80	29	284	40	398	Marginal	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.

K=3 experimental treatments will be included in the trial. A significance level of alpha=0.05 will be used, in combination with Bonferroni’s correction. The standard deviations of the responses will be assumed to be 1 in each arm. The minimum marginal power will be controlled to level 1−beta=0.8 under each of their respective least favorable configurations. The target allocation to each of the experimental arms will be the same as the control arm with the realized allocation ratios to the experimental arms at 1:1:1. Accounting for a 40% drop-out rate with an anticipating treatment effect range of 0.3-05 detectable with 80% power, the estimated sample size needed is 398-1098 overall and 10-28 participants needed per site. A Bonferroni p-value threshold 0.017 will be used to determine significance among arm comparisons.

Disjunctive Power Type

ma_est_bon_ex_dis <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "disjunctive", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Disjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_ex_dis <- ma_est_bon_ex_dis()


t_ma_est_bon_ex_dis |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 5: Example Estimation - Disjunctive Power Type**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 5: Example Estimation - Disjunctive Power Type
0.2	0.2	80	104	1036	146	1451	Disjunctive	Bonferroni	0.017
0.3	0.2	80	46	460	65	644	Disjunctive	Bonferroni	0.017
0.4	0.2	80	26	260	37	364	Disjunctive	Bonferroni	0.017
0.5	0.2	80	17	168	24	236	Disjunctive	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.

Conjunctive Power Type

ma_est_bon_ex_con <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "conjunctive", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Conjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_ex_con <- ma_est_bon_ex_con()


t_ma_est_bon_ex_con |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 6: Example Estimation - Conjunctive Power Type**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 6: Example Estimation - Conjunctive Power Type
0.2	0.2	80	241	2404	337	3366	Conjunctive	Bonferroni	0.017
0.3	0.2	80	108	1072	151	1501	Conjunctive	Bonferroni	0.017
0.4	0.2	80	61	604	85	846	Conjunctive	Bonferroni	0.017
0.5	0.2	80	39	388	55	544	Conjunctive	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.

combined <- rbind(t_ma_est_bon_ex, t_ma_est_bon_ex_dis, t_ma_est_bon_ex_con) 

combined |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 7: Example Estimation - All Power Types**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 7: Example Estimation - All Power Types
0.2	0.2	80	177	1764	247	2470	Marginal	Bonferroni	0.017
0.3	0.2	80	79	784	110	1098	Marginal	Bonferroni	0.017
0.4	0.2	80	45	444	63	622	Marginal	Bonferroni	0.017
0.5	0.2	80	29	284	40	398	Marginal	Bonferroni	0.017
0.2	0.2	80	104	1036	146	1451	Disjunctive	Bonferroni	0.017
0.3	0.2	80	46	460	65	644	Disjunctive	Bonferroni	0.017
0.4	0.2	80	26	260	37	364	Disjunctive	Bonferroni	0.017
0.5	0.2	80	17	168	24	236	Disjunctive	Bonferroni	0.017
0.2	0.2	80	241	2404	337	3366	Conjunctive	Bonferroni	0.017
0.3	0.2	80	108	1072	151	1501	Conjunctive	Bonferroni	0.017
0.4	0.2	80	61	604	85	846	Conjunctive	Bonferroni	0.017
0.5	0.2	80	39	388	55	544	Conjunctive	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.

Example Estimation: All Power Types - Dunnett NTC

ma_est_dun_ex <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "dunnett", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N/10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Dunnett",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_dun_ex <- ma_est_dun_ex()

ma_est_dun_ex_dis <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "dunnett", # multiple testing correction
        power  = "disjunctive", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Disjunctive",
    Correction = "Dunnett",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_dun_ex_dis <- ma_est_dun_ex_dis()

ma_est_dun_ex_con <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_norm(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        sigma = rep(1, 4), # standard deviations of the responses in each arm
        ratio = rep(1, 3), # allocation ratio
        correction = "dunnett", # multiple testing correction
        power  = "conjunctive", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Conjunctive",
    Correction = "Dunnett",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_dun_ex_con <- ma_est_dun_ex_con()

combined_dun <- rbind(t_ma_est_dun_ex, t_ma_est_dun_ex_dis, t_ma_est_dun_ex_con) 

combined_dun |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 8: Example Estimation - All Power Types**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Dunnett critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 8: Example Estimation - All Power Types
0.2	0.2	80	169	1688	237	2364	Marginal	Dunnett	0.02
0.3	0.2	80	76	752	106	1053	Marginal	Dunnett	0.02
0.4	0.2	80	43	424	60	594	Marginal	Dunnett	0.02
0.5	0.2	80	28	272	39	381	Marginal	Dunnett	0.02
0.2	0.2	80	98	976	137	1367	Disjunctive	Dunnett	0.02
0.3	0.2	80	44	436	62	611	Disjunctive	Dunnett	0.02
0.4	0.2	80	25	244	35	342	Disjunctive	Dunnett	0.02
0.5	0.2	80	16	160	23	224	Disjunctive	Dunnett	0.02
0.2	0.2	80	232	2316	325	3243	Conjunctive	Dunnett	0.02
0.3	0.2	80	104	1032	145	1445	Conjunctive	Dunnett	0.02
0.4	0.2	80	58	580	82	812	Conjunctive	Dunnett	0.02
0.5	0.2	80	38	372	53	521	Conjunctive	Dunnett	0.02
¹ Corrected assuming 40% drop-out rate.
² Dunnett critical p-value for significance.

Other distributions

Bernoilli

ma_est_bon_bern <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_bern(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        pi0=0.3, # control arm response rate
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N/10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_bern <- ma_est_bon_bern()

ma_est_bon_bern_dis <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
       d_multi <- des_ss_bern(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        pi0=0.3, # control arm response rate
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Disjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_bern_dis <- ma_est_bon_bern_dis()

ma_est_bon_bern_con <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_bern(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        pi0=0.3, # control arm response rate
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Conjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_bon_bern_con <- ma_est_bon_bern_con()

combined_bon_bern <- rbind(t_ma_est_bon_bern, t_ma_est_bon_bern_dis, t_ma_est_bon_bern_con) 

combined_bon_bern |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 9: Example Estimation - Bernoulli**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 9: Example Estimation - Bernoulli
0.2	0.2	80	41	408	58	572	Marginal	Bonferroni	0.017
0.3	0.2	80	18	180	26	252	Marginal	Bonferroni	0.017
0.4	0.2	80	10	96	14	135	Marginal	Bonferroni	0.017
0.5	0.2	80	6	56	8	79	Marginal	Bonferroni	0.017
0.2	0.2	80	41	408	58	572	Disjunctive	Bonferroni	0.017
0.3	0.2	80	18	180	26	252	Disjunctive	Bonferroni	0.017
0.4	0.2	80	10	96	14	135	Disjunctive	Bonferroni	0.017
0.5	0.2	80	6	56	8	79	Disjunctive	Bonferroni	0.017
0.2	0.2	80	41	408	58	572	Conjunctive	Bonferroni	0.017
0.3	0.2	80	18	180	26	252	Conjunctive	Bonferroni	0.017
0.4	0.2	80	10	96	14	135	Conjunctive	Bonferroni	0.017
0.5	0.2	80	6	56	8	79	Conjunctive	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.

Poisson

ma_est_poi_bern <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_pois(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        lambda0 = 1, # event rate in the control arm
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N/10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Marginal",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_poi_bern <- ma_est_poi_bern()

ma_est_poi_bern_dis <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_pois(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        lambda0 = 1, # event rate in the control arm
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Disjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_poi_bern_dis <- ma_est_poi_bern_dis()

ma_est_poi_bern_con <- function() {
  delta1 <- seq(0.2, 0.5, by = 0.1)
  beta <- 0.2
  N <- rep(NA, length(beta) * length(delta1))
  beta_res <- rep(NA, length(beta) * length(delta1))
  delta1_res <- rep(NA, length(beta) * length(delta1))
  mtcpv <- rep(NA, length(beta) * length(delta1))

  i <- 1
  for (d in delta1) {
    for (b in beta) {
      d_multi <- des_ss_pois(
        K = 3, # no. experimental treatments
        alpha = 0.05, # significance level
        beta = b, # beta, power = 1-beta
        lambda0 = 1, # event rate in the control arm
        delta1 = d, # interesting treatment effect on the traditional scale
        delta0 = 0.1, # uninteresting treatment effect on the traditional scale
        ratio = rep(1, 3), # allocation ratio
        correction = "bonferroni", # multiple testing correction
        power  = "marginal", # type of power to control
        integer = TRUE
      ) # return integer sample size
      N[i] <- d_multi$N
      beta_res[i] <- b
      delta1_res[i] <- d
      mtcpv[i] <- d_multi$opchar$P1[1]
      i <- i + 1
    }
  }
  df <- data.frame(
    Delta1 = delta1_res,
    Beta = beta_res,
    Power_pct = (1 - beta_res) * 100,
    n_site = ceiling(N / 10),
    N = ceiling(N),
    n_site_cor = ceiling((N * 1.4)/ 10),
    N_cor = ceiling(N * 1.4),
    PowerType = "Conjunctive",
    Correction = "Bonferroni",
    MTC_pv = round(mtcpv, digits = 3)
  )
  return(df)
}

t_ma_est_poi_bern_con <- ma_est_poi_bern_con()

combined_poi_bern <- rbind(t_ma_est_poi_bern, t_ma_est_poi_bern_dis, t_ma_est_poi_bern_con) 

combined_poi_bern |>
  gt() |>
  gt_theme() |>
  tab_header(title = md("**Table 10: Example Estimation - Poisson**")) |>
  tab_footnote(footnote = "Corrected assuming 40% drop-out rate.",
               locations = cells_column_labels(columns = c(n_site_cor, N_cor))) |>
  tab_footnote(footnote = "Bonferroni critical p-value for significance.",
               locations = cells_column_labels(columns = c(MTC_pv)))

Delta1	Beta	Power_pct	n_site	N	n_site_cor¹	N_cor¹	PowerType	Correction	MTC_pv²
Table 10: Example Estimation - Poisson
0.2	0.2	80	195	1944	273	2722	Marginal	Bonferroni	0.017
0.3	0.2	80	91	904	127	1266	Marginal	Bonferroni	0.017
0.4	0.2	80	54	532	75	745	Marginal	Bonferroni	0.017
0.5	0.2	80	36	356	50	499	Marginal	Bonferroni	0.017
0.2	0.2	80	195	1944	273	2722	Disjunctive	Bonferroni	0.017
0.3	0.2	80	91	904	127	1266	Disjunctive	Bonferroni	0.017
0.4	0.2	80	54	532	75	745	Disjunctive	Bonferroni	0.017
0.5	0.2	80	36	356	50	499	Disjunctive	Bonferroni	0.017
0.2	0.2	80	195	1944	273	2722	Conjunctive	Bonferroni	0.017
0.3	0.2	80	91	904	127	1266	Conjunctive	Bonferroni	0.017
0.4	0.2	80	54	532	75	745	Conjunctive	Bonferroni	0.017
0.5	0.2	80	36	356	50	499	Conjunctive	Bonferroni	0.017
¹ Corrected assuming 40% drop-out rate.
² Bonferroni critical p-value for significance.