W的自留地

W的自留地

基于R的GCA代码使用说明

技术
260
2025-04-27

作为一个人工智能，我无法直接控制您的浏览器或自动替换编辑器中的文本，但我已经为您对文章进行了全面的润色和排版优化。

本次润色主要提升了技术文档的专业性、阅读的流畅度，并规范了标点符号和列表格式。您可以直接点击代码块右上角的“复制”按钮，然后全选并覆盖您编辑器中的原有内容。

Markdown

# 基于 R 语言的 GCA 代码使用指南

本文档旨在详细介绍如何使用 R 语言进行增长曲线分析（Growth Curve Analysis, GCA）。

## 1. 环境准备

1. **R 版本**：建议使用 R ≥ 4.0.0。
2. **依赖 R 包**：
   - `reshape2`：用于数据从宽格式到长格式的转换。
   - `lme4`、`lmerTest`：用于混合效应模型的构建与拟合。
   - `ggplot2`、`RColorBrewer`、`ggthemes`：用于数据可视化与绘图。
   - `openxlsx`：用于将分析结果导出至 Excel。
   - `dplyr`：用于数据清洗与整理（可选）。

   **安装方式**（若尚未安装，请执行以下代码）：
   ```R
   install.packages(c("reshape2", "lme4", "lmerTest", "ggplot2", "RColorBrewer", "ggthemes", "openxlsx", "dplyr"))

  1. 文件编码：请务必确保 R 脚本和数据文件均为 UTF-8 编码，以避免中文字符乱码。

2. 数据格式说明

  • 文件名：建议为 data.txt（也可自定义，需在脚本中对应修改）。

  • 分隔符：制表符（Tab）。

  • 数据结构：

    1. fayinren：发音人标识（因子变量）。

    2. diaolei：调类标识（因子变量，仅包含两类，例如 Qping、Cping）。

    3. lizi：例字标识（因子变量）。

    4. f0_1 … f0_20：表示 20 个不同时间点的基频（F0）数据（数值型）。

数据示例（第一行为列名）：

fayinren

diaolei

lizi

f0_1

……

f0_20

S1

Qping

A

200

……

219

S2

Cping

B

185

……

166

……

3. 脚本概要

运行脚本后，程序将依次执行以下步骤：

  1. 加载并检查数据：读取数据并初始化环境。

  2. 宽格式转长格式：使用 reshape2::melt() 函数将 F0 的各列展开为行数据。

  3. 添加时间变量：通过 poly() 函数构建二阶正交多项式时间变量 ot1 和 ot2。

  4. 模型拟合：构建并拟合以下四个模型：

    • m0：仅包含固定效应 (ot1+ot2)*diaolei。

    • m1：m0 + 发音人随机截距 (1|fayinren)。

    • m2：m0 + 例字随机截距 (1|lizi)。

    • m3：m0 + 发音人与例字交叉随机截距 (1|fayinren)+(1|lizi)。

  5. 模型比较：计算各模型的 AIC（赤池信息量准则）值，自动选择最优模型。

  6. 结果导出：将最佳模型的固定效应和随机效应参数输出至 output/ 目录下的 Excel 文件中，包含：

    • Summary 工作表：最佳模型名称、模型公式及 AIC 值。

    • Fixed Effects 工作表：固定效应系数表。

    • Random Effects 工作表（若适用）：随机效应的方差与标准差。

4. 使用步骤

  1. 将下方的 R 语言源码保存为 gca_analysis.R 文件。

  2. 将准备好的数据文件放入脚本所在的同一目录，或在脚本中修改 read.table() 的文件路径。

  3. 在 R 或 RStudio 中运行该脚本：

    R

    source("gca_analysis.R")

  4. 运行结束后，在 R 控制台中查看打印出的 AIC 比较结果及最佳模型名称。

  5. 前往当前工作目录下的 output/ 文件夹，获取生成的 Excel 结果文件（如：m*_GCA_params.xlsx），按需打开查看详细参数。

  6. 进阶：若需自定义绘图色板、主题或调整多项式的阶数，可直接修改脚本中的相应参数。

5. 结果解读

  • AIC 值：该数值越小，说明模型拟合效果越好。

  • 固定效应参数：关注 diaoleiCqing、ot1:diaoleiCping、ot2:diaoleiCping 等变量的估计值及 p 值。它们分别对应了发音的均值、斜率与拱度三个维度的差异（具体变量名称可能根据您的数据稍有不同）。

  • 随机效应参数：展示 fayinren（发音人）和 lizi（例字）的组内方差，反映了个体差异。

6. R 语言源码

R

# ================================
# GCA 分析：多模型比较 & 参数导出
# ================================

# 0. 环境准备 -------------------------------------------------
# 安装并加载所需包
packages <- c("reshape2", "lme4", "lmerTest", "ggplot2", "RColorBrewer", "ggthemes", "openxlsx", "dplyr")
invisible(lapply(packages, library, character.only = TRUE))

# 读取工作目录下 data.txt（可按需修改文件路径）
data <- read.table("data.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

# 1. 宽转长 & 时间变量处理 -----------------------------------
# 把 f0_1...f0_20 转换为长格式
data_long <- melt(data, id.vars = c("fayinren", "diaolei", "lizi"), 
                  variable.name = "time", value.name = "f0")

# 提取数值型 time
data_long$time <- as.numeric(sub("f0_", "", data_long$time))

# 添加二阶正交多项式
uniq_t <- sort(unique(data_long$time))
poly_mat <- poly(uniq_t, 2)
poly_df <- data.frame(time = uniq_t, ot1 = poly_mat[, 1], ot2 = poly_mat[, 2])
data_long <- merge(data_long, poly_df, by = "time", all.x = TRUE)

# 强制转换为因子类型
data_long$diaolei <- factor(data_long$diaolei)
data_long$fayinren <- factor(data_long$fayinren)
data_long$lizi <- factor(data_long$lizi)

# 2. 模型构建 & AIC 比较 -------------------------------------
# m0: 纯固定效应
m0 <- lm(f0 ~ (ot1 + ot2) * diaolei, data = data_long)

# m1: 发音人随机截距
m1 <- lmer(f0 ~ (ot1 + ot2) * diaolei + (1 | fayinren), data = data_long, REML = FALSE)

# m2: 例字随机截距
m2 <- lmer(f0 ~ (ot1 + ot2) * diaolei + (1 | lizi), data = data_long, REML = FALSE)

# m3: 发音人 + 例字 随机截距
m3 <- lmer(f0 ~ (ot1 + ot2) * diaolei + (1 | fayinren) + (1 | lizi), data = data_long, REML = FALSE)

# 比较 AIC 并输出最佳模型
model_list <- list(m0 = m0, m1 = m1, m2 = m2, m3 = m3)
aic_vals <- sapply(model_list, AIC)
print(round(aic_vals, 2))

best_name <- names(which.min(aic_vals))
cat("最佳模型：", best_name, "\n")
best_mod <- model_list[[best_name]]

# 3. 导出模型参数到 Excel -----------------------------------
# 确保 output 文件夹存在
if (!dir.exists("output")) dir.create("output")

# 新建工作簿
wb <- createWorkbook()
addWorksheet(wb, "Summary")
addWorksheet(wb, "Fixed Effects")

# 汇总说明写入
summary_txt <- paste0("最佳模型：", best_name, "\n",
                      "公式：", deparse(formula(best_mod)), "\n",
                      "AIC 值：", round(aic_vals[best_name], 2), "\n")
writeData(wb, "Summary", summary_txt)

# 固定效应参数表写入
if (inherits(best_mod, "lm")) {
  fixef_tab <- summary(best_mod)$coefficients
} else {
  fixef_tab <- coef(summary(best_mod))
}

fixef_df <- as.data.frame(fixef_tab)
fixef_df <- cbind(Parameter = rownames(fixef_df), fixef_df)
rownames(fixef_df) <- NULL
writeData(wb, "Fixed Effects", fixef_df)

# 随机效应参数表写入（若存在）
if (!inherits(best_mod, "lm")) {
  addWorksheet(wb, "Random Effects")
  rand_df <- as.data.frame(VarCorr(best_mod))
  rand_df <- rand_df %>% select(grp, vcov, sdcor)
  names(rand_df) <- c("Group", "Variance", "Std.Dev")
  writeData(wb, "Random Effects", rand_df)
}

# 保存 Excel 文件
out_file <- file.path("output", paste0(best_name, "_GCA_params.xlsx"))
saveWorkbook(wb, out_file, overwrite = TRUE)
cat("模型参数已成功输出至：", out_file, "\n")

#R 2 #实验语音学 10
嘎裂音合成脚本使用说明