Commit ae3d0b92 authored by wuzekai's avatar wuzekai

update:方差分析与KW分析

parent d85895fb
......@@ -48,6 +48,7 @@ export(compare_props)
export(cor2df)
export(correlation)
export(cross_tabs)
export(get_single_norm)
export(goodness)
export(homo_variance_test)
export(mda)
......@@ -65,6 +66,10 @@ export(prob_unif)
export(radiant.basics)
export(radiant.basics_viewer)
export(radiant.basics_window)
export(run_anova)
export(run_homo_test)
export(run_kw)
export(run_norm_test)
export(single_mean)
export(single_prop)
import(ggplot2)
......
......@@ -6,18 +6,18 @@ mda <- function(dataset,
var,
group,
normality_type = c("overall", "by_group"),
data_filter = "",
test_method = c("anova", "kw"),
envir = parent.frame()) {
# 1. 基础参数处理
normality_type <- match.arg(normality_type, choices = c("overall", "by_group"))
test_method <- match.arg(test_method, choices = c("anova", "kw"))
df_name <- if (is_string(dataset)) dataset else deparse(substitute(dataset))
# 2. 数据提取:只保留“因变量+分组变量”
dataset <- get_data(
dataset,
vars = c(var, group), # 强制只取2个核心变量,剔除冗余列
filt = data_filter,
na.rm = FALSE, # 先不删缺失值,后续统一过滤
envir = envir
)
......@@ -51,6 +51,13 @@ mda <- function(dataset,
homo_res <- run_homo_test(valid_data, var, group) # 方差齐性检验
norm_res <- run_norm_test(valid_data, var, group, normality_type) # 正态性检验
core_test_res <- NULL
if (test_method == "anova") {
core_test_res <- run_anova(valid_data, var, group)
} else if (test_method == "kw") {
core_test_res <- run_kw(valid_data, var, group)
}
# 7. 绘图数据准备
plot_obj <- list(
norm = list(
......@@ -77,10 +84,11 @@ mda <- function(dataset,
var = var,
group = group,
normality_type = normality_type,
data_filter = if (data_filter == "") "None" else data_filter,
valid_n = nrow(valid_data), # 有效样本量
homo_res = homo_res, # 方差齐性检验结果
norm_res = norm_res, # 正态性检验结果
test_method = test_method,
core_test_res = core_test_res,
plot_obj = plot_obj
),
class = "mda"
......@@ -91,6 +99,7 @@ mda <- function(dataset,
# ------------------------------
# 辅助函数1:方差齐性检验
# ------------------------------
#' @export
run_homo_test <- function(valid_data, var, group) {
x <- valid_data[[var]]
g <- valid_data[[group]]
......@@ -165,6 +174,7 @@ run_homo_test <- function(valid_data, var, group) {
# ------------------------------
# 辅助函数2:正态性检验
# ------------------------------
#' @export
run_norm_test <- function(valid_data, var, group, normality_type) {
x <- valid_data[[var]]
g <- valid_data[[group]]
......@@ -189,6 +199,7 @@ run_norm_test <- function(valid_data, var, group, normality_type) {
# ------------------------------
# 辅助函数3:单组正态性检验
# ------------------------------
#' @export
get_single_norm <- function(x, group_label) {
res <- tibble::tibble(Group = group_label, Test = character(), Statistic = numeric(), p.value = numeric())
n <- length(x)
......@@ -282,12 +293,110 @@ get_single_norm <- function(x, group_label) {
res
}
# ------------------------------
# 辅助函数4:方差分析(ANOVA)
# ------------------------------
#' @export
run_anova <- function(valid_data, var, group) {
tryCatch({
# 1. 使用列名创建公式
formula <- reformulate(group, var)
aov_model <- stats::aov(formula, data = valid_data)
aov_summary <- summary(aov_model)[[1]]
# 2. 检查行数
if (nrow(aov_summary) < 2) {
stop("ANOVA结果不完整,无法提取统计量")
}
# 3. 稳健提取:假设最后一行是残差,其余是模型项
residual_row_idx <- nrow(aov_summary)
model_row_idx <- 1:(nrow(aov_summary) - 1)
# 4. 安全提取函数(避免空值)
safe_extract <- function(df, row_idx, col, default = NA_real_) {
if (length(row_idx) == 0 || is.na(row_idx)) return(default)
val <- df[row_idx, col]
if (is.null(val) || length(val) == 0 || all(is.na(val))) return(default)
return(as.numeric(val[1]))
}
# 5. 构建结果(确保每列长度严格为2)
res <- tibble::tibble(
Source = c("Between Groups", "Within Groups (Residuals)"),
Df = c(
as.integer(safe_extract(aov_summary, model_row_idx, "Df", NA_integer_)),
as.integer(safe_extract(aov_summary, residual_row_idx, "Df", NA_integer_))
),
Sum_Sq = c(
round(safe_extract(aov_summary, model_row_idx, "Sum Sq"), 3),
round(safe_extract(aov_summary, residual_row_idx, "Sum Sq"), 3)
),
Mean_Sq = c(
round(safe_extract(aov_summary, model_row_idx, "Mean Sq"), 3),
round(safe_extract(aov_summary, residual_row_idx, "Mean Sq"), 3)
),
F_Value = c(
round(safe_extract(aov_summary, model_row_idx, "F value"), 3),
NA_real_
),
p_Value = c(
round(safe_extract(aov_summary, model_row_idx, "Pr(>F)"), 4),
NA_real_
),
Sig = c(
sig_stars(safe_extract(aov_summary, model_row_idx, "Pr(>F)")),
""
)
)
return(res)
}, error = function(e) {
message(paste("ANOVA执行失败:", e$message))
return(tibble::tibble(
Source = c("Between Groups", "Within Groups (Residuals)"),
Df = NA_integer_, Sum_Sq = NA_real_, Mean_Sq = NA_real_,
F_Value = NA_real_, p_Value = NA_real_, Sig = ""
))
})
}
# ------------------------------
# 辅助函数5:Kruskal-Wallis检验
# ------------------------------
#' @export
run_kw <- function(valid_data, var, group) {
x <- valid_data[[var]]
g <- valid_data[[group]]
tryCatch({
kw_model <- stats::kruskal.test(x ~ g, data = valid_data)
# 整理KW结果表格
res <- tibble::tibble(
Test = "Kruskal-Wallis",
Chi_Sq = round(as.numeric(kw_model$statistic), 3),
Df = as.integer(kw_model$parameter),
p_Value = round(as.numeric(kw_model$p.value), 4),
Sig = sig_stars(kw_model$p.value)
)
return(res)
}, error = function(e) {
message(paste("Kruskal-Wallis执行失败:", e$message))
return(tibble::tibble(
Test = "Kruskal-Wallis",
Chi_Sq = NA_real_, Df = NA_integer_,
p_Value = NA_real_, Sig = ""
))
})
}
# ------------------------------
# Summary方法
# ------------------------------
#' @export
summary.mda <- function(object, dec = 3, ...) {
# 1. 基础信息
# 0. 基础信息
cat("Multigroup Difference Analysis (ANOVA/KW)\n")
cat("Data :", object$df_name, "\n")
cat("Dependent var:", object$var, "(numeric)\n")
......@@ -295,7 +404,7 @@ summary.mda <- function(object, dec = 3, ...) {
cat("Normality test:", object$normality_type, "\n")
cat("Valid samples:", object$valid_n, "\n\n")
# 2. 正态性检验结果
# 1. 正态性检验结果
cat("=== 1. Normality Test Results ===\n")
if (nrow(object$norm_res) == 0) {
cat(" No valid normality test results.\n\n")
......@@ -316,7 +425,7 @@ summary.mda <- function(object, dec = 3, ...) {
cat("\n")
}
# 3. 方差齐性检验结果
# 2. 方差齐性检验结果
cat("=== 2. Homogeneity of Variance Results ===\n")
if (nrow(object$homo_res) == 0) {
cat(" No valid homogeneity test results.\n\n")
......@@ -337,12 +446,55 @@ summary.mda <- function(object, dec = 3, ...) {
cat("\n")
}
# 4. 结论提示
# 3. 结论提示
cat("=== 3. Interpretation Tips ===\n")
cat("• 正态性:p ≥ 0.05 → 满足正态性假设\n")
cat("• 方差齐性:p ≥ 0.05 → 满足方差齐性假设\n")
cat("• 若同时满足这两个假设 → 使用方差分析(ANOVA)\n")
cat("• 若任一假设不满足 → 使用Kruskal-Wallis检验\n")
cat("• 若任一假设不满足 → 使用Kruskal-Wallis检验\n\n")
# 4.结果提示
if (is.null(object$core_test_res) || nrow(object$core_test_res) == 0) {
cat(" No valid core test result.\n\n")
} else {
# 根据选择的方法,格式化并显示结果
if (object$test_method == "anova") {
cat("=== 4. Analysis of Variance (ANOVA)===\n")
# 格式化ANOVA结果
core_formatted <- object$core_test_res %>%
dplyr::mutate(
p_Value = dplyr::case_when(
is.na(p_Value) ~ "",
p_Value < 0.001 ~ "<0.001",
p_Value < 0.01 ~ as.character(round(p_Value, 3)),
TRUE ~ as.character(round(p_Value, 4))
),
F_Value = ifelse(is.na(F_Value), "", as.character(F_Value)),
Df = as.character(Df),
Sum_Sq = as.character(Sum_Sq),
Mean_Sq = as.character(Mean_Sq)
) %>%
as.data.frame(stringsAsFactors = FALSE)
print(core_formatted, row.names = FALSE, right = FALSE)
} else if (object$test_method == "kw") {
cat("=== 4. Kruskal-Wallis Test===\n")
# 格式化KW结果
core_formatted <- object$core_test_res %>%
dplyr::mutate(
p_Value = dplyr::case_when(
is.na(p_Value) ~ "",
p_Value < 0.001 ~ "<0.001",
p_Value < 0.01 ~ as.character(round(p_Value, 3)),
TRUE ~ as.character(round(p_Value, 4))
),
Chi_Sq = as.character(Chi_Sq),
Df = as.character(Df)
) %>%
as.data.frame(stringsAsFactors = FALSE)
print(core_formatted, row.names = FALSE, right = FALSE)
}
cat("\n")
}
invisible(object)
}
......
......@@ -12,9 +12,13 @@ names(mda_plots) <- c(i18n$t("Normality: Q-Q Plot"),
i18n$t("Normality: Histogram"),
i18n$t("Homogeneity: Boxplot by Group"))
mda_test_methods <- c("anova", "kw")
names(mda_test_methods) <- c(i18n$t("Analysis of Variance (ANOVA)"),
i18n$t("Kruskal-Wallis Test"))
## 2. 函数形参
mda_args <- as.list(formals(mda))
mda_args <- mda_args[names(mda_args) %in% c("dataset", "var", "group", "normality_type", "data_filter")]
mda_args <- mda_args[names(mda_args) %in% c("dataset", "var", "group", "normality_type","test_method")]
## 3. 输入收集
mda_inputs <- reactive({
......@@ -25,12 +29,12 @@ mda_inputs <- reactive({
var = input$mda_var,
group = input$mda_group,
normality_type = input$mda_normality_type,
data_filter = if (input$show_filter) input$data_filter else "None",
test_method = input$mda_test_method,
envir = r_data
)
# 校验参数完整性
for (arg in names(mda_args)) {
if (is.null(inputs[[arg]])) inputs[[arg]] <- mda_args[[arg]]
if (is.null(inputs[[arg]]) || inputs[[arg]] == "") inputs[[arg]] <- mda_args[[arg]]
}
inputs
})
......@@ -114,7 +118,14 @@ output$ui_mda <- renderUI({
condition = "input.tabs_mda == 'Summary'",
uiOutput("ui_mda_var"),
uiOutput("ui_mda_group"),
uiOutput("ui_mda_normality_type")
uiOutput("ui_mda_normality_type"),
radioButtons(
inputId = "mda_test_method",
label = i18n$t("Select test method:"),
choices = mda_test_methods,
selected = state_single("mda_test_method", mda_test_methods, "anova"),
inline = FALSE
)
),
# Plot标签页
conditionalPanel(
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment