Files
admin 33e9543b15 Upload to Server
Uploading to server
2025-08-02 05:15:23 +07:00

216 lines
9.4 KiB
R

# Tải các thư viện cần thiết
library(shiny)
library(bslib)
library(ggplot2)
library(shinycssloaders)
# --- Giao diện người dùng (UI) ---
ui <- fluidPage(
theme = bs_theme(version = 5, bootswatch = "cerulean"),
withMathJax(),
titlePanel("Tính toán Cỡ mẫu cho Kiểm định F (so sánh 2 phương sai)"),
sidebarLayout(
sidebarPanel(
h4("Tham số đầu vào"),
numericInput("v_ratio",
label = "Tỷ lệ phương sai (\\(\\sigma_1^2 / \\sigma_2^2\\)):",
value = 2, min = 1.1, step = 0.1),
helpText("Đây là effect size. Giá trị 2 có nghĩa là phương sai của nhóm 1 được giả định lớn gấp đôi phương sai của nhóm 2."),
selectInput("alternative", "Loại kiểm định:",
choices = c("Hai phía (Two-sided)" = "two.sided",
"Nhỏ hơn (Less)" = "less",
"Lớn hơn (Greater)" = "greater")),
sliderInput("sig_level",
label = "Mức ý nghĩa (\\(\\alpha\\)):",
min = 0.01, max = 0.10, value = 0.05, step = 0.01),
sliderInput("power",
label = "Power mong muốn (\\(1 - \\beta\\)):",
min = 0.50, max = 0.99, value = 0.80, step = 0.01),
hr(),
helpText("Kết quả và đồ thị sẽ tự động cập nhật ngay lập tức.")
),
mainPanel(
tabsetPanel(
id = "results_tabs",
type = "pills",
tabPanel(
"Kết quả & Diễn giải",
h4("Kết quả tính toán"),
withSpinner(uiOutput("sample_size_output"), type = 6, color = "#007bff")
),
tabPanel(
"Đồ thị Power vs. Cỡ mẫu",
h4("Mối quan hệ giữa Power và Cỡ mẫu"),
withSpinner(plotOutput("power_plot"), type = 6, color = "#007bff")
),
tabPanel(
"Phân tích Effect Size",
h4("Mối quan hệ giữa Cỡ mẫu và Effect Size"),
p("Đồ thị này cho thấy cỡ mẫu cần thiết thay đổi như thế nào khi mức độ khác biệt về phương sai thay đổi."),
withSpinner(plotOutput("effect_analysis_plot"), type = 6, color = "#007bff")
),
tabPanel(
"Giả thuyết và Công thức (Helper)",
h4("Giả thuyết của Kiểm định F"),
p("Kiểm định F được sử dụng để so sánh phương sai của hai mẫu độc lập."),
p("$$H_0: \\sigma_1^2 = \\sigma_2^2 \\quad (\\text{tỷ lệ phương sai bằng 1})$$"),
p("$$H_a: \\sigma_1^2 \\neq \\sigma_2^2 \\quad (\\text{hoặc } > \\text{ hoặc } <\\text{)}$$"),
tags$div(class = "alert alert-warning",
tags$b("Giả định quan trọng:"), " Kiểm định F rất nhạy với giả định rằng dữ liệu trong cả hai nhóm phải tuân theo phân phối chuẩn. Nếu giả định này bị vi phạm, hãy cân nhắc sử dụng kiểm định Levene hoặc Brown-Forsythe."),
hr(),
h4("Công thức tính toán"),
tags$b("1. Thống kê kiểm định (Test Statistic):"),
p("Giá trị thống kê F được tính bằng tỷ lệ của hai phương sai mẫu:"),
p("$$ F = \\frac{s_1^2}{s_2^2} $$"),
p("Trong đó \\(s_1^2\\) và \\(s_2^2\\) là phương sai của mẫu 1 và mẫu 2, với cỡ mẫu tương ứng là \\(n_1\\) và \\(n_2\\)."),
tags$b("2. Phân phối của thống kê kiểm định:"),
p("Dưới giả thuyết \\(H_0\\) (khi \\(\\sigma_1^2 = \\sigma_2^2\\)), thống kê F tuân theo phân phối F với bậc tự do \\(df_1 = n_1 - 1\\) và \\(df_2 = n_2 - 1\\)."),
p("$$ F \\sim F(n_1 - 1, n_2 - 1) $$"),
p("Dưới giả thuyết \\(H_a\\) (khi \\(\\frac{\\sigma_1^2}{\\sigma_2^2} = \\lambda > 1\\)), giá trị \\(F/\\lambda\\) tuân theo phân phối F nói trên. Điều này rất quan trọng để tính power."),
tags$b("3. Tính toán Power:"),
p("Power là xác suất bác bỏ \\(H_0\\) một cách chính xác khi \\(H_a\\) là đúng. Công thức phụ thuộc vào loại kiểm định:"),
tags$i("a) Kiểm định hai phía (Two-sided):"),
p("$$ \\text{Power} = P(F < F_{crit, lower} | H_a) + P(F > F_{crit, upper} | H_a) $$"),
p("Trong đó \\( F_{crit, lower} = F_{\\alpha/2, n_1-1, n_2-1} \\) và \\( F_{crit, upper} = F_{1-\\alpha/2, n_1-1, n_2-1} \\)."),
tags$i("b) Kiểm định lớn hơn (Greater):"),
p("$$ \\text{Power} = P(F > F_{crit, upper} | H_a) $$"),
p("Trong đó \\( F_{crit, upper} = F_{1-\\alpha, n_1-1, n_2-1} \\)."),
tags$i("c) Kiểm định nhỏ hơn (Less):"),
p("$$ \\text{Power} = P(F < F_{crit, lower} | H_a) $$"),
p("Trong đó \\( F_{crit, lower} = F_{\\alpha, n_1-1, n_2-1} \\).")
)
)
)
)
)
# --- Logic của máy chủ (Server) ---
server <- function(input, output, session) {
# --- HÀM TÍNH POWER LÕI ---
calculate_f_power <- function(n1, n2, v_ratio, sig_level, alternative) {
if (alternative == "two.sided") {
alpha <- sig_level / 2
crit_lower <- qf(alpha, n1 - 1, n2 - 1)
crit_upper <- qf(1 - alpha, n1 - 1, n2 - 1)
power <- pf(crit_lower / v_ratio, n1 - 1, n2 - 1) + (1 - pf(crit_upper / v_ratio, n1 - 1, n2 - 1))
} else if (alternative == "greater") {
crit_upper <- qf(1 - sig_level, n1 - 1, n2 - 1)
power <- 1 - pf(crit_upper / v_ratio, n1 - 1, n2 - 1)
} else { # less
crit_lower <- qf(sig_level, n1 - 1, n2 - 1)
power <- pf(crit_lower / v_ratio, n1 - 1, n2 - 1)
}
return(power)
}
# --- PHẦN TÍNH TOÁN CHÍNH ---
# Hoàn toàn reactive, không cần debounce
main_results <- reactive({
req(input$v_ratio, input$sig_level, input$power, input$alternative)
# Tìm cỡ mẫu cần thiết
n <- 5 # Bắt đầu từ n nhỏ
power_val <- 0
while (power_val < input$power && n <= 2000) {
n <- n + 1
power_val <- calculate_f_power(n, n, input$v_ratio, input$sig_level, input$alternative)
}
required_n <- ifelse(n > 2000, NA, n)
# Tạo dữ liệu cho đồ thị Power
# Đảm bảo dãy sample_sizes hợp lệ ngay cả khi không tìm thấy required_n
upper_bound <- if (!is.na(required_n)) required_n + 50 else 200
sample_sizes <- seq(5, upper_bound, by = 1)
powers <- sapply(sample_sizes, function(s_n) {
calculate_f_power(s_n, s_n, input$v_ratio, input$sig_level, input$alternative)
})
list(
required_n = required_n,
power_plot_data = data.frame(SampleSize = sample_sizes, Power = powers)
)
})
output$sample_size_output <- renderUI({
res <- main_results()
if (!is.na(res$required_n)) {
tagList(
tags$p("Để phát hiện sự khác biệt về phương sai (tỷ lệ =", tags$b(input$v_ratio), ") với power là", tags$b(input$power), "và mức ý nghĩa", tags$b(input$sig_level), ", bạn cần một cỡ mẫu ước tính là:"),
tags$h3(style = "color: #007bff; text-align: center;", paste(res$required_n, "cho mỗi nhóm")),
tags$p(style = "text-align: center; font-style: italic;", "Tổng cỡ mẫu là ", res$required_n * 2, ".")
)
} else {
tags$div(class = "alert alert-warning", "Không thể đạt được power mong muốn với cỡ mẫu tối đa (2000). Vui lòng xem xét lại các tham số.")
}
})
output$power_plot <- renderPlot({
res <- main_results()
ggplot(res$power_plot_data, aes(x = SampleSize, y = Power)) +
geom_line(color = "#007bff", size = 1.2) +
geom_hline(yintercept = input$power, linetype = "dashed", color = "red") +
{
if(!is.na(res$required_n)) {
geom_vline(xintercept = res$required_n, linetype = "dashed", color = "darkgreen")
}
} +
labs(title = "Power vs. Cỡ mẫu (cho mỗi nhóm)", x = "Cỡ mẫu cho mỗi nhóm (n)", y = "Power (1 - β)") +
scale_y_continuous(limits = c(0, 1)) + theme_minimal(base_size = 14)
})
# --- PHẦN PHÂN TÍCH EFFECT SIZE ---
effect_analysis_plot_data <- reactive({
req(input$sig_level, input$power, input$alternative)
effect_sizes <- seq(1.2, 3.0, by = 0.1)
required_n_values <- sapply(effect_sizes, function(ratio) {
n <- 5
power_val <- 0
while (power_val < input$power && n <= 2000) {
n <- n + 1
power_val <- calculate_f_power(n, n, ratio, input$sig_level, input$alternative)
}
ifelse(n > 2000, NA, n)
})
data.frame(EffectSize = effect_sizes, RequiredN = required_n_values)
})
output$effect_analysis_plot <- renderPlot({
plot_data <- effect_analysis_plot_data()
ggplot(plot_data, aes(x = EffectSize, y = RequiredN)) +
geom_line(color = "#28a745", size = 1.2) +
geom_point(color = "#28a745", size = 3, na.rm = TRUE) +
labs(
title = paste("Cỡ mẫu cần thiết vs. Effect Size (Power cố định =", input$power, ")"),
x = "Tỷ lệ Phương sai (Effect Size)",
y = "Cỡ mẫu cần thiết cho mỗi nhóm (n)"
) +
theme_minimal(base_size = 14)
})
}
# Chạy ứng dụng Shiny
shinyApp(ui = ui, server = server)