목차
데이터 준비
한 집단의 평균값과 다른 집단의 평균값의 차이가 명백하다는 점을 최대한 과학적인 방법으로 누군가에게 납득시켜야 할 때가 있습니다. 아래에서는 빈도주의자(Frequentist)의 관점과 베이지안(Bayesian)의 과점에서 t-검정을 수행하는 방법을 살펴보고자 합니다.
데이터는 Kruschke, John. (2012) Bayesian estimation supersedes the t-test의 자료를 사용하도록 하겠습니다. 참고로 아래 데이터 준비 코드와 베이지안 t-검정 항목의 코드는 pymc3 튜토리얼의 코드를 참고했습니다.
drug = [101,100,102,104,102,97,105,105,98,101,100,123,105,103,100,95,102,106,
109,102,82,102,100,102,102,101,102,102,103,103,97,97,103,101,97,104,
96,103,124,101,101,100,101,101,104,100,101]
placebo = [99,101,100,101,102,100,97,101,104,101,102,102,100,105,88,101,100,
104,100,100,100,101,102,103,97,101,101,100,101,99,101,100,100,
101,100,99,101,100,102,99,100,99]
y1 = np.array(drug)
y2 = np.array(placebo)
y = pd.DataFrame(dict(value=np.r_[y1, y2], group=np.r_[['drug']*len(drug), ['placebo']*len(placebo)]))
2
3
4
5
6
7
8
9
10
빈도주의 t-검정
통계학 공부를 하며 접했던 다양한 자료들을 오랜만에 다시 찾아보고 나름대로 아래와 같이 정리를 해보았습니다.
t-검정(t-test)
종속변수가 연속형이고 독립변수는 단독 범주형일 때 2개의 서로 다른 수준에 대해 짝을 이루지 않거나(독립표본 t-검정), 짝을 이루는(대응표본 t-검정) 표본들이 기본 가정(독립표본 t-검정의 경우 독립성/정규성/등분산성, 대응표본 t-검정의 경우 정규성)을 만족한다는 조건 하에 각각의 수준들이 유의적인 차이를 보이는지를 검정.
독립표본 t-검정
독립성
독립변수에 따른 종속변수가 각 수준별로 짝을 이루고 있지 않아야 한다는 조건입니다(예를 들어 A 투약전/투약후, B 투약전/투약후...가 아닌 투약군 AB.../미투약군 CD...인 경우).
정규성
독립변수에 따른 종속변수는 정규분포를 만족해야 한다는 조건입니다. 주로 콜모고로프-스미노프(Kolmogorov-Smirnov) 검정을 통해 이를 밝혀내며, 파이썬에서는 statsmodels의 statsmodels.stats.diagnostic.kstest_normal 메소드를 사용할 수 있습니다.
sm.stats.diagnostic.kstest_normal(drug + placebo, dist='norm', pvalmethod='approx')
# (0.22204128505832343, 9.700938640868548e-12)
2
3
TIP
콜모고로프-스미노프 검정은 정규분포 뿐만 아니라 수준별 표본들이 동일 분포를 따르는지를 검정하는데 사용할 수 있습니다.
TIP
여기선 statsmodels의 convenience method를 이용하였지만 scipy의 scipy.stats.kstest, scipy.stats.ks_2samp와 같은 메소드들을 사용할 수도 있습니다. 이 경우 정규분포를 fitting한 다음 cdf를 해당 메소드들의 argument롤 넘겨야 하는 등의 추가작업이 필요할 수 있습니다.
등분산성
독립변수에 따른 종속변수는 각 수준별로 동일해야 한다는 조건(Homoscedasticity)입니다. 레빈(Levene) 검정을 이용하며, statsmodels.stats.weightstats.test_equal_var 메소드를 지원할 예정이지만, 현재는 2d에 대한 weight 구현 문제로 scipy의 scipy.stats.levene 메소드를 사용해야 합니다.
sp.stats.levene(drug, placebo)
# LeveneResult(statistic=4.542872436860705, pvalue=0.03587638864722)
2
3
WARNING
레빈 검정 결과 등분산성을 만족하지 못하며 따라서 최종적으로 독립표본 t-검정을 수행할 때 아래와 같이 uservar argument에 'unequal' 속성을 부여해야 합니다. 참고로 '크다' 또는 '작다'의 판단이 아닌 '다르다'를 판단해야 하기 때문에 양측검정을 해야 하고 따라서 alternative argument는 'two-sided'가 되어야 합니다.
따라서 최종 독립표본 t-검정 코드는 아래와 같습니다. statsmodels.stats.weightstats.ttest_ind를 참고하시기 바랍니다.
# cf. Levene's test p-value < 0.05, H0: Equal Variance rejected. usevar should be 'pooled'.
sm.stats.ttest_ind(drug, placebo, alternative='two-sided', usevar='unequal', value=0)
# (1.6221904572902277, 0.10975381983712841, 63.03889717357649)
2
3
4
p-value가 0.1098로 어떤 기준을 따르느냐에(0.05? 0.10?) 따라 조금 애매하다면 애매할 수 있는 수치입니다.
대응표본 t-검정
정규성 검정
정규성 검정은 위에서 말씀드린 독립표본 t-검정의 경우와 동일합니다.
대응표본 t-검정의 경우 scipy의 scipy.stats.ttest_rel 메소드를 사용하면 됩니다.
베이지안 t-검정
베이지안 데이터 분석
베이지안 데이터 분석이 익숙하지 않으신 분들께는 베이지안 데이터 분석 포스팅에 소개한 세 개의 영상을 먼저 보실 것을 추천드립니다.
독립표본 t-검정
기본 가정은 앞서 살펴본 내용과 동일하며, 아래 코드를 보시면 스튜던트 t-분포의 파라메터가 되는 group1과 group2의 평균, 표준편차, 그리고 v에 대한 사전분포(믿음)들과 이에 따른 결정론적 사항들에 대해 정의한 내용을 보실 수 있습니다.
참고로 아래 코드의 작성자는 group1, group2 각각의 평균은 샘플 전체의 평균과 표준편차의 2배를 따르는 정규분포의 어딘가에 있을 것이라고 믿고 있고(정규분포이므로 봉우리 부분에 더 큰 믿음), 각각의 표준편차는 1부터 10까지 중 어딘가에 있을 것 같지만 그 중 어느 부분에 있을 확률이 더 높은지는 전혀 모르겠다고(균등분포) 믿고 있습니다.
μ_m = y.value.mean()
μ_s = y.value.std() * 2
σ_low = 1
σ_high = 10
with pm.Model() as model:
group1_mean = pm.Normal('group1_mean', μ_m, sd=μ_s)
group2_mean = pm.Normal('group2_mean', μ_m, sd=μ_s)
group1_std = pm.Uniform('group1_std', lower=σ_low, upper=σ_high)
group2_std = pm.Uniform('group2_std', lower=σ_low, upper=σ_high)
λ1 = group1_std**-2
λ2 = group2_std**-2
ν = pm.Exponential('ν_minus_one', 1/29.) + 1
diff_of_means = pm.Deterministic('difference of means', group1_mean - group2_mean)
diff_of_stds = pm.Deterministic('difference of stds', group1_std - group2_std)
effect_size = pm.Deterministic('effect size', diff_of_means / np.sqrt((group1_std**2 + group2_std**2) / 2))
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
그 다음 Likelihood를 정의합니다. 여기서 MCMC(Markov Chain Monte Carlo) 방법을 통해 사전분포를 업데이트 하게 됩니다. 단순화시켜 말씀드리면 y1, y2에서 각각 관측된 표본들이(observed=yn) group1, group2 각각의 파라미터들에 대한 사전 믿음을 바탕으로 모델링한 분포와(여기서는 스튜던트-t) 일관된다면 이를 Accept하고 아니라면 Reject하는 과정을 거치며 각자의 목표 확률분포를 근사하게 되는 것입니다.
...(constructing a Markov chain to do Monte Carlo approximation).
MCMC는 메트로폴리스-헤이스팅스, 깁스 샘플링 등의 랜덤워크 몬테카를로 알고리즘들을 사용합니다. 더 자세한 내용은 별도의 포스팅을 통해 다루어 보도록 하겠습니다. 다만 궁금하신 분들을 위해 차분하게 잘 설명된 자료를 남겨두자면 퀀토피안의 데이터 사이언티스트인 Thomas Wiecki의 MCMC sampling for dummies를 강력 추천드리고 싶습니다.
with model:
group1 = pm.StudentT('drug', mu=group1_mean, lam=λ1, nu=ν, observed=y1)
group2 = pm.StudentT('placebo', mu=group2_mean, lam=λ2, nu=ν, observed=y2)
2
3
이제 실제로 2000회의 샘플링을 수행합니다.
with model:
trace = pm.sample(2000, cores=2)
2
아래는 Deterministic으로 정의한 부분에 대한 결과입니다. 물론 group1_mean, group2_mean, group1_std, group2_std, v_minus_one도 varnames 부분을 바꿔주기만 하면 아래와 유사한 방법으로 확인이 가능합니다.
pm.plot_posterior(trace, varnames=['difference of means','difference of stds', 'effect size'],
ref_val=0,
color='#87ceeb');
2
3
이제 우리는 경영진의 의문에 대해 이를 보다 효과적으로 설명하고, 쉽게 이해시키며, 어쩌면 설득 할 수도 있습니다. p-value가 뭔지도 모르는 사람에게 단순히 p-value가 얼마 이하라서 다르다, test 결과가 그렇다고 말 하는 것 보다, estimation의 결과를 보여주면서 우리의 사전 믿음이 어떠했고 사후 믿음은 얼마나 달라졌는지에 대한 정보도 함께 제공할 수 있는 보다 informative한 기법이기 때문입니다(인식론적 불확실성과 무작위적 불확실성의 구분).
Why
https://docs.pymc.io/notebooks/BEST.html
A more informative and effective approach for comparing groups is one based on estimation rather than testing, and is driven by Bayesian probability rather than frequentist. That is, rather than testing whether two groups are different, we instead pursue an estimate of how different they are, which is fundamentally more informative. Moreover, we include an estimate of uncertainty associated with that difference which includes uncertainty due to our lack of knowledge of the model parameters (epistemic uncertainty) and uncertainty due to the inherent stochasticity of the system (aleatory uncertainty).