#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 800
library(shiny)
library(ggplot2)
library(bslib)
library(bsicons)
ui <- page_sidebar(
theme = bs_theme(bootswatch = "materia"),
titlePanel("T-test Visualization"),
sidebar = sidebar(
numericInput("n", "Sample size per group (n):", value = 735, min = 2),
numericInput("mu0", HTML("Null hypothesized mean difference (μ<sub>0</sub>):"), value = 0),
numericInput("mu_true", HTML("True mean difference (μ<sub>true</sub>):"), value = 0), # True mean input
numericInput("sigma", "Population standard deviation (σ):", value = 2.15),
numericInput("alpha", "Significance level (α):", value = 0.05),
radioButtons("test_type", "Test Type:", choices = list("One-sided (left)" = "one left","One-sided (right)" = "one right", "Two-sided" = "two"),
selected = "two"),
actionButton("toggleAltPlot", "Show/Hide Alternative Distribution"), # Button to toggle alternative plot
actionButton("toggleCI", "Show/Hide C.I. & Rejection Region") # Button to toggle CI and rejection region
)
, actionButton("draw", "Draw Sample"),
card(
card_header("Sample"),
height = "850px",
fluidRow(
column(6,
style = "justify-content: center; align-items: center", # Center vertically and control height
uiOutput("sample_stats") # Summary stats on the left
),
column(6,
plotOutput("samplePlot", height = "175px") # Sample plot on the right, fixed height
)
)
),
card(card_header("Null Distributions"),
height = "700px",
fluidRow(
column(6,
plotOutput("meanDensityPlot", height = "175px") ), # New density plot for sample mean
column(6,
plotOutput("densityPlot", height = "175px")))), # Larger density plot below
card(card_header("Power"),
height = "350px",
card_body(
textOutput("power")))
)
server <- function(input, output, session) {
# Reactive values to track the sample, alternative plot visibility, and plot data
sample_data <- reactiveVal(NULL)
plot_data <- reactiveValues(
t_stat = NULL,
sample_mean = NULL,
sample_sd = NULL,
x_limits_null = NULL,
mean_x_limits = NULL
)
showAlt <- reactiveVal(FALSE)
observeEvent(input$toggleAltPlot, {
showAlt(!showAlt())
})
showCI <- reactiveVal(FALSE)
observeEvent(input$toggleCI, {
showCI(!showCI()) # Toggle CI and rejection region visibility
})
observeEvent(input$draw, {
# Create a bimodal sample using the true mean (µ_true)
pop1a <- rnorm(10000, mean = - 1.75, sd = (4 * input$sigma - 3.5) / 4) # First mode
pop1b <- rnorm(10000, mean = 1.75, sd = (4 * input$sigma - 3.5) / 4) # Second mode
pop2a <- rnorm(10000, mean = input$mu_true - 1.75, sd = (4 * input$sigma - 3.5) / 4) # First mode
pop2b <- rnorm(10000, mean = input$mu_true + 1.75, sd = (4 * input$sigma - 3.5) / 4) # Second mode
pop1 = c(pop1a, pop1b)
pop2 = c(pop2a, pop2b)
sample1 = sample(pop1, input$n, replace = T)
sample2 = sample(pop2, input$n, replace = T)
sample_data(data.frame(group = rep(c("Group 1", "Group 2"), each = input$n),
x = c(sample1, sample2))) # Store sample data
# Store plot-specific values in reactiveValues
plot_data$sample_mean <- mean(sample2) - mean(sample1)
plot_data$sample_sd <- sqrt((var(sample1)+var(sample2))/2)
plot_data$t_stat <- (plot_data$sample_mean - input$mu0) / (sqrt(2*plot_data$sample_sd^2 / input$n)) # t-statistic
# Dynamically set x-axis limits for density plots, based on current input values
plot_data$x_limits_null <- c(
qt(0.0001, df = 2*input$n - 2),
qt(0.9999, df = 2*input$n - 2) + (input$mu_true - input$mu0) / (sqrt(2 * plot_data$sample_sd^2 / input$n))
)
plot_data$mean_x_limits <- c(
qnorm(0.0001, mean = input$mu0, sd = (sqrt(2*plot_data$sample_sd^2 / input$n))),
qnorm(0.9999, mean = input$mu_true, sd = (sqrt(2*plot_data$sample_sd^2 / input$n)))
)
# Render the summary statistics in the left panel
output$sample_stats <- renderUI({
summary_stats <- data.frame(
Statistic = c("Sample Mean", "Sample Standard Deviation", "Estimated Standard Error", "T-test Statistic"),
Value = c(
round(plot_data$sample_mean, 3),
round(plot_data$sample_sd, 3),
round((sqrt(2*plot_data$sample_sd^2 / input$n)), 3),
paste0("<span style='color:springgreen4;'>", round(plot_data$t_stat, 3), "</span>")
),
stringsAsFactors = FALSE
)
HTML(
paste0("<table style='border-collapse: collapse; width: 100%;font-size: 16px;'>",
"<tr><th style='text-align: center; border: 1px solid black;'>Statistic</th>",
"<th style='text-align: center; border: 1px solid black;'>Value</th></tr>",
paste0("<tr><td style='text-align: center; border: 1px solid black;'>", summary_stats$Statistic, "</td>",
"<td style='text-align: center;border: 1px solid black;'>", summary_stats$Value, "</td></tr>",
collapse = ""),
"</table>")
)
})
# Plot for the t-test statistic distribution
output$densityPlot <- renderPlot({
# Base plot for the null hypothesis density
p <- ggplot(data.frame(x = plot_data$x_limits_null), aes(x)) +
stat_function(fun = stats::dt, args = list(df = input$n - 1), color = "firebrick2") +
geom_vline(xintercept = plot_data$t_stat, color = "springgreen4", linetype = "dashed") +
geom_text(aes(x = plot_data$t_stat, label = paste("T =", round(plot_data$t_stat, 2)), y = Inf),
color = "springgreen4", angle = 90, vjust = 0, hjust = 1, size = 4.5) +
labs(title = "Test Statistic Distribution under Null Hypothesis",
x = "T statistic", y = "Density") + theme_classic()
# Add alternative distribution if showAlt is TRUE
if (showAlt()) {
p <- p + stat_function(
fun = dt,
args = list(
# 'ncp' is the non-centrality parameter
ncp = input$mu_true / (sqrt(2*plot_data$sample_sd^2 / input$n)),
# 'df' is the degrees of freedom
df = 2*input$n - 2
),
color = "lightblue", alpha = 0.7) +
labs(title = "Test Statistic Distribution: Null vs Alternative")
}
# Show rejection region if the toggle is on
if (showCI()) {
alpha <- input$alpha
test_type <- input$test_type
if (test_type == "two") {
# Calculate critical t-values for two-sided test
critical_value_low <- qt(alpha / 2, df = 2*input$n - 2)
critical_value_high <- qt(1 - alpha / 2, df = 2*input$n - 2)
# Add shaded rejection regions using geom_rect
p <- p + geom_rect(aes(xmin = -Inf, xmax = critical_value_low, ymin = 0, ymax = Inf),
fill = "red", alpha = 0.05) +
geom_rect(aes(xmin = critical_value_high, xmax = Inf, ymin = 0, ymax = Inf),
fill = "red", alpha = 0.05)
} else if (test_type == "one left") {
critical_value_high <- qt(1-alpha, df = 2*input$n - 2)
p <- p + geom_rect(aes(xmin = -Inf, xmax = critical_value_high, ymin = 0, ymax = Inf),
fill = "red", alpha = 0.05)
} else {
# One-sided test
critical_value_low <- qt(1-alpha, df = 2*input$n - 2)
# Add a shaded rejection region for one-sided test
p <- p + geom_rect(aes(xmin = critical_value_low, xmax = Inf, ymin = 0, ymax = Inf),
fill = "red", alpha = 0.05)
} }
p
})
# Plot for the sample histogram
output$samplePlot <- renderPlot({
ggplot(sample_data(), aes(x = x, fill = group)) +
geom_histogram(bins = 10, alpha = 0.6, color = "black", position = "identity") +
scale_fill_manual(values = c("steelblue4", "orange3")) +
labs(title = paste0("Sample (µ = ", input$mu_true,")"),
x = "Sample Values", y = "Frequency") + theme_classic()
})
# Plot for the expected sample mean distribution
output$meanDensityPlot <- renderPlot({
# Set appropriate limits for the expected sample mean density plot
mean_x_limits <- c(qnorm(0.0001, mean = input$mu0, sd = (sqrt(2*plot_data$sample_sd^2 / input$n))),
qnorm(0.9999, mean = input$mu_true, sd = (sqrt(2*plot_data$sample_sd^2 / input$n)) ))
se <- (sqrt(2*plot_data$sample_sd^2 / input$n))
alpha <- input$alpha
test_type <- input$test_type
if (test_type == "one left") {
ci_high<- plot_data$sample_mean + qt(1-alpha, df = 2*input$n-2)*se
ci_low <- -Inf # For one-sided, CI extends to infinity on the right side
} else if (test_type == "one right"){
ci_low <- plot_data$sample_mean - qt(1-alpha, df = 2*input$n-2)*se
ci_high <- Inf # For one-sided, CI extends to infinity on the right side
} else {
ci_low <- plot_data$sample_mean - qt(1-alpha/2, df = 2*input$n-2)*se
ci_high <- plot_data$sample_mean + qt(1-alpha/2, df = 2*input$n-2)*se
}
p <- ggplot(data.frame(x = plot_data$mean_x_limits), aes(x)) +
stat_function(fun = stats::dnorm, args = list(mean = input$mu0, sd = (sqrt(2*plot_data$sample_sd^2 / input$n))), color = "firebrick2") +
geom_vline(xintercept = plot_data$sample_mean, color = "orange3", linetype = "dashed") +
geom_text(aes(x = plot_data$sample_mean, label = paste("X̄ =", round(plot_data$sample_mean, 2)), y = Inf),
color = "orange3", angle = 90, vjust = 0, hjust = 1, size = 4.5) +
labs(title = "Distribution of Sample Mean",
x = "Sample Mean (X̄)", y = "Density") + theme_classic()
# Show confidence interval if the toggle is on
if (showCI()) {
p <- p + geom_segment(aes(x = ci_low, xend = ci_high, y = 0, yend = 0),
color = "blue", size = 2) +
geom_rect(aes(xmin = ci_low, xmax = ci_high, ymin = 0, ymax = Inf),
fill = "blue", alpha = 0.05) +
geom_vline(aes(xintercept = input$mu0), color = "firebrick1", linetype = 2)
}
# Add alternative distribution if showAlt is TRUE
if (showAlt()) {
p <- p + stat_function(fun = dnorm,
args = list(mean = input$mu_true,
sd = (sqrt(2*plot_data$sample_sd^2 / input$n))),
color = "lightblue", alpha = 0.7) +
labs(title = "Sample Mean: Null vs Alternative")
}
p # Return the plot
})
})
calc_power <- reactive({
if(input$mu0 == input$mu_true) {
return("Power is not well-defined under the null hypothesis")
} else {
# Map test type to appropriate alternative hypothesis for power.t.test()
sides <- switch(input$test_type,
"two" = "two.sided",
"one left" = "one.sided",
"one right" = "one.sided")
# Calculate power
pwr <- power.t.test(delta = abs(input$mu_true - input$mu0),
n = input$n,
sd = input$sigma,
sig.level = input$alpha,
type = "two.sample",
alternative = sides)$power
# Return power value as a string
return(paste0("Power: ", round(pwr, digits = 4)))
}
})
# Output power as text
output$power <- renderText({
calc_power() # Use the reactive expression to get the power
})
# Trigger power calculation on app startup (using an observer)
observe({
output$power <- renderText({
calc_power()
})
})
}
shinyApp(ui = ui, server = server)