简介：SmokingGeneExpression 是一份面向生物信息学与统计分析学习者的 R 语言项目资料，围绕吸烟对健康肺组织基因表达的短期因果效应展开。内容将因果推断思路落到基因表达场景，涉及工具变量、倾向得分匹配、有向无环图等混杂控制方法，以及 RNA-seq 差异表达分析、PCA 可视化与生物学解释，适合希望理解观察性数据因果分析流程并借助 R 复现实战的读者。压缩包共 45 个文件，约 3.1MB，以 png 图表、txt 结果说明、xlsx 与 csv 数据表为主，另含 R 脚本和 pdf 文档，覆盖分析代码、数据样例与结果呈现等用途。已有 148 人学习，可作为课程作业、课题预研或方法对照的参考。读者可从中获得较完整的分析素材，按目录梳理研究逻辑、复现关键步骤并积累 R 语言实践。

1. 项目定位：用 R 把吸烟对肺组织基因表达的因果效应算清楚

一个做肺组织转录组的同行跟我吐槽过：差异表达分析跑完，吸烟组和从不吸烟组之间挂着三千多个显著基因，图表也漂亮，结果报告交上去，第一个问题就是“你怎么知道这些差异不是吸烟者普遍年龄更大、BMI 更高造成的”。相关性分析能告诉你两组人不一样，但不能告诉你是吸烟把某个基因的表达抬上去的。SmokingGeneExpression 这个 R 项目要解决的正是这一层：把“吸烟”当作处理变量，把“某个基因的表达量”当作结局，用因果推断的估计量把混杂因素压下去，给出每个基因在健康肺组织里对吸烟暴露的短期因果效应估计。

适合谁来读：手上有队列或公开表达资源、会写基本 R 脚本、但还没系统跑过倾向评分加权和双重稳健估计的人。项目的技术栈是纯 R——数据层用 data.table 和基础的 Bioconductor 表达处理，估计层主要靠 WeightIt、grf、sandwich 这几个包，可视化走 ggplot2。它不训练模型、不做深度学习，整套流程跑完就是一张带置信区间和 q 值的基因效应表，外加几张诊断图。对做生信但被审稿意见追着问因果的人，这个骨架的参考价值比多跑十个差异分析大。

2. 工程结构：R 脚本怎么分工、依赖怎么装、第一次怎么跑通

这个项目不是包（package），是一套按分析阶段切开的脚本工程。这意味着它没有R CMD check那套约束，好处是改起来自由，坏处是环境不一致时排错特别费劲——这是后面避坑章节的重要伏笔。

2.1 依赖清单：哪些包是必需的，哪些是可选的

依赖分成三档。第一档是数据处理，data.table、dplyr、tidyr、readr；第二档是因果估计核心，WeightIt（倾向评分与权重）、grf（因果森林与 AIPW 的机器学习 nuisance 估计）、sandwich加lmtest（聚类稳健标准误）；第三档是多重检验与出图，qvalue、ggplot2、patchwork。并行那块用future+furrr，规模化跑全基因组时会开关它。

环境锁定用renv，项目根目录有一份renv.lock。恢复环境的标准动作：

# 在项目根目录执行，renv 会按 lock 文件把包版本拉回一致状态 Rscript -e 'if (!requireNamespace("renv", quietly = TRUE)) install.packages("renv"); renv::restore()'

这条命令的重点不是“安装依赖”，而是“安装指定版本的依赖”。grf和WeightIt这两个包在版本升级时函数签名改过，尤其是WeightIt的estimand参数和grf的honesty默认值，跨版本跑出来的数会对不上。我第一次拿到这套脚本时偷懒没跑 restore，直接用系统里现装的包，结果 IPW 和 AIPW 的效应估计差了将近 8%，排查了两天才发现是权重归一化方式变了。血泪经验：这类因果估计工程，包版本必须锁死。

如果机器上装不了qvalue（Bioconductor 包），退路是用p.adjust(method = "BH")，两者在大多数情况下结论一致，只是 q 值的具体数值略有差异，报告里注明用的是哪一种即可。

2.2 目录分层：data、R、scripts、output 各管什么

工程目录大致是这样切的：

这个分层是有讲究的。R/里的函数不碰全局变量，参数从外面传进去，方便单独做单元测试；scripts/里的脚本只负责“读什么、调哪个函数、写到哪里”。很多生信分析脚本写着写着就变成一个两千行的巨型文件，中间变量到处飞，改一个参数要重跑全流程，这种结构就是为了防这个。常见的做法是每个scripts/开头统一 source 一遍函数层：

# scripts/01_prepare_data.R 的开头，固定写这一块 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 统一从项目根目录解析路径，避免依赖 setwd proj_root <- here::here() source(file.path(proj_root, "R", "io.R")) source(file.path(proj_root, "R", "preprocess.R")) source(file.path(proj_root, "R", "estimators.R")) # 全局参数集中在这一处，后面所有脚本读同一个配置 cfg <- list( expr_file = file.path(proj_root, "data/raw/expression_matrix.tsv"), pheno_file = file.path(proj_root, "data/raw/phenotype.tsv"), covar_file = file.path(proj_root, "data/raw/covariates.tsv"), min_expr = 1, # 表达过滤阈值，log2 尺度下取 1 大致对应低表达噪声层 min_samples = 0.3, # 一个基因至少在 30% 样本中超过阈值才保留 seed = 20240101 )

cfg这个列表是整条流水线的黑匣子开关。min_expr和min_samples直接决定后面有多少基因进入多重检验，改动了结论就得重出；seed因为因果森林内部有随机抽样，不固定的话两次跑出来的效应估计在小数点后几位会飘。路径用here::here()而不是相对路径，是为了让脚本从 RStudio、命令行、还是调度系统提交都能找到文件，这个小习惯能省掉大量“在我机器上能跑”的扯皮。

2.3 第一次跑通：从 01 到 05 的执行顺序

五个脚本是一条串行链，但中间结果都落盘，所以可以从任意一步重跑：

# 全流程，按编号顺序执行 Rscript scripts/01_prepare_data.R Rscript scripts/02_fit_ipw.R Rscript scripts/03_fit_aipw.R Rscript scripts/04_causal_forest.R Rscript scripts/05_report.R

01做数据对齐和表达预处理，输出data/derived/analysis_set.rds；02和03分别在同一个分析集上跑 IPW 和 AIPW，各自产出一张基因级效应表；04只在候选基因子集上跑因果森林，用来观察异质性；05把三张表按基因 ID 合并、统一多重检验校正、出图和写报告。

第一次跑建议先在01里把基因数砍到几百个做冒烟测试，确认整条链没断再放开全量。基因数上万的表达矩阵走完全套，单核大概要几十分钟到几小时，取决于样本量和是否开并行。冒烟测试这一步别省，04的因果森林在小样本上最容易因为某个基因方差为零直接抛错，早发现早改。

3. 数据层拆解：暴露、结局与协变量怎么对齐到同一张表

因果推断里最容易翻车的从来不是估计量写错，而是数据对齐错了。这个项目的核心数据结构只有一张宽表：每一行是一个样本，列包括样本 ID、吸烟暴露、几十上百个协变量、以及每个基因的表达值。所有估计函数都只认这张表，保证口径一致。

3.1 吸烟暴露变量的编码：二值化之前的三个判断

原始表型里吸烟状态通常是分好几档的：从不、曾经、当前、偶尔、已戒若干年，甚至还有每天支数。直接二值化会丢信息，但项目走的是“当前吸烟 vs 从不吸烟”这条二值路线，理由有三：短期因果效应下，“曾经吸烟”人群的暴露已经停止，把它塞进处理组会稀释效应；每天支数在多数队列里缺失率太高，做连续暴露需要另一套剂量-反应模型；而二值暴露能直接套用倾向评分框架，诊断工具也成熟。

编码这一步要把“曾经吸烟”和缺失显式剔出去，而不是当成对照组：

# 把原始吸烟编码映射成 0/1 的处理变量，并显式记录剔除原因 build_exposure <- function(pheno) { stopifnot(all(c("sample_id", "smoke_status") %in% names(pheno))) out <- pheno %>% mutate( smoke_status = tolower(trimws(smoke_status)), # 只保留两种明确状态，其余一律置 NA W = case_when( smoke_status %in% c("current", "current smoker", "daily") ~ 1L, smoke_status %in% c("never", "never smoker") ~ 0L, TRUE ~ NA_integer_ ), # 单独打个标记，方便统计各剔除原因的数量 drop_reason = case_when( is.na(W) & smoke_status %in% c("former", "past") ~ "former_smoker", is.na(W) & is.na(smoke_status) ~ "missing_status", is.na(W) & smoke_status == "" ~ "blank_status", is.na(W) ~ "other_status", TRUE ~ NA_character_ ) ) attr(out, "exposure_counts") <- table(out$W, useNA = "always") attr(out, "drop_counts") <- table(out$drop_reason, useNA = "no") out }

case_when的匹配顺序很重要，%in%里放了多种写法是为了兼容不同来源的表型表，真实数据里 “Current Smoker” 和 “current” 混着出现是常态。drop_reason这一列看起来多余，但它是后面写方法学部分的证据——审稿人一定会问你为什么最终样本是 N，你得能答出“剔除了多少既往吸烟者、多少状态缺失”。attr挂上去的计数会被01脚本打印到日志里，养成这个习惯，能省掉回头重算的麻烦。

还有一层容易忽略的：暴露的时间窗。项目说的是“短期”，那就得确认表型采集时间和组织取样时间接近，如果表型是问卷回顾性的“过去一年是否吸烟”，而样本是五年前存的，这个短期假设就不成立。数据里如果有采样日期字段，务必在预处理阶段做个时间差统计。

3.2 表达矩阵的预处理：过滤、标准化、以及要不要去批次

表达矩阵进模型前要过三道关。第一道是低表达过滤，用min_expr和min_samples两个参数控制，目的是减少检验数量、降低多重校正的惩罚。第二道是分布变换，如果是原始 count 就走log2(CPM + 1)或vst，如果上游已经给了 log2 尺度就跳过。第三道是批次，这一步是分水岭：批次如果和吸烟状态强相关，去批次会把效应一起减掉；如果批次独立于暴露，去不去影响都不大。

保守做法是先检查批次与暴露的关联，再做决定：

# 过滤 + 标准化，并在返回前检查批次是否与暴露混杂 prep_expression <- function(expr, meta, cfg) { mat <- as.matrix(expr[, -1, drop = FALSE]) rownames(mat) <- expr[[1]] # 1) 低表达过滤：两个条件同时满足才保留 keep <- rowMeans(mat > cfg$min_expr) >= cfg$min_samples mat <- mat[keep, , drop = FALSE] # 2) log2 变换，已经在 log 尺度上的矩阵跳过 if (max(mat, na.rm = TRUE) > 50) mat <- log2(mat + 1) # 3) 分位数标准化，抹平样本间整体分布差异 mat <- limma::normalizeBetweenArrays(mat, method = "quantile") # 4) 批次与暴露的关联检查：卡方检验，p 小说明批次和吸烟纠缠在一起 if (!is.null(meta$batch) && !all(is.na(meta$batch))) { tab <- table(meta$batch, meta$W) p_batch <- suppressWarnings(chisq.test(tab)$p.value) message(sprintf("批次-暴露关联检验 p = %.4g", p_batch)) if (!is.na(p_batch) && p_batch < 0.05) { warning("批次与吸烟状态显著关联，去批次可能吸收掉部分处理效应，慎用 ComBat") } } list(expr = mat, n_gene = nrow(mat), n_sample = ncol(mat)) }

rowMeans(mat > cfg$min_expr)这一行是过滤的核心表达式，返回每个基因超过阈值的样本比例，再和min_samples比。normalizeBetweenArrays走 quantile 方法，对偏态分布的表达值比中位数标准化更稳。最值得看的是第 4 段——批次卡方检验。如果 p 值很小，说明这批数据里“哪台机器测的”和“这个人吸不吸烟”高度相关，此时上 ComBat 或limma::removeBatchEffect，等于把处理效应的一部分当成批次削掉了，后面估计出来的效应会系统性偏小。这种情况我一般不去批次，而是把批次当作协变量放进倾向评分模型，让加权去平衡它，而不是从数据里删掉它。

矩阵的样本顺序和表型的样本顺序必须严格一致，这是最容易出错的地方。prep_expression返回前应该强制按meta$sample_id重排，或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次断言，别指望两边碰巧一致。

3.3 协变量集与设计矩阵共线性检查

协变量集的选择原则是：只放暴露之前就已经确定的变量（年龄、性别、遗传背景主成分、测序深度、批次），不要放可能被暴露影响的中介变量（比如某些炎症指标、肺功能）。放中介进去会把总效应部分遮掉，估出来的是直接效应，解释起来是另一回事。

确定变量集之后，必须检查设计矩阵的秩。协变量之间完全共线（比如两个高度相关的技术指标，或者某个分类变量某一水平只有一个人）会让倾向评分模型直接崩掉，glm可能不报错但给出荒谬的系数。

# 构建协变量设计矩阵并做秩与稀疏水平检查 check_design <- function(meta, covars) { f <- reformulate(covars) mm <- model.matrix(f, data = meta) # 秩亏说明存在完全共线，glm 的结果不可信 r <- qr(mm)$rank message(sprintf("设计矩阵维度 %d x %d，秩 = %d", nrow(mm), ncol(mm), r)) if (r < ncol(mm)) warning("设计矩阵秩亏，存在完全共线协变量，需手动剔除") # 分类变量稀疏水平：某个水平样本过少会让该水平的倾向评分极端 for (v in covars) { if (is.character(meta[[v]]) || is.factor(meta[[v]])) { tb <- table(meta[[v]]) if (any(tb < 20)) warning(sprintf("变量 %s 存在样本量 < 20 的水平：%s", v, paste(names(tb)[tb < 20], collapse = ", "))) } } invisible(mm) }

qr(mm)$rank小于列数就是秩亏，最常见的成因是分类变量做了独热编码后某一列全是零，或者两个连续变量做了线性组合。tb < 20这个阈值不是硬标准，20 只是个提醒线，真正要看的是该水平在两组里的分布——如果某个批次只有 3 个人且全在吸烟组，倾向评分在这个批次上会趋近 1，加权时必然出极端权重，这就是第 5 章要展开的问题。提前在这里发现，比跑到加权阶段再看权重分布图省事得多。

4. 因果估计核心：三种估计量在单个基因上的完整链路

数据对齐之后，估计层的输入就是一个向量：结局Y（某个基因的表达值）、处理W（0/1）、协变量矩阵X。项目给三种估计量，不是让你三选一，而是让它们互相印证：IPW 看加权后的可比人群，AIPW 用两个模型互相兜底，因果森林看效应是不是在亚群里不一致。三者结论方向一致，才敢往报告里写。

4.1 倾向评分与 IPW：把不可比的人群“造”成可比

倾向评分e(X) = P(W=1|X)的直觉是：给每个样本算一个“被吸烟的概率”，然后用这个概率的倒数做权重，让协变量分布在两组之间拉平。吸烟概率高的人如果实际没吸烟，权重就大，代表他“替”了很多吸烟者。

实现上用WeightIt，它比手写glm多了一层诊断：

# 单个基因的 IPW 估计，输入是结局向量 + 处理向量 + 协变量数据框 estimate_ipw <- function(y, w, x, trim = 0.01) { df <- data.frame(y = y, W = w, x) # 用 logistic 回归估倾向评分，estimand = "ATE" 表示估的是全人群平均效应 wt <- WeightIt::weightit(W ~ ., data = df, method = "glm", estimand = "ATE", stabilize = TRUE) ps <- wt$ps # 重叠性诊断：倾向评分贴边说明两组人根本不可比 overlap <- c(min = min(ps), max = max(ps), n_lt01 = sum(ps < 0.01), n_gt99 = sum(ps > 0.99)) # 权重截尾，压掉极端权重的杠杆 wts <- wt$weights if (!is.null(trim) && trim > 0) { q <- quantile(wts, probs = c(trim, 1 - trim)) wts <- pmin(pmax(wts, q[1]), q[2]) } # 加权最小二乘给 ATE，robust = TRUE 拿稳健标准误 fit <- lm(y ~ W, data = df, weights = wts) se <- sqrt(sandwich::vcovHC(fit, type = "HC1")["W", "W"]) list(est = coef(fit)[["W"]], se = se, ps = ps, overlap = overlap, ess = sum(wts)^2 / sum(wts^2)) }

estimand = "ATE"和"ATT"的区别要拎清：ATE 估的是“如果全人群都吸烟 vs 全人群都不吸烟”的平均差异，ATT 只估“吸烟者如果没吸烟会怎样”。做暴露危害评估一般用 ATE。stabilize = TRUE把权重除以处理组的边际概率，能显著降低权重的方差。overlap里那两个n_lt01/n_gt99计数是诊断的核心——如果有一批人倾向评分大于 0.99，说明模型认为他们“几乎肯定吸烟”，而那批人恰好没吸烟，权重会飙到 100 以上，一个样本就能主导整个估计。trim = 0.01是在 1% 和 99% 分位处截尾，代价是引入轻微偏差，但换来方差可控，实践中这笔交易通常划算。

ess（有效样本量）是另一个必看指标，sum(w)^2 / sum(w^2)给出加权后相当于多少“等效样本”。如果原始 500 例加权后 ESS 只剩 60，说明权重分布极其不均，这个估计不该信。我的经验阈值是 ESS 低于原始的 30% 就要回头查协变量集。

4.2 AIPW 双重稳健：两个模型互相兜底

IPW 只依赖倾向评分模型正确，一旦模型设定错了（比如该放交互项没放），估计就有偏。AIPW 加了结局模型，两个模型只要有一个对，估计就一致，这叫双重稳健。代价是计算量大一些，但在这个场景里值得。

结局模型可以用简单的线性回归，也可以用grf的回归森林。项目里默认走线性，因为基因表达和协变量的关系大多近似线性，森林在小样本上反而方差大：

# AIPW 估计：需要倾向评分 e(X) 和两个结局预测 mu1(X), mu0(X) estimate_aipw <- function(y, w, x, use_forest = FALSE) { df <- data.frame(y = y, W = w, x) # 倾向评分模型 ps_fit <- glm(W ~ ., data = df, family = binomial()) e <- fitted(ps_fit) # 夹一下，避免 e 正好等于 0 或 1 导致除零 e <- pmin(pmax(e, 1e-6), 1 - 1e-6) # 结局模型：分两组各拟合一个，预测反事实结局 if (use_forest) { cf <- grf::regression_forest(cbind(df$x, W = df$W), df$y) mu1 <- predict(cf, cbind(df$x, W = 1))$predictions mu0 <- predict(cf, cbind(df$x, W = 0))$predictions } else { m1 <- lm(y ~ ., data = subset(df, W == 1)) m0 <- lm(y ~ ., data = subset(df, W == 0)) mu1 <- predict(m1, newdata = df) mu0 <- predict(m0, newdata = df) } # AIPW 影响函数逐样本取值 psi <- (mu1 - mu0) + w * (y - mu1) / e - (1 - w) * (y - mu0) / (1 - e) est <- mean(psi) se <- sd(psi) / sqrt(length(psi)) list(est = est, se = se, mu1 = mean(mu1), mu0 = mean(mu0), # 结局模型在两组上的外推幅度，过说明模型在硬撑 extrap = c(max_mu1 = max(mu1), min_mu1 = min(mu1), max_mu0 = max(mu0), min_mu0 = min(mu0))) }

psi那四行是整套方法的核心。前两项mu1 - mu0是结局模型的预测效应，后面两个减项是“模型残差乘以逆倾向评分”的修正项。直觉上：如果结局模型预测得准，修正项接近零，估计主要由mu1 - mu0决定；如果倾向评分模型更准，修正项会把结局模型的偏差拉回来。

e夹到[1e-6, 1-1e-6]是必须的，现实数据里fitted(ps_fit)出现 1.0000 是常见的，除法直接炸出 Inf。extrap那几个极值是我加的诊断，看mu1的最大值有没有远超实际观测到的表达量范围——如果超过了，说明结局模型在做荒谬的外推，这条估计不能采信，先回头检查协变量。

use_forest = TRUE那条分支用grf::regression_forest，把W当成一个普通特征塞进X，预测时分别把W置 1 和 0。这种写法比训练两个独立森林更高效，但对样本量要求更高，小于 200 例时慎用。

4.3 因果森林：平均效应之外的异质性

IPW 和 AIPW 给的是全人群平均效应，但吸烟对基因表达的影响很可能在亚群里不一致——比如年龄大的人反应更强，或者只在不吸烟的遗传背景子群里显现。因果森林直接估每个样本的条件平均处理效应tau(X)，然后看tau有没有系统性的分化。

# 因果森林：估条件平均处理效应，看异质性 fit_causal_forest <- function(y, w, x, num_trees = 2000) { X <- cbind(as.data.frame(x)) cf <- grf::causal_forest( X = X, Y = y, W = w, num.trees = num_trees, honesty = TRUE, # 诚实划分，避免用同一批数据既分裂又估值 min.node.size = 10, # 叶子最小样本，太小会过拟合 sample.fraction = 0.5, # 子采样比例 seed = 20240101 ) # 全人群平均效应，target.sample = "overlap" 只在有重叠的人群上平均 ate <- grf::average_treatment_effect(cf, target.sample = "overlap") # 变量重要性，看哪些协变量在驱动异质性 vimp <- grf::variable_importance(cf) list(cf = cf, ate = ate, vimp = vimp) }

honesty = TRUE是grf的默认值，作用是样本分裂和效应估计用不同的子样本，防止过拟合把噪声当成异质性。这个参数不要关掉，关了之后tau的分布会变得极其“好看”但全是假的。min.node.size = 10是叶子节点的最小样本量，基因表达数据的噪声大，设太小会切出一堆只有几个样本的叶子，估出来的tau方差爆炸。

target.sample = "overlap"这一项值得单独说。默认是"all"，在全员上平均；改成"overlap"只在倾向评分有重叠的区域上平均，更稳健但也更保守。做暴露效应的基因筛查时用"overlap"，结论不容易被少数极端个体带跑。

variable_importance的输出是一个矩阵，第一列是原始重要性，第二列是归一化到 0-1 的版本。看前五名就够——如果年龄或某个技术协变量排第一，说明异质性主要来自技术因素而非生物学，这时候要警惕，别急着往机制上解释。

在基因层面，因果森林只在候选基因子集上跑，不是全基因组。原因很实际：2000 棵树乘以几百个基因，单机跑一晚上都不一定完，而且大多数基因的异质性检验并不显著，跑全量是浪费算力。

4.4 从 p 值到 q 值：结果表的组装

三种估计跑完，每个基因会拿到三套est和se。组装结果表要解决两个问题：一是统一效应方向（基因表达上调还是下调），二是多重检验。基因数量上千，直接拿 p < 0.05 筛出来的基本都是假阳性。

# 合并三种估计量的结果，做 BH 校正，输出可读的结果表 assemble_results <- function(ipw_list, aipw_list, cf_list = NULL) { build <- function(lst, tag) { data.frame( gene_id = names(lst), est = vapply(lst, `[[`, numeric(1), "est"), se = vapply(lst, `[[`, numeric(1), "se"), stringsAsFactors = FALSE ) %>% setNames(c("gene_id", paste0(tag, c("_est", "_se")))) } res <- Reduce(function(a, b) merge(a, b, by = "gene_id"), list(build(ipw_list, "ipw"), build(aipw_list, "aipw"))) # 用 AIPW 作为主估计，它双重稳健 res$z <- res$aipw_est / res$aipw_se res$p <- 2 * pnorm(-abs(res$z)) res$q <- p.adjust(res$p, method = "BH") # 方向标记，方便画火山图 res$dir <- ifelse(res$q < 0.05 & res$aipw_est > 0, "up", ifelse(res$q < 0.05 & res$aipw_est < 0, "down", "ns")) # 一致性旗标：两种方法方向相反时标出来，这类基因不写进结论 res$flag_inconsistent <- sign(res$ipw_est) != sign(res$aipw_est) res[order(res$q), ] }

p.adjust(method = "BH")是 Benjamini-Hochberg，控制 FDR，比 Bonferroni 温和得多，在基因表达场景是标准做法。dir那一列把显著基因分成上调和下调，画图直接映射颜色。

flag_inconsistent这一列是经验产物。IPW 和 AIPW 的方向如果相反，说明两个模型的假设冲突得很厉害，可能是倾向评分模型设定有问题，也可能是结局模型的线性假设在这个基因上不成立。这类基因即使 q 值很小也不该写进结论，最多放进补充材料，标注方法间不一致，等着用敏感性分析去解释。

主估计选 AIPW 而不是 IPW，理由是双重稳健：只要倾向评分和结局模型有一个对，估计就一致，比单靠一个模型的 IPW 更扛得住设定误差。

5. 避坑清单：因果推断流水线里最容易翻车的五处

这套流程跑通不难，难的是跑出来的数能不能信。下面五条是我在这类项目上反复踩到的，每一条都按“现象 → 原因 → 解决”整理，配上排查代码。

5.1 权重爆炸与重叠性缺失

现象：IPW 估出来的效应比 AIPW 大好几倍，标准误也大得离谱，ESS只剩原始样本的 10% 不到。

原因：倾向评分模型太“自信”，给某些样本赋了接近 0 或 1 的概率。常见于某个协变量几乎能完全预测处理状态——比如某一批次的所有样本都来自吸烟组，模型在这个批次上学到“批次 = 吸烟”，权重自然爆掉。另一个成因是分类变量的稀疏水平，某一类只有几个人，逻辑回归在该维度上给出极端的系数。

解决：先看重叠性诊断，再决定是修剪样本还是截尾权重。

# 权重诊断：分布、ESS、以及极端权重来自哪些样本 diagnose_weights <- function(wts, meta) { ess <- sum(wts)^2 / sum(wts^2) cat(sprintf("有效样本量 ESS = %.1f / %d (%.1f%%)\n", ess, length(wts), 100 * ess / length(wts))) # 找出权重最大的 10 个样本，看看它们有什么共同点 idx <- order(wts, decreasing = TRUE)[1:min(10, length(wts))] print(cbind(meta[idx, , drop = FALSE], weight = wts[idx])) invisible(ess) }

打印极端权重样本的协变量，是定位问题最快的路径。如果这十个样本全在同一个批次，那就是批次导致的；如果全是高龄，那是年龄分布尾部的问题。ESS低于 30% 时，我一般先尝试trim = 0.01，还不行就改用重叠权重（WeightIt里method = "glm", estimand = "ATO"），重叠权重对极端倾向评分不敏感，代价是估计目标变成“重叠人群的效应”，解释时要写清楚。

5.2 批次、共线性与“静默”的样本错位

现象：脚本一声不吭跑完，结果表也有，但效应估计整体偏小，或者某个已知的强效应基因死活不显著。

原因：三种可能。批次与暴露强相关，去批次时把处理效应一起减掉了；协变量共线导致glm系数不稳定，倾向评分本身就不准；表达矩阵和表型的样本顺序没对齐，merge或cbind之后列和行的对应关系已经乱了，但 R 不会报错。

解决：三个检查必须都在预处理阶段做完——批次卡方检验（3.2 节）、设计矩阵秩检查（3.3 节）、样本 ID 断言。

# 样本对齐的最后一道保险，放在所有合并操作之后 assert_alignment <- function(expr, meta) { stopifnot( "表达矩阵列名缺失" = !is.null(colnames(expr)), "表型样本 ID 缺失" = "sample_id" %in% names(meta), "样本数不一致" = ncol(expr) == nrow(meta), "样本顺序错位" = identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }

identical这一条是硬断言，不做容错，顺序不一致直接抛错。我宁可让脚本在第一步就炸，也不想在写了两周报告之后发现样本对错了。

5.3 因果森林在小样本上的“看起来很准”

现象：因果森林跑出来的tau分布特别漂亮，亚群间差异明显，但换一个随机种子结果就变了。

原因：honesty = FALSE或者min.node.size设得太小，树在噪声上过拟合。基因表达数据本身噪声大，样本量几百例时，2000 棵树的平均还压不住方差，tau的分布就是随机的。

解决：三个参数一起调——保持honesty = TRUE，min.node.size不低于 10，sample.fraction取 0.4 到 0.6 之间；然后用两个不同种子各跑一遍，看tau的相关系数。

# 稳定性检验：换种子重跑，看条件效应是否可复现 check_cf_stability <- function(y, w, x, seeds = c(1, 2024)) { taus <- lapply(seeds, function(s) { cf <- grf::causal_forest(X = as.data.frame(x), Y = y, W = w, num.trees = 2000, honesty = TRUE, min.node.size = 10, seed = s) predict(cf)$predictions }) r <- cor(taus[[1]], taus[[2]]) cat(sprintf("跨种子 tau 相关系数 = %.3f\n", r)) if (r < 0.7) warning("条件效应不稳定，异质性结论不可靠") invisible(r) }

相关系数低于 0.7 就别写异质性结论了，老老实实报平均效应。这个检验成本不高，值得每个候选基因都过一遍。

5.4 多重检验的两种误用

现象：q 值算出来了，显著基因上千个，但和已知的吸烟相关基因对不上。

原因：要么用了 Bonferroni（在基因表达场景过于严苛，真正的信号被压掉），要么在筛选和验证之间反复调阈值，最后报的是“调出来的”显著性。还有一种隐蔽的：IPW 和 AIPW 各算一遍 q 值，挑显著的那个报，这是变相的多重比较。

解决：主分析只跑一套（AIPW），q 值只算一次，阈值预先定死（比如 q < 0.05），敏感性分析的结果单独列，不参与筛选。方法学部分要写清楚“主估计用的是哪一个”。

5.5 用同一份数据调参又评估

现象：所有诊断都完美，重叠性好、ESS 高、q 值分布漂亮，但换一个独立队列就复现不出来。

原因：倾向评分模型的协变量集、trim阈值、因果森林的参数，如果在看到结果之后反复调整，等于把结局信息泄漏进了模型选择。这类偏差不会体现在任何一张诊断图上。

解决：把样本按 7:3 随机切开，前 70% 用来定协变量集和超参数，后 30% 只跑一次，不看结果不回头。如果样本量不够切，就在分析计划里把参数全部预先写死，中途改了要记录在日志里。

# 预注册式的参数冻结：把最终参数写进文件，后续脚本只读不猜 freeze_params <- function(cfg, path = "output/tables/frozen_params.rds") { saveRDS(cfg, path) message("参数已冻结至 ", path, "，后续分析请勿修改") }

这条看着像流程洁癖，但它是唯一能防住“调参调到显著”的办法。玄学的地方在于，你自己调的时候完全不觉得在作弊，每一步都有“合理”的理由。

6. 二次开发：把单基因流程改成批量筛选与可复现报告

拿到这套脚本的人，十有八九第一件事是想跑全基因组。单基因函数写好了，批量跑不是简单套个for循环——直接循环会撞上三个问题：内存、算力和失败恢复。我给的做法是把估计函数包一层，加上并行和缓存。

先看批量层怎么写：

# 批量估计：并行 + 单基因失败不影响整体 + 结果落盘缓存 run_batch <- function(expr_mat, w, x, method = c("ipw", "aipw"), n_workers = 4) { method <- match.arg(method) future::plan(future::multisession, workers = n_workers) on.exit(future::plan(future::sequential), add = TRUE) gene_ids <- rownames(expr_mat) est_fun <- if (method == "ipw") estimate_ipw else estimate_aipw res <- furrr::future_map(gene_ids, function(g) { y <- as.numeric(expr_mat[g, ]) # 零方差基因直接跳过，否则 glm 会抛错并中断整个批次 if (stats::var(y, na.rm = TRUE) < 1e-8) return(NULL) tryCatch( est_fun(y = y, w = w, x = x), error = function(e) { structure(list(err = conditionMessage(e)), class = "fail") } ) }, .options = furrr::furrr_options(seed = TRUE)) names(res) <- gene_ids fails <- vapply(res, inherits, logical(1), "fail") if (any(fails)) message(sprintf("%d 个基因估计失败，已跳过", sum(fails))) res[!fails] }

furrr_options(seed = TRUE)保证每个 worker 的随机数流可复现，不加这一句，因果森林那部分每次跑出来都不一样。tryCatch把单基因的失败隔离掉，glm在完全分离的数据上会抛 warning 甚至 error，一个基因炸掉不该拖垮整批。零方差判断放在最前面，是最常见的崩溃点——某些基因在所有样本里表达值都一样，lm直接返回 NA 系数。n_workers不建议超过物理核心数减一，grf内部本身也会占内存，开太多反而会因为内存交换变慢。

第二个改造点是报告层。原始脚本输出的是 CSV，二次开发可以把它变成参数化的 Rmd 或 Quarto 文档，把 q 值分布、火山图、前若干个基因的森林图一次性生成，换一批数据只改cfg里的路径。这里的关键是把“结果表”和“图”解耦：图只读表，不重算估计。

第三个改造点是加一层敏感性分析。做法很简单：把trim从 0.01 换成 0.05 和 0.1，各跑一遍 AIPW，看显著基因的重叠比例。

这张表建议直接贴进补充材料，比任何文字描述都有说服力。审稿人问“你的结论稳健吗”，把重叠率甩过去就行。

最后一个技巧是日志。把sessionInfo()、cfg、每个脚本的运行时间和输入文件的行列数统一写进output/logs/，一个脚本一个文件。这个习惯我是被坑出来的：有次报告里的数和三个月前跑的对不上，翻遍代码也没改过，最后发现是上游给的表达矩阵悄悄更新了一版，样本数从 512 变成了 508。从那以后我每次跑因果估计流水线，都强制在第一步把sessionInfo()和输入文件的 md5 一起写进日志，跑完再核对一遍才敢往报告里放数。这类工程的可靠性不在估计量写得多漂亮，而在每一步都能回溯。希望帮到你。

本文还有配套的精品资源，点击获取