228 lines
9.1 KiB
R
228 lines
9.1 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 Z hai mẫu độc lập"),
|
|
|
|
sidebarLayout(
|
|
sidebarPanel(
|
|
h4("Tham số đầu vào"),
|
|
|
|
numericInput("d",
|
|
label = "Effect Size (Cohen's d):",
|
|
value = 0.5, min = 0.1, step = 0.1),
|
|
helpText("Cohen's d đo lường độ lớn của sự khác biệt giữa hai trung bình, được chuẩn hóa bằng độ lệch chuẩn TỔNG THỂ (đã biết)."),
|
|
|
|
selectInput("alternative", "Loại kiểm định:",
|
|
choices = c("Hai phía (Two-sided)" = "two.sided",
|
|
"Lớn hơn (Greater)" = "greater",
|
|
"Nhỏ hơn (Less)" = "less")),
|
|
|
|
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(),
|
|
h5("Công thức tính cỡ mẫu (mỗi nhóm):"),
|
|
# Hiển thị công thức tính mẫu trực tiếp
|
|
uiOutput("sample_size_formula_display")
|
|
),
|
|
|
|
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 Effect Size thay đổi."),
|
|
withSpinner(plotOutput("effect_analysis_plot"), type = 6, color = "#007bff")
|
|
),
|
|
|
|
tabPanel(
|
|
"Giả thuyết và Công thức (Helper)",
|
|
h4("Sự khác biệt chính: Z-test vs. T-test (cho hai mẫu)"),
|
|
tags$div(class = "alert alert-info",
|
|
tags$b("Điểm mấu chốt:"),
|
|
p("Sự lựa chọn giữa Z-test và T-test phụ thuộc vào việc bạn có biết phương sai của tổng thể (\\(\\sigma^2\\)) hay không."),
|
|
tags$ul(
|
|
tags$li(tags$b("Sử dụng Z-test khi:"), "Phương sai của CẢ HAI tổng thể (\\(\\sigma_1^2\\) và \\(\\sigma_2^2\\)) đều đã biết. Đây là một giả định rất chặt và hiếm khi xảy ra trong thực tế."),
|
|
tags$li(tags$b("Sử dụng T-test khi:"), "Phương sai của tổng thể không xác định và phải được ước tính từ dữ liệu mẫu. Đây là trường hợp gần như luôn luôn xảy ra trong nghiên cứu thực tế.")
|
|
)
|
|
),
|
|
hr(),
|
|
h4("Giả thuyết của Kiểm định Z hai mẫu"),
|
|
p("Kiểm định này so sánh trung bình của hai nhóm độc lập (\\(\\mu_1\\) và \\(\\mu_2\\)), khi phương sai của chúng đã biết."),
|
|
p("$$H_0: \\mu_1 = \\mu_2$$"),
|
|
p("$$H_a: \\mu_1 \\neq \\mu_2 \\quad (\\text{hoặc } > \\text{ hoặc } <\\text{)}$$"),
|
|
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ê Z được tính như sau:"),
|
|
p("$$ Z = \\frac{(\\bar{x}_1 - \\bar{x}_2) - 0}{\\sqrt{\\frac{\\sigma_1^2}{n_1} + \\frac{\\sigma_2^2}{n_2}}} $$"),
|
|
p("Trong đó \\(\\sigma_1^2\\) và \\(\\sigma_2^2\\) là các phương sai TỔNG THỂ đã biết.")
|
|
)
|
|
)
|
|
)
|
|
)
|
|
)
|
|
|
|
# --- Logic của máy chủ (Server) ---
|
|
server <- function(input, output, session) {
|
|
|
|
# Hiển thị công thức tính cỡ mẫu
|
|
output$sample_size_formula_display <- renderUI({
|
|
alpha <- input$sig_level
|
|
beta <- 1 - input$power
|
|
d <- input$d
|
|
|
|
if (input$alternative == "two.sided") {
|
|
withMathJax(paste0("$$ n = \\frac{2(Z_{1-\\alpha/2} + Z_{1-\\beta})^2}{d^2} $$"))
|
|
} else {
|
|
withMathJax(paste0("$$ n = \\frac{2(Z_{1-\\alpha} + Z_{1-\\beta})^2}{d^2} $$"))
|
|
}
|
|
})
|
|
|
|
# --- HÀM TÍNH POWER LÕI CHO Z-TEST 2 MẪU ---
|
|
calculate_z_power <- function(n, d, sig_level, alternative) {
|
|
# Đối với 2 mẫu, mean của Z-statistic dưới Ha là d * sqrt(n/2)
|
|
non_centrality_param <- d * sqrt(n / 2)
|
|
|
|
if (alternative == "two.sided") {
|
|
alpha <- sig_level / 2
|
|
crit_val_upper <- qnorm(1 - alpha)
|
|
power <- pnorm(crit_val_upper - non_centrality_param, lower.tail = FALSE) +
|
|
pnorm(-crit_val_upper - non_centrality_param, lower.tail = TRUE)
|
|
} else { # one.sided
|
|
crit_val <- qnorm(1 - sig_level)
|
|
power <- pnorm(crit_val - non_centrality_param, lower.tail = FALSE)
|
|
}
|
|
return(power)
|
|
}
|
|
|
|
# --- PHẦN TÍNH TOÁN CHÍNH ---
|
|
main_results <- reactive({
|
|
req(input$d, input$sig_level, input$power, input$alternative)
|
|
|
|
# Tìm cỡ mẫu cần thiết
|
|
n <- 2 # Bắt đầu từ n nhỏ
|
|
power_val <- 0
|
|
while (power_val < input$power && n <= 5000) {
|
|
n <- n + 1
|
|
power_val <- calculate_z_power(n, input$d, input$sig_level, input$alternative)
|
|
}
|
|
required_n <- ifelse(n > 5000, NA, n)
|
|
|
|
# Tạo dữ liệu cho đồ thị Power
|
|
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_z_power(s_n, input$d, 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.null(res) || is.na(res$required_n)) {
|
|
return(tags$div(class = "alert alert-warning", "Không thể đạt được power mong muốn với cỡ mẫu tối đa (5000). Vui lòng xem xét lại các tham số."))
|
|
}
|
|
|
|
tagList(
|
|
tags$p("Để phát hiện một effect size (Cohen's d) là", tags$b(input$d), "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, ".")
|
|
)
|
|
})
|
|
|
|
output$power_plot <- renderPlot({
|
|
res <- main_results()
|
|
req(res$power_plot_data)
|
|
|
|
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") +
|
|
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(0.1, 1.5, by = 0.05)
|
|
|
|
required_n_values <- sapply(effect_sizes, function(d_val) {
|
|
n <- 2
|
|
power_val <- 0
|
|
while (power_val < input$power && n <= 5000) {
|
|
n <- n + 1
|
|
power_val <- calculate_z_power(n, d_val, input$sig_level, input$alternative)
|
|
}
|
|
ifelse(n > 5000, NA, n)
|
|
})
|
|
|
|
data.frame(EffectSize = effect_sizes, RequiredN = required_n_values)
|
|
})
|
|
|
|
output$effect_analysis_plot <- renderPlot({
|
|
plot_data <- effect_analysis_plot_data()
|
|
res <- main_results()
|
|
|
|
p <- 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 = "Effect Size (Cohen's d)",
|
|
y = "Cỡ mẫu cần thiết cho mỗi nhóm (n)"
|
|
) +
|
|
theme_minimal(base_size = 14)
|
|
|
|
if (!is.null(res) && !is.na(res$required_n)) {
|
|
p <- p +
|
|
geom_vline(xintercept = input$d, linetype = "dotted", color = "blue", size = 1) +
|
|
geom_point(aes(x = input$d, y = res$required_n), color = "blue", size = 5, shape = 18) +
|
|
annotate("text", x = input$d, y = res$required_n,
|
|
label = paste("n =", res$required_n), vjust = -1.5, color = "blue", fontface = "bold")
|
|
}
|
|
|
|
p
|
|
})
|
|
}
|
|
|
|
# Chạy ứng dụng Shiny
|
|
shinyApp(ui = ui, server = server)
|