이원분산분석(two-way ANOVA)
관련된 내용
- 분산분석 (Analysis of variance)의 개요
- 일원분산분석(One-way ANOVA)
- 사후분석(Post-hoc test)
- 이원분산분석(two-way ANOVA)
일원분산분석(one-way anova)의 경우 명목형인 설명변수(요인)의 각 효과와 연속형 반응변수의 관계를 추정하는 방법입니다. 요인이 2개 이상일 경우 고려해야 할 변동이 증가합니다. 이러한 다요인 분산분석에서 2개의 요인을 포함하는 자료의 분석을 이원분산분석(two-way anova)라고 하며 데이터 구조는 표 1과 같습니다.
| 요인1(α) | 요인2(β) | 합 | 평균 | |||
|---|---|---|---|---|---|---|
| 처리1 | 처리2 | … | 처리k | |||
| 1 | y11 | y12 | … | y1k | T1. | y1. |
| 2 | y21 | y22 | … | y2k | T2. | y2. |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| n | yn1 | yn2 | … | ynk | Tn. | yn. |
| 합 | T.1 | T.2 | … | T.k | T | |
| 평균 | y.1 | y.2 | … | y.k | y | |
표 1은 다음과 같이 구성되어있습니다.
- 요인(α, β, 설명변수): 2
- 반응변수(y)
- 수준(처리): 인자 1의 n(α), 인자 2의 k(β)
표 1의 각 값은 식 1과 같이 나타낼 수 있습니다.
| yijl = μ + αi + βj + αi βj + eijl | (식 1) |
| i: 요인 1의 인덱스, 1, 2, ..., n | |
| j: 요인 2의 인덱스, 1, 2, ..., k | |
| l: i, j에 대응하는 반응변수의 인덱스 |
위 모형은 n개 수준의 α와 k개 수준의 β에 의한 효과인 주효과(main effect)와 두 요인에 의한 교호 효과(interaction effect)는 αβ, 그리고 각 관찰치와의 오차항으로 구성됩니다. 식 1의 모형은 식 2와 같이 각 값의 변동으로 나타낼 수 있습니다.
| \begin{align}\text{SST} &= \text{SSA}+\text{SSB}+\text{SSAB}+\text{SSE}\\ \sum(\text{y}_{ijl}-\bar{\text{y}}_{...})^2&= \sum(\text{y}_{i..}-\bar{\text{y}}_{...})^2\\&+ \sum(\text{y}_{.j.}-\bar{\text{y}}_{...})^2\\&+ \sum(\bar{\text{y}}_{ij.}-\bar{\text{y}}_{i..}-\bar{\text{y}}_{.j.}+\bar{\text{y}}_{...})^2\\&+ \sum(\text{y}_{ijl}-\bar{\text{y}}_{i..}-\bar{\text{y}}_{.j.}-\bar{\text{y}}_{ij.})^2\end{align} | (식 2) |
| SST: 총제곱합 | |
| SSA: 요인1 제곱합 | |
| SSB: 요인2 제곱합 | |
| SSAB: 교호효과 제곱합 | |
| SSE: 오차 제곱합 |
식 1의 모형을 근거로 다음의 가설을 검정합니다.
| H0: 요인들에 대한 효과가 존재하지 않습니다. |
| H1: 요인들에 대한 효과가 존재합니다. |
분산분석은 각 모집단이 독립이라는 가정을 합니다. 그러므로 요인간에 상호작용이 존재하지 않아야 합니다. 즉, 특정 요인들 사이에 결합된 결과가 각 값들의 합과 같아야 한다는 것입니다. 이 가정이 위배되는 경우 두 인자들 사이의 교호작용을 고려해 주어야 합니다. 일반적으로 데이터 중에서 최대 평균과 최소 평균의 비가 0.5 이상이 아닐 경우 각 변수들이 독립이라고 가정합니다.
이원분산분석의 결과는 표 2로 요약될 수 있습니다.
| Effect | sample 크기 | df | SS | MSS | F |
|---|---|---|---|---|---|
| 요인A 주효과 | a | a-1 | SSA | MSA | MSA/MSE ~ F |
| 요인B 주효과 | b | b-1 | SSB | MSB | MSB/MSE ~ F |
| 교호효과(AB) | ab | (a-1)(b-1) | SSAB | MSAB | MSAB/MSE ~ F |
| 오차변동 | n | 각 자유도의 차이 | SSE | MSE | |
| 총변동 | n | n-1 | SST |
이원분산분석은 일원분산분석과 같이 anova_lm() 함수를 사용합니다. 단지 식 1과 같이 모형에 각 요인효과와 두 요인간의 교호효과를 나타내는 항들이 존재해야 하는 복잡성에서 일원분산분석과의 차이가 있을 뿐입니다.
예 1)
다음은 5개의 그룹, 각각에 3개의 교육방법에 따른 결과입니다. 이 자료에 대해 그룹간, 교육방법간에 효과가 존재하는지를 검정합니다.
| group | Type | ||
|---|---|---|---|
| A | B | C | |
| 1 | 7,13,8,9,8,7 | 11,8,8,7,14,10 | 13,8,11,10,8,12 |
| 2 | 11,14,8,8,14,7 | 12,9,9,10,8,13 | 12,141,2,11,13,10 |
| 3 | 8,8,7,9,9,10 | 14,9, 8,13,8,12 | 9,7,11,12,13,8 |
| 4 | 14,11,11,9,10,7 | 14,11,14,11,13,13 | 9,9,8,7,13,9 |
| 5 | 8,12,12,12,11,7 | 10,7,8,13,14,11, | 7,11,7,9,12,9 |
anova_lm() 함수를 적용하기 위해서는 전달하는 인수는 다음과 같이 구성된 모형입니다.
반응변수 ~ 설명변수
위 표에서 요인 group과 type이 설명변수, 그들에 대응하는 값들이 반응변수입니다. 위 모형을 작성하기 위해서는 자료는 group, type, 그리고 반응변수로 그룹화하여 서로 적합하게 대응되도록 조정해야 합니다. 즉, 분산분석을 위해 적합한 데이터 형태로 전환이 필요합니다. 변환된 자료은 객체의 요소를 반복하는 np.repeat() 함수와 객체 자체를 반복하는 np.tile() 함수를 사용하여 생성할 수 있습니다.
data=np.array([ 7, 13, 8, 9, 8, 7, 11, 14, 8, 8, 14, 7, 8, 8, 7, 9, 9,
10, 14, 11, 11, 9, 10, 7, 8, 12, 12, 12, 11, 7, 11, 8, 8, 7,
14, 10, 12, 9, 9, 10, 8, 13, 14, 9, 8, 13, 8, 12, 14, 11, 14,
11, 13, 13, 10, 7, 8, 13, 14, 11, 13, 8, 11, 10, 8, 12, 12, 14,
12, 11, 13, 10, 9, 7, 11, 12, 13, 8, 9, 9, 8, 7, 13, 9, 7,
11, 7, 9, 12, 9])
d={'type':np.repeat(['a','b','c'], 30),'group': np.tile(np.repeat([1,2,3,4, 5], 6),3), 'data':data}
df=pd.DataFrame.from_dict(d)
df.head(3)
| type | group | data | |
|---|---|---|---|
| 0 | a | 1 | 7 |
| 1 | a | 1 | 13 |
| 2 | a | 1 | 8 |
anova_lm() 함수에 전달되는 모형은 ols() 클래스에 의해 생성됩니다. 모형에서 각 요인들 사이의 교호작용은 *으로 나타냅니다. 클래스 ols()에 전달되는 모델이 교호작용 항만으로 구성되는 경우 각 변수의 주효과는 암시적으로 입력된 형태입니다(식 6.2.3).
| value ~ C(group) + C(type) + C(group) * C(type) | (식 6.2.3) |
| ⇔ value ~ C(group) * C(type) |
mod=ols('data~C(group)*C(type)', data=df).fit()
sm.stats.anova_lm(mod)
| df | sum_sq | mean_sq | F | PR(>F) | |
|---|---|---|---|---|---|
| C(group) | 4.0 | 24.333333 | 6.083333 | 1.239249 | 0.301614 |
| C(type) | 2.0 | 18.200000 | 9.100000 | 1.853780 | 0.163750 |
| C(group):C(type) | 8.0 | 61.800000 | 7.725000 | 1.573676 | 0.147128 |
| Residual | 75.0 | 368.166667 | 4.908889 | NaN | NaN |
위 결과는 두 요인 group과 type에 대한 주효과에 대한 유의확률이 유의수준 0.05보다 크며 교호작용 역시 0.05보다 큰 유의확률을 보이므로 다음의 모든 귀무가설을 기각할 수 없습니다.
| H01: group의 효과는 존재하지 않습니다. |
| H02: type의 효과는 존재하지 않습니다. |
| H03: type과 group의 교호작용에 의한 효과는 존재하지 않습니다. |
예 2)
코스피주가지수(kos), 환율(wu), kodex 반도체(ksemi)의 일간 변화량사이에 유의 한 차이가 존재하는지를 파악합니다.
| kos | uk | ksemi | |
| 0 | 2 | 2 | 3.016632 |
| 1 | 1 | 2 | 2.667315 |
| 2 | 2 | 1 | 2.722964 |
| ⋮ | ⋮ | ⋮ | ⋮ |
위 자료를 호출하고 작성하기 위한 코드는 다음과 같습니다. 위 자료에서 kos와 uk는 명목변수입니다. 즉, 호출된 연속변수를 명목변수로 변환한 것으로 함수 pd.cut()을 사용합니다.
st=pd.Timestamp(2020, 5, 1)
et=pd.Timestamp(2024, 5, 30)
code=['KS11', 'USD/KRW', "091160"]
nme=['kos','uk', 'ksemi']
da={}
for i, j in zip(code, nme):
x=fdr.DataReader(i,st, et)['Close']
da[j]=((x.pct_change())*100).dropna()
mn, mx=da['kos'].min(), da['kos'].max()
kos1=pd.cut(da["kos"], bins=np.linspace(mn, mx+0.1, 5), labels=range(4))
mn, mx=da['uk'].min(), da['uk'].max()
uk1=pd.cut(da["uk"], bins=np.linspace(mn, mx+0.1, 5), labels=range(4))
data=pd.concat([kos1, uk1, da['ksemi']], join="inner", axis=1)
data.index=range(len(data))
data.columns=['kos','uk', 'ksemi']
da=data.dropna()
da.head(3)
이원분산분석을 위한 모델은 두 요인인 kos와 uk를 설명변수로하여 반응변수인 ksemi에 미치는 효과를 결정하기 위한 것입니다.
mod=ols('ksemi~C(kos)*C(uk)', data=data).fit()
np.around(sm.stats.anova_lm(mod), 4)
| df | sum_sq | mean_sq | F | PR(>F) | |
|---|---|---|---|---|---|
| C(kos) | 3.0 | 1096.42 | 365.47 | 184.03 | 0.00 |
| C(uk) | 3.0 | 9.75 | 3.25 | 1.64 | 0.18 |
| C(kos):C(uk) | 9.0 | 11.51 | 1.28 | 0.64 | 0.76 |
| Residual | 988.0 | 1962.13 | 1.99 | NaN | NaN |
위 결과에 의하면 반응변수 ksemi에 대한 변수 kos에 의한 주효과만이 존재합니다. 이 결과의 신뢰성은 분석의 전제조건의 충족 여부에 좌우 됩니다. 먼저 정규성 부합 여부를 확인하기 위해 생성된 모델의 잔차를 대상으로 분석을 진행합니다.
res=mod.resid
nor=stats.shapiro(res)
print("통계량: %.3f, p-value: %.3f" %(nor[0], nor[1]))
통계량: 0.991, p-value: 0.000
위 결과는 정규성에 부합한다는 귀무가설을 기각할 수 있습니다. 이 결과를 시각화하면 그림 1과 같습니다.
plt.figure(figsize=(5,3)) stats.probplot(res, plot=plt) plt.show()
그림 1에 의하면 두 개 이상의 가능한 이상치가 발견되고 양쪽 꼬리부분에서 정규성의 이탈이 발견되었습니다. 정규성 가정에 대한 개선을 위해 각 요인들의 표본분포를 적용합니다. 위 데이터는 시간에 따른 것으로 모든 변수의 데이터들의 표집은 동 시간대에 이루어져야 합니다. 이 점을 고려하여 다음과 같이 표본들의 평균으로 구성된 데이터를 생성하였습니다.
st=pd.Timestamp(2020, 5, 1)
et=pd.Timestamp(2024, 5, 30)
code=['KS11', 'USD/KRW', "091160"]
nme=['kos','uk', 'ksemi']
da={}
for i, j in zip(code, nme):
x=fdr.DataReader(i,st, et)['Close']
da[j]=((x.pct_change())*100).dropna()
st=pd.Timestamp(2020, 5, 1)
et=pd.Timestamp(2024, 5, 30)
code=['KS11', 'USD/KRW', "091160"]
nme=['kos','uk', 'ksemi']
da={}
for i, j in zip(code, nme):
x=fdr.DataReader(i,st, et)['Close']
da[j]=((x.pct_change())*100).dropna()
data=(pd.DataFrame.from_dict(da)).dropna()
data_sam=pd.DataFrame()
for i in range(1000):
tst=pd.DataFrame(data.sample(10, random_state=i).mean(axis=0)).T
data_sam=pd.concat([data_sam, tst])
data_sam.index=range(1,1001)
data_sam.head(3).round(2)
| kos | uk | ksemi | |
| 1 | 0.05 | 0.16 | 0.32 |
| 2 | -0.29 | 0.05 | 0.51 |
| 3 | 0.01 | 0.10 | 0.15 |
변수 kos, uk를 목록화합니다.
da_sam=data_sam.copy() da_sam["kos"]=pd.qcut(data_sam["kos"], q=2, labels=range(1, 3)) da_sam["uk"]=pd.qcut(data_sam["uk"], q=2, labels=range(1, 3)) da_sam.head(3).round(2)
| kos | uk | ksemi | |
| 1 | 1 | 2 | 0.32 |
| 2 | 1 | 2 | 0.51 |
| 3 | 1 | 2 | 0.15 |
분산분석을 시행합니다.
mod2=ols("ksemi~C(kos)*C(uk)", data=data_sam).fit()
sm.stats.anova_lm(mod2, type=2).round(2)
| df | sum_sq | mean_sq | F | PR(>F) | |
| C(kos) | 1.0 | 138.02 | 138.02 | 588.58 | 0.00 |
| C(uk) | 1.0 | 0.62 | 0.62 | 2.65 | 0.10 |
| C(kos):C(uk) | 1.0 | 0.61 | 0.61 | 2.59 | 0.11 |
| Residual | 996.0 | 233.55 | 0.23 | NaN | NaN |
분산분석 결과는 이전 결과와 유사합니다. 즉, 요인 kos에 의한 효과만 인정됩니다. 분석을 위한 자료의 정규성과 등분산성을 검정합니다. 다음 코드는 정규성 검정 결과입니다.
res2=mod2.resid
nor2=stats.shapiro(res2)
print("통계량: %.3f, p-value: %.3f" %(nor2[0], nor2[1]))
통계량: 0.998, p-value: 0.228
위 결과는 정규분포에 부합한다는 귀무가설을 기각 할 수 없습니다. 등분산성을 테스트하기 위해 함수 scipy.stats.levene()을 적용합니다. 이 함수의 인수는 클래스로 구성된 각 그룹이 되므로 먼저 각 그룹의 자료수를 파악하기 위해 pandas.pivot_table() 함수를 사용하여 피봇테이블을 작성합니다.
pivot=pd.pivot_table(da_sam, values="ksemi", index="kos", columns="uk", aggfunc=lambda x: len(x)) pivot
| uk | 1 | 2 |
|---|---|---|
| kos | ||
| 1 | 222 | 278 |
| 2 | 278 | 222 |
위 결과는 각 그룹에 포함된 샘플의 규모를 나타내지만 실질적으로 필요한 부분은 각 포함된 샘플의 값입니다. 이것은 함수 pandas.groupby()를 사용하여 확인할 수 있습니다.
gr=da_sam.groupby(["kos","uk"]) gr.groups
{(1, 1): [5, 7, 9, 10, 11, 13, 20, 21, 22, 24, 27, 37, 39, 42, 43, 44, …
위 결과는 변수 "kos"와 "uk"의 각 클래스로 분류된 "ksemi"의 인덱스를 사전형식(dictionary)으로 나타냅니다. 즉, 위 결과의 key는 두 설명변수의 클래스입니다.
gr.groups.keys()
dict_keys([(1, 1), (1, 2), (2, 1), (2, 2)])
예를 들어 위 결과를 사용하여 "kos"와 "uk"의 클래스 (1,1)에 해당하는 ksemi 자료는 다음과 같습니다.
da_sam.ksemi[gr.groups[(1,1)]].head(3)
5 -0.824925 7 -0.490587 9 -0.439029 Name: ksemi, dtype: float64
위와 같은 형식으로 각 그룹에 포함된 ksemi 자료를 호출하여 levene 검정을 실시하면 다음과 같습니다.
leveneRe=stats.levene(da_sam.ksemi[gr.groups[(1,1)]],da_sam.ksemi[gr.groups[(1,2)]],da_sam.ksemi[gr.groups[(2,1)]], da_sam.ksemi[gr.groups[(2,2)]])
print("통계량: %.3f, p-value: %.3f" %(leveneRe[0], leveneRe[1]))
통계량: 1.674, p-value: 0.171
그룹별로 호출된 자료에 대해 lenvene 검정을 실시합니다.
leveneRe=stats.levene(new[0], new[1],new[2],new[3])
print("통계량: %.3f, p-value: %.3f" %(leveneRe[0], leveneRe[1]))
통계량: 2.190, p-value: 0.088
위 결과는 자료는 등분산성을 만족한다는 귀무가설을 기각 할 수 없습니다.
분산분석에 대한 사후분석 tukeyHSD를 실시합니다. 요인변수(factor)가 두개 이상일 경우 statsmodels.stats.multicomp.pairwise_tukeyhsd(x, group, alpha=0.05) 함수를 적용할 수 없습니다. 대신에 쌍(pair)분석을 위한 대상에 대해 t-검정을 실시합니다. 두 변수의 클래스로 분류된 총 4 그룹들에서 분석을 위한 두 개의 그룹의 추출은 조합으로 계산할 수 있으며 itertools.combinations() 함수를 사용합니다.
import itertools idx=list(itertools.combinations(gr.groups.keys(), 2)) idx
[((1, 1), (1, 2)), ((1, 1), (2, 1)), ((1, 1), (2, 2)), ((1, 2), (2, 1)), ((1, 2), (2, 2)), ((2, 1), (2, 2))]
t-test는 stats.ttest_ind() 함수를 적용합니다. 또한 t-test는 등분산의 성립 여부에 따라 합동 표준편차의 계산방식이 다르므로 비교 대상의 등분산 여부를 먼저 결정합니다. stats.bartlett() 함수를 적용합니다.
["pValue: %.3f" %stats.bartlett(da_sam.ksemi[gr.groups[i]], da_sam.ksemi[gr.groups[j]])[1] for i, j in idx]
['pValue: 0.380', 'pValue: 0.048', 'pValue: 0.728', 'pValue: 0.241', 'pValue: 0.214', 'pValue: 0.020']
["%s, %s, pValue: %.3f" %(i, j, stats.ttest_ind(da_sam.ksemi[gr.groups[i]], da_sam.ksemi[gr.groups[j]])[1]) for i, j in idx]
['(1, 1), (1, 2), pValue: 0.069', '(1, 1), (2, 1), pValue: 0.000', '(1, 1), (2, 2), pValue: 0.000', '(1, 2), (2, 1), pValue: 0.000', '(1, 2), (2, 2), pValue: 0.000', '(2, 1), (2, 2), pValue: 0.036']
위 결과는 (1, 1):(1, 2)와 (2, 1):(2, 2)의 분석에서만 유의한 차이는 보이지 않습니다. 이 두 룹의 경우 uk 클래스만 다른 경우입니다. 그러므로 ksemi에 효과를 나타내는 변수는 kos라고 할 수 있으며 이 결과는 위의 모델 mod2에 대한 분산분석 결과와 같습니다.
댓글
댓글 쓰기