## -----------------------------------------------------------------------------
#| echo: false
library(rpact)
library(dplyr)
library(tidyr)
library(magrittr)
# packageVersion("rpact")


## -----------------------------------------------------------------------------
#| echo: true
# Load the rpact package
library(rpact)

# Define the design parameters for sample size calculation
design <- getDesignGroupSequential(
    kMax = 1, # Only one analysis (classic fixed design)
    alpha = 0.025, # Significance level
    beta = 0.20, # 80% power
    sided = 1 # One-sided test
)

# Estimate sample size
sampleSizeResult <- design |>
    getSampleSizeMeans(
        groups = 2, # Two groups: Treatment vs. Placebo
        alternative = 12, # Expected effect size
        stDev = 15 # Common standard deviation
    )


## -----------------------------------------------------------------------------
#| echo: true
# Print the summary of the results
sampleSizeResult |> summary()


## -----------------------------------------------------------------------------
#| echo: true
# Fetch the number of subjects required
sampleSizeResult |> fetch("Number of subjects fixed")


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
# Define parameters based on initial assumptions
calculatedSampleSize <- ceiling(sampleSizeResult$numberOfSubjects)

powerResult <- design |>
    getPowerMeans(
        groups = 2, # Two groups: Treatment and Placebo
        alternative = 12, # Expected effect size
        stDev = 15, # Common standard deviation
        maxNumberOfSubjects = calculatedSampleSize
    )


## -----------------------------------------------------------------------------
# Print the summary of the the results
powerResult |> summary()


## -----------------------------------------------------------------------------
#| echo: true
# Define scenarios for adjustments
scenarios <- list(
    list(alternative = 10, stDev = 15), # 1: Reduced effect
    list(alternative = 10, stDev = 14), # 2: Reduced effect + sd
    list(alternative = 11, stDev = 15), # 3: Reduced effect
    list(alternative = 11, stDev = 14), # 4: Reduced effect + sd
    list(alternative = 12, stDev = 15), # 5: Base scenario
    list(alternative = 13, stDev = 16), # 6: Incr. sd + effect
    list(alternative = 13, stDev = 17), # 7: Incr. sd + effect
    list(alternative = 12, stDev = 16), # 8: Increased sd
    list(alternative = 12, stDev = 17) # 9: Increased sd
)

# Run calculations for each scenario
results <- scenarios |>
    lapply(function(scenario) {
        getPowerMeans(
            design = design,
            groups = 2,
            alternative = scenario$alternative,
            stDev = scenario$stDev,
            maxNumberOfSubjects = 52
        )
    })

# Fetch only the power from the result objects
x <- sapply(results, function(result) {
    result |> fetch("Overall reject")
})


## -----------------------------------------------------------------------------
#| echo: true
# Display results for each scenario
scenarios |>
    bind_rows() |>
    mutate(power = x) |>
    as.data.frame() |>
    set_rownames(paste("Scenario", 1:length(x))) |>
    kable()


## -----------------------------------------------------------------------------
#| echo: true
# Define the initial design
design <- getDesignGroupSequential(
    kMax = 1,
    alpha = 0.025,
    beta = 0.2,
    sided = 1
)

# Define scenarios for different allocation ratios
scenarios <- list(
    list(allocationRatio = 1), # 1:1 allocation
    list(allocationRatio = 2), # 2:1 allocation
    list(allocationRatio = 3) # 3:1 allocation
)

# Run calculations for each scenario with specified allocation ratios
results <- scenarios |>
    lapply(function(scenario) {
        getPowerMeans(
            design = design,
            groups = 2,
            alternative = 12,
            stDev = 15,
            maxNumberOfSubjects = 50,
            allocationRatioPlanned = scenario$allocationRatio
        )
    })

# Fetch only the power from the result objects
x <- sapply(results, function(result) {
    result |> fetch("Overall reject")
})


## -----------------------------------------------------------------------------
#| echo: true
# Display results for each scenario
scenarios |>
    bind_rows() |>
    mutate(power = x) |>
    as.data.frame() |>
    set_rownames(paste("Scenario", 1:length(x))) |>
    kable()


## -----------------------------------------------------------------------------
#| echo: true
numberOfSubjectsFixed <- getDesignGroupSequential(
    kMax = 1,
    alpha = 0.025,
    beta = 0.2,
    sided = 1
) |>
    getSampleSizeMeans(
        groups = 2,
        alternative = 12,
        stDev = 15,
        allocationRatioPlanned = 2
    ) |>
    fetch("Number of subjects fixed") |>
    ceiling()
numberOfSubjectsFixed


## -----------------------------------------------------------------------------
#| echo: true
# Define scenarios
scenarios <- list(
    list(informationRate = 0.5),
    list(informationRate = 0.6),
    list(informationRate = 0.7),
    list(informationRate = 0.8),
    list(informationRate = 0.9)
)


## -----------------------------------------------------------------------------
#| echo: true
# Run calculations for each scenario
results <- scenarios |>
    lapply(function(scenario) {
        getDesignGroupSequential(
            informationRates = c(scenario$informationRate, 1),
            alpha = 0.025,
            sided = 1,
            typeOfDesign = "OF"
        ) |>
            getPowerMeans(
                groups = 2,
                alternative = 12,
                stDev = 15,
                maxNumberOfSubjects = numberOfSubjectsFixed, # 58
                allocationRatioPlanned = 2
            )
    })


## -----------------------------------------------------------------------------
#| echo: true
# Display results for each scenario
x <- lapply(results, function(result) {
    result |> fetch("Expected number of subjects", "Early stop")
}) |>
    bind_rows()

#| echo: true
scenarios |>
    bind_rows() |>
    cbind(x) |>
    as.data.frame() |>
    set_rownames(paste("Scenario", 1:length(scenarios))) |>
    kable()


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
design <- getDesignGroupSequential(
    informationRates = c(0.7, 1),
    alpha = 0.025, # Overall significance level
    beta = 0.2, # Power 80%
    sided = 1, # One-sided test
    typeOfDesign = "OF" # O'Brien & Fleming design
)

sampleSizeResult <- getSampleSizeMeans(
    design = design,
    groups = 2,
    alternative = 12, # Target effect
    stDev = 15, # Common standard deviation
    allocationRatioPlanned = 2 # 2:1 allocation
)

# Print the summary of the results
sampleSizeResult |>
    summary()


## -----------------------------------------------------------------------------
#| echo: true
sampleSizeResult |>
    fetch(
        "Expected number of subjects under H1",
        "Early stop"
    )


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
powerResult <- getPowerMeans(
    design = design,
    groups = 2,
    alternative = 12, # Expected effect size
    stDev = 15, # Common standard deviation
    # Sample size per stage from previous calculation
    maxNumberOfSubjects = ceiling(sampleSizeResult$numberOfSubjects)[2, 1],
    allocationRatioPlanned = 2 # Allocation ratio 2:1
)

# Print the summary of the results
powerResult |>
    summary()


## -----------------------------------------------------------------------------
#| echo: true

# Define scenarios
scenarios <- list(
    list(futilityBounds = -1),
    list(futilityBounds = -0.5),
    list(futilityBounds = 0),
    list(futilityBounds = 0.5),
    list(futilityBounds = 1),
    list(futilityBounds = 1.5)
)


## -----------------------------------------------------------------------------
#| echo: true
# Run calculations for each scenario
results <- scenarios |>
    lapply(function(scenario) {
        getDesignGroupSequential(
            informationRates = c(0.7, 1),
            futilityBounds = scenario$futilityBounds,
            alpha = 0.025,
            sided = 1,
            typeOfDesign = "OF"
        ) |>
            getSampleSizeMeans(
                groups = 2,
                alternative = 12,
                stDev = 15,
                allocationRatioPlanned = 2
            )
    })


## -----------------------------------------------------------------------------
#| echo: true
# Display results for each scenario
x <- sapply(results, function(result) {
    result |> fetch("Futility bounds (treatment effect scale)")
})

#| echo: true
scenarios |>
    bind_rows() |>
    mutate("Futility bounds (treatment effect scale)" = x) |>
    as.data.frame() |>
    set_rownames(paste("Scenario", 1:length(x))) |>
    kable()


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
?uniroot


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
uniroot(
    f,
    interval,
    ...,
    lower = min(interval),
    upper = max(interval),
    f.lower = f(lower, ...),
    f.upper = f(upper, ...),
    extendInt = c("no", "yes", "downX", "upX"),
    check.conv = FALSE,
    tol = .Machine$double.eps^0.25,
    maxiter = 1000,
    trace = 0
)


## -----------------------------------------------------------------------------
#| echo: true
soughtBoundaryTreatmentEffectScale <- 2
futilityBound <- uniroot(
    function(x) {
        soughtBoundaryTreatmentEffectScale -
            getDesignGroupSequential(
                informationRates = c(0.7, 1),
                futilityBounds = x,
                alpha = 0.025,
                sided = 1,
                typeOfDesign = "OF"
            ) |>
            getSampleSizeMeans(
                groups = 2,
                alternative = 12,
                stDev = 15,
                allocationRatioPlanned = 2
            ) |>
            fetch("Futility bounds (treatment effect scale)") |>
            as.numeric()
    },
    lower = 0,
    upper = 2
)$root
futilityBound


## -----------------------------------------------------------------------------
#| echo: true

n1 <- sampleSizeResult$numberOfSubjects[2, 1]
n11 <- n1 * 2 / 3
n21 <- n1 * 1 / 3

futilityBoundClosed <- soughtBoundaryTreatmentEffectScale * sqrt(0.7 * n11 * n21 / (n11 + n21)) / 15
futilityBoundClosed


## -----------------------------------------------------------------------------
#| echo: true
sampleSizeResult <- getDesignGroupSequential(
    informationRates = c(0.7, 1),
    futilityBounds = futilityBoundClosed,
    alpha = 0.025,
    sided = 1,
    typeOfDesign = "OF"
) |>
    getSampleSizeMeans(
        groups = 2,
        alternative = 12,
        stDev = 15,
        allocationRatioPlanned = 2
    )


## -----------------------------------------------------------------------------
#| echo: true
sampleSizeResult |> summary()


## -----------------------------------------------------------------------------
#| echo: true
# nMax <- ceiling(sampleSizeResult$numberOfSubjects)[2,1]
getDesignGroupSequential(
    informationRates = c(0.7, 1),
    futilityBounds = futilityBound,
    alpha = 0.025,
    sided = 1,
    typeOfDesign = "OF"
) |>
    plot(type = 1, nMax = 60)


## -----------------------------------------------------------------------------
#| echo: true
#| output-location: slide
getDesignGroupSequential(
    informationRates = c(0.7, 1),
    futilityBounds = futilityBound,
    alpha = 0.025,
    sided = 1,
    typeOfDesign = "OF"
) |>
    getSampleSizeMeans(
        groups = 2,
        alternative = 12,
        stDev = 15,
        allocationRatioPlanned = 2
    ) |>
    plot(type = 1)


## -----------------------------------------------------------------------------
#| echo: true
simResults <- getDesignGroupSequential(
    informationRates = c(0.7, 1),
    futilityBounds = futilityBound,
    alpha = 0.025,
    sided = 1,
    typeOfDesign = "OF"
) |> 
    getSimulationMeans(
        groups = 2,
        alternative = c(0, 12), # H0 and H1
        stDev = 15,
        allocationRatioPlanned = 2,
        plannedSubjects = c(41, 59),
        maxNumberOfIterations = 1000,
        seed = 12345
    )


## -----------------------------------------------------------------------------
#| echo: true
simResults |> summary()


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
trialDataStage1 <- read.csv(here::here("data/trial_data_stage1.csv"))
trialDataStage1$group <- factor(trialDataStage1$group)
trialDataStage1 |>
    as_tibble()


## -----------------------------------------------------------------------------
#| echo: true
trialData <- trialDataStage1 |>
    mutate(bloodPressureDiff = bloodPressure - bloodPressureBaseline)
trialData |>
    as_tibble()


## -----------------------------------------------------------------------------
#| echo: true
dataSummary <- trialData |>
    group_by(group) |>
    summarise(
        n = n(),
        meanBaseline = mean(bloodPressureBaseline),
        stDevBaseline = sd(bloodPressureBaseline),
        meanTherapy = mean(bloodPressure),
        stDevTherapy = sd(bloodPressure),
        mean = mean(bloodPressureDiff),
        stDev = sd(bloodPressureDiff)
    )
dataSummary |>
    kable()


## -----------------------------------------------------------------------------
#| echo: true
library(ggplot2)
ggplot(
    aes(x = group, y = bloodPressureDiff, fill = group),
    data = trialData
) +
    geom_boxplot(na.rm = TRUE)


## -----------------------------------------------------------------------------
#| echo: true
library(tidyr)
datasetStage1 <- dataSummary |>
    select(group, n, mean, stDev) |>
    mutate(stage = rep(1, nrow(dataSummary))) |>
    pivot_wider(
        names_from = "group",
        values_from = c("n", "mean", "stDev"),
        names_sep = ""
    ) |>
    getDataset()


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
getAnalysisResults(
    design = design,
    dataInput = datasetStage1,
    directionUpper = FALSE,
    nPlanned = 18,
    allocationRatioPlanned = 2,
    stDevH1 = 15
) |>
    summary()


## -----------------------------------------------------------------------------
#| echo: true
#| eval: true
trialData <- read.csv(here::here("data/trial_data.csv"))
trialData$group <- factor(trialData$group)
trialData |>
    as_tibble()


## -----------------------------------------------------------------------------
#| echo: true
library(emmeans)
trialData$group <- relevel(trialData$group, ref = "2")
trialData$bloodPressureDiff <-
    trialData$bloodPressure - trialData$bloodPressureBaseline
dataset <- getDataset(
    lm(bloodPressureDiff ~ group, data = trialData, subset = (stage == 1)) |>
        emmeans("group"),
    lm(bloodPressureDiff ~ group, data = trialData, subset = (stage == 2)) |>
        emmeans("group")
)


## -----------------------------------------------------------------------------
#| echo: true
dataset |> summary()


## -----------------------------------------------------------------------------
#| echo: true
getAnalysisResults(
    design = design,
    dataInput = dataset,
    directionUpper = FALSE
) |>
    summary()