사후분석(Post-hoc test)
관련된 내용
- 분산분석의 개요
- 일원분산분석(one-way ANOVA)
- 사후분석(Post-hoc test)
- 이원분산분석(two-way ANOVA)
통계 분석은 데이터의 정규성, 등분산성, 그리고 독립성등의 기본가정의 충족 여부에 결과의 합리성을 확보할 수 있습니다. 또한 분산분석의 경우 다중 그룹의 비교이므로 각 그룹간의 효과를 분석할 필요가 존재합니다. 이러한 분석은 분산분석 모형 구축과는 별도로 진행하므로 사후분석이라 합니다. 사후분석에서 가정의 적합성을 검정하기 위해 적용할 수 있는 방법이나 함수는 다음과 같습니다.
- 정규성 : 반응변수는 정규분포를 따릅니다.
- 등분산성
정규성
분석 결과에 대한 신뢰는 데이터가 통계적 테스트의 기본 가정을 만족시키는 정도에 달려 있습니다. 일원 분산 분석에서 각 그룹의 모집단은 정규 분포를 따르고 동일한 분산을 가정합니다. 그러나 모집단의 정규성을 검정하는 것은 어려울 수 있습니다. 대신에 모델에 의한 잔차의 정규성 검정으로 대신합니다. Q-Q 플롯을 사용하여 정규성 가정을 시각적으로 평가할 수 있습니다.
예1)
일정기간의 코스피(kos), 코스탁(kq), 다우존스 주가지수(dj) 그리고 원-달러(WonDol)의 일일 변화량에 대한 분산분석에 대한 사후분석을 진행합니다.
다음의 코드에 의해 호출할 수 있으며 각 자료의 Open과 Close 사이의 일일 변화율을 계산한 자료를 독립변수와 의존변수(일원 분산분석을 위한 자료구조 참조)로 구분합니다.
st=pd.Timestamp(2024,1, 1)
et=pd.Timestamp(2024, 5, 30)
code=["KS11", "KQ11", "DJI", "USD/KRW"]
nme=['kos','kq','dj','WonDol']
da=pd.DataFrame()
for i in code:
x=fdr.DataReader(i,st, et)['Close']
x1=x.pct_change()
da=pd.concat([da, x1], axis=1)
da.columns=nme
da.index=range(len(da))
da1=da.dropna()
da1.head(2).round(3)
da2=pd.melt(da1, value_vars=['kos', 'kq', 'dj', 'WonDol'], var_name="idx", value_name="val")
da2.head(3).round(3)
| idx | val | |
|---|---|---|
| 0 | kos | -0.023 |
| 1 | kos | -0.008 |
| 2 | kos | -0.003 |
위 자료에 대한 분산분석 모델은 다음과 같습니다.
model=ols("val~C(idx)", data=da2).fit()
sm.stats.anova_lm(model).round(3)
| df | sum_sq | mean_sq | F | PR(>F) | |
|---|---|---|---|---|---|
| C(idx) | 3.0 | 0.000 | 0.0 | 0.222 | 0.881 |
| Residual | 376.0 | 0.027 | 0.0 | NaN | NaN |
model=ols("val~C(idx)", data=da2).fit()
구축한 모델에 의한 잔차의 정규성 검정 결과입니다. 그림 1은 scipy.stats.probplot() 함수를 적용하여 시각적으로 나타낼 수 있습니다.
res=model.resid res[:5]
0 -0.0232 1 -0.0076 2 -0.0033 3 -0.0038 4 -0.0024 dtype: float64
plt.figure(figsize=(4, 3)) varQQplot1=stats.probplot(res, plot=plt)
그림 1로 모델의 잔차는 정규분포를 나타내는 이론적 분위수로부터의 큰 이탈을 보이지 않습니다. 이점을 정량적으로 나타내기 위해 AD검증을 실행합니다.
stat, critic, plev=stats.anderson(res)
print(f"통계량:{stat.round(4)}\n임계값:{np.around(critic,4)}\n유의수준:{np.around(plev, 4)}")
통계량:1.8528 임계값:[0.57 0.649 0.779 0.909 1.081] 유의수준:[15. 10. 5. 2.5 1. ]
생성된 모델에 의한 예측값과 실측값의 차이인 잔차들의 분포는 정규성에 부합한다는 귀무가설을 기각할 수 있습니다.
이상값 검정
분산 방법론의 분석은 이상값의 존재에 민감할 수 있습니다. 위에서 소개한 것과 같이 IQR을 사용하여 이상값을 결정할 수 있습니다. 또한 statsmodels.stats.outliers_influence 클래스의 outlier_test() 함수를 사용하여 Bonferroni 방법으로 이상치를 결정합니다. 이 함수의 인수는 분산분석 또는 회귀분석 모델입니다. Studentized residuals, the unadjusted p-value, 그리고 Bonferroni 방법에 의해 수정된 corrected p-value을 반환합니다. 이 수정된 p-value 값이 1이하 값이 되면 이상치로 결정할 수 있습니다.
Bonferroni correction
가설 검정에서 유의 수준은 귀무가설의 기각 여부를 판단하는 기준입니다. 이 기준은 가설이 참인 경우 거짓으로 판단할 오류 즉, α오류를 포함합니다. 분석 대상이 단일한 경우 예를 들어 그룹 A, B의 평균이 같다는 가설인 경우 유의수준 0.05를 설정한 경우 단일한 분석으로 α오류는 1 - 0.95 = 0.05가 됩니다. 그러나 A = B = C인 경우를 분산분석이 아닌 개별적 분석을 할 경우 A, B와 B, C 그리고 C, A의 3번의 개별적인 비교를 시행해야 합니다. 이 경우의 α오류는 1 - (0.95)3 = 0.143으로 증가하게 됩니다. 이러한 오류를 보정하기 위해 유의 수준을 factor의수로 나누어 주는 방법이 사용됩니다. 이러한 유의수준의 조정방법을 Bonferroni 교정이라 합니다.
round(1-(1-0.05/3)**3, 3)
0.049
from statsmodels.stats.outliers_influence import outlier_test outl=outlier_test(model) outl.head(3)
| student_resid | unadj_p | bonf(p) | |
|---|---|---|---|
| 0 | -2.765 | 0.006 | 1.0 |
| 1 | -0.896 | 0.371 | 1.0 |
| 2 | -0.386 | 0.700 | 1.0 |
out2=np.where(outl.iloc[:,2]<1)[0] print(out2)
[ 21 164]
위 결과의 감지된 이상값를 제외한 데이터에 대해 새로운 모델을 작성하였습니다.
da3=da2.drop(out2) da3.shape
(378, 2)
model1=ols("val~C(idx)", data=da3).fit()
sm.stats.anova_lm(model1).round(3)
| df | sum_sq | mean_sq | F | PR(>F) | |
|---|---|---|---|---|---|
| C(idx) | 3.0 | 0.000 | 0.0 | 0.447 | 0.719 |
| Residual | 374.0 | 0.026 | 0.0 | NaN | NaN |
res1=model1.resid
stat1, critic1, plev1=stats.anderson(res1)
print(f"statis:{stat1.round(4)}\nCritic point:{np.around(critic1,4)}\nsigni.level:{np.around(plev1, 4)}")
statis:1.5851 Critic point:[0.57 0.649 0.779 0.909 1.081] signi.level:[15. 10. 5. 2.5 1. ]
위 결과에 의하면 이상치 제거에 의해 자료의 정규성 가정에 대한 개선 효과는 나타나지 않았습니다. 그러나 새로 생성된 모델에 대한 분산분석의 귀무가설에 대한 P-value가 작아짐을 보입니다. 이는 데이터 제거가 전체적인 성격의 변화를 초래할 수 있음을 나타내는 것으로 주의해야 할 부분입니다. 데이터 수의 증가 또는 표분분포의 적용으로 정규성 가정의 개선을 모색할 수 있습니다.
plt.figure(figsize=(4, 3)) varQQplot1=stats.probplot(res1, plot=plt) plt.show()
등분산 검정
scipy.stats.levene() 함수를 적용합니다.
equalVal=stats.levene(da1['kos'], da1['kq'], da1['dj'], da1['WonDol'])
print("statistic: %.3f, pValue: %.3f" %(equalVal[0], equalVal[1]))
statistic: 16.539, pValue: 0.000
위 결과에 의하면 이 모델의 등분산성 가정은 만족되지 않습니다.
다중비교(Multi comparison)
분산분석은 여러 그룹간의 차이를 분석하는 것으로 두 그룹간의 차이를 계산하지 않습니다. 위에서 소개한 것과 같이 여러개 그룹들에 대한 개별 비교는 α 오류를 증가시킵니다. 그러므로 이러한 점을 고려한 Tukey HSD 방법을 적용합니다. 이 방법은 그룹내에 모든 데이터에 대한 비교 대신 각 그룹의 평균들을 사용하여 식 1과 같은 검정통계량을 계산합니다.
| $$q=\frac{\vert {\mu_1 - \mu_2}\vert}{\sqrt{\frac{\text{MSE}}{n}}}$$ | (식 1) |
두 그룹의 평균 차에 대한 절대값을 모델의 MSE를 데이터 수(n)으로 나눈수의 제곱으로 계산한 q를 studentized range distribution(SRD) 분포를 기준으로 검정합니다. 이 과정은 statsmodels.stats.multicomp.pairwise_tukeyhsd(x, group, alpha=0.05) 또는 scipy.stats.tukey_hsd(x, y, ...)로 실시할 수 있습니다.
import statsmodels as stm posthoc=stm.stats.multicomp.pairwise_tukeyhsd(da2.iloc[:,1], da2.iloc[:,0]) posthoc.summary()
| group1 | group2 | meandiff | p-adj | lower | upper | reject |
|---|---|---|---|---|---|---|
| WonDol | dj | -0.0004 | 0.9863 | -0.0036 | 0.0028 | False |
| WonDol | kos | -0.0007 | 0.9342 | -0.0039 | 0.0025 | False |
| WonDol | kq | -0.0009 | 0.8688 | -0.0041 | 0.0022 | False |
| dj | kos | -0.0003 | 0.9944 | -0.0035 | 0.0029 | False |
| dj | kq | -0.0005 | 0.9739 | -0.0037 | 0.0027 | False |
| kos | kq | -0.0002 | 0.9981 | -0.0034 | 0.003 | False |
위 tukeyHSD 검정 결과는 각각 두 그룹의 비교에 대한 결과는 모두 유의한 차이를 보이지 않습니다. 그러므로 비교되는 두 그룹의 평균이 같다는 귀무가설을 기각할 수 없습니다. 이 결과는 위 분산분석의 결과와 같습니다.
위 결과를 시각화하여 나타낼 수 있습니다(그림 3).
posthoc.plot_simultaneous(figsize=(5,4))
위 분석을 stats.tukey_hsd()로 분석한 결과는 아래와 같습니다.
hsd=scipy.stats.tukey_hsd(da1['kos'], da1['kq'], da1['dj'], da1['WonDol']) print(hsd)
Tukey's HSD Pairwise Group Comparisons (95.0% Confidence Interval) Comparison Statistic p-value Lower CI Upper CI (0 - 1) 0.000 0.998 -0.003 0.003 (0 - 2) -0.000 0.994 -0.004 0.003 (0 - 3) -0.001 0.934 -0.004 0.002 (1 - 0) -0.000 0.998 -0.003 0.003 (1 - 2) -0.001 0.974 -0.004 0.003 (1 - 3) -0.001 0.869 -0.004 0.002 (2 - 0) 0.000 0.994 -0.003 0.004 (2 - 1) 0.001 0.974 -0.003 0.004 (2 - 3) -0.000 0.986 -0.004 0.003 (3 - 0) 0.001 0.934 -0.002 0.004 (3 - 1) 0.001 0.869 -0.002 0.004 (3 - 2) 0.000 0.986 -0.003 0.004
댓글
댓글 쓰기