[파이썬 기초] 통계 기초 (scipy)

SciPy는 과학 계산을 위한 파이썬 라이브러리입니다.

그 중 scipy.stats 모듈은 통계 검정, 확률 분포, 상관 분석 등을 제공합니다.

생물정보학에서는 발현량 차이 검정, p-value 계산, 다중 검정 보정 등에 자주 사용합니다.

pip install scipy

---

기술 통계

분포의 기본적인 특성을 수치로 요약합니다.

import numpy as np
from scipy import stats

np.random.seed(42)
coverage = np.random.normal(30, 8, 100)

# 기본 기술 통계
print(coverage.mean())      # 평균
print(np.median(coverage))  # 중앙값
print(coverage.std())       # 표준편차
print(coverage.var())       # 분산

# scipy로 한꺼번에 확인
result = stats.describe(coverage)
print(result)
# DescribeResult(nobs=100, minmax=(...), mean=..., variance=..., skewness=..., kurtosis=...)
print(f"왜도: {result.skewness:.4f}")   # 분포의 비대칭 정도
print(f"첨도: {result.kurtosis:.4f}")   # 분포의 뾰족한 정도

# 사분위수
q1, q2, q3 = np.percentile(coverage, [25, 50, 75])
iqr = stats.iqr(coverage)    # 사분위 범위 (Q3 - Q1)
print(f"Q1: {q1:.2f}, Q2: {q2:.2f}, Q3: {q3:.2f}, IQR: {iqr:.2f}")

---

정규성 검정

많은 통계 검정이 데이터가 정규분포를 따른다는 가정을 합니다. 검정 전에 정규성을 확인하는 것이 좋습니다.

Shapiro-Wilk 검정

샘플 크기가 작을 때(n < 50) 적합합니다.

stat, p_value = stats.shapiro(coverage[:50])
print(f"통계량: {stat:.4f}, p-value: {p_value:.4f}")

if p_value > 0.05:
    print("정규분포를 따른다고 볼 수 있음 (귀무가설 기각 실패)")
else:
    print("정규분포를 따르지 않음 (귀무가설 기각)")

Kolmogorov-Smirnov 검정

샘플 크기가 클 때 사용합니다.

stat, p_value = stats.kstest(coverage, "norm",
                              args=(coverage.mean(), coverage.std()))
print(f"KS 통계량: {stat:.4f}, p-value: {p_value:.4f}")

D'Agostino-Pearson 검정

stat, p_value = stats.normaltest(coverage)
print(f"통계량: {stat:.4f}, p-value: {p_value:.4f}")

---

가설 검정

귀무가설(H₀): 차이가 없다 / 대립가설(H₁): 차이가 있다
p-value < 0.05이면 귀무가설을 기각하고 통계적으로 유의하다고 판단하는 것이 일반적입니다.

단일 표본 t-검정

표본 평균이 특정 값과 차이가 있는지 검정합니다.

np.random.seed(42)
coverage = np.random.normal(28, 8, 50)

# 커버리지 평균이 30과 다른가?
stat, p_value = stats.ttest_1samp(coverage, popmean=30)
print(f"t-통계량: {stat:.4f}, p-value: {p_value:.4f}")

독립 표본 t-검정 (two-sample t-test)

두 독립 집단의 평균 차이를 검정합니다.

np.random.seed(42)
group_a = np.random.normal(30, 6, 30)    # 처리군
group_b = np.random.normal(25, 7, 30)    # 대조군

# 등분산 가정
stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=True)
print(f"t-통계량: {stat:.4f}, p-value: {p_value:.4f}")

# Welch's t-test: 등분산 가정하지 않음 (더 일반적으로 권장)
stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"Welch t-통계량: {stat:.4f}, p-value: {p_value:.4f}")

대응 표본 t-검정 (paired t-test)

같은 대상에서 처리 전후를 비교할 때 사용합니다.

np.random.seed(42)
before = np.random.normal(20, 5, 20)
after = before + np.random.normal(3, 2, 20)    # 처리 후 평균 3 증가

stat, p_value = stats.ttest_rel(before, after)
print(f"t-통계량: {stat:.4f}, p-value: {p_value:.4f}")

등분산 검정 (Levene's test)

두 집단의 분산이 같은지 확인합니다. t-검정 전에 수행합니다.

stat, p_value = stats.levene(group_a, group_b)
print(f"Levene 통계량: {stat:.4f}, p-value: {p_value:.4f}")

if p_value > 0.05:
    print("등분산 가정 성립 → ttest_ind(equal_var=True) 사용")
else:
    print("등분산 가정 불성립 → Welch's t-test 사용")

---

비모수 검정

데이터가 정규분포를 따르지 않거나 샘플 크기가 작을 때 사용합니다.

Mann-Whitney U 검정

독립 표본 t-검정의 비모수 대안입니다.

stat, p_value = stats.mannwhitneyu(group_a, group_b, alternative="two-sided")
print(f"U 통계량: {stat:.4f}, p-value: {p_value:.4f}")

Wilcoxon 부호순위 검정

대응 표본 t-검정의 비모수 대안입니다.

stat, p_value = stats.wilcoxon(before, after)
print(f"통계량: {stat:.4f}, p-value: {p_value:.4f}")

Kruskal-Wallis 검정

세 집단 이상의 비교에 사용합니다. 일원분산분석(ANOVA)의 비모수 대안입니다.

np.random.seed(42)
group_c = np.random.normal(35, 5, 30)

stat, p_value = stats.kruskal(group_a, group_b, group_c)
print(f"H 통계량: {stat:.4f}, p-value: {p_value:.4f}")

---

상관 분석

두 변수 간의 선형 관계를 수치로 표현합니다.

Pearson 상관계수

정규분포를 따르는 연속형 변수에 적합합니다.

np.random.seed(42)
coverage = np.random.normal(30, 8, 50)
gc_content = 0.3 * coverage / 100 + np.random.normal(0.5, 0.03, 50)

r, p_value = stats.pearsonr(coverage, gc_content)
print(f"Pearson r: {r:.4f}, p-value: {p_value:.4f}")

Spearman 상관계수

순위 기반 상관계수로 비모수적 방법입니다. 정규분포 가정 없이 사용할 수 있습니다.

rho, p_value = stats.spearmanr(coverage, gc_content)
print(f"Spearman ρ: {rho:.4f}, p-value: {p_value:.4f}")

---

다중 검정 보정

여러 유전자의 발현량을 동시에 검정할 때처럼, 검정을 많이 수행할수록 우연히 유의하게 나오는 경우(위양성)가 증가합니다. 이를 보정하기 위해 다중 검정 보정을 적용합니다.

from statsmodels.stats.multitest import multipletests

# 예시: 100개 유전자의 p-value
np.random.seed(42)
p_values = np.random.uniform(0, 1, 100)
p_values[:10] = np.random.uniform(0, 0.01, 10)    # 일부는 실제 유의미한 값

# Benjamini-Hochberg FDR 보정 (생물정보학에서 가장 많이 사용)
reject, p_adjusted, _, _ = multipletests(p_values, method="fdr_bh")
print(f"보정 전 유의한 수: {(p_values < 0.05).sum()}")
print(f"보정 후 유의한 수: {reject.sum()}")

# Bonferroni 보정 (더 보수적)
reject_bf, p_adj_bf, _, _ = multipletests(p_values, method="bonferroni")
print(f"Bonferroni 보정 후 유의한 수: {reject_bf.sum()}")

multipletests는 statsmodels 패키지에 포함되어 있습니다.

pip install statsmodels

---

확률 분포

scipy.stats는 다양한 확률 분포 객체를 제공합니다.

# 정규분포
normal = stats.norm(loc=0, scale=1)    # 평균 0, 표준편차 1
print(normal.pdf(0))      # 확률밀도함수 (PDF)
print(normal.cdf(1.96))   # 누적분포함수 (CDF): P(X <= 1.96)
print(normal.ppf(0.975))  # 분위수 (퍼센타일)

# p-value에서 z-score 임계값 확인
z_crit = stats.norm.ppf(0.975)    # 양측 5% 유의수준의 z값
print(f"z 임계값: {z_crit:.4f}")   # 1.96

# t-분포 (자유도 지정)
t_dist = stats.t(df=29)
print(t_dist.ppf(0.975))    # 자유도 29에서 양측 5% 유의수준의 t값

# 이항분포
binom = stats.binom(n=100, p=0.5)
print(binom.pmf(50))     # P(X = 50)
print(binom.cdf(50))     # P(X <= 50)

---

검정 선택 가이드

상황	정규분포 따를 때	정규분포 아닐 때
단일 집단 vs 기준값	단일 표본 t-검정	Wilcoxon 부호순위 검정
두 독립 집단 비교	독립 표본 t-검정	Mann-Whitney U 검정
대응 표본 비교	대응 표본 t-검정	Wilcoxon 부호순위 검정
세 집단 이상 비교	일원분산분석 (ANOVA)	Kruskal-Wallis 검정
두 변수 간 관계	Pearson 상관	Spearman 상관