통계는 실전이야: ANOVA(feat. Tukey HSD)

오늘 해볼 건 어노바, 그러니까 분산분석입니다. 이건 뭐 없고 얘네'들' 중 뭔가 다른게 있나를 보는겁니다. 얘네'들'? 예. 보통 집단이 세개 이상이면 어노바 하고요, 비모수로 빠질거면 크러스칼 월리스(윌리스 아니고) 검정이라고 있어요. 둘 다 구체적으로 뭐가 다른지 보려면 후속 절차도 들어가야 하는데 어노바는 튜키(Tukey HSD)랑 놀고 크러스칼은 듄 테스트(Dunn's test)랑 놉니다.

---

일원분산분석

그 어노바가 독립변수 개수에 따라서 일원, 이원, 삼원(3개), 다변량(...) 이렇게 있습니다. 근데 나도 3원 이상은 안해봄... 여기서는 일원이랑 이원만 할겁니다. 예. 아무튼 일원은 뭐냐... 요인이 하나다 이겁니다. 예를 들자면 컨트롤, ACC(에틸렌 전구체), 에틸렌을 투여했을 때 식물체의 길이같은 거 말이다. 호르몬이 두 종류 아니냐고? 하나가 전구체임다. 그니까 식물이 빨아먹고 에틸렌으로 만드는거. 아니면 옥신을 농도별로 본다던가?

 

# 예시 데이터
fert_a = np.random.normal(10, 2, 30)
fert_b = np.random.normal(10, 2, 30)
fert_c = np.random.normal(13, 2, 30)

일원분산분석 예시는 비료다. 비료를 주고 식물이 얼마나 크는지를 봤다 이거죠.

 

어노바의 귀무가설은 '얘네들 다 똑같다(=또이또이, 쌤쌤똔똔)'고 대립가설은 '얘네들 중 적어도 하나는 다르다'이다. 이거 잘 기억하세요. 그래야 튜키 왜 하는지도 이해가 됨.

# ANOVA
f_stat, p_val = stats.f_oneway(fert_a, fert_b, fert_c)

print(f'p-value: {p_val:.4e}')

if p_val > 0.05:
    print('귀무가설이 기각되지 않았습니다. ')
else: 
    print('튜키 드가자!!!')

p-value: 8.0855e-07
튜키 드가자!!!

P-value가 말하고 있습니다. 얘네들 중 뭔가 다른 놈이 있다고.

 

# 데이터프레임 생성
df_anova = pd.DataFrame({
    # 성장 데이터 합치기
    'Growth' : np.concatenate([fert_a, fert_b, fert_c]),
    # 그룹
    'Fertilizer' : ['A'] * 30 + ['B'] * 30 + ['C'] * 30
})

튜키 할라고 묶었음... 음 사실 그것때문에 묶은건 아니고, 빡스플롯 만들라고 묶은거긴 한데 어쨌든 튜키도 이걸로 볼겁니다.

 

Tukey HSD

tukey_fert = pairwise_tukeyhsd(endog=df_anova['Growth'],      # 결과값 (수치형)
                            groups=df_anova['Fertilizer'],  # 비교 집단 (범주형)
                            alpha=0.05)

print(tukey_fert)

Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     A      B   0.4804 0.6782 -0.8804 1.8412  False
     A      C   3.0533    0.0  1.6926 4.4141   True
     B      C    2.573 0.0001  1.2122 3.9338   True
---------------------------------------------------

다른거 보지 말고 reject가 True인 것만 보십쇼. reject가 True인 애들이 '까봤더니 얘네들이 다르네?'다.

 

fig = tukey_fert.plot_simultaneous(figsize=(10, 6))
plt.title("Tukey HSD: Fertilizer comparison")
plt.axvline(np.mean(fert_c), linestyle='--')
plt.xlabel("Growth")
plt.show()

이거는 뭐 하실 분들만 하십쇼. 튜키도 시각화가 됩니다.

 

이원분산분석

# 데이터 생성
data = pd.DataFrame({
    'Drug': ['A']*20 + ['B']*20,
    'Exercise': (['High']*10 + ['Low']*10) * 2,
    'Weight_Loss': [
        # Drug A: High(7~9), Low(3~5)
        *np.random.normal(8, 1, 10), *np.random.normal(4, 1, 10),
        # Drug B: High(6~8), Low(5~7)
        *np.random.normal(7, 1, 10), *np.random.normal(6, 1, 10)
    ]
})

요인이 두개지요. 예시 데이터에서는 운동과 보조제로 나뉘어져있는데, 뭐 풀때기를 옥신이랑 비료로 하던가, 여러가지 있다. 예전에 캐글 데이터로 팀플했던 것들 중에는 고객 가입률을 연락 빈도+연령대 두개로 묶어서 본 것도 있었다.

 

model = ols('Weight_Loss ~ C(Drug) * C(Exercise)', data=data).fit()
anova_table = sm.stats.anova_lm(model, typ=2)

print(anova_table)

                        sum_sq    df          F        PR(>F)
C(Drug)               4.048614   1.0   4.530293  4.021375e-02
C(Exercise)          50.983335   1.0  57.049015  6.261649e-09
C(Drug):C(Exercise)  23.696062   1.0  26.515272  9.513253e-06
Residual             32.172335  36.0        NaN           NaN

아 위에꺼랑 다르다고요? 둘 중 하나 암거나 쓰시면 됩니다. 아마 scipy.stats에서도 이원분석 될거예요. 아무튼 저기러는 뭘 봐야 하나요? PR(>F)가 피밸류다. 4쩜얼마에 -02 붙었으면 어쨌든 0.05보다는 작은거니까 후속분석 들어가자.

 

Tukey HSD

# 1. '약 + 운동' 조합 컬럼 만들기 (사후검정을 위해 그룹을 하나로 합침)
data['Combination'] = data['Drug'] + " / " + data['Exercise']

# 2. Tukey HSD 실행
tukey = pairwise_tukeyhsd(endog=data['Weight_Loss'],     # 종속변수
                            groups=data['Combination'],   # 그룹화 변수
                            alpha=0.05)                   # 유의수준

print("=== Tukey HSD 사후검정 결과 ===")
print(tukey)

=== Tukey HSD 사후검정 결과 ===
  Multiple Comparison of Means - Tukey HSD, FWER=0.05   
========================================================
 group1   group2  meandiff p-adj   lower   upper  reject
--------------------------------------------------------
A / High  A / Low  -3.7973    0.0 -4.9359 -2.6587   True
A / High B / High  -0.9031 0.1612 -2.0417  0.2356  False
A / High  B / Low  -1.6217 0.0026 -2.7603  -0.483   True
 A / Low B / High   2.8942    0.0  1.7556  4.0329   True
 A / Low  B / Low   2.1756 0.0001   1.037  3.3143   True
B / High  B / Low  -0.7186 0.3387 -1.8572    0.42  False
--------------------------------------------------------

얘는 요인 두 개를 묶어서 본 거고, 요인별로 따로 보는 방법도 있긴 있다.

 

# 운동
print("=== 운동량별 차이 (Tukey HSD) ===")
tukey_exec = pairwise_tukeyhsd(endog=data['Weight_Loss'], groups=data['Exercise'], alpha=0.05)
print(tukey_exec)

# 보조제
print("\n=== 보조제별 차이 (Tukey HSD) ===")
tukey_drug = pairwise_tukeyhsd(endog=data['Weight_Loss'], groups=data['Drug'], alpha=0.05)
print(tukey_drug)

=== 운동량별 차이 (Tukey HSD) ===
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj  lower   upper  reject
---------------------------------------------------
  High    Low  -2.2579   0.0 -3.0618 -1.4541   True
---------------------------------------------------

=== 보조제별 차이 (Tukey HSD) ===
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     A      B   0.6363 0.2376 -0.4372 1.7098  False

롸? 둘이 뭔 차이예요? 위는 교호작용(상호작용)도 보는거고 아래는 다른 요인 떼놓고 개별로 보는거다. 그러니까 위는 보조제와 운동의 상호작용을 보는거고, 아래는 운동만+보조제만 보는거다. 이 케이스는 상호작용이랑 두 요인 다 P-value가 유의수준(0.05)보다 작게 나왔으니 둘 다 봐야 하는데, 캐글 은행 데이터는 아래 것만 봤었다.

 

롸? 왜요? 캐글 그 은행 데이터에서 이원분석 요소가 고객 연령대랑 연락 빈도였거든요? 근데 그래프를 그려보니까 연령대랑 상관없이 연락을 많이 하면 가입률이 떨어지는겁니다. 그니까 우리는 그 양상을 시각화를 통해 확인했고, 통계적으로 그 차이가 유의한지를 보기 위해서 어노바를 한거거든요. 그러니까 시각화 해보니까 양상이 이런데, 진짜 그런지를 확인해본거다 이겁니다.

 

# Tukey 결과 시각화
fig = tukey.plot_simultaneous()
plt.title('Tukey HSD Multiple Comparison')
plt.xlabel('Weight Loss (kg)')
plt.show()

이것도 옵션이긴 한데 시각화고요

 

Interaction plot

# 폰트는 이미 NanumSquare로 설정되어 있으니 바로 그립니다.
fig = interaction_plot(x=data['Exercise'], 
                        trace=data['Drug'], 
                        response=data['Weight_Loss'],
                        colors=['red', 'blue'], 
                        markers=['D', '^'])

plt.title('Weight Loss Interaction: Drug & Exercise')
plt.show()

이건 또 머시여… 나도 처음봐요 이거. 그 상호작용 여부를 나타낸건데 x축이 운동 강도고 세로축이 체중 감량, 두 선은 보조제거든요? 자, 선을 보니까 교차하고 있습니다. 네. 이건 운동량과 보조제가 상호작용을 하고 있다는 의미입니다. 아니 그럼 상호작용을 안 하면 어떻게 되는데요? 등호마냥 평행선 되죠.

 

아까 어노바랑 튜키를 통해서 확인한 결과가 저기 다 들어있습니다. 운동 빡시게 하는 집단은 어떤 보조제를 쓰더라도 비슷비슷해서 reject가 False였고, 운동 저강도로 하는 집단은 보조제에 따라 유의한 체중 감량 효과가 있었기때문에 reject가 True였죠? 보조제 A는 빡시게 하는 그룹, 저강도로 하는 그룹간에 차이가 있었지만 보조제 B는 그렇지 못해서 reject가 False가 떴고요. 그런겁니다.

 

근데 이원분석 해보셨다면서 저건 왜 처음봐요? 위에도 썼지만 그건 시각화를 통해서 이미 아 그렇군을 확인하고 진짜 통계적으로 그 차이가 유의한지를 본 거라 그렇습니다.