超全的卡方分析流程！

最新推荐文章于 2025-01-08 23:53:57 发布

爱做科研的桶

最新推荐文章于 2025-01-08 23:53:57 发布

阅读量1.5k

点赞数 35

分类专栏： R语言数据分析文章标签： python 开发语言数据分析 r语言科技

本文链接：https://blog.youkuaiyun.com/llthxx/article/details/144066524

版权

R语言数据分析专栏收录该内容

8 篇文章

订阅专栏

卡方分析介绍：

卡方分析（Chi-Square Analysis）是一种用于检验两个或更多分类变量之间关系的统计方法，广泛应用于社会科学、市场研究、生物学等领域。它主要用于判断观察到的数据与期望数据之间的差异是否显著，从而帮助我们理解变量之间的关系或独立性。

1. 卡方分析的基本概念

卡方分析是基于卡方分布（Chi-Square Distribution）的一种假设检验方法，常用于处理类别数据（也叫定性数据或名义数据）。卡方分析的基本思路是，通过比较观察到的频数（实际数据）和期望的频数（在零假设下的理论数据）之间的差异，来判断变量之间是否存在显著关系。

关键点：

观察频数（Observed Frequency）：数据实际发生的频数。
期望频数（Expected Frequency）：如果零假设成立的情况下，每个类别的预期频数。

2. 卡方分析的常见类型

卡方分析主要有两种形式：

a. 卡方独立性检验（Chi-Square Test of Independence）

用于检验两个分类变量之间是否存在关联或独立性。通常，卡方独立性检验用于交叉列联表（Contingency Table）中，分析不同类别变量是否独立无关。

例子：调查性别和是否吸烟之间是否存在关联。例如，性别（男性、女性）和吸烟习惯（吸烟、不吸烟）是否有关系。

b. 卡方拟合优度检验（Chi-Square Goodness-of-Fit Test）

用于检验一个单一分类变量的观察频数是否符合预期的分布。它比较观察到的样本频率与某一已知分布（如均匀分布、正态分布等）之间的差异。

例子：检查骰子投掷的结果是否均匀，即各个点数（1-6）的出现频率是否相等。

3. 卡方检验的应用

卡方独立性检验：常用于市场研究、社会学研究、医学研究等，检查不同类别之间是否相关。例如，分析性别和购买意图是否相关。
卡方拟合优度检验：常用于生物学实验、物理学研究等，检验实际数据是否符合理论模型或分布。例如，验证某种遗传模式是否符合预期的基因分布。

4. 卡方检验的注意事项

样本量要求：卡方检验对样本量有一定要求，尤其是期望频数不能过小（通常要求期望频数大于5）。如果样本量太小，可以考虑合并类别或者使用其他检验方法（如Fisher精确检验）。
数据类型：卡方检验仅适用于分类数据（名义数据或顺序数据），对于连续变量需要进行其他检验（如t检验、ANOVA等）。
独立性假设：卡方检验假设每个观察值是独立的。如果数据存在依赖关系（如重复测量、配对数据），卡方检验可能不适用。

5. 卡方分析的限制

期望频数问题：如果期望频数太小（小于5），卡方检验可能不可靠。在这种情况下，使用Fisher精确检验等方法更为合适。
无法衡量关联强度：卡方检验只能判断变量之间是否有显著关系，但无法反映关联的强度或方向。可以通过计算Cramér's V等其他指标来衡量关联的强弱。
假设检验的限制：卡方检验的假设通常是简单的独立性或拟合优度检验，对于更复杂的模型可能需要其他统计方法（如多变量分析）。

代码：

第一部分：设置环境

print("-------绘图-------")
setwd("C:/Users/Public/李立统/学术学习/论文工程/研究/GPT研究/文章思路/分析思路/选择比例分析1")
getwd()

closeAllConnections()
rm(list = ls())

设置工作目录：
- setwd()：指定工作目录，所有的文件操作都会在这个目录下进行。
- getwd()：确认当前工作目录。
清理环境：
- closeAllConnections()：关闭所有打开的文件连接。
- rm(list = ls())：删除当前工作环境中的所有对象，确保代码运行时环境干净。

第二部分：加载必要的包

library(openxlsx)
library(gtsummary)
library(dplyr)
library(vcd)
library(purrr)
library(ggstatsplot)
library(ggplot2)
library(jmv)

加载分析和可视化所需的包：
- openxlsx：读取 Excel 文件。
- gtsummary 和 dplyr：数据处理和统计分析。
- vcd：关联分析和统计量计算。
- purrr：函数式编程（用于批量操作）。
- ggstatsplot 和 ggplot2：生成统计图形。
- jmv：交叉表分析。

第三部分：数据导入与预处理

mydata <- read.xlsx("gpt.xlsx")
mydata <- subset(mydata, mydata$final != 3)

数据读取：
- read.xlsx("gpt.xlsx")：读取 Excel 文件 gpt.xlsx。
数据过滤：
- subset(mydata, mydata$final != 3)：保留 final 列中值不为 3 的数据。

mydata <- mydata %>%
  mutate(
    gpt = factor(gpt, levels = c("gpt3.5", "gpt4")),
    system = factor(system, levels = c("system2", "system1")),
    query = factor(query, levels = c("query1", "query2")),
    final_new = factor(final, labels = c("不参加", "参加")),
    sort_new = factor(sort, labels = paste0(seq(10, 90, by = 10), "%"))
  )

数据类型转换和新变量创建：
- 使用 mutate 对列重新编码：
  - gpt、system 和 query 被转为有序因子，指定排序规则。
  - final_new 和 sort_new 是新的因子变量，赋予更直观的标签。

第四部分：数据子集生成

mydata1 <- mydata %>% filter(system == "system1", query == "query1")
mydata2 <- mydata %>% filter(system == "system1", query == "query2")
mydata3 <- mydata %>% filter(system == "system2", query == "query1")
mydata4 <- mydata %>% filter(system == "system2", query == "query2")

基于 system 和 query 列的值筛选数据，生成 4 个子集（每个子集对应一种 system-query 组合）。

第五部分：计算频数表和行百分比

compute_table_and_props <- function(data, row_var, col_var) {
  freq_table <- table(data[[row_var]], data[[col_var]])
  row_props <- prop.table(freq_table, margin = 1)
  return(list(freq_table = freq_table, row_props = row_props))
}

函数compute_table_and_props：
- 输入数据和行、列变量。
- 生成二维频数表（freq_table）和按行百分比（row_props）。

mydata1_analysis <- compute_table_and_props(mydata1, "gpt", "final_new")

调用函数计算各子集的频数表和行百分比。

第六部分：卡方检验和关联统计量

run_chi_square_analysis <- function(freq_table) {
  chi_test <- chisq.test(freq_table, correct = FALSE)
  assoc_stats <- assocstats(freq_table)
  return(list(chi_test = chi_test, assoc_stats = assoc_stats))
}

函数run_chi_square_analysis：
- 进行卡方检验（chisq.test）。
- 使用 assocstats 计算关联统计量（如 Cramer V 系数）。

第七部分：将列联表转换为数据框

table_to_dataframe <- function(freq_table) {
  as.data.frame(as.table(freq_table)) %>%
    rename(gpt = Var1, final_new = Var2, frequency = Freq)
}

table_to_dataframe 函数：
- 将频数表转换为易于处理的数据框。

第八部分：可视化分析

generate_ggbarstats <- function(data, title) {
  ggbarstats(
    data = data,
    x = gpt,
    y = final_new,
    counts = frequency,
    title = title,
    xlab = "Model",
    ylab = "Count",
    ggtheme = ggplot2::theme_minimal()
  )
}

generate_ggbarstats 函数：
- 用 ggbarstats 创建分组条形图，展示 gpt 和 final_new 之间的频数关系。

第九部分：保存图形

for (i in seq_along(plots)) {
  ggsave(
    filename = paste0("plot_", names(tables)[i], ".png"),
    plot = plots[[i]],
    width = 8,
    height = 6
  )
}

使用 ggsave 将每个图形保存为 PNG 文件。

总结：

环境设置与包加载： 初始化工作目录，清理工作环境，并加载所需的 R 包。
数据导入与预处理： 从 Excel 文件读取数据，进行预处理（包括因子变量转换）。
数据子集生成： 基于不同的 system 和 query 筛选数据，生成 4 个子集。
频数表与百分比计算： 使用 table 和 prop.table 计算频数表和按行百分比。
卡方检验与关联统计： 对每个子集进行卡方检验，并计算关联统计量（如 Cramer V 系数）。
转换为数据框： 将频数表转换为数据框，便于进一步处理和可视化。
计算交叉表相关系数： 使用 contTables 计算交叉表的相关系数。
图形生成与保存： 使用 ggbarstats 生成条形图，并将图形保存为 PNG 文件。

# ------- 设置环境 -------
print("-------绘图-------")
setwd("C:/Users/Public/李立统/学术学习/论文工程/研究/GPT研究/文章思路/分析思路/选择比例分析1")
getwd()

closeAllConnections()
rm(list = ls())

# 加载必要包
library(openxlsx)
library(gtsummary)
library(dplyr)
library(vcd)
library(purrr)
library(ggstatsplot)
library(ggplot2)
library(jmv)

# ------- 数据导入与预处理 -------
mydata <- read.xlsx("gpt.xlsx")
mydata <- subset(mydata, mydata$final != 3)

# 数据预处理：转换为因子并赋予适当的标签
mydata <- mydata %>%
  mutate(
    gpt = factor(gpt, levels = c("gpt3.5", "gpt4")),
    system = factor(system, levels = c("system2", "system1")),
    query = factor(query, levels = c("query1", "query2")),
    final_new = factor(final, labels = c("不参加", "参加")),
    sort_new = factor(sort, labels = paste0(seq(10, 90, by = 10), "%"))
  )

# ------- 生成数据子集 -------
# 根据系统和查询条件筛选数据子集
mydata1 <- mydata %>% filter(system == "system1", query == "query1")
mydata2 <- mydata %>% filter(system == "system1", query == "query2")
mydata3 <- mydata %>% filter(system == "system2", query == "query1")
mydata4 <- mydata %>% filter(system == "system2", query == "query2")

# ------- 计算频数表和行百分比 -------
# 定义函数：计算频数表和按行百分比
compute_table_and_props <- function(data, row_var, col_var) {
  freq_table <- table(data[[row_var]], data[[col_var]])
  row_props <- prop.table(freq_table, margin = 1)  # 按行计算百分比
  return(list(freq_table = freq_table, row_props = row_props))
}

# 计算每个子集的频数表和行百分比
mydata1_analysis <- compute_table_and_props(mydata1, "gpt", "final_new")
mydata2_analysis <- compute_table_and_props(mydata2, "gpt", "final_new")
mydata3_analysis <- compute_table_and_props(mydata3, "gpt", "final_new")
mydata4_analysis <- compute_table_and_props(mydata4, "gpt", "final_new")

# 打印每个子集的行百分比
print("---- mydata1 的行百分比 ----")
print(mydata1_analysis$row_props)

print("---- mydata2 的行百分比 ----")
print(mydata2_analysis$row_props)

print("---- mydata3 的行百分比 ----")
print(mydata3_analysis$row_props)

print("---- mydata4 的行百分比 ----")
print(mydata4_analysis$row_props)

# ------- 执行卡方检验和计算关联统计量 -------
# 定义函数：执行卡方检验并计算关联统计量
run_chi_square_analysis <- function(freq_table) {
  chi_test <- chisq.test(freq_table, correct = FALSE)  # 不进行连续性校正
  
  # 输出期望频数和卡方检验结果
  print("期望频数:")
  print(chi_test$expected)
  print("卡方检验结果:")
  print(chi_test)
  
  # 计算 Cramer V 系数等关联统计量
  assoc_stats <- assocstats(freq_table)
  print("关联统计量 (Cramer V等):")
  print(assoc_stats)
  
  return(list(chi_test = chi_test, assoc_stats = assoc_stats))
}

# 对每个数据子集执行卡方检验
print("---- mydata1 的卡方检验 ----")
result1 <- run_chi_square_analysis(mydata1_analysis$freq_table)

print("---- mydata2 的卡方检验 ----")
result2 <- run_chi_square_analysis(mydata2_analysis$freq_table)

print("---- mydata3 的卡方检验 ----")
result3 <- run_chi_square_analysis(mydata3_analysis$freq_table)

print("---- mydata4 的卡方检验 ----")
result4 <- run_chi_square_analysis(mydata4_analysis$freq_table)

# ------- 将列联表转换为数据框 -------
# 定义函数：将频数表转换为数据框
table_to_dataframe <- function(freq_table) {
  as.data.frame(as.table(freq_table)) %>%
    rename(gpt = Var1, final_new = Var2, frequency = Freq)
}

# 将每个频数表转换为数据框
mydata1_table <- table_to_dataframe(mydata1_analysis$freq_table)
mydata2_table <- table_to_dataframe(mydata2_analysis$freq_table)
mydata3_table <- table_to_dataframe(mydata3_analysis$freq_table)
mydata4_table <- table_to_dataframe(mydata4_analysis$freq_table)

# 打印每个数据框
print("---- mydata1_table ----")
print(mydata1_table)

print("---- mydata2_table ----")
print(mydata2_table)

print("---- mydata3_table ----")
print(mydata3_table)

print("---- mydata4_table ----")
print(mydata4_table)

# ------- 计算交叉表相关系数 -------
contTables(formula = frequency ~ gpt:final_new, 
           data = mydata1_table, contCoef = TRUE, phiCra = TRUE)
contTables(formula = frequency ~ gpt:final_new, 
           data = mydata2_table, contCoef = TRUE, phiCra = TRUE)
contTables(formula = frequency ~ gpt:final_new, 
           data = mydata3_table, contCoef = TRUE, phiCra = TRUE)
contTables(formula = frequency ~ gpt:final_new, 
           data = mydata4_table, contCoef = TRUE, phiCra = TRUE)

# ------- 可视化：生成 ggbarstats 图形 -------
# 定义函数：生成 ggbarstats 图形
generate_ggbarstats <- function(data, title) {
  ggbarstats(
    data = data,
    x = gpt,
    y = final_new,
    counts = frequency,
    title = title,
    xlab = "Model",
    ylab = "Count",
    ggtheme = ggplot2::theme_minimal()
  )
}

# 创建一个包含所有数据框及标题的列表
tables <- list(
  mydata1_table = mydata1_table,
  mydata2_table = mydata2_table,
  mydata3_table = mydata3_table,
  mydata4_table = mydata4_table
)
titles <- c("System1-Query1", "System1-Query2", "System2-Query1", "System2-Query2")

# 批量生成图形并保存
plots <- purrr::map2(tables, titles, ~ generate_ggbarstats(.x, .y))

# 打印每个图形
for (plot in plots) {
  print(plot)
}

# 可选：将图形保存为文件
for (i in seq_along(plots)) {
  ggsave(
    filename = paste0("plot_", names(tables)[i], ".png"),
    plot = plots[[i]],
    width = 8,
    height = 6
  )
}