貝氏層次模型（Bayesian Hierarchical Model, BHM）的發展歷史，是一段從「統計學界的邊緣理論」走向「現代數據科學與生物資訊核心」的逆襲史。它的誕生與普及，緊密伴隨著計算機算力的革命以及高維度數據（如轉錄組學、空間統計）的爆發。
以下為您梳理貝氏層次模型的發展軌跡與現代多元的應用範疇：
⏳ 貝氏層次模型的發展歷史
貝氏層次模型的演進可以大致分為四個里程碑階段：
1. 奠基期：詹姆斯-斯坦估計與收縮理論（1950 - 1960年代）
雖然托馬斯·貝氏（Thomas Bayes）在 18 世紀就提出了貝氏定理，但「層次模型」的統計學基礎直到 20 世紀中葉才出現。
關鍵突破：1961 年，統計學家查爾斯·斯坦（Charles Stein）與威拉德·詹姆斯（Willard James）提出了一個震驚統計學界的定理（詹姆斯-斯坦悖論）。他們證明：在估計三個或更多獨立的常態分佈均值時，將所有獨立樣本向全局平均值「收縮（Shrinkage）」後得到的估計值，其總體誤差必定小於傳統的最大概似估計（MLE）。
歷史意義：這顛覆了「各人自掃門前雪（獨立估計）」的傳統觀念，為層次模型中「訊息共享」的核心思想奠定了數學正當性。
2. 理論成型期：貝氏架構的全面引入（1970 - 1980年代）
隨後，統計學家發現詹姆斯-斯坦估計可以用貝氏機率的框架給予完美的解釋。
關鍵突破：1970 年代，布萊德利·艾弗隆（Bradley Efron，Bootstrap 的發明者）與卡爾·莫里斯（Carl Morris）發表了一系列論文，正式將收縮理論與貝氏事前分佈結合，展示了如何用「超參數（Hyperparameters）」來建構多層次的統計結構。
歷史意義：此時，層次模型（Hierarchical Model）與經驗貝氏（Empirical Bayes）的理論框架正式成型，統計學家開始意識到這種結構在處理群聚數據（Clustered data）時的巨大潛力。
3. 計算革命期：MCMC 演算法與 BUGS 的誕生（1990年代）
儘管 80 年代框架已成，但當時的貝氏層次模型面臨一個致命瓶頸：無法計算後驗分佈。高維度的多層積分在數學上沒有解析解，使得模型只能停留在理論階段。
計算機革命：1990 年，艾倫·蓋爾芬德（Alan Gelfand）和阿德里安·史密斯（Adrian Smith）將馬可夫鏈蒙地卡羅法（MCMC，特別是吉布斯採樣 Gibbs Sampling）引入貝氏統計，將複雜的高維積分問題轉化為計算機隨機抽樣。
軟體普及：隨後，BUGS（Bayesian inference Using Gibbs Sampling）軟體專案啟動，科學家終於有了工程工具可以真正「執行」貝氏層次模型。BHM 自此進入實用化爆炸期。
4. 現代爆發期：高維數據與機率程式語言（2000年代至今）
進入 21 世紀，隨著基因體學、大數據與人工智慧的興起，數據特徵數（如幾萬個基因）遠大於樣本數的現象成為常態。
現代技術：傳統的吉布斯採樣在高維空間容易卡頓，科學家開發了基於哈密頓力學的 HMC（Hamiltonian Monte Carlo） 與 NUTS（No-U-Turn Sampler） 演算法。
生態圈：這催生了新一代的機率程式語言（PPL），如您代碼中使用的 PyMC，以及 Stan、Pyro。現在，科學家只需要幾行代碼，就能在幾分鐘內完成過去需要幾天才能算完的高維貝氏層次推論。
🌍 貝氏層次模型的現代應用領域
由於 BHM 擅長處理「數據具有層級結構（如：細胞內有基因、學校內有學生、區域內有個體）」以及「小樣本、高雜訊」的問題，它在現代科學中得到了極其廣泛的應用：
1. 生物資訊學與基因體學（Bioinformatics）
正如您目前執行的轉錄組學分析，BHM 是多體學（Multi-omics）數據分析的黃金標準之一。
差異表達分析（Differential Expression）：在小樣本（如 $n=3$）下，利用全局基因背景約束個體變異數。經典工具如
limma、DESeq2的核心數學原理，本質上都是經驗貝氏層次模型（Empirical Bayes）。單細胞 RNA 測序（scRNA-seq）：單細胞數據具有極高比例的「零值雜訊（Dropout）」，BHM 可以透過細胞群體與基因群體的雙重層次，對缺失值進行穩健的填補（Imputation）與降噪。
2. 流行病學與空間疾病地圖（Epidemiology & Spatial Mapping）
在公共衛生領域，科學家需要評估不同行政區的疾病發病率。
小區域估計（Small Area Estimation）：某些偏遠小鎮人口極少，只要出現 1 個癌症病例，表面上的發病率就會飆高。
解法：流行病學常用 BYM 模型（Besag-York-Mollié，一種空間貝氏層次模型），讓小鎮的發病率向「全縣平均（全局）」以及「鄰近鄉鎮（空間相關）」進行收縮。這能有效剃除因人口過少產生的統計雜訊，繪製出真正精準的疾病風險地圖。
3. 生態學與環境科學（Ecology）
生態學的野生動物調查往往面臨極端惡劣的觀測條件。
動植物種群動態：科學家在野外架設紅外線相機捕捉動物，觀測數據受到天氣、地形、相機故障等重重干擾。
解法：透過層次模型，第一層建模「動物真實存在的空間分佈（狀態）」，第二層建模「相機捕捉到動物的捕獲機率（觀測雜訊）」。BHM 能夠成功將「生物真實狀態」從「環境干擾雜訊」中剝離出來。
4. 臨床試驗與統合分析（Meta-Analysis）
在醫學藥物開發中，科學家需要綜合評估全球數十個不同中心、不同樣本量發表的臨床報告。
統合分析：每個臨床試驗（Study）都有其特異性（如不同國家、不同年齡層）。
解法：BHM 將每個試驗視為第二層個體，全球總效果視為第一層超事前分佈。它能完美評估「研究間的異質性（Heterogeneity）」，即使某些小規模試驗流於偏差，模型也能給予合理的權重收縮，給出最客觀的藥效評估。
5. 金融、行銷與商業決策
在商業零售中，企業需要預測成千上萬種商品的未來銷量，或是評估不同地區廣告的投放效果。
多層次行銷模型（Hierarchical MMM）：將消費者依據「城市」、「年齡層」分層，利用 BHM 既能學到全國消費者的共同行為趨勢（超參數），又能捕捉到特定城市（如台北 vs 高雄）的特異性偏好。這在數據稀疏的細分市場預測中表現極其優異。
📝 歷史與應用的總結
從詹姆斯-斯坦悖論引發的觀念革命，到 MCMC 演算法引爆的計算革命，貝氏層次模型之所以能在當代各大領域成為大宗，是因為它承載了人類對數據認知的轉變：我們不再孤立地看待每一個觀測對象，而是將其納入系統性的全局網絡中。 這種「和而不同、資訊共享」的哲學，使其成為當今應對複雜、高雜訊、小樣本數據最精密的統計武器。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 3 09:15:55 2026
核心建立在貝氏層次模型（Bayesian Hierarchical Model）之上，用來處理生物統計學中典型「小樣本、
高維度（高基因數）」的轉錄組特徵選取（Feature Selection）挑戰。
使用貝氏機率理論（Bayesian Inference）進行特徵選取（Feature Selection）是最佳策略。
透過貝氏方法，我們可以利用收縮估計（Shrinkage Estimation）或層次模型（Hierarchical Models）
共享基因間的變異數資訊，從而做出穩健的推論。
執行上述程式後，你會得到一個統計表格，篩選邏輯如下：
1. 核心指標：後驗機率（Posterior Probability）Prob_UP (AT>MT)：代表處理後基因表現量確實
上升的客觀機率。如果該值為 0.97，代表有 97% 的把握該基因在處理後被活化。Prob_DOWN (AT<MT)：
代表處理後表現量下降的機率。
2. 效果量（Bayesian Effect Size）由於樣本數只有 3 組，有時候雖然平均差值大，但可能只是單一
鰻魚的極端值（離群值）。Bayesian_Effect_Size 納入了變異度（$\mu / \sigma$）。數值絕對值
越大（例如 $> 1.5$ 或 $< -1.5$），代表該基因的表達量變化相對於其波動而言非常顯著，發育程度的
差異與該基因高度相關。
3. 如何結合「卵巢發育程度的差異」？如果你有進一步紀錄這三隻鰻魚處理後的卵巢發育成熟度分數（例如：
Fish 1 熟度 80%, Fish 2 熟度 50%, Fish 3 熟度 20%），你可以將模型升級為貝氏線性迴歸模型：
透過檢定斜率 beta_1 的後驗分佈是否遠離 0，就能精準篩選出「其表達量變化與卵巢發育程度呈正/負相關」
的核心關鍵基因。
程式碼在執行時遇到了轉錄組高維度數據分析最常見的兩個嚴重瓶頸：嚴重的發散收斂問題（Divergences）
與 $\hat{R} > 1.01$： 因為每組基因只有 3 個成對樣本（樣本數極少），當你對每個基因「獨立」建
立貝氏模型時，模型很難單憑 3 個點準確估計出該基因的變異度 sigma。這在貝氏統計中被稱為「漏斗效應
（Neal's Funnel）」，會導致 MCMC 採樣器瘋狂報錯、結果不可信。效率極低（維度災難）： 你的轉錄
組矩陣包含 21,651 個基因。如果用 for gene in genes: 迴圈一個一個跑 pm.sample()，跑完兩萬
多個基因可能需要耗費幾十個小時甚至好幾天。💡 完美的解決方案：貝氏層次模型（Bayesian
Hierarchical Model）我們應該利用貝氏的核心優勢——層次模型（Hierarchical Modeling），將所
有基因放進同一個模型中進行矩陣化（Vectorized）平行運算。透過讓所有基因共享同一個上位的事前分佈
（Hyper-priors），模型會自動進行「收縮估計（Shrinkage Estimation）」。白話來說，就是讓表達
量穩定的基因去「幫助」變異大的基因修正其變異數。如此一來：徹底解決少樣本導致的發散（Divergences）
問題。速度提升數千倍（20,000+ 基因可在幾分鐘內透過矩陣運算一次採樣完成）。
@author: yshuang
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
import time
from pathlib import Path
from typing import Tuple, Optional
import warnings
# ============================================================
# 配置區塊（集中管理所有路徑與參數）
# ============================================================
class Config:
"""集中管理所有配置參數，避免硬編碼分散在代碼中"""
# 路徑配置
INPUT_FILE = Path("/Users/yshuang/Documents/Python/FC_GSEA.xlsx")
OUTPUT_DIR = Path("/Users/yshuang/Documents/Python")
OUTPUT_CSV = OUTPUT_DIR / "Bayesian_Phenotype_Selection_Results.csv"
OUTPUT_RNK = OUTPUT_DIR / "Eel_Ovary_AT_vs_MT_Phenotype_Bayesian.rnk"
# 數據配置
SHEET_NAME = "data (TPM)"
USE_COLS = "A, E:J"
N_SAMPLES_PER_GROUP = 3 # 每組樣本數
# 💡【外表型差型定量變數】
# 根據 1440 > 2613 >> 2003 邏輯賦予發育推進權重 (AT - MT 的表型變化量)
# 這裡的 Key 必須與您的 Excel 樣本名稱中的個體編號完全對應
PHENOTYPE_DELTA_MAP = {
"1440": 3.0, # 強
"2613": 2.0, # 中
"2003": 0.5 # 弱
}
# MCMC 配置
DRAWS = 1000
TUNE = 1000
CHAINS = 2
TARGET_ACCEPT = 0.95
RANDOM_SEED = 42
# 特徵選取門檻 (改看斜率為正/負的後驗機率)
PROB_THRESHOLD = 0.95
TOP_N_DISPLAY = 10
def load_and_preprocess_data(config: Config) -> Tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, np.ndarray, np.ndarray, list]:
"""
讀取轉錄組數據，並提取與個體嚴格對齊的外表型差值變數
"""
if not config.INPUT_FILE.exists():
raise FileNotFoundError(f"❌ 輸入檔案不存在: {config.INPUT_FILE}")
print(f"📂 讀取數據: {config.INPUT_FILE}")
df = pd.read_excel(
config.INPUT_FILE,
sheet_name=config.SHEET_NAME,
index_col=0,
usecols=config.USE_COLS
)
df.columns = df.columns.str.strip()
if df.shape[1] == config.N_SAMPLES_PER_GROUP * 2:
Y_raw = df.T.astype(float)
else:
Y_raw = df.astype(float)
print(f"✅ 基因矩陣載入完成: {Y_raw.shape[0]} 樣本 × {Y_raw.shape[1]} 基因")
Y_asinh = np.arcsinh(Y_raw)
mt_samples = sorted([col for col in Y_asinh.index if 'MT' in col.upper()])
at_samples = sorted([col for col in Y_asinh.index if 'AT' in col.upper()])
if len(mt_samples) != len(at_samples) or len(mt_samples) == 0:
raise ValueError("❌ 樣本配對失敗，請檢查 MT 與 AT 命名。")
print(f"📊 已成功配對樣本: 對照組 {mt_samples} ↔️ 處理組 {at_samples}")
# 計算成對差值
data_mt = Y_asinh.loc[mt_samples].values
data_at = Y_asinh.loc[at_samples].values
diff_matrix = data_at - data_mt
X_data = diff_matrix.T # (n_genes, n_samples)
genes = Y_raw.columns.tolist()
# 💡【核心新增】：動態解析樣本名稱，並匹配外表型定量值
phenotype_list = []
for sample_name in mt_samples:
matched = False
for key, val in config.PHENOTYPE_DELTA_MAP.items():
if key in sample_name:
phenotype_list.append(val)
matched = True
break
if not matched:
raise ValueError(f"❌ 樣本 {sample_name} 無法匹配到任何外表型設定，請檢查 Config.PHENOTYPE_DELTA_MAP")
X_phenotype = np.array(phenotype_list)
print(f"📈 樣本對齊的外表型推進值 (ΔPhenotype): {dict(zip(mt_samples, X_phenotype))}")
diff_df = pd.DataFrame(diff_matrix, index=mt_samples, columns=genes)
return Y_raw, Y_asinh, diff_df, X_data, X_phenotype, genes
def build_bayesian_regression_model(X_data: np.ndarray, X_phenotype: np.ndarray, config: Config) -> az.InferenceData:
"""
建立並執行貝氏層次線性迴歸模型
使得基因表達量差值 (Y) 與外表型推進值 (X) 進行耦合
"""
n_genes, n_samples = X_data.shape
print(f"🚀 啟動貝氏層次迴歸 MCMC 採樣...")
start_time = time.time()
with pm.Model() as hierarchical_regression_model:
# --- 1. 迴歸截距 (Intercept) 層次分佈 ---
alpha_global = pm.Normal("alpha_global", mu=0, sigma=2)
sigma_alpha = pm.HalfNormal("sigma_alpha", sigma=2)
alpha_offset = pm.Normal("alpha_offset", mu=0, sigma=1, shape=n_genes)
alpha_genes = pm.Deterministic("alpha_genes", alpha_global + alpha_offset * sigma_alpha)
# --- 2. 迴歸斜率 (Slope, 對外表型的敏感度) 層次分佈 ---
# beta_global 代表全轉錄組受到該外表型驅動的平均斜率趨勢
beta_global = pm.Normal("beta_global", mu=0, sigma=2)
sigma_beta = pm.HalfNormal("sigma_beta", sigma=2)
beta_offset = pm.Normal("beta_offset", mu=0, sigma=1, shape=n_genes)
# beta_genes 代表每個基因各自對外表型發育強弱的「響應斜率」
beta_genes = pm.Deterministic("beta_genes", beta_global + beta_offset * sigma_beta)
# --- 3. 群體共享的殘差變異度 ---
sigma_noise = pm.HalfNormal("sigma_noise", sigma=2)
# --- 4. 線性模型方程式： μ = 基礎差值 + 斜率 * 外表型變化量 ---
mu_predicted = alpha_genes[:, None] + beta_genes[:, None] * X_phenotype[None, :]
# 似然函數
likelihood = pm.Normal(
"y_obs",
mu=mu_predicted,
sigma=sigma_noise,
observed=X_data
)
# 執行採樣
trace = pm.sample(
draws=config.DRAWS,
tune=config.TUNE,
chains=config.CHAINS,
target_accept=config.TARGET_ACCEPT,
return_inferencedata=True,
random_seed=config.RANDOM_SEED,
progressbar=True
)
elapsed = time.time() - start_time
print(f"🎉 MCMC 完成！耗時: {elapsed:.1f} 秒")
return trace
def extract_posterior_statistics(
trace: az.InferenceData,
X_data: np.ndarray,
genes: list,
config: Config
) -> pd.DataFrame:
"""從迴歸模型的後驗分佈提取指標（改以斜率 Beta 為核心）"""
print("📊 計算外表型相關後驗統計量...")
# 提取響應外表型的斜率後驗樣本 (n_genes, n_draws)
beta_samples = az.extract(trace, var_names="beta_genes").values
alpha_samples = az.extract(trace, var_names="alpha_genes").values
sigma_noise_samples = az.extract(trace, var_names="sigma_noise").values
# 計算統計量
mean_diff = X_data.mean(axis=1)
post_beta_mean = beta_samples.mean(axis=1)
post_alpha_mean = alpha_samples.mean(axis=1)
# 💡 計算該基因與「外表型推進」呈正相關或負相關的後驗機率
prob_up = np.mean(beta_samples > 0, axis=1)
prob_down = np.mean(beta_samples < 0, axis=1)
# 貝氏效果量：以斜率（信噪比）作為核心
sigma_broadcast = np.tile(sigma_noise_samples, (len(genes), 1))
epsilon = 1e-10
effect_size = np.mean(beta_samples / (sigma_broadcast + epsilon), axis=1)
df_results = pd.DataFrame({
"Raw_Mean_Diff": mean_diff,
"Bayesian_Intercept(Alpha)": post_alpha_mean,
"Bayesian_Slope(Beta)": post_beta_mean,
"Prob_Phenotype_Pos(Beta>0)": prob_up,
"Prob_Phenotype_Neg(Beta<0)": prob_down,
"Bayesian_Effect_Size": effect_size
}, index=genes)
# 特徵選取：斜率方向信心大於門檻者
df_results["Selected"] = (
(df_results["Prob_Phenotype_Pos(Beta>0)"] > config.PROB_THRESHOLD) |
(df_results["Prob_Phenotype_Neg(Beta<0)"] > config.PROB_THRESHOLD)
)
# 按表型驅動效果量絕對值排序
df_results = df_results.sort_values(
by="Bayesian_Effect_Size",
key=lambda x: np.abs(x),
ascending=False
)
print(f"\n--- 前 {config.TOP_N_DISPLAY} 個強烈受外表型驅動基因 ---")
print(df_results.head(config.TOP_N_DISPLAY).to_string())
n_selected = df_results["Selected"].sum()
print(f"\n📈 特徵選取摘要: {n_selected}/{len(genes)} 基因與外表型發育軌跡強烈耦合 (門檻 {config.PROB_THRESHOLD*100:.0f}%)")
return df_results
def export_gsea_ranklist(df_results: pd.DataFrame, config: Config) -> None:
"""匯出基於外表型迴歸斜率排序的 GSEA .rnk 檔案"""
print("\n🚚 生成外表型驅動型 GSEA Rank List...")
config.OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
df_gsea = df_results[["Bayesian_Effect_Size"]].copy().dropna()
if df_gsea.index.duplicated().any():
df_gsea = df_gsea[~df_gsea.index.duplicated(keep="first")]
df_gsea = df_gsea.reset_index()
df_gsea.columns = ["Gene", "Score"]
df_gsea = df_gsea.sort_values(by="Score", ascending=False).reset_index(drop=True)
df_gsea.to_csv(config.OUTPUT_RNK, sep="\t", header=False, index=False)
print(f"✅ GSEA .rnk 匯出成功: {config.OUTPUT_RNK}")
def run_diagnostics(trace: az.InferenceData) -> None:
"""收斂診斷"""
print("\n🔍 收斂診斷...")
summary = az.summary(trace, var_names=["alpha_global", "beta_global", "sigma_noise"])
print("\n--- 全局超參數迴歸收斂摘要 ---")
print(summary[["mean", "sd", "r_hat"]].to_string())
def main():
config = Config()
try:
# 1. 數據載入與預處理 (包含外表型對齊)
Y_raw, Y_asinh, diff_df, X_data, X_phenotype, genes = load_and_preprocess_data(config)
# 2. 建立貝氏層次迴歸模型
trace = build_bayesian_regression_model(X_data, X_phenotype, config)
# 3. 收斂診斷
run_diagnostics(trace)
# 4. 提取迴歸統計量
df_results = extract_posterior_statistics(trace, X_data, genes, config)
# 5. 匯出完整結果
df_results.to_csv(config.OUTPUT_CSV)
print(f"\n💾 迴歸分析結果已儲存: {config.OUTPUT_CSV}")
# 6. 匯出 GSEA 檔案
export_gsea_ranklist(df_results, config)
print("\n🎊 外表型耦合分析全數完成！")
except Exception as e:
print(f"\n❌ 執行錯誤: {e}")
raise
if __name__ == "__main__":
main()
