基于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：用于数据整理（可选）。

   安装方式（若尚未安装）：

   install.packages(c(
     "reshape2", "lme4", "lmerTest", "ggplot2",
     "RColorBrewer", "ggthemes", "openxlsx", "dplyr"
   ))

3. 文件编码：请确保脚本和数据文件为 UTF-8 编码。

2.数据格式说明

• 文件名：data.txt（或自定义文件名）。

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

• 列结构：

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

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

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

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

示例（第一行为列名）：

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 中运行：

   source("gca_analysis.R")

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

5. 在当前目录下的 output/ 文件夹中，找到生成的 Excel 文件：

   m*_GCA_params.xlsx

   按需打开查看参数详情。

6. 若要自定义色板、主题或多项式阶数，可在脚本对应部分修改。

5.结果解读

• AIC 值：数值越小，模型拟合越优。

• 固定效应参数：diaoleiCqing、 ot1:diaoleiCping、ot2:diaoleiCping 的估计值及 p 值，对应均值、斜率与拱度三个维度。根据文件的变量命名，显示可能不同。

• 随机效应参数：展示 fayinren、lizi 的组内方差。

6.R语言源码

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

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

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

# 1. 宽转长 & 时间变量处理 -----------------------------------
# 把 f0_1...f0_20 转到 long 格式
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)
}

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