분산분석(ANOVA) 2

집단간의 차이가 존재하는지 통계적으로 검정하기 위해서는, 실험 모집단에 대한 모형이 있어야 합니다. 이를 위하여, $i$번째 집단에서의 반응치가 평균${\mu }_i$, 분산이 ${\sigma }^2$인 정규모집단을 따른다고 가정하겠습니다.

그렇다면, $i$번째 집단에서 $j$번째 반응값을 $Y_{ij}$라고 하면 $Y_{ij}$는 아래와 같이 표현할 수 있습니다.

이 때, 오차항 ${\epsilon }_{ij}$은 모두 독립이고 평균이 0, 분산이 ${\sigma }^2$인 정규분포를 따릅니다.

$k$개의 집단에서 모평균 차이가 없다는 귀무가설은 아래와 같습니다.

만약, 귀무가설이 맞다면 처리제곱합(SStr)의 값이 작아질 것이고 평균처리제곱(MStr)도 작아질 것입니다. 반면에 귀무가설이 틀리다면, 평균처리제곱(MStr)은 커질 것입니다. 따라서 평균처리제곱의 값으로 귀무가설 기각여부를 판단해야 하는데, 그 기준으로 평균오차제곱(MSE)가 활용됩니다.(처리제곱합 등에 대한 개념: direction-f.tistory.com/41)

최종적으로, 평균처리제곱과 평균오차제곱의 비율 $\frac{MStr}{MSE}$이 검정통계량이 됩니다. 이 검정통계량의 분포는 자유도가($k-1$, $n-k$)인 F분포를 따르게 됩니다. 이 때 $n$은 전체 반응치 수(=$\sum n_i$)입니다.

F분포는 양수의 구간에서만 확률 값을 가지는 분포이고 모양도 대칭이 아닙니다. F분포에서 분자의 자유도는 $v_1$, 분모의 자유도는 $v_2$로 표기합니다.  

모양은 아래와 같습니다.

https://ko.wikipedia.org/wiki/F_%EB%B6%84%ED%8F%AC

F분포를 활용한 모평균의 동일성 검정은 아래와 같이 정의 될 수 있습니다.

 

k개의 집단의 모평균에 대한 가설 $H_0$: ${\mu }_1{\mu }_2=\cdots ={\mu }_k$를 검정하기 위한 검정통계량은 다음과 같습니다.

이는 자유도가 $(k-1,n-k)$인 F분포를 따르고 유의수준 $\alpha$의 기각역은 다음과 같습니다.

$F_{\alpha }(k-1,n-k)$는 F분포에서 상위 $\alpha$의 확률을 주는 경계값입니다.

 

Python을 활용해서 실습해보겠습니다.

from scipy.stats import f
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import researchpy as rp

df = pd.read_csv("https://raw.githubusercontent.com/researchpy/Data-sets/master/difficile.csv")
## dose : 약 복욕량, libido: 약물에 대한 효과 -> 약 복용량에 따른 약물에 대한 효과를 검증함

df = df.drop(["person"], axis =1)
df['dose'].replace({1: 'placebo', 2: 'low', 3: 'high'}, inplace= True)

results =rp.summary_cont(df['libido'].groupby(df['dose']))
print(results)

## 결과
#
#          N  Mean      SD      SE  95% Conf.  Interval
# dose                                                 
# high     5   5.0  1.5811  0.7071     3.0368    6.9632
# low      5   3.2  1.3038  0.5831     1.5811    4.8189
# placebo  5   2.2  1.3038  0.5831     0.5811    3.8189

##

## F-Statstic 도출
## MStr 도출
SStr = (results.loc["high"]["N"]* (results.loc["high"]["Mean"]-df["libido"].mean())**2)\
+(results.loc["low"]["N"]* (results.loc["low"]["Mean"]-df["libido"].mean())**2)\
+(results.loc["placebo"]["N"]* (results.loc["placebo"]["Mean"]-df["libido"].mean())**2)

print("SStr:{}".format(SStr))

MStr = SStr/(len(results.index)-1)

print("MStr:{}".format(MStr))

## MSE 도출
SSE = (df.loc[df["dose"]=="high"]["libido"]-results.loc["high"]["Mean"]).pow(2).sum()\
+ (df.loc[df["dose"]=="low"]["libido"]-results.loc["low"]["Mean"]).pow(2).sum()\
+ (df.loc[df["dose"]=="placebo"]["libido"]-results.loc["placebo"]["Mean"]).pow(2).sum()

print("SSE:{}".format(round(SSE,2)))

MSE = SSE/(len(df)-len(results.index))

print("MStr:{}".format(round(MSE,2)))


#F-statistic
F_ = MStr/MSE

print("F-Statistc:{}".format(F_))

## 결과
# SStr:20.133333333333333
# MStr:10.066666666666666
# SSE:23.6
# MStr:1.97
# F-Statistc:5.11864406779661
##

## F-distribution
df1 = (len(results.index)-1) ## 3
df2 = (len(df)-len(results.index)) ## 12

f_dist = f(df1, df2)
f_05 = f_dist.ppf(0.95)

print(f_05) # 3.8852938346523933 

if F_ >= f_05:
    print("H0 기각") 
else :
    print("H0 채택")

## H0 기각 -> 집단간 평균이 동일하지 않음

p_val  = 1-f_dist.cdf(F_)
print(p_val) # 0.024694289538222614


##시각화
x = np.linspace(f_dist.ppf(0.01), f_dist.ppf(0.99),1000)
plt.plot(x, f_dist.pdf(x),'g-', lw=2, alpha=0.6, label='F-pdf')
plt.title("F-distribution")
plt.legend()
plt.show()

## package를 활용한 ANOVA 분석

import scipy.stats as stats

stats.f_oneway(df['libido'][df['dose'] == 'high'],
               df['libido'][df['dose'] == 'low'],
               df['libido'][df['dose'] == 'placebo'])
               
## F_onewayResult(statistic=5.11864406779661, pvalue=0.024694289538222603)
## 직접계산한 값과 statstic,P-value 동일