278 lines
15 KiB
R
278 lines
15 KiB
R
# ==============================================================================
|
|
# ỨNG DỤNG SHINY TÍNH CỠ MẪU CHO HỒI QUY LOGISTIC (PHIÊN BẢN ỔN ĐỊNH)
|
|
# - Sửa lỗi hiển thị công thức bằng cách escape ký tự `\` trong R string.
|
|
# Author: Gemini & User Collaboration
|
|
# Date: 2025-10-17
|
|
# ==============================================================================
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# THIẾT LẬP: Tải các thư viện cần thiết
|
|
# ------------------------------------------------------------------------------
|
|
library(shiny)
|
|
library(bslib)
|
|
library(ggplot2)
|
|
library(shinycssloaders)
|
|
library(dplyr)
|
|
library(purrr)
|
|
library(scales)
|
|
|
|
# ==============================================================================
|
|
# PHẦN 1: CÁC HÀM HỖ TRỢ (HELPER FUNCTIONS)
|
|
# ==============================================================================
|
|
|
|
#' Tính cỡ mẫu cho hồi quy logistic theo công thức xấp xỉ của Hsieh (1989).
|
|
calc_n_logistic_hsieh <- function(beta, var_x, p_mean, power = 0.8, sig.level = 0.05) {
|
|
z_a <- qnorm(1 - sig.level / 2)
|
|
z_b <- qnorm(power)
|
|
denom <- (beta^2) * var_x * p_mean * (1 - p_mean)
|
|
if (denom <= 0) return(NA_real_)
|
|
n <- (z_a + z_b)^2 / denom
|
|
ceiling(n)
|
|
}
|
|
|
|
#' Ước tính công suất thực nghiệm cho hồi quy logistic bằng mô phỏng Monte Carlo.
|
|
simulate_logistic_power <- function(n, beta, intercept = 0, x_dist = list(type = "normal", mean = 0, sd = 1),
|
|
n_sim = 1000, alpha = 0.05) {
|
|
sim_results <- replicate(n_sim, {
|
|
x <- switch(x_dist$type,
|
|
"normal" = rnorm(n, mean = x_dist$mean, sd = x_dist$sd),
|
|
"binary" = rbinom(n, size = 1, prob = x_dist$prob),
|
|
"uniform" = runif(n, min = x_dist$min, max = x_dist$max))
|
|
|
|
linpred <- intercept + beta * x
|
|
p <- 1 / (1 + exp(-linpred))
|
|
y <- rbinom(n, 1, p)
|
|
|
|
fit <- try(glm(y ~ x, family = binomial), silent = TRUE)
|
|
if (inherits(fit, "try-error")) return(NA_real_)
|
|
|
|
coefs <- summary(fit)$coefficients
|
|
if ("x" %in% rownames(coefs)) {
|
|
as.numeric(coefs["x", "Pr(>|z|)"] < alpha)
|
|
} else {
|
|
0
|
|
}
|
|
})
|
|
mean(sim_results, na.rm = TRUE)
|
|
}
|
|
|
|
# ==============================================================================
|
|
# PHẦN 2: GIAO DIỆN NGƯỜI DÙNG (USER INTERFACE - UI)
|
|
# ==============================================================================
|
|
|
|
ui <- fluidPage(
|
|
theme = bs_theme(version = 5, bootswatch = "cerulean"),
|
|
withMathJax(),
|
|
|
|
titlePanel("Công Cụ Tính Cỡ Mẫu Cho Hồi Quy Logistic"),
|
|
|
|
sidebarLayout(
|
|
sidebarPanel(
|
|
width = 4,
|
|
h4("Nhập Tham Số"),
|
|
br(),
|
|
h5("Tham số mô hình"),
|
|
sliderInput("beta", "Giá trị β (log-odds) kỳ vọng:", value = 0.693, min = -2, max = 2, step = 0.01),
|
|
sliderInput("p_mean_event", "Tỉ lệ 'event' trung bình trong quần thể:", value = 0.25, min = 0.01, max = 0.99),
|
|
selectInput("x_type", "Phân phối của biến độc lập X:",
|
|
choices = c("Nhị phân (0/1)" = "binary", "Chuẩn (Normal)" = "normal", "Đồng nhất (Uniform)" = "uniform")),
|
|
conditionalPanel("input.x_type == 'normal'", numericInput("x_sd", "Độ lệch chuẩn của X:", value = 1, min = 0.1)),
|
|
conditionalPanel("input.x_type == 'binary'", sliderInput("x_prob", "Tỉ lệ P(X=1):", value = 0.3, min = 0.01, max = 0.99)),
|
|
conditionalPanel("input.x_type == 'uniform'", numericInput("x_min", "Min của X:", value = -1), numericInput("x_max", "Max của X:", value = 1)),
|
|
|
|
hr(),
|
|
h5("Tham số kiểm định"),
|
|
sliderInput("power_log", "Công suất mong muốn (1 - \\(\\beta\\)):", value = 0.8, min = 0.5, max = 0.99),
|
|
sliderInput("alpha_log", "Mức ý nghĩa (\\(\\alpha\\)):", value = 0.05, min = 0.01, max = 0.1),
|
|
hr(),
|
|
actionButton("go_log", "Tính toán & Phân tích", class = "btn-primary w-100", icon = icon("calculator"))
|
|
),
|
|
|
|
mainPanel(
|
|
width = 8,
|
|
tabsetPanel(
|
|
id = "log_results_tabs",
|
|
type = "pills",
|
|
tabPanel("Kết quả & Diễn giải",
|
|
withSpinner(uiOutput("log_result_text"), type = 6, color = "#007bff")),
|
|
tabPanel("Đồ thị Power & Mô phỏng",
|
|
withSpinner(plotOutput("log_power_plot"), type = 6, color = "#007bff"),
|
|
hr(),
|
|
h4("Phân tích Mô phỏng Chi tiết"),
|
|
p("Chạy mô phỏng sâu hơn với một cỡ mẫu cụ thể để xem phân phối của p-value và ước tính power chính xác hơn."),
|
|
fluidRow(
|
|
column(6, numericInput("n_for_sim_detail", "Nhập cỡ mẫu (n) để mô phỏng:", value = 200, min = 10)),
|
|
column(6, actionButton("run_sim_detail", "Chạy mô phỏng chi tiết", icon = icon("play-circle"), class="btn-success", style="margin-top: 25px;"))
|
|
),
|
|
withSpinner(plotOutput("log_sim_detail_plot"), type = 6, color = "#007bff")),
|
|
tabPanel("Phân tích Effect Size",
|
|
withSpinner(uiOutput("log_effect_size_ui"), type = 6, color = "#007bff")),
|
|
|
|
# --- TAB CÔNG THỨC ĐƯỢC ĐỊNH NGHĨA TĨNH TRONG UI (SỬA LỖI) ---
|
|
tabPanel("Giả thuyết, Công thức & Ví dụ",
|
|
h4("Giả thuyết thống kê"),
|
|
p("Kiểm định trong hồi quy logistic thường tập trung vào việc xem hệ số β có khác 0 một cách có ý nghĩa thống kê hay không."),
|
|
p("$$H_0: \\beta = 0$$"),
|
|
p(em("Diễn giải: Biến độc lập X không có ảnh hưởng đến log-odds của kết quả Y (tương đương Odds Ratio = 1).")),
|
|
p("$$H_a: \\beta \\neq 0$$"),
|
|
p(em("Diễn giải: Biến độc lập X có ảnh hưởng đến log-odds của kết quả Y (tương đương Odds Ratio ≠ 1).")),
|
|
hr(),
|
|
|
|
h4("Công thức tính cỡ mẫu (Xấp xỉ Hsieh, 1989)"),
|
|
p("$$n \\approx \\frac{(Z_{1-\\alpha/2} + Z_{\\text{power}})^2}{\\beta^2 \\, \\text{Var}(X) \\, \\bar{p}(1-\\bar{p})}$$"),
|
|
p(strong("Trong đó:")),
|
|
tags$ul(
|
|
tags$li("\\(Z_{1-\\alpha/2}\\) là giá trị Z-score ứng với mức ý nghĩa \\(\\alpha\\) (ví dụ: 1.96 cho \\(\\alpha=0.05\\))."),
|
|
tags$li("\\(Z_{\\text{power}}\\) là giá trị Z-score ứng với công suất mong muốn (ví dụ: 0.84 cho power=0.8)."),
|
|
tags$li("\\(\\beta\\) là hệ số hồi quy log-odds kỳ vọng (effect size)."),
|
|
tags$li("\\(\\text{Var}(X)\\) là phương sai của biến độc lập X."),
|
|
tags$li("\\(\\bar{p}\\) là tỉ lệ trung bình của kết quả (event) trong quần thể.")
|
|
),
|
|
hr(),
|
|
|
|
h4("Ví dụ trong Y tế công cộng"),
|
|
p(strong("Bối cảnh nghiên cứu:")),
|
|
p("Một nhà nghiên cứu muốn tính cỡ mẫu để xem liệu việc ", strong("hút thuốc lá (biến X)"), " có phải là yếu tố nguy cơ cho ", strong("bệnh tăng huyết áp (biến Y)"), " hay không."),
|
|
p(strong("Diễn giải các tham số:")),
|
|
tags$ul(
|
|
tags$li(strong("Kết quả (Event):"), " Bị tăng huyết áp (Y=1)."),
|
|
tags$li(strong("Biến độc lập (X):"), " Hút thuốc lá (X=1) vs. không hút (X=0). Đây là biến nhị phân."),
|
|
tags$li(strong("\\(\\beta\\) (Effect size):"), " Dựa trên y văn, nhà nghiên cứu kỳ vọng Odds Ratio (OR) của việc hút thuốc gây tăng huyết áp là 2.0. Do đó, effect size mong muốn là \\(\\beta = \\ln(OR) = \\ln(2.0) \\approx 0.693\\)."),
|
|
tags$li(strong("\\(\\bar{p}\\) (Tỉ lệ event TB):"), " Tỉ lệ mắc tăng huyết áp chung trong quần thể ước tính là 25% (\\(\\bar{p} = 0.25\\))."),
|
|
tags$li(strong("\\(\\text{Var}(X)\\):"), " Tỉ lệ người hút thuốc trong quần thể là 30% (P(X=1)=0.3). Phương sai của X là \\(\\text{Var}(X) = P(X=1) \\times (1-P(X=1)) = 0.3 \\times 0.7 = 0.21\\)."),
|
|
tags$li(strong("Power và \\(\\alpha\\):"), " Nghiên cứu được thiết kế để có 80% cơ hội phát hiện mối liên quan (power=0.8) với mức ý nghĩa 5% (\\(\\alpha=0.05\\)).")
|
|
),
|
|
p("Nhà nghiên cứu sẽ nhập các giá trị trên vào công cụ để ước tính cỡ mẫu cần thiết. Các giá trị này cũng đã được đặt làm mặc định trong ứng dụng để bạn dễ hình dung.")
|
|
)
|
|
)
|
|
)
|
|
)
|
|
)
|
|
|
|
# ==============================================================================
|
|
# PHẦN 3: LOGIC MÁY CHỦ (SERVER)
|
|
# ==============================================================================
|
|
server <- function(input, output, session) {
|
|
|
|
rv <- reactiveValues()
|
|
|
|
observeEvent(input$go_log, {
|
|
varx <- switch(input$x_type,
|
|
"normal" = input$x_sd^2,
|
|
"binary" = input$x_prob * (1 - input$x_prob),
|
|
"uniform" = (input$x_max - input$x_min)^2 / 12)
|
|
n_hsieh <- calc_n_logistic_hsieh(input$beta, varx, input$p_mean_event, input$power_log, input$alpha_log)
|
|
rv$log_n_hsieh <- n_hsieh
|
|
rv$log_varx <- varx
|
|
|
|
if (!is.na(n_hsieh)) {
|
|
n_range_log <- unique(round(seq(max(30, n_hsieh * 0.5), n_hsieh * 1.5, length.out = 15)))
|
|
xdist_log <- switch(input$x_type,
|
|
"normal" = list(type = "normal", mean = 0, sd = input$x_sd),
|
|
"binary" = list(type = "binary", prob = input$x_prob),
|
|
"uniform" = list(type = "uniform", min = input$x_min, max = input$x_max))
|
|
|
|
power_values_log <- map_dbl(n_range_log, ~simulate_logistic_power(n = .x, beta = input$beta, x_dist = xdist_log, n_sim = 500, alpha = input$alpha_log))
|
|
rv$log_plot_data <- tibble(SampleSize = n_range_log, Power = power_values_log)
|
|
} else {
|
|
rv$log_plot_data <- NULL
|
|
}
|
|
})
|
|
|
|
output$log_result_text <- renderUI({
|
|
if (input$go_log == 0) {
|
|
return(tags$div(class="alert alert-info", "Nhập các tham số và nhấn 'Tính toán & Phân tích' để xem kết quả."))
|
|
}
|
|
|
|
req(rv$log_n_hsieh)
|
|
n <- rv$log_n_hsieh
|
|
|
|
if (is.na(n)) {
|
|
return(tags$div(class = "alert alert-danger", "Không thể tính toán. Kiểm tra lại các tham số (ví dụ: mẫu số của công thức có thể bằng 0)."))
|
|
}
|
|
|
|
tagList(
|
|
h4("Kết quả tính toán (Xấp xỉ Hsieh)"),
|
|
p("Với Var(X) ≈", tags$b(round(rv$log_varx, 3)), ", để phát hiện một \\(\\beta\\) là", tags$b(input$beta), "với power", tags$b(input$power_log), "và \\(\\alpha\\) là", tags$b(input$alpha_log), ", cỡ mẫu ước tính là:"),
|
|
tags$h3(style = "color: #007bff; text-align: center;", n, " quan sát"),
|
|
hr(),
|
|
tags$div(class = "alert alert-light",
|
|
tags$b("Ghi chú:"), " Đây là kết quả xấp xỉ. Biểu đồ trong tab 'Đồ thị Power & Mô phỏng' được tạo ra bằng mô phỏng và có thể cho kết quả chính xác hơn.")
|
|
)
|
|
})
|
|
|
|
output$log_power_plot <- renderPlot({
|
|
req(rv$log_plot_data, rv$log_n_hsieh)
|
|
ggplot(rv$log_plot_data, aes(x = SampleSize, y = Power)) +
|
|
geom_line(color = "#007bff", size = 1.2) +
|
|
geom_point(color = "#007bff", size = 3) +
|
|
geom_hline(yintercept = input$power_log, linetype = "dashed", color = "red") +
|
|
geom_vline(xintercept = rv$log_n_hsieh, linetype = "dotted", color = "darkorange", size=1.2) +
|
|
labs(title = "Power thực nghiệm (Mô phỏng) vs. Cỡ mẫu", x = "Cỡ mẫu (n)", y = "Power (1 - \\(\\beta\\))") +
|
|
annotate("text", x = rv$log_n_hsieh * 1.05, y = 0.1, label = paste("Hsieh Approx.\nn =", rv$log_n_hsieh), color = "darkorange", hjust = 0) +
|
|
scale_y_continuous(labels = percent, limits = c(0, 1)) +
|
|
theme_minimal(base_size = 14)
|
|
})
|
|
|
|
output$log_effect_size_ui <- renderUI({
|
|
or <- exp(input$beta)
|
|
interpretation <- if (or > 1) {
|
|
paste0("một sự gia tăng ", round((or - 1) * 100, 1), "% trong \"odds\" của kết quả (event) xảy ra.")
|
|
} else {
|
|
paste0("một sự sụt giảm ", round((1-or) * 100, 1), "% trong \"odds\" của kết quả (event) xảy ra.")
|
|
}
|
|
|
|
tagList(
|
|
h4("Phân tích Effect Size: Tỷ số chênh (Odds Ratio - OR)"),
|
|
p("Trong hồi quy logistic, effect size thường được biểu diễn bằng Tỷ số chênh (OR), được tính bằng công thức \\(OR = e^\\beta\\)."),
|
|
p("Nó cho biết odds của một \"event\" (kết quả Y=1) thay đổi như thế nào khi biến độc lập X tăng lên một đơn vị."),
|
|
hr(),
|
|
p("Với giá trị \\(\\beta\\) bạn đã chọn là ", tags$b(input$beta), ", Tỷ số chênh tương ứng là:"),
|
|
tags$h3(style = "color: #007bff; text-align: center;", round(or, 3)),
|
|
tags$div(class = "alert alert-light",
|
|
tags$b("Diễn giải:"),
|
|
p("Khi biến X tăng lên một đơn vị, điều này tương ứng với ", tags$b(interpretation))
|
|
)
|
|
)
|
|
})
|
|
|
|
observeEvent(input$run_sim_detail, {
|
|
xdist_detail <- switch(input$x_type,
|
|
"normal" = list(type = "normal", mean = 0, sd = input$x_sd),
|
|
"binary" = list(type = "binary", prob = input$x_prob),
|
|
"uniform" = list(type = "uniform", min = input$x_min, max = input$x_max))
|
|
|
|
pvals <- replicate(1000, {
|
|
x <- switch(xdist_detail$type, "normal" = rnorm(input$n_for_sim_detail, 0, xdist_detail$sd), "binary" = rbinom(input$n_for_sim_detail, 1, xdist_detail$prob), "uniform" = runif(input$n_for_sim_detail, xdist_detail$min, xdist_detail$max))
|
|
p <- 1 / (1 + exp(-(0 + input$beta * x)))
|
|
y <- rbinom(input$n_for_sim_detail, 1, p)
|
|
fit <- try(glm(y ~ x, family = binomial), silent = TRUE)
|
|
if (!inherits(fit, "try-error") && "x" %in% rownames(summary(fit)$coefficients)) {
|
|
summary(fit)$coefficients["x", "Pr(>|z|)"]
|
|
} else {
|
|
NA_real_
|
|
}
|
|
})
|
|
rv$log_sim_detail_data <- tibble(p_value = pvals)
|
|
})
|
|
|
|
output$log_sim_detail_plot <- renderPlot({
|
|
req(rv$log_sim_detail_data)
|
|
power_est <- mean(rv$log_sim_detail_data$p_value < input$alpha_log, na.rm=TRUE)
|
|
|
|
ggplot(rv$log_sim_detail_data, aes(x = p_value)) +
|
|
geom_histogram(bins = 30, fill = "#28a745", color = "black", boundary=0) +
|
|
geom_vline(xintercept = input$alpha_log, linetype = "dashed", color = "red", size = 1) +
|
|
labs(
|
|
title = paste0("Phân phối P-value từ mô phỏng (n=", input$n_for_sim_detail, ")"),
|
|
subtitle = paste0("Công suất thực nghiệm ước tính: ", percent(power_est, accuracy = 0.1)),
|
|
x = "P-value", y = "Tần suất"
|
|
) +
|
|
theme_minimal(base_size = 14)
|
|
})
|
|
}
|
|
|
|
# ==============================================================================
|
|
# PHẦN 4: CHẠY ỨNG DỤNG
|
|
# ==============================================================================
|
|
shinyApp(ui, server) |