组均值及 95% CI 趋势图

# 组均值及 95% CI 趋势图：三个随访时间点的描述性趋势
# 项目原创脚本，MIT 许可；模拟数据 CC0-1.0；图形 CC BY 4.0。
# 在全新 R 会话运行。每次运行会更新下方 output_dir 中的本示例文件。
# 图形标签使用英文以避免不同电脑缺少中文字体；可在参数和 labs() 中修改。

# 1. 依赖：首次使用时手动运行以下安装命令，脚本不会自动安装。
# install.packages(c("ggplot2", "jsonlite"))
required_packages <- c("ggplot2", "jsonlite")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages)) stop("Please install: ", paste(missing_packages, collapse = ", "))
library(ggplot2)

# 2. 可修改参数：颜色、字体、导出尺寸和路径。输出路径相对于当前工作目录。
set.seed(20261014)
colors <- c("#007F96", "#B66A46")
font_family <- "sans"
figure_width <- 8
figure_height <- 5.5
output_dir <- file.path(getwd(), "rfig-time-trend-output")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
theme_set(theme_minimal(base_size = 12, base_family = font_family) +
 theme(plot.background = element_rect(fill = "white", color = NA), panel.grid.minor = element_blank(),
       plot.title = element_text(face = "bold", size = 15), legend.position = "bottom",
       plot.margin = margin(14, 18, 14, 14)))

# 3. 示例数据：每行一位受试者在一个时间点的测量；每人三行。 保留 id、group、week、value 长表结构；真实失访需说明各时间点分母和处理方式。
# 使用自己的 CSV 时，将以下数据生成段替换为：
# data <- read.csv("your-data.csv", stringsAsFactors = FALSE)
# 再按本例要求设置分类顺序、检查字段、单位及缺失值。
subjects <- data.frame(id = 1:50, group = rep(c("Control", "Treatment"), each = 25), random_intercept = rnorm(50, 0, 2))
data <- merge(subjects, data.frame(week = c(0, 4, 12)))
data$value <- 15 + data$random_intercept - .05 * data$week - .18 * data$week * (data$group == "Treatment") + rnorm(nrow(data), 0, 1.5)
data$random_intercept <- NULL
data$group <- factor(data$group, levels = c("Control", "Treatment"))
stopifnot(!anyNA(data))

# 4. 统计计算与绘图。每组 25 人，每人测量 3 次；每个组×时间点单独计算均值及 t 95% CI。各时间点内不同受试者独立，但跨时间点相关；这些描述性区间不能替代纵向模型或组×时间交互检验。演示无失访。
cells <- split(data, interaction(data$group, data$week, drop = TRUE))
summary_data <- do.call(rbind, lapply(cells, function(d) {
 n <- nrow(d); m <- mean(d$value); margin <- qt(.975, n-1) * sd(d$value) / sqrt(n)
 data.frame(group = d$group[1], week = d$week[1], n = n, mean = m, lower = m-margin, upper = m+margin)
}))
p <- ggplot(summary_data, aes(week, mean, group = group, color = group)) +
 geom_line(aes(linetype = group), linewidth = .8) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .4) +
 geom_point(aes(shape = group), size = 3) + scale_color_manual(values = colors) + scale_x_continuous(breaks = c(0,4,12)) +
 labs(x = "Follow-up (weeks)", y = "Mean measurement (95% CI)", color = NULL, shape = NULL, linetype = NULL, title = "Descriptive trends across follow-up")
metrics <- list(cell_n = summary_data$n, mean = summary_data$mean, lower = summary_data$lower, upper = summary_data$upper)
print(p)

# 5. 导出：CSV 与绘图所用 data 完全一致，preview 与科研图片均来自同一个 p。
write.csv(data, file.path(output_dir, "data.csv"), row.names = FALSE, fileEncoding = "UTF-8")
ggsave(file.path(output_dir, "preview.png"), p, width = figure_width, height = figure_height, dpi = 160, bg = "white")
ggsave(file.path(output_dir, "figure.png"), p, width = figure_width, height = figure_height, dpi = 300, bg = "white")
ggsave(file.path(output_dir, "figure.pdf"), p, width = figure_width, height = figure_height, device = grDevices::pdf, bg = "white")
ggsave(file.path(output_dir, "figure.tiff"), p, width = figure_width, height = figure_height, dpi = 300, compression = "lzw", bg = "white")
jsonlite::write_json(metrics, file.path(output_dir, "metrics.json"), auto_unbox = TRUE, digits = 16, pretty = TRUE)