[data analysis] 일원분산분석(One-way ANOVA)

일원분산분석(One-way ANOVA)

관련된 내용

• 분산분석 (Analysis of variance)의 개요
• 일원분산분석(one-way ANOVA)
• 사후분석(Post-hoc test)
• 이원분산분석(two-way ANOVA)

분산분석의 귀무가설은 다음과 같습니다.

H0 : µ₁ = µ₂ = · · · = µ_n

분산분석을 위해 다음을 가정합니다.

• 각 모집단은 정규분포를 따릅니다.
• 모든 모집단의 분산은 동일합니다.
• 관측치들은 독립적이어야 합니다.

위의 정규성 가정은 각 그룹에 대응하는 모집단을 검정하는 것은 어렵기 때문에 모델의 잔차에 대한 검정으로 대신합니다. 또한 독립성은 자료의 수집단계의 정보에 의해 판단되는 것으로 분석 중에 그 검정은 쉽지 않습니다.

표 1에서 나타낸 것과 같이 one-way ANOVA는 각 factor에 포함되는 수준 즉 처리(treatment, 요인수준)가 없습니다. 그러므로 일원분산분석에서는 요인과 treatment가 같으며 각 요인에 포함된 값들(반응변수)을 그룹화합니다. 이 구조에서 각 그룹내의 변동과 각 요인들 사이의 변동을 비교합니다.

표 1 일원분산분석을 위한 자료구조

요인(처리)	1	2	…	t
반응(값)	x11	x12	…	x1t
	x21	x22	…	x2t
	⋮	⋮	⋮	⋮
	xn1	xn2	…	xnt
평균	X1	X2	…	X.t
총평균	X..

표 1로부터 각 값들은 식 1과 같이 일반화한 모형으로 나타낼 수 있습니다.

xij = μj + eij	(식 1)
xij: 각 값
μj: 그룹 j의 평균
eij: xij에 대응하는 오차
i:1, 2,…, n(그룹내 값의 수)
j:1, 2, …, t(그룹의 수)

이 모형에서 각 변수는 독립적이고 정규분포에 부합한다고 가정했으므로 오차항(e) 역시 평균이 0이고 일정한 분산을 가진 정규분포를 따른다고 가정할 수 있습니다. 즉, 식 2와 같이 나타낼 수 있습니다.

e ∼ N(0, σ2)

(식 2)

각 값의 오차는 식 3과 같이 그룹간 오차와 그룹내 오차의 합으로 계산됩니다.

$$x_{ij}-\bar{x}=(\bar{x_j}-\bar{x})+(x_{ij}-\bar{x_j})$$

(식 3)

식 3에서 $\bar{x}$는 전체 평균, $\bar{x_j}$는 그룹 j의 평균을 나타냅니다.

평균을 기준으로 각 값의 편차의 총합은 0 또는 0에 근접하므로 그것은 자료의 특성을 나타낼 수 없습니다. 그러므로 편차 대신에 편차의 제곱합(sum of squares, SS)이 분포의 특성을 위한 지표로 사용됩니다. 식 3의 각 오차항의 제곱합은 식 4와 같이 양변을 제곱하고 합하는 것으로 생성할 수 있습니다. 그 식의 좌항은 그룹(요인)의 구분없이 모든 데이터에 대한 편차의 제곱이 되며 우항은 그룹간과 그룹내의 편차의 제곱합을 나타냅니다. 여러 그룹이 관계된 모델에서는 공분산(covariance) 고려하기 위해 우항에 $\sum^t_{j=1}\text{Cov}_j$ 즉, 그룹간 공분산의 총합을 첨가 합니다.

\begin{align}\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x})^2&=\sum^t_{j=1}(\bar{x_j}-\bar{x})^2+\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x_j})^2+\sum^t_{j=1}\text{Cov}_j\\&\text{t: 그룹의 수}\\&\text{n: 각 그룹내의 값의 수} \end{align}

(식 4)

그러나 분산분석은 각 그룹이 독립임을 가정하므로 식 4의 공분산은 0이 되므로 식 5과 같이 정리됩니다.

\begin{align}&\begin{aligned}\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x})^2&=\sum^t_{j=1}(\bar{x_j}-\bar{x})^2+\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x_j})^2 \\ \text{SST}&=\quad\text{SSTr} \qquad + \;\text{SSE} \end{aligned}\\ &\text{SST}:\text{ sum of squares total, 전체 제곱합} \\&\text{SSTr}:\text{ sum of squares treatment, 처리제곱합(between group)} \\&\text{SSE}:\text{ sum of squares error, 오차제곱합(within group)}\end{align}

(식 5)

식 5의 각 제곱합을 대응되는 자유도로 나눈 값을 제곱평균(mean of squares)라고 하며 식 6과 같이 정의합니다.

\begin{align}&\begin{aligned}\text{MST}&=\frac{\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x})^2}{n\cdot t - 1}\\&=\frac{\text{SST}}{n\cdot t - 1}\\ \text{MSTr}&=\frac{\sum^t_{j=1}(\bar{x_j}-\bar{x})^2}{t-1}\\&=\frac{\text{SSTr}}{t-1}\\ \text{MSE}&=\frac{\sum^t_{j=1}\sum^n_{i=1}(x_{ij}-\bar{x_j})^2 }{n\cdot t-t}\\&=\frac{\text{SSE}}{n\cdot t-t}\end{aligned}\\ &\text{MST: mean of squares total, 전체 제곱평균}\\ &\text{MSTr: mean of squares treatment, 처리제곱평균}\\ &\text{MSE: mean of squares error, 오차제곱평균}\end{align}

(식 6)

예 1)

다음 코드로 생성한 자료(da)에 대해 sstr, sse, mstr, mse를 계산해 봅니다.

np.random.seed(1)
da=np.random.randint(1, 11, size=(10, 4))

print(da[:3,:])

[[ 6  9 10  6]
 [ 1  1  2  8]
 [ 7 10  3  5]]

자료 da 에 대한 각 통계량을 계산하면 다음과 같습니다.

다음 코드의 .shape 속성은 객체의 모양 즉, 행과 열의 수를 나타냅니다.

n, t=shape(da)
N=n*t
print(f'Total #:{N}, factor #:{t}, # per factor:{n}')

Total #:40, factor #:4, # per factor:10

da.shape

(10, 4)

Mu=np.mean(da)
mu=np.mean(da, axis=0)
print(f'전체평균: {Mu}, 각 그룹 평균: {mu}')

전체평균: 6.1, 각 그룹 평균: [6.5 6.8 4.6 6.5]

sst=np.sum((np.ravel(da)-Mu)**2)
mst=sst/(N-1)
print(f'SST: {round(sst,3)}, MST: {round(mst, 3)}')

SST: 377.6, MST: 9.682

sstr=np.sum((Mu-mu)**2*np.repeat(10, 4))
mstr=sstr/(t-1)
print(f'SSTr: {round(sstr,3)}, MSTr: {round(mstr, 3)}')

SSTr: 30.6, MSTr: 10.2

sse=np.sum((da-mu)**2)
mse=sse/(N-t)
print(f'SSE: {round(sse,3)}, MSE: {round(mse, 3)}')

SSE: 347.0, MSE: 9.639

위 통계량은 파이썬의 패키지인 statsmodels.api의 anova_lm(model) 함수를 사용하여 계산할 수 있습니다. 이 함수에 전달하는 인수는 statsmodels.api.ols()에 의해 생성되는 회귀모델입니다. 회귀모델은 설명변수와 반응변수 사이에 관계를 수식화하는 모델을 생성합니다. 이 함수를 사용하기 위해서는 다음 코드와 같이 패키지를 임포트 해야 합니다.

import statsmodels.api as sm
from statsmodels.formula.api import ols

ols() 함수에 전달하는 인수는 설명변수와 반응변수로 구분하여야 합니다. 그러므로 자료 da의 구조를 각 그룹과 값으로 표 2와 같이 변형합니다.

표 2 자료 구조의 변화

		⇒	설명변수	반응변수
group 1	group 2		group 1	x11
x11	x12		group 1	x21
x21	x21		group 2	x12
			group 2	x22

위 변환을 위해 다음 코드에서는 np.repeat()와 np.ravel()함수를 적용하였습니다.

da1=np.c_[np.repeat(range(1, 5), 10), np.ravel(da, order='F')]
da1=pd.DataFrame(da1, columns=["ind",'de'])
da1.head()

	ind	de
0	1	6
1	1	1
2	1	7
3	1	6
4	1	5

model=ols('de~C(ind)', data=da1).fit()
sm.stats.anova_lm(model)

	df	sum_sq	mean_sq	F	PR(>F)
C(ind)	3.0	30.6	10.200000	1.058213	0.378876
Residual	36.0	347.0	9.638889	NaN	NaN

위 결과에서 설명변수들 사이의 오차는 1행, 그룹내 오차는 2행에 대응하는 값입니다. 즉, 1행과 2행은 각각 그룹간, 그룹내의 통계량을 나타낸 것입니다.

정규성과 독립성에 부합하는 자료의 선형변환 역시 정규분포를 따릅니다. 각 자료들의 편차는 원시자료의 선형변환 결과이므로 변환된 자료와 그들의 제곱은 각각 정규분포와 χ²분포에 부합할 것입니다. 또한 χ²에 부합하는 자료들의 분산비는 F-분포를 형성합니다. 분산분석은 이 가정들을 전제조건으로 시행하는 것으로 그 분석의 통계량인 MSTr과 MSE의 비(식 7)에 대한 검정은 F분포를 기반으로 시행합니다.

\begin{align}F&=\frac{\text{MSTr}}{\text{MSE}}\\&=\frac{\frac{SSTr}{t-1}}{\frac{SSE}{n-t}}\end{align}

(식 7)

다음은 자료에 대한 정규성 검정 결과입니다.

pval={}
for i in range(0, da.shape[1]):
    pval[i]=round(stats.shapiro(da[:,i])[1],3)
pval

{0: 0.262, 1: 0.015, 2: 0.324, 3: 0.535}

자료의 두 번째 그룹을 제외하고 모두 정규분포에 부합한다고 할 수 있습니다. 두번쨰 그룹은 정규성 부합에 위배됩니다. 이러한 경우 자료의 조정이나 분산분석 대상에서 제외를 고려해야 하지만 이 예제에서 그대로 분석을 진행합니다.

독립성 검정을 위해 scipy.stats.chi2_contingency() 함수를 적용합니다. 이 함수의 경우 필요한 데이터 형태는 교차표이므로 위의 da1을 사용하여 교차표를 생성합니다.

ta=pd.crosstab(da1.iloc[:,0], da1.iloc[:,1])
indeRe=stats.chi2_contingency(ta)
p_val=indeRe[1]; p_val.round(3)

0.441

위 두 검정의 결과는 분산분석의 전제조건에 부합함을 나타냅니다. 그러므로 anova_lm()의 결과를 적용할 수 있습니다. 이 결과의 F값은 검정통계량되며 그에 대응하는 유의확률 (p-value)은 다음의 귀무가설을 기각할 수 없음을 의미합니다. 즉, 4개 그룹의 평균은 같다고 할 수 있습니다.

H0: μ₁ = μ₂ = … =μ_n

다음 예에서 ols() 클래스 적용 없이 분산분석을 실행합니다.

예 2)

위 예제의 자료 da에 대해 분산분석을 실시합니다.

Fratio=mstr/mse
print(f"F 비: {round(Fratio, 3)}")

F 비: 1.058

ci=stats.f.interval(0.95, t-1, N-t)
print(f"하한: %.3f, 상한: %.3f" %(ci[0], ci[1]))

하한: 0.071, 상한: 3.505

pV=stats.f.sf(Fratio, t-1, N-t)
print(f'p-value: {round(pV, 3)}')

p-value: 0.379

위 결과에 의하면 F 비(F ratio)는 신뢰구간 내에 존재하며 이 통계량의 유의확률(p value) 역시 유의수준(0.05)에 비해 크므로 귀무가설을 기각할 수 없습니다. 즉, 각 그룹의 모평균을 같다고 할 수 있습니다. 이 결과는 위 statsmodels.stats.anoval_lm() 함수를 적용한 결과와 같습니다.

예 3)

일정기간의 코스피(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)

	kos	kq	dj	WonDol
1	-0.023	-0.008	-0.008	0.011
2	-0.008	-0.006	0.000	0.001

위 데이터에 대한 분산분석 모형을 생성하기 위해서는 표 2와 같이 데이터를 변수이름과 각 이름에 대응하는 값들의 형태로 재구조화 해야 합니다. 이러한 변형은 pd.melt() 함수를 사용합니다.

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

위 자료에 대한 anova 분석 결과는 다음과 같습니다.

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은 다양한 속성을 가지고 있습니다. 이 속성을 사용하여 위 표의 각 값들을 호출할 수 있습니다.

print("f value: %.4f\np-value:%.4f\nmsetr:%.4f\nmse: %.4f"%(model.fvalue,model.f_pvalue, model.mse_model, model.mse_resid))

f value: 0.2223
p-value:0.8809
msetr:0.0000
mse: 0.0001

scipy.stats.f_oneway(x1, x2,...) 함수를 적용하여 같은 결과를 확인할 수 있습니다. 이 함수는 검정통계량과 p value 만을 반환합니다.

static, pV=stats.f_oneway(da1["kos"],da1["kq"],da1["dj"],da1["WonDol"])
print("static: %.3f, p-value: %.3f" %(static, pV))

static: 0.222, p-value: 0.881