跳到正文
数据统计代码教程

可重复研究：R 项目结构与依赖管理

投稿返修时最难受的场景，是重跑三个月前的分析。包升级后 aggregate() 行为没变但绘图默认值变了，原始数据被同门「顺手」改过一格，上一版的中间结果不知道存在哪。这三件事都跟分析水平无关，全是工程问题。可重复研究（reproducible research）要达成的目标很朴素：任何人，包括半年后的你自己，拿到项目目录，能重跑出图表里同样的数字。

先把目录结构定下来

Section titled “先把目录结构定下来”

项目根目录长这样：

toothgrowth/
├── toothgrowth.Rproj
├── R/
│   ├── 01-clean.R
│   └── 02-analysis.R
├── data/
│   ├── raw/
│   └── processed/
├── output/
├── reports/
│   └── report.Rmd
├── renv.lock
└── .gitignore

名字都是约定俗成的，不要自创：R/ 放脚本和自定义函数，data/raw/ 放原始数据，data/processed/ 放清洗后的中间数据，output/ 放图和表的成品，reports/ 放 Rmd 源文件。别人接手时不用问「这个文件是干嘛的」。

两条硬规则。第一，data/raw/ 只读：原始数据一旦在 Excel 里被手工改过一格，所有下游结果就失去了追溯性，你没法说清某个异常值是仪器测出来的还是手滑改的。清洗逻辑必须写成代码，即使只是改一个列名。第二，data/processed/ 和 output/ 是产物，写进 .gitignore 不提交，因为它们能由代码重建。留下生成它们的脚本就够了：

.Rproj.user/
renv/library/
data/processed/
output/

脚本用数字前缀标明执行顺序（01-clean.R、02-analysis.R），再补一个入口把整条链路串起来：

R/run-all.R
source(here::here("R", "01-clean.R"))
source(here::here("R", "02-analysis.R"))

不要 setwd()

Section titled “不要 setwd()”

setwd("/Users/zhangsan/Desktop/analysis") 在写下的那一刻是有效的，换台机器、换个人、换个目录名都会失效，而且失效时往往报的是「文件不存在」，让人怀疑数据丢了。更麻烦的是它污染全局状态：你在控制台手动 setwd() 过一次，随后所有的相对路径都变成了相对那个目录，脚本和交互式会话的行为不再一致。

正确做法是用 RStudio 项目。File → New Project 创建的目录里会生成一个 .Rproj 文件，双击它打开 RStudio，工作目录自动设为项目根，不用写一行代码。

命令行下没有这个机制，Rscript 的工作目录是调用它的 shell 所在目录，所以脚本内部要用 here 包定位路径：

here::here("data", "processed", "growth.rds")
[1] "/home/user/projects/toothgrowth/data/processed/growth.rds"

here::here() 的定位方式是：从当前工作目录往上找 .Rproj、.git 之类的锚点文件，找到哪个就把哪一级当项目根。所以不管你是从 RStudio 渲染，还是从任意子目录用 Rscript 调用，结果都一致。这个包没有任何依赖，值得每个项目都装。

Rmd 渲染是路径错误的重灾区：knitr 默认把工作目录设成 Rmd 文件所在的目录，而不是项目根。Rmd 放在 reports/ 里时，read.csv("data/raw/x.csv") 必然失败。在 setup 块里统一一次就好：

# 放在 Rmd 最前面的 setup 块中
knitr::opts_knit$set(root.dir = here::here())

用 renv 锁定包版本

Section titled “用 renv 锁定包版本”

路径问题解决后，剩下最大的变量是包版本。同一句 filter() 在 dplyr 1.0 和 1.1 下的边缘行为可能不同，stringsAsFactors 的默认值在 R 4.0 前后就变过一次。记不住这些，就交给工具记。

renv::init()

renv::init() 做三件事：扫描代码里用到的包、在项目内建一个私有库（renv/library/，不污染系统库）、生成 renv.lock 记录每个包的确切版本和来源，同时在 .Rprofile 里加一行激活代码。之后的工作流只有三条命令：

renv::snapshot()   # 把当前用到的包版本写进 renv.lock
renv::status()     # 检查 lock 文件与当前库是否一致
renv::restore()    # 按 renv.lock 把包装回来

renv.lock 是纯文本 JSON，要提交到 Git；renv/library/ 不用提交，别人 restore() 一下就有了。

lock 文件的粒度值得留意。它列出的是项目库里每个包的确切版本，其中大部分是间接依赖：代码里可能只写了 library(dplyr)，lock 文件里却会连带出现 rlang、vctrs、tibble、pillar 等 dplyr 自己依赖的包。这和手写一份「只列直接依赖」的清单不是一回事，两种信息都有用，但用途不同：想知道代码里真正用到哪些包，用 renv::dependencies() 扫一遍，它按文件列出每一处 library() 和 :: 调用；想知道复现这个环境需要装什么，看 renv.lock。

{
  "R": {
    "Version": "4.5.1",
    "Repositories": [
      { "Name": "CRAN", "URL": "https://cloud.r-project.org" }
    ]
  },
  "Packages": {
    "knitr": {
      "Package": "knitr",
      "Version": "1.49",
      "Source": "Repository",
      "Repository": "CRAN"
    }
  }
}

这是节选，真实文件里会列出每个包及其依赖。要注意 renv 的能力边界：它锁的是 R 包的版本，锁不住 R 本体版本、系统库、外部的生信工具。所以论文的方法部分除了 renv.lock，还要写清楚 R 版本。

数据版本控制

Section titled “数据版本控制”

原始数据只读这件事，靠自觉不如靠工具。三条做法按项目规模选：

  • 数据量小（几十 MB 以内）：把 data/raw/ 一并提交到 Git。文件一旦入库就不要再改，需要修正就加一个新文件，在脚本里显式说明用了哪一版。
  • 数据量大：用 Git LFS 或者 DVC 管理大文件，仓库里只存指针。
  • 数据绝对不能出实验室：至少要记录校验和。tools::md5sum() 或者 digest::digest() 算一次原始文件的哈希，写进 README，事后能验证数据有没有被动过。

文件名带上日期也是有效的土办法：raw/2026-01-06-clinical.xlsx。看出问题的那一刻，你就知道该找哪一版。

有一类操作要绝对避免：手工编辑 data/processed/ 里的文件。一旦手工改过，你就有了两个事实来源——文件里的事实和脚本能重跑出来的事实，它们迟早会分叉，而且没有任何记录告诉你差在哪。

随机种子：固定的是整条流

Section titled “随机种子：固定的是整条流”

前面几节管的是路径、包和数据。还有一类变量不影响文件，只影响数字：随机数。抽样、聚类初始化、交叉验证划分、bootstrap 都从随机数流里取数，不固定种子，重跑一次结果就变。

set.seed(42)
run1 <- rnorm(5)
set.seed(42)
run2 <- rnorm(5)
print(round(run1, 4))
print(round(run2, 4))
print(identical(run1, run2))
[1]  1.3710 -0.5647  0.3631  0.6329  0.4043
[1]  1.3710 -0.5647  0.3631  0.6329  0.4043
[1] TRUE

set.seed() 设定的是 R 全局随机数流的状态，rnorm()、runif()、sample()、kmeans() 都从这一条流里取数。既然是「一条流」，它就是按调用顺序推进的：中途插进任何一次抽样，后面的数全部错位。

set.seed(42)
x <- rnorm(3)


set.seed(42)
noise <- runif(1)   # 中途插进来的一次抽样
y <- rnorm(3)


print(round(x, 4))
print(round(y, 4))
[1] 1.3710 -0.5647  0.3631
[1] 1.5307 0.9559 0.0479

两段代码开头都写了 set.seed(42)，x 和 y 却完全不同。那次 runif(1) 从流里取走一个数，y 拿到的是流的第二段。这个错误在真实项目里很隐蔽：调试时临时加一行 sample() 看数据，忘了删，此后所有需要复现的结果都变了，而代码本身看不出任何异常。

把种子放在调用处而不是函数体里，能避免函数被调用的顺序影响结果：

kmeans_fit <- function(x, k, seed) {
  set.seed(seed)
  kmeans(x, centers = k)
}


set.seed(1)
a <- kmeans_fit(iris[, 1:4], 3, 42)
set.seed(999)     # 调用前的全局状态被覆盖
b <- kmeans_fit(iris[, 1:4], 3, 42)
print(identical(a$cluster, b$cluster))
[1] TRUE

对照一个不在函数里设种子的版本：

kmeans_loose <- function(x, k) kmeans(x, centers = k)


set.seed(1)
c1 <- kmeans_loose(iris[, 1:4], 3)
set.seed(2)
c2 <- kmeans_loose(iris[, 1:4], 3)
print(identical(c1$cluster, c2$cluster))
[1] FALSE

同样的数据、同样的 k = 3，换个种子就得到不同分组。聚类结果对初始中心敏感，这不是实现缺陷，是算法本身的性质，所以种子和 k 一样属于必须记录的分析参数。

withr 包提供了更干净的写法，它在局部临时设种，退出时把全局流恢复原状：

library(withr)


set.seed(1)
a <- rnorm(2)
with_seed(99, rnorm(2))   # 只在这一行生效
b <- rnorm(2)


print(round(a, 4))
print(round(b, 4))
[1] 0.2140 0.4797
[1] -0.6265  0.1836

b 的取值和完全没有 with_seed() 那一行时一致，说明这次局部设种没有推进全局流。需要在一段代码里用不同种子做多次独立抽样时，with_seed() 比反复 set.seed() 更不容易出错。

sample() 同样依赖种子，划分训练集和测试集时要固定它：

set.seed(7)
print(sample(10, 3))
set.seed(7)
print(sample(10, 3))
[1] 10  3  7
[1] 10  3  7

方法部分应当写明固定种子的位置和数值，例如「用 set.seed(42) 配合 sample() 按 7:3 划分训练集」。种子数字本身没有意义，写出来的目的是让别人能复核你的划分方式。

记录环境

Section titled “记录环境”

每条结论都要能追溯到具体的软件版本，sessionInfo() 一行就能打出来：

sessionInfo()
R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 24.04.1 LTS


Matrix products: default


locale:
[1] LC_CTYPE=zh_CN.UTF-8       LC_NUMERIC=C
[3] LC_TIME=zh_CN.UTF-8        LC_COLLATE=zh_CN.UTF-8


attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

后面还会列出本次会话加载的全部包及其版本。写论文时把前几行抄进方法部分就够用：R 版本、平台、操作系统。审稿人问「你们用的什么软件」，这比事后回忆可靠得多。

更省事的做法是让报告自带环境快照：在 Rmd 末尾加一个块，把块选项设成 echo=FALSE，代码不显示、输出照常出现，每份报告末尾就自动带上一份当时的环境记录。

# 末尾的环境附录块，块选项写 {r, echo=FALSE}
sessionInfo()

想让快照成为可归档的产物，把它写进 output/，和图表放在一起。这样每次完整跑一遍分析，环境记录就自动刷新一次，不依赖事后回忆：

writeLines(capture.output(sessionInfo()), here::here("output", "session-info.txt"))

capture.output() 把本该打印到控制台的内容收成字符向量，writeLines() 落盘。会话里加载的包越多，文件越长，用 length(readLines()) 能看到行数。这个文件和 renv.lock 的分工是：renv.lock 声明「应该装哪些包」，session-info.txt 记录「这次实际用的是哪些包、什么版本、什么平台」。

脚本从干净状态跑起

Section titled “脚本从干净状态跑起”

R 会话里残留的对象同样会破坏可重复性。你在控制台定义过 df，脚本里少写一行读数据的代码照样能跑通，因为它用的是内存里那个 df。换台机器重跑，报错 object 'df' not found。

rm(list = ls())   # 清空全局环境里的对象

rm(list = ls()) 能清掉全局环境里的对象，但清不掉已经加载的包和改过的 options()，所以它只是近似。真正的检验是重开一个 R 会话从头执行。Rscript 每次启动一个全新的进程，没有任何残留状态：

终端窗口
Rscript -e 'cat(paste(search(), collapse = " | "), "\n")'
.GlobalEnv | package:stats | package:graphics | package:grDevices | package:utils | package:datasets | package:methods | Autoloads | package:base

搜索路径里只有默认挂载的基础包，前面会话里定义的任何变量都不存在。Rscript 在中途报错时返回非零退出码，可以直接挂进 make 或 CI。脚本改完先用 Rscript R/run-all.R 跑一遍，比在 RStudio 里点几次运行可靠。

R Markdown 渲染是同一回事。RStudio 的 Knit 按钮默认在当前会话里执行，如果 Rmd 依赖了控制台里事先定义的变量，本机渲染成功，同门拿到后必然失败。在渲染设置里勾选「在独立会话中运行」，或者改用命令行：

终端窗口
Rscript -e 'rmarkdown::render("reports/report.Rmd")'

命令行渲染不需要图形界面，也更容易在服务器或 CI 上复现。

一个能重跑的完整脚本

Section titled “一个能重跑的完整脚本”

把前面的做法落进 R/02-analysis.R。真实项目里第一行是从 data/processed/ 读入，这里用内置数据演示同一段逻辑：

R/02-analysis.R
dir.create(here::here("output"), showWarnings = FALSE)


data(PlantGrowth)   # 真实项目替换为 readRDS(here::here("data", "processed", "growth.rds"))


growth <- within(PlantGrowth, {
  group <- factor(group, levels = c("ctrl", "trt1", "trt2"))
})


means <- aggregate(weight ~ group, data = growth, FUN = mean)
sds   <- aggregate(weight ~ group, data = growth, FUN = sd)


means$mean_weight <- round(means$weight, 2)
means$sd_weight   <- round(sds$weight, 2)
means$weight      <- NULL


print(means)
write.csv(means, here::here("output", "growth-means.csv"), row.names = FALSE)
  group mean_weight sd_weight
1  ctrl        5.03      0.58
2  trt1        4.66      0.45
3  trt2        5.53      0.57

这段代码里没有一处绝对路径、没有 setwd()、没有依赖控制台里残留的变量。换台机器，renv::restore() 装回包，Rscript R/run-all.R 就能得到同一张表。PlantGrowth 是 R 内置的盆栽重量数据，三个处理组各 10 盆，weight 是干重；表里 ctrl 组 5.03、trt2 组 5.53，组间差异的检验放到 /r/modeling/anova/ 里讲。

几个反复出现的坑

Section titled “几个反复出现的坑”

用 .RData 当数据交换格式。 它保存的是 R 对象本身，跨大版本可能读不出来，别人用 Python 也打不开。中间结果用 .rds 或 CSV，交付数据一律用 CSV 或 Parquet。

把 renv::init() 当成一次性动作。 装完新包不 snapshot()，renv.lock 就停在原地，别人 restore() 出来的环境缺包。养成装完包顺手 renv::status() 的习惯。

论文里的数字是从控制台复制过去的。 只要有一处手工复制，改数据后正文和图表就不同步了。数字一律用行内代码从数据里取，Rmd 的做法见 /r/reporting/rmarkdown-basics/。

只提交了代码，没提交数据版本信息。 半年后你自己也说不清当时的输入是哪一版，结论就无法复核。

忘了固定随机种子。 抽样、聚类初始化、交叉验证划分都依赖随机数。脚本里不写 set.seed()，同一个脚本两次运行会给出不同的分组和不同的 p 值，而且不报任何错，你只会觉得「这批数据的结果不太稳定」。

把包装进全局库调试。 临时 install.packages() 到系统库或用户库，当时能跑，换台机器或换个账号就找不到那个包。项目用到的包装进 renv 的项目库，装完顺手 renv::snapshot() 更新 lock 文件。

同样的要求对 Python 项目一样成立：Notebook 要能从头跑到尾、依赖写进 requirements.txt 或 environment.yml、路径不写死。对照 /python/basics/jupyter/ 可以看另一套工具链下的同类做法。