q-값 (q-Value)

Storey의 FDR 추정 — BH보다 강력한 다중 검정 보정과 local FDR

Benjamini-Hochberg(BH) 절차는 FDR을 제어하지만 참 귀무가설 비율 π₀를 보수적으로 1로 가정한다. Storey(2002)의 q-값은 π₀를 데이터로 추정하여 BH보다 더 많은 가설을 기각하면서 FDR을 더 정밀하게 추정한다. q-값의 정의, π₀ 추정 방법, BH p-값과의 관계, local FDR과의 비교를 직관과 수식으로 설명한다.

Statistics
저자

Kwangmin Kim

공개

2026년 04월 09일

1 개요

다중 검정에서 Benjamini-Hochberg(BH) 절차는 FDR을 \(q\) 이하로 제어한다. 그런데 BH는 암묵적으로 모든 \(m\) 개의 귀무가설이 참일 수 있다고 가정한다 — 즉, 참 귀무가설 비율 \(\pi_0 = m_0 / m = 1\) 로 취급한다.

실제로는 \(m_0 < m\) 이다 (일부 귀무가설은 거짓). \(\pi_0\) 를 1보다 작게 추정하면: - FDR 추정이 더 정확해지고 - 더 많은 귀무가설을 기각할 수 있다 (검정력 향상)

Storey(2002)와 Storey & Tibshirani(2003)는 \(\pi_0\) 를 데이터로 추정하여 각 귀무가설에 q-값(q-value) 을 부여하는 방법을 제안했다.

q-값은 “이 귀무가설을 기각할 때의 예상 FDR”이며, p-값이 “이 검정 통계량 이상이 귀무 분포에서 나올 확률”인 것과 대응된다.

2 왜 \(\pi_0\) 추정이 필요한가

2.1 BH의 보수성 원인

BH 절차가 보장하는 것은 (Benjamini & Hochberg, 1995):

\[\text{FDR}_{\text{BH}} \leq \frac{m_0}{m} \cdot q = \pi_0 \cdot q \leq q\]

BH는 \(\pi_0 = 1\) 을 가정하므로 실제 FDR은 \(\pi_0 \cdot q\) 인데, \(\pi_0 < 1\) 이면 실제 FDR이 목표 \(q\) 보다 훨씬 낮다 → BH가 지나치게 보수적이다.

직관: 1,000개의 유전자 중 실제로 200개만 발현 차이가 있다고 하자 (\(\pi_0 = 0.8\)). BH는 마치 1,000개 모두 “무관한 유전자일 수 있다”고 가정하므로, 실제로는 더 많은 유전자를 유의하게 선택해도 FDR이 목표치를 넘지 않는다.

2.2 \(\pi_0\) 를 추정하면

FDR을 다음과 같이 더 정밀하게 추정할 수 있다:

\[\widehat{\text{FDR}}(t) = \frac{\hat{\pi}_0 \cdot m \cdot t}{R(t)}\]

여기서: - \(t\) : p-값 임계값 - \(R(t) = \#\{j : p_j \leq t\}\) : \(t\) 이하 p-값의 수 (기각 수) - \(\hat{\pi}_0\) : 추정된 참 귀무가설 비율 - \(\hat{\pi}_0 \cdot m \cdot t\) : 임계값 \(t\) 이하의 위양성 수 기대값 (참 귀무가설에서 균등 분포 가정)

3 q-값의 정의

q-값 (Storey, 2002)

귀무가설 \(H_{0j}\) 에 대한 q-값은:

\[q(p_j) = \min_{t \geq p_j} \widehat{\text{FDR}}(t) = \min_{t \geq p_j} \frac{\hat{\pi}_0 \cdot m \cdot t}{R(t)}\]

즉, “관측된 p-값이 \(p_j\) 이하인 가설들을 모두 기각할 때의 최소 예상 FDR”이다.

p-값과 q-값의 대응 관계:

p-값 q-값
정의 \(H_0\) 하에서 이 이상 극단적인 통계량의 확률 이 가설을 기각할 때의 예상 FDR
제어 대상 Type I Error (단일 검정) FDR (다중 검정)
귀무 분포 이론적 분포 또는 순열 p-값들의 경험 분포 활용
해석 “유의성 증거의 강도” “발견 품질의 기대치”
사용 맥락 단일/소규모 검정 대규모 다중 검정

q-값은 단조 증가하도록 강제된다 (\(\min_{t \geq p_j}\)): p-값이 작을수록 q-값도 작다.

4 \(\pi_0\) 추정: Storey의 방법

4.1 핵심 아이디어

참 귀무가설 하에서 p-값은 \([0,1]\) 의 균등 분포를 따른다. 큰 p-값 (예: \(p > \lambda\), \(\lambda = 0.5\)) 은 거의 대부분 참 귀무가설에서 온 것이다.

따라서 \(p > \lambda\) 인 p-값의 수를 세어 \(\pi_0\) 를 추정한다:

\[\hat{\pi}_0(\lambda) = \frac{\#\{j : p_j > \lambda\}}{m(1-\lambda)} \tag{Storey 1}\]

직관: \(H_0\) 가 참이면 p-값이 \((\lambda, 1]\) 에 떨어질 확률은 \(1-\lambda\) 다. \(m_0\) 개의 참 귀무가설이 있으면 이 구간에 평균 \(m_0 (1-\lambda)\) 개가 있다. 거짓 귀무가설의 p-값은 0 근처에 몰리므로, \(\lambda\) 가 충분히 크면 \(p > \lambda\) 인 대부분은 참 귀무가설이다.

\[\hat{\pi}_0(\lambda) \approx \frac{m_0 (1-\lambda)}{m(1-\lambda)} = \frac{m_0}{m} = \pi_0\]

\(\lambda\) 선택의 딜레마
  • \(\lambda\) 가 작으면: 거짓 귀무가설의 p-값도 포함 → \(\hat{\pi}_0\) 과대추정 (보수적)
  • \(\lambda\) 가 크면: 추정에 쓰이는 p-값 수가 적어 분산 증가

Storey & Tibshirani(2003)는 여러 \(\lambda\) 값에서 \(\hat{\pi}_0(\lambda)\) 를 계산하고 자연 스플라인(natural spline) 으로 \(\lambda \to 1\) 의 극한을 추정한다:

\[\hat{\pi}_0 = \hat{f}(1)\]

여기서 \(\hat{f}\) 는 \(\{(\lambda, \hat{\pi}_0(\lambda))\}\) 에 적합한 스플라인이다.

5 q-값 계산 알고리즘

Storey q-값 계산 절차

입력: p-값 \(p_1, \ldots, p_m\), 선택적으로 \(\lambda\) 격자

  1. p-값을 오름차순 정렬: \(p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}\)
  2. \(\hat{\pi}_0\) 를 추정한다 (여러 \(\lambda\) 에서 스플라인 추정)
  3. 각 임계값 \(t = p_{(j)}\) 에서 FDR을 추정한다: \[\widehat{\text{FDR}}(p_{(j)}) = \frac{\hat{\pi}_0 \cdot m \cdot p_{(j)}}{j}\]
  4. q-값을 단조 증가 제약 하에 계산한다: \[\hat{q}(p_{(j)}) = \min_{k \geq j} \widehat{\text{FDR}}(p_{(k)})\]
  5. 원래 순서로 복원한다.

수치 예시 (\(m = 5\), \(\hat{\pi}_0 = 0.8\)):

순위 \(j\) \(p_{(j)}\) \(\widehat{\text{FDR}}\) q-값
1 0.006 \(0.8 \times 5 \times 0.006 / 1 = 0.024\) 0.024
2 0.012 \(0.8 \times 5 \times 0.012 / 2 = 0.024\) 0.024
3 0.601 \(0.8 \times 5 \times 0.601 / 3 = 0.801\) 0.801
4 0.756 \(0.8 \times 5 \times 0.756 / 4 = 0.756\) 0.756
5 0.918 \(0.8 \times 5 \times 0.918 / 5 = 0.734\) 0.734

BH 절차 (\(\hat{\pi}_0 = 1\)) 와 비교 시, q-값은 \(\hat{\pi}_0 < 1\) 을 반영하여 더 작다 → 같은 FDR 기준에서 더 많은 가설을 기각.

6 BH p-값과 q-값의 수치적 관계

BH 절차의 조정된 p-값(adjusted p-value, BH p-값)과 Storey q-값은 개념적으로 유사하다:

  • BH 조정 p-값: \(p_j^{\text{BH}} = \min_{k \geq j} \frac{m \cdot p_{(k)}}{k}\) (단조 증가)
  • Storey q-값: \(q_j = \min_{k \geq j} \frac{\hat{\pi}_0 \cdot m \cdot p_{(k)}}{k}\)

관계: \(q_j = \hat{\pi}_0 \cdot p_j^{\text{BH}}\)

\(\hat{\pi}_0 \leq 1\) 이므로 q-값 \(\leq\) BH 조정 p-값 → q-값이 항상 더 작거나 같다.

ISLR의 R 출력에서:

q.values.BH[1:10]
[1] 0.08989 0.99149 0.12212 0.92343 0.95604 ...

이 “q-값”은 실제로 BH 조정 p-값이다 (\(\hat{\pi}_0 = 1\) 가정). 진정한 Storey q-값은 qvalue 패키지를 사용해야 한다.

7 Local FDR (국소 FDR)

7.1 개념

q-값은 “p-값 \(\leq p_j\) 인 모든 가설을 기각했을 때의 FDR”이다. 이는 global (또는 tail-area) FDR이다.

반면 local FDR (Efron et al., 2001)은 “p-값이 정확히 \(p_j\) 근방에 있는 가설들만의 위양성 비율”이다:

\[\text{fdr}(p) = \frac{\pi_0 \cdot f_0(p)}{f(p)}\]

여기서: - \(f_0(p)\): 참 귀무가설에서 p-값의 밀도 (균등 분포 = 1) - \(f(p)\): 관측된 p-값의 혼합 밀도

Global FDR (q-값) Local FDR
정의 \(E(V/R)\) 추정 \(\pi_0 f_0(p) / f(p)\)
해석 이 임계값까지 기각 시 위양성 비율 이 p-값을 가진 가설 자체가 참일 확률
역할 기각 집합 전체 품질 개별 가설의 신빙성
계산 p-값 정렬 + 비율 밀도 추정 필요

8 실무 지침

8.1 q-값 사용 시나리오

  • 유전체/전사체 분석: DESeq2, edgeR의 기본 출력이 BH 조정 p-값; Storey q-값으로 교체 시 발견 수 증가
  • GWAS: LD 구조가 있어 q-값 계산이 복잡 → 순열 기반 FDR 선호
  • 일반 대규모 탐색: 가설 수 \(m \geq 1{,}000\), p-값 분포의 귀무 성분이 명확할 때

8.2 BH를 쓸지 Storey q-값을 쓸지

상황 권장 이유
규제 제출, 확증 연구 BH 보수적이고 증명 쉬움
탐색적, 대규모 Storey q-값 더 높은 검정력
\(m\) 작음 (\(<100\)) BH \(\hat{\pi}_0\) 추정 불안정
p-값 분포가 불균등 Storey q-값 \(\pi_0 < 1\) 을 활용

9 코드 예시

9.1 Step 1: 순수 Python — Storey q-값 직접 구현

import numpy as np
from scipy.interpolate import UnivariateSpline

def estimate_pi0(pvalues: np.ndarray, lambdas: np.ndarray = None) -> float:
    """
    Storey & Tibshirani(2003) 방법으로 π₀ 추정.
    여러 λ에서 π₀(λ)를 계산 후 자연 스플라인으로 λ→1 극한 추정.
    """
    m = len(pvalues)
    if lambdas is None:
        lambdas = np.arange(0.05, 0.96, 0.05)
    
    # 각 λ에서 π₀(λ) 계산
    pi0_lambda = np.array([
        np.sum(pvalues > lam) / (m * (1 - lam))
        for lam in lambdas
    ])
    
    # 자연 스플라인으로 λ=1에서의 극한 추정
    spline = UnivariateSpline(lambdas, pi0_lambda, k=3, s=len(lambdas))
    pi0_hat = float(spline(1.0))
    
    return min(pi0_hat, 1.0)  # π₀ ≤ 1 제약


def storey_qvalues(pvalues: np.ndarray, pi0: float = None) -> np.ndarray:
    """
    Storey(2002) q-값 계산.
    pvalues: (m,) p-값 배열
    pi0: 미리 추정된 π₀ (None이면 자동 추정)
    반환: (m,) q-값 배열
    """
    m = len(pvalues)
    
    if pi0 is None:
        pi0 = estimate_pi0(pvalues)
    
    order = np.argsort(pvalues)
    sorted_p = pvalues[order]
    
    # FDR 추정: π₀ * m * p_{(j)} / j
    fdr_est = pi0 * m * sorted_p / np.arange(1, m + 1)
    
    # 단조 증가 제약: q_{(j)} = min_{k≥j} FDR_est(k)
    # 뒤에서부터 누적 최솟값
    qvalues_sorted = np.minimum.accumulate(fdr_est[::-1])[::-1]
    qvalues_sorted = np.clip(qvalues_sorted, 0, 1)
    
    # 원래 순서로 복원
    qvalues = np.empty(m)
    qvalues[order] = qvalues_sorted
    
    return qvalues, pi0


# ── 시뮬레이션 예시 ──────────────────────────────────────────
rng = np.random.default_rng(42)
m = 1000
m0 = 800   # 참 귀무가설 80%
m1 = m - m0

# 참 귀무가설: 균등 분포
p_null = rng.uniform(0, 1, size=m0)
# 거짓 귀무가설: 작은 p-값 집중
p_alt = rng.beta(0.3, 5, size=m1)
pvalues = np.concatenate([p_null, p_alt])
true_null = np.array([True]*m0 + [False]*m1)

# q-값 계산
qvalues, pi0_hat = storey_qvalues(pvalues)

print(f"추정된 π₀: {pi0_hat:.4f} (실제: {m0/m:.2f})")
print(f"q-값 ≤ 0.05인 가설 수: {(qvalues <= 0.05).sum()}")

# BH 조정 p-값과 비교
from statsmodels.stats.multitest import multipletests
_, bh_pvals, _, _ = multipletests(pvalues, alpha=0.05, method='fdr_bh')
print(f"BH 조정 p-값 ≤ 0.05인 가설 수: {(bh_pvals <= 0.05).sum()}")

9.2 Step 2: R의 qvalue 패키지 (실무 활용)

# install.packages("BiocManager")
# BiocManager::install("qvalue")
library(qvalue)

# 시뮬레이션 데이터
set.seed(42)
m <- 1000; m0 <- 800
pvalues <- c(runif(m0), rbeta(m - m0, 0.3, 5))

# q-값 계산
qobj <- qvalue(p = pvalues)

cat("추정된 π₀:", qobj$pi0, "\n")
cat("q-값 ≤ 0.05인 가설 수:", sum(qobj$qvalues <= 0.05), "\n")

# π₀ 추정 그래프
plot(qobj)

# 요약
summary(qobj)

9.3 Step 3: BH vs Storey q-값 비교 시각화

import matplotlib.pyplot as plt
import numpy as np

rng = np.random.default_rng(0)
m = 500
pi0_true = 0.7
m0 = int(m * pi0_true)
pvalues = np.sort(np.concatenate([rng.uniform(0, 1, m0), rng.beta(0.2, 5, m - m0)]))

qvals, pi0_hat = storey_qvalues(pvalues)

# BH 조정 p-값 (π₀=1 가정)
bh_vals = np.minimum.accumulate((m * pvalues / np.arange(1, m+1))[::-1])[::-1]
bh_vals = np.clip(bh_vals, 0, 1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# q-값 vs BH 비교
ax1.plot(pvalues, bh_vals, label=f"BH 조정 p-값 (π₀=1)", color="darkorange")
ax1.plot(pvalues, qvals, label=f"Storey q-값 (π₀̂={pi0_hat:.2f})", color="steelblue")
ax1.axhline(0.05, color="red", linestyle="--", label="q=0.05 임계값")
ax1.set_xlabel("정렬된 p-값")
ax1.set_ylabel("q-값 / BH 조정 p-값")
ax1.set_title("BH vs Storey q-값 비교")
ax1.legend(fontsize=8)
ax1.set_xlim(0, 0.15)

# π₀ 추정 과정
lambdas = np.arange(0.05, 0.96, 0.05)
pi0_lam = np.array([np.sum(pvalues > lam) / (m * (1-lam)) for lam in lambdas])
ax2.scatter(lambdas, pi0_lam, color="steelblue", label=r"$\hat{\pi}_0(\lambda)$")
ax2.axhline(pi0_true, color="green", linestyle="--", label=f"실제 π₀={pi0_true}")
ax2.axhline(pi0_hat, color="red", linestyle=":", label=f"추정 π₀={pi0_hat:.3f}")
ax2.set_xlabel(r"$\lambda$")
ax2.set_ylabel(r"$\hat{\pi}_0(\lambda)$")
ax2.set_title(r"$\pi_0$ 추정 과정")
ax2.legend(fontsize=8)

plt.tight_layout()
plt.savefig("qvalue_comparison.png", dpi=120)
plt.show()

10 관련 주제

선행 지식

후속 주제

Subscribe

Enjoy this blog? Get notified of new posts by email: