簡(jiǎn)介SmokingGeneExpression 是一份面向生物信息學(xué)與統(tǒng)計(jì)分析學(xué)習(xí)者的 R 語(yǔ)言項(xiàng)目資料圍繞吸煙對(duì)健康肺組織基因表達(dá)的短期因果效應(yīng)展開。內(nèi)容將因果推斷思路落到基因表達(dá)場(chǎng)景涉及工具變量、傾向得分匹配、有向無(wú)環(huán)圖等混雜控制方法以及 RNA-seq 差異表達(dá)分析、PCA 可視化與生物學(xué)解釋適合希望理解觀察性數(shù)據(jù)因果分析流程并借助 R 復(fù)現(xiàn)實(shí)戰(zhàn)的讀者。壓縮包共 45 個(gè)文件約 3.1MB以 png 圖表、txt 結(jié)果說明、xlsx 與 csv 數(shù)據(jù)表為主另含 R 腳本和 pdf 文檔覆蓋分析代碼、數(shù)據(jù)樣例與結(jié)果呈現(xiàn)等用途。已有 148 人學(xué)習(xí)可作為課程作業(yè)、課題預(yù)研或方法對(duì)照的參考。讀者可從中獲得較完整的分析素材按目錄梳理研究邏輯、復(fù)現(xiàn)關(guān)鍵步驟并積累 R 語(yǔ)言實(shí)踐。1. 項(xiàng)目定位用 R 把吸煙對(duì)肺組織基因表達(dá)的因果效應(yīng)算清楚一個(gè)做肺組織轉(zhuǎn)錄組的同行跟我吐槽過差異表達(dá)分析跑完吸煙組和從不吸煙組之間掛著三千多個(gè)顯著基因圖表也漂亮結(jié)果報(bào)告交上去第一個(gè)問題就是“你怎么知道這些差異不是吸煙者普遍年齡更大、BMI 更高造成的”。相關(guān)性分析能告訴你兩組人不一樣但不能告訴你是吸煙把某個(gè)基因的表達(dá)抬上去的。SmokingGeneExpression 這個(gè) R 項(xiàng)目要解決的正是這一層把“吸煙”當(dāng)作處理變量把“某個(gè)基因的表達(dá)量”當(dāng)作結(jié)局用因果推斷的估計(jì)量把混雜因素壓下去給出每個(gè)基因在健康肺組織里對(duì)吸煙暴露的短期因果效應(yīng)估計(jì)。適合誰(shuí)來(lái)讀手上有隊(duì)列或公開表達(dá)資源、會(huì)寫基本 R 腳本、但還沒系統(tǒng)跑過傾向評(píng)分加權(quán)和雙重穩(wěn)健估計(jì)的人。項(xiàng)目的技術(shù)棧是純 R——數(shù)據(jù)層用 data.table 和基礎(chǔ)的 Bioconductor 表達(dá)處理估計(jì)層主要靠 WeightIt、grf、sandwich 這幾個(gè)包可視化走 ggplot2。它不訓(xùn)練模型、不做深度學(xué)習(xí)整套流程跑完就是一張帶置信區(qū)間和 q 值的基因效應(yīng)表外加幾張?jiān)\斷圖。對(duì)做生信但被審稿意見追著問因果的人這個(gè)骨架的參考價(jià)值比多跑十個(gè)差異分析大。2. 工程結(jié)構(gòu)R 腳本怎么分工、依賴怎么裝、第一次怎么跑通這個(gè)項(xiàng)目不是包package是一套按分析階段切開的腳本工程。這意味著它沒有R CMD check那套約束好處是改起來(lái)自由壞處是環(huán)境不一致時(shí)排錯(cuò)特別費(fèi)勁——這是后面避坑章節(jié)的重要伏筆。2.1 依賴清單哪些包是必需的哪些是可選的依賴分成三檔。第一檔是數(shù)據(jù)處理data.table、dplyr、tidyr、readr第二檔是因果估計(jì)核心WeightIt傾向評(píng)分與權(quán)重、grf因果森林與 AIPW 的機(jī)器學(xué)習(xí) nuisance 估計(jì)、sandwich加lmtest聚類穩(wěn)健標(biāo)準(zhǔn)誤第三檔是多重檢驗(yàn)與出圖qvalue、ggplot2、patchwork。并行那塊用futurefurrr規(guī)?；苋蚪M時(shí)會(huì)開關(guān)它。環(huán)境鎖定用renv項(xiàng)目根目錄有一份renv.lock?；謴?fù)環(huán)境的標(biāo)準(zhǔn)動(dòng)作# 在項(xiàng)目根目錄執(zhí)行renv 會(huì)按 lock 文件把包版本拉回一致狀態(tài) Rscript -e if (!requireNamespace(renv, quietly TRUE)) install.packages(renv); renv::restore()這條命令的重點(diǎn)不是“安裝依賴”而是“安裝指定版本的依賴”。grf和WeightIt這兩個(gè)包在版本升級(jí)時(shí)函數(shù)簽名改過尤其是WeightIt的estimand參數(shù)和grf的honesty默認(rèn)值跨版本跑出來(lái)的數(shù)會(huì)對(duì)不上。我第一次拿到這套腳本時(shí)偷懶沒跑 restore直接用系統(tǒng)里現(xiàn)裝的包結(jié)果 IPW 和 AIPW 的效應(yīng)估計(jì)差了將近 8%排查了兩天才發(fā)現(xiàn)是權(quán)重歸一化方式變了。血淚經(jīng)驗(yàn)這類因果估計(jì)工程包版本必須鎖死。如果機(jī)器上裝不了qvalueBioconductor 包退路是用p.adjust(method BH)兩者在大多數(shù)情況下結(jié)論一致只是 q 值的具體數(shù)值略有差異報(bào)告里注明用的是哪一種即可。2.2 目錄分層data、R、scripts、output 各管什么工程目錄大致是這樣切的路徑內(nèi)容說明data/raw/原始表達(dá)矩陣、表型表、協(xié)變量表只讀任何腳本不得寫回data/derived/標(biāo)準(zhǔn)化后的表達(dá)、構(gòu)建好的分析集中間產(chǎn)物可全部刪除重跑R/函數(shù)層io.R、preprocess.R、estimators.R、sensitivity.R、utils.R只定義函數(shù)不執(zhí)行分析scripts/執(zhí)行層01_prepare_data.R到05_report.R有編號(hào)按順序跑output/tables/基因效應(yīng)表、診斷統(tǒng)計(jì)表交付物output/figures/森林圖、權(quán)重分布圖、QQ 圖交付物這個(gè)分層是有講究的。R/里的函數(shù)不碰全局變量參數(shù)從外面?zhèn)鬟M(jìn)去方便單獨(dú)做單元測(cè)試scripts/里的腳本只負(fù)責(zé)“讀什么、調(diào)哪個(gè)函數(shù)、寫到哪里”。很多生信分析腳本寫著寫著就變成一個(gè)兩千行的巨型文件中間變量到處飛改一個(gè)參數(shù)要重跑全流程這種結(jié)構(gòu)就是為了防這個(gè)。常見的做法是每個(gè)scripts/開頭統(tǒng)一 source 一遍函數(shù)層# scripts/01_prepare_data.R 的開頭固定寫這一塊 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 統(tǒng)一從項(xiàng)目根目錄解析路徑避免依賴 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)) # 全局參數(shù)集中在這一處后面所有腳本讀同一個(gè)配置 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, # 表達(dá)過濾閾值log2 尺度下取 1 大致對(duì)應(yīng)低表達(dá)噪聲層 min_samples 0.3, # 一個(gè)基因至少在 30% 樣本中超過閾值才保留 seed 20240101 )cfg這個(gè)列表是整條流水線的黑匣子開關(guān)。min_expr和min_samples直接決定后面有多少基因進(jìn)入多重檢驗(yàn)改動(dòng)了結(jié)論就得重出seed因?yàn)橐蚬謨?nèi)部有隨機(jī)抽樣不固定的話兩次跑出來(lái)的效應(yīng)估計(jì)在小數(shù)點(diǎn)后幾位會(huì)飄。路徑用here::here()而不是相對(duì)路徑是為了讓腳本從 RStudio、命令行、還是調(diào)度系統(tǒng)提交都能找到文件這個(gè)小習(xí)慣能省掉大量“在我機(jī)器上能跑”的扯皮。2.3 第一次跑通從 01 到 05 的執(zhí)行順序五個(gè)腳本是一條串行鏈但中間結(jié)果都落盤所以可以從任意一步重跑# 全流程按編號(hào)順序執(zhí)行 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.R01做數(shù)據(jù)對(duì)齊和表達(dá)預(yù)處理輸出data/derived/analysis_set.rds02和03分別在同一個(gè)分析集上跑 IPW 和 AIPW各自產(chǎn)出一張基因級(jí)效應(yīng)表04只在候選基因子集上跑因果森林用來(lái)觀察異質(zhì)性05把三張表按基因 ID 合并、統(tǒng)一多重檢驗(yàn)校正、出圖和寫報(bào)告。第一次跑建議先在01里把基因數(shù)砍到幾百個(gè)做冒煙測(cè)試確認(rèn)整條鏈沒斷再放開全量?；驍?shù)上萬(wàn)的表達(dá)矩陣走完全套單核大概要幾十分鐘到幾小時(shí)取決于樣本量和是否開并行。冒煙測(cè)試這一步別省04的因果森林在小樣本上最容易因?yàn)槟硞€(gè)基因方差為零直接拋錯(cuò)早發(fā)現(xiàn)早改。3. 數(shù)據(jù)層拆解暴露、結(jié)局與協(xié)變量怎么對(duì)齊到同一張表因果推斷里最容易翻車的從來(lái)不是估計(jì)量寫錯(cuò)而是數(shù)據(jù)對(duì)齊錯(cuò)了。這個(gè)項(xiàng)目的核心數(shù)據(jù)結(jié)構(gòu)只有一張寬表每一行是一個(gè)樣本列包括樣本 ID、吸煙暴露、幾十上百個(gè)協(xié)變量、以及每個(gè)基因的表達(dá)值。所有估計(jì)函數(shù)都只認(rèn)這張表保證口徑一致。3.1 吸煙暴露變量的編碼二值化之前的三個(gè)判斷原始表型里吸煙狀態(tài)通常是分好幾檔的從不、曾經(jīng)、當(dāng)前、偶爾、已戒若干年甚至還有每天支數(shù)。直接二值化會(huì)丟信息但項(xiàng)目走的是“當(dāng)前吸煙 vs 從不吸煙”這條二值路線理由有三短期因果效應(yīng)下“曾經(jīng)吸煙”人群的暴露已經(jīng)停止把它塞進(jìn)處理組會(huì)稀釋效應(yīng)每天支數(shù)在多數(shù)隊(duì)列里缺失率太高做連續(xù)暴露需要另一套劑量-反應(yīng)模型而二值暴露能直接套用傾向評(píng)分框架診斷工具也成熟。編碼這一步要把“曾經(jīng)吸煙”和缺失顯式剔出去而不是當(dāng)成對(duì)照組# 把原始吸煙編碼映射成 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)), # 只保留兩種明確狀態(tài)其余一律置 NA W case_when( smoke_status %in% c(current, current smoker, daily) ~ 1L, smoke_status %in% c(never, never smoker) ~ 0L, TRUE ~ NA_integer_ ), # 單獨(dú)打個(gè)標(biāo)記方便統(tǒng)計(jì)各剔除原因的數(shù)量 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%里放了多種寫法是為了兼容不同來(lái)源的表型表真實(shí)數(shù)據(jù)里 “Current Smoker” 和 “current” 混著出現(xiàn)是常態(tài)。drop_reason這一列看起來(lái)多余但它是后面寫方法學(xué)部分的證據(jù)——審稿人一定會(huì)問你為什么最終樣本是 N你得能答出“剔除了多少既往吸煙者、多少狀態(tài)缺失”。attr掛上去的計(jì)數(shù)會(huì)被01腳本打印到日志里養(yǎng)成這個(gè)習(xí)慣能省掉回頭重算的麻煩。還有一層容易忽略的暴露的時(shí)間窗。項(xiàng)目說的是“短期”那就得確認(rèn)表型采集時(shí)間和組織取樣時(shí)間接近如果表型是問卷回顧性的“過去一年是否吸煙”而樣本是五年前存的這個(gè)短期假設(shè)就不成立。數(shù)據(jù)里如果有采樣日期字段務(wù)必在預(yù)處理階段做個(gè)時(shí)間差統(tǒng)計(jì)。3.2 表達(dá)矩陣的預(yù)處理過濾、標(biāo)準(zhǔn)化、以及要不要去批次表達(dá)矩陣進(jìn)模型前要過三道關(guān)。第一道是低表達(dá)過濾用min_expr和min_samples兩個(gè)參數(shù)控制目的是減少檢驗(yàn)數(shù)量、降低多重校正的懲罰。第二道是分布變換如果是原始 count 就走log2(CPM 1)或vst如果上游已經(jīng)給了 log2 尺度就跳過。第三道是批次這一步是分水嶺批次如果和吸煙狀態(tài)強(qiáng)相關(guān)去批次會(huì)把效應(yīng)一起減掉如果批次獨(dú)立于暴露去不去影響都不大。保守做法是先檢查批次與暴露的關(guān)聯(lián)再做決定# 過濾 標(biāo)準(zhǔn)化并在返回前檢查批次是否與暴露混雜 prep_expression - function(expr, meta, cfg) { mat - as.matrix(expr[, -1, drop FALSE]) rownames(mat) - expr[[1]] # 1) 低表達(dá)過濾兩個(gè)條件同時(shí)滿足才保留 keep - rowMeans(mat cfg$min_expr) cfg$min_samples mat - mat[keep, , drop FALSE] # 2) log2 變換已經(jīng)在 log 尺度上的矩陣跳過 if (max(mat, na.rm TRUE) 50) mat - log2(mat 1) # 3) 分位數(shù)標(biāo)準(zhǔn)化抹平樣本間整體分布差異 mat - limma::normalizeBetweenArrays(mat, method quantile) # 4) 批次與暴露的關(guān)聯(lián)檢查卡方檢驗(yàn)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(批次-暴露關(guān)聯(lián)檢驗(yàn) p %.4g, p_batch)) if (!is.na(p_batch) p_batch 0.05) { warning(批次與吸煙狀態(tài)顯著關(guān)聯(lián)去批次可能吸收掉部分處理效應(yīng)慎用 ComBat) } } list(expr mat, n_gene nrow(mat), n_sample ncol(mat)) }rowMeans(mat cfg$min_expr)這一行是過濾的核心表達(dá)式返回每個(gè)基因超過閾值的樣本比例再和min_samples比。normalizeBetweenArrays走 quantile 方法對(duì)偏態(tài)分布的表達(dá)值比中位數(shù)標(biāo)準(zhǔn)化更穩(wěn)。最值得看的是第 4 段——批次卡方檢驗(yàn)。如果 p 值很小說明這批數(shù)據(jù)里“哪臺(tái)機(jī)器測(cè)的”和“這個(gè)人吸不吸煙”高度相關(guān)此時(shí)上 ComBat 或limma::removeBatchEffect等于把處理效應(yīng)的一部分當(dāng)成批次削掉了后面估計(jì)出來(lái)的效應(yīng)會(huì)系統(tǒng)性偏小。這種情況我一般不去批次而是把批次當(dāng)作協(xié)變量放進(jìn)傾向評(píng)分模型讓加權(quán)去平衡它而不是從數(shù)據(jù)里刪掉它。矩陣的樣本順序和表型的樣本順序必須嚴(yán)格一致這是最容易出錯(cuò)的地方。prep_expression返回前應(yīng)該強(qiáng)制按meta$sample_id重排或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次斷言別指望兩邊碰巧一致。3.3 協(xié)變量集與設(shè)計(jì)矩陣共線性檢查協(xié)變量集的選擇原則是只放暴露之前就已經(jīng)確定的變量年齡、性別、遺傳背景主成分、測(cè)序深度、批次不要放可能被暴露影響的中介變量比如某些炎癥指標(biāo)、肺功能。放中介進(jìn)去會(huì)把總效應(yīng)部分遮掉估出來(lái)的是直接效應(yīng)解釋起來(lái)是另一回事。確定變量集之后必須檢查設(shè)計(jì)矩陣的秩。協(xié)變量之間完全共線比如兩個(gè)高度相關(guān)的技術(shù)指標(biāo)或者某個(gè)分類變量某一水平只有一個(gè)人會(huì)讓傾向評(píng)分模型直接崩掉glm可能不報(bào)錯(cuò)但給出荒謬的系數(shù)。# 構(gòu)建協(xié)變量設(shè)計(jì)矩陣并做秩與稀疏水平檢查 check_design - function(meta, covars) { f - reformulate(covars) mm - model.matrix(f, data meta) # 秩虧說明存在完全共線glm 的結(jié)果不可信 r - qr(mm)$rank message(sprintf(設(shè)計(jì)矩陣維度 %d x %d秩 %d, nrow(mm), ncol(mm), r)) if (r ncol(mm)) warning(設(shè)計(jì)矩陣秩虧存在完全共線協(xié)變量需手動(dòng)剔除) # 分類變量稀疏水平某個(gè)水平樣本過少會(huì)讓該水平的傾向評(píng)分極端 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小于列數(shù)就是秩虧最常見的成因是分類變量做了獨(dú)熱編碼后某一列全是零或者兩個(gè)連續(xù)變量做了線性組合。tb 20這個(gè)閾值不是硬標(biāo)準(zhǔn)20 只是個(gè)提醒線真正要看的是該水平在兩組里的分布——如果某個(gè)批次只有 3 個(gè)人且全在吸煙組傾向評(píng)分在這個(gè)批次上會(huì)趨近 1加權(quán)時(shí)必然出極端權(quán)重這就是第 5 章要展開的問題。提前在這里發(fā)現(xiàn)比跑到加權(quán)階段再看權(quán)重分布圖省事得多。4. 因果估計(jì)核心三種估計(jì)量在單個(gè)基因上的完整鏈路數(shù)據(jù)對(duì)齊之后估計(jì)層的輸入就是一個(gè)向量結(jié)局Y某個(gè)基因的表達(dá)值、處理W0/1、協(xié)變量矩陣X。項(xiàng)目給三種估計(jì)量不是讓你三選一而是讓它們互相印證IPW 看加權(quán)后的可比人群AIPW 用兩個(gè)模型互相兜底因果森林看效應(yīng)是不是在亞群里不一致。三者結(jié)論方向一致才敢往報(bào)告里寫。4.1 傾向評(píng)分與 IPW把不可比的人群“造”成可比傾向評(píng)分e(X) P(W1|X)的直覺是給每個(gè)樣本算一個(gè)“被吸煙的概率”然后用這個(gè)概率的倒數(shù)做權(quán)重讓協(xié)變量分布在兩組之間拉平。吸煙概率高的人如果實(shí)際沒吸煙權(quán)重就大代表他“替”了很多吸煙者。實(shí)現(xiàn)上用WeightIt它比手寫glm多了一層診斷# 單個(gè)基因的 IPW 估計(jì)輸入是結(jié)局向量 處理向量 協(xié)變量數(shù)據(jù)框 estimate_ipw - function(y, w, x, trim 0.01) { df - data.frame(y y, W w, x) # 用 logistic 回歸估傾向評(píng)分estimand ATE 表示估的是全人群平均效應(yīng) wt - WeightIt::weightit(W ~ ., data df, method glm, estimand ATE, stabilize TRUE) ps - wt$ps # 重疊性診斷傾向評(píng)分貼邊說明兩組人根本不可比 overlap - c(min min(ps), max max(ps), n_lt01 sum(ps 0.01), n_gt99 sum(ps 0.99)) # 權(quán)重截尾壓掉極端權(quán)重的杠桿 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]) } # 加權(quán)最小二乘給 ATErobust TRUE 拿穩(wěn)健標(biāo)準(zhǔn)誤 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的區(qū)別要拎清ATE 估的是“如果全人群都吸煙 vs 全人群都不吸煙”的平均差異ATT 只估“吸煙者如果沒吸煙會(huì)怎樣”。做暴露危害評(píng)估一般用 ATE。stabilize TRUE把權(quán)重除以處理組的邊際概率能顯著降低權(quán)重的方差。overlap里那兩個(gè)n_lt01/n_gt99計(jì)數(shù)是診斷的核心——如果有一批人傾向評(píng)分大于 0.99說明模型認(rèn)為他們“幾乎肯定吸煙”而那批人恰好沒吸煙權(quán)重會(huì)飆到 100 以上一個(gè)樣本就能主導(dǎo)整個(gè)估計(jì)。trim 0.01是在 1% 和 99% 分位處截尾代價(jià)是引入輕微偏差但換來(lái)方差可控實(shí)踐中這筆交易通常劃算。ess有效樣本量是另一個(gè)必看指標(biāo)sum(w)^2 / sum(w^2)給出加權(quán)后相當(dāng)于多少“等效樣本”。如果原始 500 例加權(quán)后 ESS 只剩 60說明權(quán)重分布極其不均這個(gè)估計(jì)不該信。我的經(jīng)驗(yàn)閾值是 ESS 低于原始的 30% 就要回頭查協(xié)變量集。4.2 AIPW 雙重穩(wěn)健兩個(gè)模型互相兜底IPW 只依賴傾向評(píng)分模型正確一旦模型設(shè)定錯(cuò)了比如該放交互項(xiàng)沒放估計(jì)就有偏。AIPW 加了結(jié)局模型兩個(gè)模型只要有一個(gè)對(duì)估計(jì)就一致這叫雙重穩(wěn)健。代價(jià)是計(jì)算量大一些但在這個(gè)場(chǎng)景里值得。結(jié)局模型可以用簡(jiǎn)單的線性回歸也可以用grf的回歸森林。項(xiàng)目里默認(rèn)走線性因?yàn)榛虮磉_(dá)和協(xié)變量的關(guān)系大多近似線性森林在小樣本上反而方差大# AIPW 估計(jì)需要傾向評(píng)分 e(X) 和兩個(gè)結(jié)局預(yù)測(cè) mu1(X), mu0(X) estimate_aipw - function(y, w, x, use_forest FALSE) { df - data.frame(y y, W w, x) # 傾向評(píng)分模型 ps_fit - glm(W ~ ., data df, family binomial()) e - fitted(ps_fit) # 夾一下避免 e 正好等于 0 或 1 導(dǎo)致除零 e - pmin(pmax(e, 1e-6), 1 - 1e-6) # 結(jié)局模型分兩組各擬合一個(gè)預(yù)測(cè)反事實(shí)結(jié)局 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 影響函數(shù)逐樣本取值 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), # 結(jié)局模型在兩組上的外推幅度過說明模型在硬撐 extrap c(max_mu1 max(mu1), min_mu1 min(mu1), max_mu0 max(mu0), min_mu0 min(mu0))) }psi那四行是整套方法的核心。前兩項(xiàng)mu1 - mu0是結(jié)局模型的預(yù)測(cè)效應(yīng)后面兩個(gè)減項(xiàng)是“模型殘差乘以逆傾向評(píng)分”的修正項(xiàng)。直覺上如果結(jié)局模型預(yù)測(cè)得準(zhǔn)修正項(xiàng)接近零估計(jì)主要由mu1 - mu0決定如果傾向評(píng)分模型更準(zhǔn)修正項(xiàng)會(huì)把結(jié)局模型的偏差拉回來(lái)。e夾到[1e-6, 1-1e-6]是必須的現(xiàn)實(shí)數(shù)據(jù)里fitted(ps_fit)出現(xiàn) 1.0000 是常見的除法直接炸出 Inf。extrap那幾個(gè)極值是我加的診斷看mu1的最大值有沒有遠(yuǎn)超實(shí)際觀測(cè)到的表達(dá)量范圍——如果超過了說明結(jié)局模型在做荒謬的外推這條估計(jì)不能采信先回頭檢查協(xié)變量。use_forest TRUE那條分支用grf::regression_forest把W當(dāng)成一個(gè)普通特征塞進(jìn)X預(yù)測(cè)時(shí)分別把W置 1 和 0。這種寫法比訓(xùn)練兩個(gè)獨(dú)立森林更高效但對(duì)樣本量要求更高小于 200 例時(shí)慎用。4.3 因果森林平均效應(yīng)之外的異質(zhì)性IPW 和 AIPW 給的是全人群平均效應(yīng)但吸煙對(duì)基因表達(dá)的影響很可能在亞群里不一致——比如年齡大的人反應(yīng)更強(qiáng)或者只在不吸煙的遺傳背景子群里顯現(xiàn)。因果森林直接估每個(gè)樣本的條件平均處理效應(yīng)tau(X)然后看tau有沒有系統(tǒng)性的分化。# 因果森林估條件平均處理效應(yīng)看異質(zhì)性 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, # 誠(chéng)實(shí)劃分避免用同一批數(shù)據(jù)既分裂又估值 min.node.size 10, # 葉子最小樣本太小會(huì)過擬合 sample.fraction 0.5, # 子采樣比例 seed 20240101 ) # 全人群平均效應(yīng)target.sample overlap 只在有重疊的人群上平均 ate - grf::average_treatment_effect(cf, target.sample overlap) # 變量重要性看哪些協(xié)變量在驅(qū)動(dòng)異質(zhì)性 vimp - grf::variable_importance(cf) list(cf cf, ate ate, vimp vimp) }honesty TRUE是grf的默認(rèn)值作用是樣本分裂和效應(yīng)估計(jì)用不同的子樣本防止過擬合把噪聲當(dāng)成異質(zhì)性。這個(gè)參數(shù)不要關(guān)掉關(guān)了之后tau的分布會(huì)變得極其“好看”但全是假的。min.node.size 10是葉子節(jié)點(diǎn)的最小樣本量基因表達(dá)數(shù)據(jù)的噪聲大設(shè)太小會(huì)切出一堆只有幾個(gè)樣本的葉子估出來(lái)的tau方差爆炸。target.sample overlap這一項(xiàng)值得單獨(dú)說。默認(rèn)是all在全員上平均改成overlap只在傾向評(píng)分有重疊的區(qū)域上平均更穩(wěn)健但也更保守。做暴露效應(yīng)的基因篩查時(shí)用overlap結(jié)論不容易被少數(shù)極端個(gè)體帶跑。variable_importance的輸出是一個(gè)矩陣第一列是原始重要性第二列是歸一化到 0-1 的版本?？辞拔迕蛪颉绻挲g或某個(gè)技術(shù)協(xié)變量排第一說明異質(zhì)性主要來(lái)自技術(shù)因素而非生物學(xué)這時(shí)候要警惕別急著往機(jī)制上解釋。在基因?qū)用嬉蚬种辉诤蜻x基因子集上跑不是全基因組。原因很實(shí)際2000 棵樹乘以幾百個(gè)基因單機(jī)跑一晚上都不一定完而且大多數(shù)基因的異質(zhì)性檢驗(yàn)并不顯著跑全量是浪費(fèi)算力。4.4 從 p 值到 q 值結(jié)果表的組裝三種估計(jì)跑完每個(gè)基因會(huì)拿到三套est和se。組裝結(jié)果表要解決兩個(gè)問題一是統(tǒng)一效應(yīng)方向基因表達(dá)上調(diào)還是下調(diào)二是多重檢驗(yàn)?；驍?shù)量上千直接拿 p 0.05 篩出來(lái)的基本都是假陽(yáng)性。# 合并三種估計(jì)量的結(jié)果做 BH 校正輸出可讀的結(jié)果表 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 作為主估計(jì)它雙重穩(wěn)健 res$z - res$aipw_est / res$aipw_se res$p - 2 * pnorm(-abs(res$z)) res$q - p.adjust(res$p, method BH) # 方向標(biāo)記方便畫火山圖 res$dir - ifelse(res$q 0.05 res$aipw_est 0, up, ifelse(res$q 0.05 res$aipw_est 0, down, ns)) # 一致性旗標(biāo)兩種方法方向相反時(shí)標(biāo)出來(lái)這類基因不寫進(jìn)結(jié)論 res$flag_inconsistent - sign(res$ipw_est) ! sign(res$aipw_est) res[order(res$q), ] }p.adjust(method BH)是 Benjamini-Hochberg控制 FDR比 Bonferroni 溫和得多在基因表達(dá)場(chǎng)景是標(biāo)準(zhǔn)做法。dir那一列把顯著基因分成上調(diào)和下調(diào)畫圖直接映射顏色。flag_inconsistent這一列是經(jīng)驗(yàn)產(chǎn)物。IPW 和 AIPW 的方向如果相反說明兩個(gè)模型的假設(shè)沖突得很厲害可能是傾向評(píng)分模型設(shè)定有問題也可能是結(jié)局模型的線性假設(shè)在這個(gè)基因上不成立。這類基因即使 q 值很小也不該寫進(jìn)結(jié)論最多放進(jìn)補(bǔ)充材料標(biāo)注方法間不一致等著用敏感性分析去解釋。主估計(jì)選 AIPW 而不是 IPW理由是雙重穩(wěn)健只要傾向評(píng)分和結(jié)局模型有一個(gè)對(duì)估計(jì)就一致比單靠一個(gè)模型的 IPW 更扛得住設(shè)定誤差。5. 避坑清單因果推斷流水線里最容易翻車的五處這套流程跑通不難難的是跑出來(lái)的數(shù)能不能信。下面五條是我在這類項(xiàng)目上反復(fù)踩到的每一條都按“現(xiàn)象 → 原因 → 解決”整理配上排查代碼。5.1 權(quán)重爆炸與重疊性缺失現(xiàn)象IPW 估出來(lái)的效應(yīng)比 AIPW 大好幾倍標(biāo)準(zhǔn)誤也大得離譜ESS只剩原始樣本的 10% 不到。原因傾向評(píng)分模型太“自信”給某些樣本賦了接近 0 或 1 的概率。常見于某個(gè)協(xié)變量幾乎能完全預(yù)測(cè)處理狀態(tài)——比如某一批次的所有樣本都來(lái)自吸煙組模型在這個(gè)批次上學(xué)到“批次 吸煙”權(quán)重自然爆掉。另一個(gè)成因是分類變量的稀疏水平某一類只有幾個(gè)人邏輯回歸在該維度上給出極端的系數(shù)。解決先看重疊性診斷再?zèng)Q定是修剪樣本還是截尾權(quán)重。# 權(quán)重診斷分布、ESS、以及極端權(quán)重來(lái)自哪些樣本 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))) # 找出權(quán)重最大的 10 個(gè)樣本看看它們有什么共同點(diǎn) idx - order(wts, decreasing TRUE)[1:min(10, length(wts))] print(cbind(meta[idx, , drop FALSE], weight wts[idx])) invisible(ess) }打印極端權(quán)重樣本的協(xié)變量是定位問題最快的路徑。如果這十個(gè)樣本全在同一個(gè)批次那就是批次導(dǎo)致的如果全是高齡那是年齡分布尾部的問題。ESS低于 30% 時(shí)我一般先嘗試trim 0.01還不行就改用重疊權(quán)重WeightIt里method glm, estimand ATO重疊權(quán)重對(duì)極端傾向評(píng)分不敏感代價(jià)是估計(jì)目標(biāo)變成“重疊人群的效應(yīng)”解釋時(shí)要寫清楚。5.2 批次、共線性與“靜默”的樣本錯(cuò)位現(xiàn)象腳本一聲不吭跑完結(jié)果表也有但效應(yīng)估計(jì)整體偏小或者某個(gè)已知的強(qiáng)效應(yīng)基因死活不顯著。原因三種可能。批次與暴露強(qiáng)相關(guān)去批次時(shí)把處理效應(yīng)一起減掉了協(xié)變量共線導(dǎo)致glm系數(shù)不穩(wěn)定傾向評(píng)分本身就不準(zhǔn)表達(dá)矩陣和表型的樣本順序沒對(duì)齊merge或cbind之后列和行的對(duì)應(yīng)關(guān)系已經(jīng)亂了但 R 不會(huì)報(bào)錯(cuò)。解決三個(gè)檢查必須都在預(yù)處理階段做完——批次卡方檢驗(yàn)3.2 節(jié)、設(shè)計(jì)矩陣秩檢查3.3 節(jié)、樣本 ID 斷言。# 樣本對(duì)齊的最后一道保險(xiǎn)放在所有合并操作之后 assert_alignment - function(expr, meta) { stopifnot( 表達(dá)矩陣列名缺失 !is.null(colnames(expr)), 表型樣本 ID 缺失 sample_id %in% names(meta), 樣本數(shù)不一致 ncol(expr) nrow(meta), 樣本順序錯(cuò)位 identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }identical這一條是硬斷言不做容錯(cuò)順序不一致直接拋錯(cuò)。我寧可讓腳本在第一步就炸也不想在寫了兩周報(bào)告之后發(fā)現(xiàn)樣本對(duì)錯(cuò)了。5.3 因果森林在小樣本上的“看起來(lái)很準(zhǔn)”現(xiàn)象因果森林跑出來(lái)的tau分布特別漂亮亞群間差異明顯但換一個(gè)隨機(jī)種子結(jié)果就變了。原因honesty FALSE或者min.node.size設(shè)得太小樹在噪聲上過擬合?；虮磉_(dá)數(shù)據(jù)本身噪聲大樣本量幾百例時(shí)2000 棵樹的平均還壓不住方差tau的分布就是隨機(jī)的。解決三個(gè)參數(shù)一起調(diào)——保持honesty TRUEmin.node.size不低于 10sample.fraction取 0.4 到 0.6 之間然后用兩個(gè)不同種子各跑一遍看tau的相關(guān)系數(shù)。# 穩(wěn)定性檢驗(yàn)換種子重跑看條件效應(yīng)是否可復(fù)現(xiàn) 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 相關(guān)系數(shù) %.3f\n, r)) if (r 0.7) warning(條件效應(yīng)不穩(wěn)定異質(zhì)性結(jié)論不可靠) invisible(r) }相關(guān)系數(shù)低于 0.7 就別寫異質(zhì)性結(jié)論了老老實(shí)實(shí)報(bào)平均效應(yīng)。這個(gè)檢驗(yàn)成本不高值得每個(gè)候選基因都過一遍。5.4 多重檢驗(yàn)的兩種誤用現(xiàn)象q 值算出來(lái)了顯著基因上千個(gè)但和已知的吸煙相關(guān)基因?qū)Σ簧稀Ｔ蛞从昧?Bonferroni在基因表達(dá)場(chǎng)景過于嚴(yán)苛真正的信號(hào)被壓掉要么在篩選和驗(yàn)證之間反復(fù)調(diào)閾值最后報(bào)的是“調(diào)出來(lái)的”顯著性。還有一種隱蔽的IPW 和 AIPW 各算一遍 q 值挑顯著的那個(gè)報(bào)這是變相的多重比較。解決主分析只跑一套AIPWq 值只算一次閾值預(yù)先定死比如 q 0.05敏感性分析的結(jié)果單獨(dú)列不參與篩選。方法學(xué)部分要寫清楚“主估計(jì)用的是哪一個(gè)”。5.5 用同一份數(shù)據(jù)調(diào)參又評(píng)估現(xiàn)象所有診斷都完美重疊性好、ESS 高、q 值分布漂亮但換一個(gè)獨(dú)立隊(duì)列就復(fù)現(xiàn)不出來(lái)。原因傾向評(píng)分模型的協(xié)變量集、trim閾值、因果森林的參數(shù)如果在看到結(jié)果之后反復(fù)調(diào)整等于把結(jié)局信息泄漏進(jìn)了模型選擇。這類偏差不會(huì)體現(xiàn)在任何一張?jiān)\斷圖上。解決把樣本按 7:3 隨機(jī)切開前 70% 用來(lái)定協(xié)變量集和超參數(shù)后 30% 只跑一次不看結(jié)果不回頭。如果樣本量不夠切就在分析計(jì)劃里把參數(shù)全部預(yù)先寫死中途改了要記錄在日志里。# 預(yù)注冊(cè)式的參數(shù)凍結(jié)把最終參數(shù)寫進(jìn)文件后續(xù)腳本只讀不猜 freeze_params - function(cfg, path output/tables/frozen_params.rds) { saveRDS(cfg, path) message(參數(shù)已凍結(jié)至 , path, 后續(xù)分析請(qǐng)勿修改) }這條看著像流程潔癖但它是唯一能防住“調(diào)參調(diào)到顯著”的辦法。玄學(xué)的地方在于你自己調(diào)的時(shí)候完全不覺得在作弊每一步都有“合理”的理由。6. 二次開發(fā)把單基因流程改成批量篩選與可復(fù)現(xiàn)報(bào)告拿到這套腳本的人十有八九第一件事是想跑全基因組。單基因函數(shù)寫好了批量跑不是簡(jiǎn)單套個(gè)for循環(huán)——直接循環(huán)會(huì)撞上三個(gè)問題內(nèi)存、算力和失敗恢復(fù)。我給的做法是把估計(jì)函數(shù)包一層加上并行和緩存。先看批量層怎么寫# 批量估計(jì)并行 單基因失敗不影響整體 結(jié)果落盤緩存 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 會(huì)拋錯(cuò)并中斷整個(gè)批次 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 個(gè)基因估計(jì)失敗已跳過, sum(fails))) res[!fails] }furrr_options(seed TRUE)保證每個(gè) worker 的隨機(jī)數(shù)流可復(fù)現(xiàn)不加這一句因果森林那部分每次跑出來(lái)都不一樣。tryCatch把單基因的失敗隔離掉glm在完全分離的數(shù)據(jù)上會(huì)拋 warning 甚至 error一個(gè)基因炸掉不該拖垮整批。零方差判斷放在最前面是最常見的崩潰點(diǎn)——某些基因在所有樣本里表達(dá)值都一樣lm直接返回 NA 系數(shù)。n_workers不建議超過物理核心數(shù)減一grf內(nèi)部本身也會(huì)占內(nèi)存開太多反而會(huì)因?yàn)閮?nèi)存交換變慢。第二個(gè)改造點(diǎn)是報(bào)告層。原始腳本輸出的是 CSV二次開發(fā)可以把它變成參數(shù)化的 Rmd 或 Quarto 文檔把 q 值分布、火山圖、前若干個(gè)基因的森林圖一次性生成換一批數(shù)據(jù)只改cfg里的路徑。這里的關(guān)鍵是把“結(jié)果表”和“圖”解耦圖只讀表不重算估計(jì)。第三個(gè)改造點(diǎn)是加一層敏感性分析。做法很簡(jiǎn)單把trim從 0.01 換成 0.05 和 0.1各跑一遍 AIPW看顯著基因的重疊比例。敏感性維度參數(shù)觀察指標(biāo)警戒線權(quán)重截尾trim∈ {0.01, 0.05, 0.1}顯著基因重疊率低于 70% 說明結(jié)論依賴閾值結(jié)局模型線性 vs 回歸森林效應(yīng)估計(jì)相關(guān)系數(shù)低于 0.8 說明模型設(shè)定敏感協(xié)變量集加/減批次變量q 0.05 基因數(shù)變化數(shù)量差一倍以上要解釋這張表建議直接貼進(jìn)補(bǔ)充材料比任何文字描述都有說服力。審稿人問“你的結(jié)論穩(wěn)健嗎”把重疊率甩過去就行。最后一個(gè)技巧是日志。把sessionInfo()、cfg、每個(gè)腳本的運(yùn)行時(shí)間和輸入文件的行列數(shù)統(tǒng)一寫進(jìn)output/logs/一個(gè)腳本一個(gè)文件。這個(gè)習(xí)慣我是被坑出來(lái)的有次報(bào)告里的數(shù)和三個(gè)月前跑的對(duì)不上翻遍代碼也沒改過最后發(fā)現(xiàn)是上游給的表達(dá)矩陣悄悄更新了一版樣本數(shù)從 512 變成了 508。從那以后我每次跑因果估計(jì)流水線都強(qiáng)制在第一步把sessionInfo()和輸入文件的 md5 一起寫進(jìn)日志跑完再核對(duì)一遍才敢往報(bào)告里放數(shù)。這類工程的可靠性不在估計(jì)量寫得多漂亮而在每一步都能回溯。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取