파이썬으로하는 일원배치 분산분석 (one-way ANOVA)

1. 일원배치 분산분석 (one-way ANOVA)¶

• 분산분석은 두 개이상의 집단에서 그룹 평균 간 차이를 그룹 내 변동에 비교하여 살펴보는 데이터 분석방법 : 여러 그룹간의 평균의 차이가 통계적으로 유의미 한지를 판단하기 위한 시험법
• 일원배치 분산분석은 반응값에 대해 하나의 범주형 변수의 영향을 알아보기 위해 사용되는 검증방법 : 한가지 변수의 변화가 결과 변수에 미치는 영향을 보기 위해 사용
• F 검정 통계량을 이용
• 각 집단의 측정치는 서로 독립적이며 정규분포를 따른다.(정규성 가정) / 각 집단 측정치의 분산은 같다. (등분산 가정)

 

 

2. Example 1¶

• 22명의 심장 수술을 받은 환자를 3가지 그룹으로 나누고 적혈구의 엽산 수치를 24시간 이후에 측정하였음.

In [1]:

import numpy as np
import urllib

url = 'https://raw.githubusercontent.com/thomas-haslwanter/statsintro_python/master/ipynb/Data/data_altman/altman_910.txt'
data = np.genfromtxt(urllib.request.urlopen(url), delimiter=',')
data

Out[1]:

array([[243.,   1.],
       [251.,   1.],
       [275.,   1.],
       [291.,   1.],
       [347.,   1.],
       [354.,   1.],
       [380.,   1.],
       [392.,   1.],
       [206.,   2.],
       [210.,   2.],
       [226.,   2.],
       [249.,   2.],
       [255.,   2.],
       [273.,   2.],
       [285.,   2.],
       [295.,   2.],
       [309.,   2.],
       [241.,   3.],
       [258.,   3.],
       [270.,   3.],
       [293.,   3.],
       [328.,   3.]])

 

In [2]:

import seaborn as sns
import matplotlib.pyplot as plt
sns.set_style('whitegrid')

group1 = data[data[:,1]==1,0]  
group2 = data[data[:,1]==2,0]
group3 = data[data[:,1]==3,0]

plot_data = [group1, group2, group3]
plt.boxplot(plot_data)
plt.show()

 

 

boxplot에서 보듯이 평균값의 차이가 실제로 의미가 있는 차이인지, 분산이 커서 그런것인지 애매한 상황이다. 이런 상황에서 분산분석을 통해 통계적 유의성을 알아볼 수 있음.

 

 

Statsmodel을 사용한 일원분산분석¶

In [3]:

import pandas as pd
import warnings
warnings.filterwarnings('ignore')
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

df = pd.DataFrame(data, columns=['value','treatment'])
print(df.head(3))

model = ols('value ~ C(treatment)', df).fit()
print(anova_lm(model))

 

   value  treatment
0  243.0        1.0
1  251.0        1.0
2  275.0        1.0
                df        sum_sq      mean_sq         F    PR(>F)
C(treatment)   2.0  15515.766414  7757.883207  3.711336  0.043589
Residual      19.0  39716.097222  2090.320906       NaN       NaN

 

 

3. Example 2¶

• iris데이터를 이용해 종별로 꽃받침의 폭(sepal.width)의 평균이 같은지 혹은 차이가 있는지를 확인
• 귀무가설 : 세가지 종에 대해 sepal.width평균은 모두 같다.
• 대립가설 : 적어도 하나의 종에 대한 sepal.width의 평균값에는 차이가 있다.

In [4]:

from sklearn.datasets import load_iris

iris = load_iris()
data = pd.DataFrame(data=iris.data, columns=iris.feature_names)
data['target'] = iris.target
data.columns = ['sepal_length','sepal_width','petal_length','petal_width','target']
data.head(3)

Out[4]:

	sepal_length	sepal_width	petal_length	petal_width	target
0	5.1	3.5	1.4	0.2	0
1	4.9	3.0	1.4	0.2	0
2	4.7	3.2	1.3	0.2	0

 

 

3.1 step1 일원분산분석 가정 확인¶

 

자료 수집이 무선표집(random sampling)되었다면 독립성은 만족한다고 봄

In [5]:

#정규성 확인 - 집단(수준)별로 실시

from  scipy.stats import shapiro
print(shapiro(data.sepal_width[data.target==0]))
print(shapiro(data.sepal_width[data.target==1]))
print(shapiro(data.sepal_width[data.target==2]))

 

ShapiroResult(statistic=0.97171950340271, pvalue=0.2715264856815338)
ShapiroResult(statistic=0.9741330742835999, pvalue=0.33798879384994507)
ShapiroResult(statistic=0.9673910140991211, pvalue=0.1809043288230896)

 

모두 p-value가 0.05 보다 크므로 각 집단의 자료가 정규성을 만족한다고 볼 수 있음.

(정규성 검정 귀무가설-> 정규성을 만족한다 이므로 위의 수치는 유의수준 0.05보다 크기 때문에 귀무가설을 기각하지 않는다. )

 

In [6]:

#등분산성 확인 - 레빈 검증
from scipy.stats import levene
print(levene(data.sepal_width[data.target==0],
      data.sepal_width[data.target==1],
      data.sepal_width[data.target==2]))

#등분산성 확인 - 바틀렛 검증
from scipy.stats import bartlett
print(bartlett(data.sepal_width[data.target==0],
      data.sepal_width[data.target==1],
      data.sepal_width[data.target==2]))

 

LeveneResult(statistic=0.5902115655853319, pvalue=0.5555178984739075)
BartlettResult(statistic=2.0910752014391774, pvalue=0.35150280041581317)

 

두 테스트 모두, 세 집단의 모분산에 유의미한 차이를 발견하지 못함. 등분산성 가정이 유지됨

 

 

3.2 step2 일원분산분석 수행¶

In [7]:

model = ols('sepal_width ~ C(target)', data).fit()
anova_lm(model)

Out[7]:

	df	sum_sq	mean_sq	F	PR(>F)
C(target)	2.0	11.344933	5.672467	49.16004	4.492017e-17
Residual	147.0	16.962000	0.115388	NaN	NaN

 

• Pr(>F)가 p-value. 이 값이 유의수준 0.05하에서 귀무가설을 기각함. 따라서 세가지 종에따른 꽃받침 폭이 모두 동일하지는 않다고 결론내릴 수 있다. 즉, 종별 꽃받침 폭의 평균값들 중에서 적어도 어느 하나의 종은 통계적으로 유의한 차이가 있다.
• SSA의 자유도는 2(집단의 수-1=3-1), SST의 자유도는 147(관측값의 수-집단의 수=150-3)

 

 

3.3 step3 사후분석¶

 

• 분산분석의 결과 귀무가설이 기각되어 적어도 한 집단에서 평균의 차이가 있음이 통계적으로 증명되었을 경우, 어떤 집단들에 대해서 평균의 차이가 존재하는지를 알아보기 위해 실시하는 분석
• 조합 가능한 모든 쌍에 대해 비교를 하므로 과잉검증으로 인한 FWER 증가
• 널리 쓰이는 봉페로니 교정과 투키의 HSD를 소개

 

In [8]:

#사후분석을 위한 준비
from statsmodels.sandbox.stats.multicomp import MultiComparison
import scipy.stats

comp = MultiComparison(data.sepal_width, data.target)

 

In [9]:

#봉페로니 교정

result = comp.allpairtest(scipy.stats.ttest_ind, method='bonf')
print(result[0])

#투키의 HSD - Tuckey's Honestly Significant Difference = "진정으로 유의미한 차이"
from statsmodels.stats.multicomp import pairwise_tukeyhsd
hsd = pairwise_tukeyhsd(data['sepal_width'], data['target'], alpha=0.05)
hsd.summary()

 

Test Multiple Comparison ttest_ind 
FWER=0.05 method=bonf
alphacSidak=0.02, alphacBonf=0.017
=============================================
group1 group2   stat   pval  pval_corr reject
---------------------------------------------
     0      1   9.455    0.0       0.0   True
     0      2  6.4503    0.0       0.0   True
     1      2 -3.2058 0.0018    0.0055   True
---------------------------------------------

 

Out[9]:

Multiple Comparison of Means - Tukey HSD, FWER=0.05

group1	group2	meandiff	p-adj	lower	upper	reject
0	1	-0.658	0.001	-0.8189	-0.4971	True
0	2	-0.454	0.001	-0.6149	-0.2931	True
1	2	0.204	0.0088	0.0431	0.3649	True

 

pval=p-value 모든 종들에 대해서 꽃받침 폭의 평균값은 각각 통계적으로 유의한 차이가 있다는 것을 알 수 있음. 종 0과 1의 meandiff가 음수이기 때문에 꽃받침의 폭은 종이 0일때보다 1일때가 통계적으로 유의하게 큰 값을 가진다고 해석할 수 있음.