Comparison of two means

# app.R

# Minimal Shiny power / sample-size app (two-sample t-test), inspired by Lenth's interface style.

# Run: shiny::runApp()

library(shiny)

library(ggplot2)

# Helper: keep slider + numeric input in sync

sync_numeric_slider <- function(input, output, session, slider_id, numeric_id) {

    observeEvent(input[[slider_id]], {

        updateNumericInput(session, numeric_id, value = input[[slider_id]])

    }, ignoreInit = TRUE)

    observeEvent(input[[numeric_id]], {

        updateSliderInput(session, slider_id, value = input[[numeric_id]])

    }, ignoreInit = TRUE)

}

ui <- fluidPage(

    titlePanel("Power & Sample Size: Two-sample t-test"),

    sidebarLayout(

        sidebarPanel(

            helpText("Prospective planning tool for a two-sample t-test (difference in means)."),

            selectInput(

                "solve_for", "Solve for",

                choices = c("Power" = "power", "n (per group)" = "n"),

                selected = "power"

            ),

            tags$hr(),

            selectInput(

                "alternative", "Alternative",

                choices = c("Two-sided" = "two.sided", "One-sided (greater)" = "one.sided"),

                selected = "two.sided"

            ),

            sliderInput("alpha", "Significance level (alpha)", min = 0.001, max = 0.20, value = 0.05, step = 0.001),

            numericInput("alpha_n", NULL, value = 0.05, min = 0.0001, max = 0.5, step = 0.001),

            sliderInput("sd", "SD (common)", min = 0.1, max = 50, value = 10, step = 0.1),

            numericInput("sd_n", NULL, value = 10, min = 0.001, step = 0.1),

            sliderInput("delta", "Effect (delta = mean1 - mean2)", min = 0.01, max = 50, value = 5, step = 0.01),

            numericInput("delta_n", NULL, value = 5, min = 0, step = 0.01),

            sliderInput("n", "n per group", min = 2, max = 1000, value = 50, step = 1),

            numericInput("n_n", NULL, value = 50, min = 2, step = 1),

            sliderInput("power", "Target power", min = 0.05, max = 0.999, value = 0.80, step = 0.001),

            numericInput("power_n", NULL, value = 0.80, min = 0.01, max = 0.999, step = 0.001),

            tags$hr(),

            sliderInput("ratio", "Allocation ratio (n2 / n1)", min = 0.25, max = 4, value = 1, step = 0.01),

            helpText("If ratio = 1, groups are balanced. Otherwise, n2 = ratio * n1."),

            tags$hr(),

            selectInput(

                "xvar", "Plot: power vs",

                choices = c("n (per group)" = "n", "delta (effect)" = "delta"),

                selected = "n"

            ),

            sliderInput("grid_n_max", "Max n for plot", min = 20, max = 5000, value = 300, step = 10),

            checkboxInput("show_points", "Show points", value = FALSE)

        ),

        mainPanel(

            h3("Result"),

            verbatimTextOutput("result_text"),

            tags$hr(),

            h3("Graph"),

            plotOutput("power_plot", height = "420px"),

            tags$hr(),

            h3("Notes"),

            tags$ul(

                tags$li("Two-sample t-test uses a common SD and the noncentral t distribution via power.t.test()."),

                tags$li("For unequal allocation, the app approximates power by using an effective n per group based on harmonic mean."),

                tags$li("This is intended for prospective planning, not post-hoc 'observed power'.")

            )

        )

    )

)

server <- function(input, output, session) {

    # Sync sliders and numeric inputs

    sync_numeric_slider(input, output, session, "alpha", "alpha_n")

    sync_numeric_slider(input, output, session, "sd", "sd_n")

    sync_numeric_slider(input, output, session, "delta", "delta_n")

    sync_numeric_slider(input, output, session, "n", "n_n")

    sync_numeric_slider(input, output, session, "power", "power_n")

    # Use numeric versions as canonical

    alpha <- reactive(input$alpha_n)

    sd    <- reactive(input$sd_n)

    delta <- reactive(input$delta_n)

    n1    <- reactive(input$n_n)

    pow_t <- reactive(input$power_n)

    # Effective n for unequal allocation:

    # If n2 = ratio*n1, then harmonic mean:

    # n_eff = 2 / (1/n1 + 1/n2) = 2*n1*n2/(n1+n2)

    n_eff <- reactive({

        n2 <- input$ratio * n1()

        2 * n1() * n2 / (n1() + n2)

    })

    # Core calculation (solve for power or n)

    calc <- reactive({

        alt <- input$alternative

        sig <- alpha()

        if (input$solve_for == "power") {

            # compute power given n

            res <- power.t.test(

                n = n_eff(), delta = delta(), sd = sd(),

                sig.level = sig, type = "two.sample", alternative = alt

            )

            list(mode = "power", obj = res)

        } else {

            # solve for n (effective) to reach target power

            res <- power.t.test(

                n = NULL, delta = delta(), sd = sd(),

                sig.level = sig, power = pow_t(),

                type = "two.sample", alternative = alt

            )

            # res$n is per-group n in balanced design. We interpret it as effective n.

            list(mode = "n", obj = res)

        }

    })

    output$result_text <- renderText({

        cobj <- calc()$obj

        # Back out n1/n2 suggestion if ratio != 1 and solving for n

        if (input$solve_for == "n") {

            n_eff_needed <- cobj$n

            r <- input$ratio

            # Choose n1 such that harmonic mean equals n_eff_needed:

            # n_eff = 2*n1*(r*n1)/(n1+r*n1) = 2*r*n1^2 / ((1+r)*n1) = 2*r*n1/(1+r)

            # => n1 = n_eff * (1+r)/(2*r)

            n1_suggest <- n_eff_needed * (1 + r) / (2 * r)

            n2_suggest <- r * n1_suggest

            paste0(

                "Two-sample t-test\n",

                "Alternative: ", input$alternative, "\n",

                "alpha: ", format(alpha(), digits = 4), "\n",

                "SD: ", format(sd(), digits = 4), "\n",

                "delta: ", format(delta(), digits = 4), "\n\n",

                "Solved effective n per group (balanced-equivalent): ", format(n_eff_needed, digits = 5), "\n",

                "With ratio n2/n1 = ", format(r, digits = 4), ":\n",

                "  Suggested n1: ", ceiling(n1_suggest), "\n",

                "  Suggested n2: ", ceiling(n2_suggest), "\n",

                "  (n2 = ratio * n1)\n\n",

                "Achieved power (in the balanced-equivalent calculation): ", format(cobj$power, digits = 5), "\n"

            )

        } else {

            # solving for power

            n2_curr <- input$ratio * n1()

            paste0(

                "Two-sample t-test\n",

                "Alternative: ", input$alternative, "\n",

                "alpha: ", format(alpha(), digits = 4), "\n",

                "SD: ", format(sd(), digits = 4), "\n",

                "delta: ", format(delta(), digits = 4), "\n",

                "n1: ", n1(), "\n",

                "n2: ", ceiling(n2_curr), " (ratio = ", format(input$ratio, digits = 4), ")\n",

                "Effective n (harmonic mean): ", format(n_eff(), digits = 6), "\n\n",

                "Power: ", format(cobj$power, digits = 6), "\n"

            )

        }

    })

    output$power_plot <- renderPlot({

        alt <- input$alternative

        sig <- alpha()

        sdd <- sd()

        del <- delta()

        r   <- input$ratio

        if (input$xvar == "n") {

            nmax <- input$grid_n_max

            ngrid <- unique(round(seq(2, nmax, length.out = 120)))

            # compute effective n for each n1 on grid

            n2grid <- r * ngrid

            neffgrid <- 2 * ngrid * n2grid / (ngrid + n2grid)

            pwr <- vapply(neffgrid, function(ne) {

                power.t.test(

                    n = ne, delta = del, sd = sdd,

                    sig.level = sig, type = "two.sample", alternative = alt

                )$power

            }, numeric(1))

            df <- data.frame(n1 = ngrid, power = pwr)

            g <- ggplot(df, aes(x = n1, y = power)) +

                geom_line() +

                labs(

                    x = "n per group (group 1)",

                    y = "Power",

                    title = "Power vs n"

                ) +

                coord_cartesian(ylim = c(0, 1))

            if (isTRUE(input$show_points)) g <- g + geom_point()

            # Show current setting as vertical line

            g + geom_vline(xintercept = n1(), linetype = 2)

        } else {

            delgrid <- seq(0.01, max(0.5, del * 2), length.out = 140)

            # current n effective

            ne <- n_eff()

            pwr <- vapply(delgrid, function(dd) {

                power.t.test(

                    n = ne, delta = dd, sd = sdd,

                    sig.level = sig, type = "two.sample", alternative = alt

                )$power

            }, numeric(1))

            df <- data.frame(delta = delgrid, power = pwr)

            g <- ggplot(df, aes(x = delta, y = power)) +

                geom_line() +

                labs(

                    x = "Effect size (delta)",

                    y = "Power",

                    title = "Power vs effect size"

                ) +

                coord_cartesian(ylim = c(0, 1))

            if (isTRUE(input$show_points)) g <- g + geom_point()

            g + geom_vline(xintercept = del, linetype = 2)

        }

    })

}

shinyApp(ui, server)