지난번 포스팅에서는 Python의 statsmodels  모듈을 이용하여 여러개의 숫자형 변수에 대해 집단 간 평균의 차이가 있는지를 for loop 순환문을 사용하여 검정하는 방법(rfriend.tistory.com/639)을 소개하였습니다. 

 

Python에서 for loop 문을 사용하면 순차적으로 처리 (sequential processing) 를 하게 되므로, 일원분산분석을 해야 하는 숫자형 변수의 개수가 많아질 수록 선형적으로 처리 시간이 증가하게 됩니다. 

 

Greenplum에서 PL/Python (또는 PL/R) 을 사용하면 일원분산분석의 대상의 되는 숫자형 변수가 매우 많고 데이터 크기가 크더라도 분산병렬처리 (distributed parallel processing) 하여 ANOVA test 를 처리할 수 있으므로 신속하게 분석을 할 수 있는 장점이 있습니다.

 

더불어서, 데이터가 저장되어 있는 DB에서 데이터의 이동 없이(no data I/O, no data movement), In-DB 처리/분석이 되므로 work-flow 가 간소화되고 batch scheduling 하기에도 편리한 장점이 있습니다.  

 

만약 데이터는 DB에 있고, 애플리케이션도 DB를 바라보고 있고, 분석은 Python 서버 (또는 R 서버)에서 하는 경우라면, 분석을 위해 DB에서 데이터를 samfile 로 떨구고, 이를 Python에서 pd.read_csv()로 읽어들여서 분석하고, 다시 결과를 text file로 떨구고, 이 text file을 ftp로 DB 서버로 이동하고, psql로 COPY 문으로 테이블에 insert 하는 workflow ... 관리 포인트가 많아서 정신 사납고 복잡하지요?!  

 

자, 이제 Greenplum database에서 PL/Python으로 일원분산분석을 병렬처리해서 집단 간 여러개의 개별 변수별 평균 차이가 있는지 검정을 해보겠습니다. 

 

 

(1) 여러 개의 변수를 가지는 샘플 데이터 만들기

 

정규분포로 부터 난수를 발생시켜서 3개 그룹별로 각 30개 씩의 샘플 데이터를 생성하였습니다. 숫자형 변수로는 'x1', 'x2', 'x3', 'x4'의 네 개의 변수를 생성하였습니다. 이중에서  'x1', 'x2'는 3개 집단이 모두 동일한 평균과 분산을 가지는 정규분포로 부터 샘플을 추출하였고, 반면에 'x3', 'x4'는 3개 집단 중 2개는 동일한 평균과 분산의 정규분포로 부터 샘플을 추출하고 나머지 1개 집단은 다른 평균을 가지는 정규분포로 부터 샘플을 추출하였습니다. (뒤에  one-way ANOVA 검정을 해보면 'x3', 'x4'에 대한 집단 간 평균 차이가 있는 것으로 결과가 나오겠지요?!)

 

import numpy as np
import pandas as pd

# generate 90 IDs 
id = np.arange(90) + 1 

# Create 3 groups with 30 observations in each group. 
from itertools import chain, repeat 
grp = list(chain.from_iterable((repeat(number, 30) for number in [1, 2, 3]))) 

# generate random numbers per each groups from normal distribution 
np.random.seed(1004) 

# for 'x1' from group 1, 2 and 3
x1_g1 = np.random.normal(0, 1, 30) 
x1_g2 = np.random.normal(0, 1, 30) 
x1_g3 = np.random.normal(0, 1, 30) 

# for 'x2' from group 1, 2 and 3 
x2_g1 = np.random.normal(10, 1, 30) 
x2_g2 = np.random.normal(10, 1, 30) 
x2_g3 = np.random.normal(10, 1, 30) 

# for 'x3' from group 1, 2 and 3 
x3_g1 = np.random.normal(30, 1, 30) 
x3_g2 = np.random.normal(30, 1, 30) 
x3_g3 = np.random.normal(50, 1, 30) 

# different mean 
x4_g1 = np.random.normal(50, 1, 30) 
x4_g2 = np.random.normal(50, 1, 30) 
x4_g3 = np.random.normal(20, 1, 30) 

# different mean # make a DataFrame with all together 
df = pd.DataFrame({
    'id': id, 'grp': grp, 
    'x1': np.concatenate([x1_g1, x1_g2, x1_g3]), 
    'x2': np.concatenate([x2_g1, x2_g2, x2_g3]), 
    'x3': np.concatenate([x3_g1, x3_g2, x3_g3]), 
    'x4': np.concatenate([x4_g1, x4_g2, x4_g3])}) 
    

df.head()

 

id

grp

x1

x2

x3

x4

1

1

0.594403

10.910982

29.431739

49.232193

2

1

0.402609

9.145831

28.548873

50.434544

3

1

-0.805162

9.714561

30.505179

49.459769

4

1

0.115126

8.885289

29.218484

50.040593

5

1

-0.753065

10.230208

30.072990

49.601211

 

 

위에서 만든 가상의 샘플 데이터를 Greenplum DB에 'sample_tbl' 이라는 이름의 테이블로 생성해보겠습니다.  Python pandas의  to_sql()  메소드를 사용하면 pandas DataFrame을 쉽게 Greenplum DB (또는 PostgreSQL DB)에 uploading 할 수 있습니다.  

 

# creating a table in Greenplum by importing pandas DataFrame
conn = "postgresql://gpadmin:changeme@localhost:5432/demo"

df.to_sql('sample_tbl', 
         conn, 
         schema = 'public', 
         if_exists = 'replace', 
         index = False)
 

 

 

 

Jupyter Notebook에서 Greenplum DB에 접속해서 SQL로 이후 작업을 진행하겠습니다.

(Jupyter Notebook에서 Greenplum DB access 하고 SQL query 실행하는 방법은 rfriend.tistory.com/572 참조하세요)

 

-- 여기서 부터는 Jupyter Notebook에서 실행한 것입니다. --

%load_ext sql

# postgresql://Username:Password@Host:Port/Database
%sql postgresql://gpadmin:changeme@localhost:5432/demo
[Out][
'Connected: gpadmin@demo'

 

 

 

위 (1) 에서 pandas 의 to_sql() 로 importing 한 sample_tbl 테이블에서 5개 행을 조회해보겠습니다. 

 

%sql select * from sample_tbl order by id limit 5;
 * postgresql://gpadmin:***@localhost:5432/demo
5 rows affected.

[Out]
id	grp	x1	                x2	                x3	                x4
1	1	0.594403067344276	10.9109819091195	29.4317394311833	49.2321928075563
2	1	0.402608708677309	9.14583073327387	28.54887315985  	50.4345438286737
3	1	-0.805162233589535	9.71456131309311	30.5051787625131	49.4597693977764
4	1	0.115125695763445	8.88528940547472	29.2184835450055	50.0405932387396
5	1	-0.753065219532709	10.230207786414 	30.0729900069999	49.6012106088522

 

 

 

(2) 데이터 구조 변경: reshaping from wide to long

 

PL/Python에서 작업하기 쉽도록 테이블 구조를  wide format에서 long format 으로 변경하겠습니다. union all 로 해서 칼럼 갯수 만큼 위/아래로 append  해나가면 되는데요, DB 에서 이런 형식의 데이터를 관리하고 있다면 아마도 이미 long format 으로 관리하고 있을 가능성이 높습니다. (새로운 데이터가 수집되면 계속  insert into 하면서 행을 밑으로 계속 쌓아갈 것이므로...)

 

%%sql
-- reshaping a table from wide to long
drop table if exists sample_tbl_long;
create table sample_tbl_long as (
    select id, grp, 'x1' as col, x1 as val from sample_tbl
    union all 
    select id, grp, 'x2' as col, x2 as val from sample_tbl
    union all 
    select id, grp, 'x3' as col, x3 as val from sample_tbl
    union all 
    select id, grp, 'x4' as col, x4 as val from sample_tbl
) distributed randomly;

 * postgresql://gpadmin:***@localhost:5432/demo
Done.
360 rows affected.



%sql select * from sample_tbl_long order by id, grp, col limit 8;

[Out]
 * postgresql://gpadmin:***@localhost:5432/demo
8 rows affected.
id	grp	col	val
1	1	x1	0.594403067344276
1	1	x2	10.9109819091195
1	1	x3	29.4317394311833
1	1	x4	49.2321928075563
2	1	x1	0.402608708677309
2	1	x2	9.14583073327387
2	1	x3	28.54887315985
2	1	x4	50.4345438286737

 

 

 

(3) 분석 결과 반환 composite type 정의

 

일원분산분석 결과를 반환받을 때 각 분석 대상 변수 별로 (a) F-통계량, (b) p-value 의 두 개 값을 float8 데이터 형태로  반환받는 composite type 을 미리 정의해놓겠습니다. 

%%sql
-- Creating a coposite return type
drop type if exists plpy_anova_type cascade;
create type plpy_anova_type as (
    f_stat float8
    , p_val float8
);

 * postgresql://gpadmin:***@localhost:5432/demo
Done.
Done.

 

 

 

(4)  일원분산분석(one-way ANOVA) PL/Python 사용자 정의함수 정의

 

집단('grp')과 측정값('val')을 input 으로 받고, statsmodels 모듈의 sm.stats.anova_lm() 메소드로 일원분산분석을 하여 결과 테이블에서 'F-통계량'과 'p-value'만 인덱싱해서 반환하는 PL/Python 사용자 정의 함수를 정의해보겠습니다. 

 

%%sql
-- Creating the PL/Python UDF of ANOVA
drop function if exists plpy_anova_func(text[], float8[]);
create or replace function plpy_anova_func(grp text[], val float8[])
returns plpy_anova_type 
as $$
    import pandas as pd  
    import statsmodels.api as sm
    from statsmodels.formula.api import ols
    
    df = pd.DataFrame({'grp': grp, 'val': val})
    model = ols('val ~ grp', data=df).fit()
    anova_result = sm.stats.anova_lm(model, typ=1)
    return {'f_stat': anova_result.loc['grp', 'F'], 
            'p_val': anova_result.loc['grp', 'PR(>F)']}
$$ language 'plpythonu';

 * postgresql://gpadmin:***@localhost:5432/demo
Done.
Done.

 

 

 

(5) 일원분산분석(one-way ANOVA) PL/Python 함수 분산병렬처리 실행

 

PL/Python 사용자 정의함수는 SQL query 문으로 실행합니다. 이때 PL/Python 이 'F-통계량'과 'p-value'를 반환하도록 UDF를 정의했으므로 아래처럼 (plpy_anova_func(grp_arr, val_arr)).* 처럼 ().* 으로 해서 모든 결과('F-통계량' & 'p-value')를 반환하도록 해줘야 합니다. (빼먹고 실수하기 쉬우므로 ().*를 빼먹지 않도록 주의가 필요합니다)

 

이렇게 하면 변수별로 segment nodes 에서 분산병렬로 각각 처리가 되므로, 변수가 수백~수천개가 있더라도 (segment nodes가 많이 있다는 가정하에) 분산병렬처리되어 신속하게 분석을 할 수 있습니다. 그리고 결과는 바로 Greenplum DB table에 적재가 되므로 이후의 application이나 API service에서 가져다 쓰기에도 무척 편리합니다. 

 

%%sql
-- Executing the PL/Python UDF of ANOVA
drop table if exists plpy_anova_table;
create table plpy_anova_table as (
    select 
        col
        , (plpy_anova_func(grp_arr, val_arr)).*
    from (
        select
            col
            , array_agg(grp::text order by id) as grp_arr
            , array_agg(val::float8 order by id) as val_arr
        from sample_tbl_long
        group by col
    ) a
) distributed randomly;

 * postgresql://gpadmin:***@localhost:5432/demo
Done.
4 rows affected.

 

 

총 4개의 각 변수별 일원분산분석 결과를 조회해보면 아래와 같습니다. 

%%sql
select * from plpy_anova_table order by col;

[Out]
 * postgresql://gpadmin:***@localhost:5432/demo
4 rows affected.
col	f_stat	p_val
x1	0.773700830155438	0.46445029458511966
x2	0.20615939957339052	0.8140997216173114
x3	4520.512608893724	1.2379278415456727e-88
x4	9080.286130418674	1.015467388498996e-101

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 샘플 크기가 다른 2개 이상의 집단에 대해 평균의 차이가 존재하는지를 검정하는 일원분산분석(one-way ANOVA)에 대해 scipy 모듈의 scipy.stats.f_oneway() 메소드를 사용해서 분석하는 방법(rfriend.tistory.com/638)을 소개하였습니다. 

 

이번 포스팅에서는 2개 이상의 집단에 대해 pandas DataFrame에 들어있는 여러 개의 숫자형 변수(one-way ANOVA for multiple numeric variables in pandas DataFrame) 별로 일원분산분석 검정(one-way ANOVA test)을 하는 방법을 소개하겠습니다. 

 

숫자형 변수와 집단 변수의 모든 가능한 조합을 MultiIndex 로 만들어서 statsmodels.api 모듈의 stats.anova_lm() 메소드의 모델에 for loop 순환문으로 변수를 바꾸어 가면서 ANOVA 검정을 하도록 작성하였습니다. 

 

 

 

먼저, 3개의 집단('grp 1', 'grp 2', 'grp 3')을 별로 'x1', 'x2', 'x3, 'x4' 의 4개의 숫자형 변수를 각각 30개씩 가지는 가상의 pandas DataFrame을 만들어보겠습니다. 이때 숫자형 변수는 모두 정규분포로 부터 난수를 발생시켜 생성하였으며, 'x3'와 'x4'에 대해서는 집단3 ('grp 3') 의 평균이 다른 2개 집단의 평균과는 다른 정규분포로 부터 난수를 발생시켜 생성하였습니다.  

 

아래의 가상 데이터셋은 결측값이 없이 만들었습니다만, 실제 기업에서 쓰는 데이터셋에는 혹시 결측값이 존재할 수도 있으므로 결측값을 없애거나 또는 결측값을 그룹 별 평균으로 대체한 후에 one-way ANOVA 를 실행하기 바랍니다. 

 

## Creating sample dataset
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# generate 90 IDs
id = np.arange(90) + 1

# Create 3 groups with 30 observations in each group.
from itertools import chain, repeat
grp = list(chain.from_iterable((repeat(number, 30) for number in [1, 2, 3])))

# generate random numbers per each groups from normal distribution
np.random.seed(1004)

# for 'x1' from group 1, 2 and 3
x1_g1 = np.random.normal(0, 1, 30)
x1_g2 = np.random.normal(0, 1, 30)
x1_g3 = np.random.normal(0, 1, 30)

# for 'x2' from group 1, 2 and 3
x2_g1 = np.random.normal(10, 1, 30)
x2_g2 = np.random.normal(10, 1, 30)
x2_g3 = np.random.normal(10, 1, 30)

# for 'x3' from group 1, 2 and 3
x3_g1 = np.random.normal(30, 1, 30)
x3_g2 = np.random.normal(30, 1, 30)
x3_g3 = np.random.normal(50, 1, 30) # different mean

x4_g1 = np.random.normal(50, 1, 30)
x4_g2 = np.random.normal(50, 1, 30)
x4_g3 = np.random.normal(20, 1, 30) # different mean

# make a DataFrame with all together
df = pd.DataFrame({'id': id, 
                   'grp': grp, 
                   'x1': np.concatenate([x1_g1, x1_g2, x1_g3]), 
                   'x2': np.concatenate([x2_g1, x2_g2, x2_g3]), 
                   'x3': np.concatenate([x3_g1, x3_g2, x3_g3]), 
                   'x4': np.concatenate([x4_g1, x4_g2, x4_g3])})
                   
df.head()
[Out] 

id	grp	x1	x2	x3	x4
0	1	1	0.594403	10.910982	29.431739	49.232193
1	2	1	0.402609	9.145831	28.548873	50.434544
2	3	1	-0.805162	9.714561	30.505179	49.459769
3	4	1	0.115126	8.885289	29.218484	50.040593
4	5	1	-0.753065	10.230208	30.072990	49.601211


df[df['grp'] == 3].head()
[Out] 

id	grp	x1	x2	x3	x4
60	61	3	-1.034244	11.751622	49.501195	20.363374
61	62	3	0.159294	10.043206	50.820755	19.800253
62	63	3	0.330536	9.967849	50.461775	20.993187
63	64	3	0.025636	9.430043	50.209187	17.892591
64	65	3	-0.092139	12.543271	51.795920	18.883919

 

 

 

가령, 'x3' 변수에 대해 집단별로 상자 그래프 (Box plot for 'x3' by groups) 를 그려보면, 아래와 같이 집단1과 집단2는 유사한 반면에 집단3은 평균이 차이가 많이 나게 가상의 샘플 데이터가 생성되었음을 알 수 있습니다. 

 

## Boxplot for 'x3' by 'grp'
plt.rcParams['figure.figsize'] = [10, 6]
sns.boxplot(x='grp', y='x3', data=df)
plt.show()

 

 

여러개의 변수에 대해 일원분산분석을 하기 전에, 먼저 이해를 돕기 위해 Python의 statsmodels.api 모듈의 stats.anova_lm() 메소드를 사용해서 'x1' 변수에 대해 집단(집단 1/2/3)별로 평균이 같은지 일원분산분석으로 검정을 해보겠습니다. 

 

    - 귀무가설(H0) : 집단1의 x1 평균 = 집단2의 x1 평균 = 집단3의 x1 평균

    - 대립가설(H1) : 적어도 1개 이상의 집단의 x1 평균이 다른 집단의 평균과 다르다. (Not H0)

 

# ANOVA for x1 and grp
import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols('x1 ~ grp', data=df).fit()
sm.stats.anova_lm(model, typ=1)
[Out]

df	sum_sq	mean_sq	F	PR(>F)
grp	1.0	0.235732	0.235732	0.221365	0.639166
Residual	88.0	93.711314	1.064901	NaN	NaN

 

일원분산분석 결과 F 통계량이 0.221365, p-value가 0.639 로서 유의수준 5% 하에서 귀무가설을 채택합니다. 즉, 3개 집단 간 x1의 평균의 차이는 없다고 판단할 수 있습니다. (정규분포 X ~ N(0, 1) 를 따르는 모집단으로 부터 무작위로 3개 집단의 샘플을 추출했으므로 차이가 없게 나오는게 맞겠습니다.)

 

 

한개의 변수에 대한 일원분산분석하는 방법을 알아보았으니, 다음으로는 3개 집단별로 여러개의 연속형 변수인 'x1', 'x2', 'x3', 'x4' 에 대해서 for loop 순환문으로 돌아가면서 일원분산분석을 하고, 그 결과를 하나의 DataFrame에 모아보도록 하겠습니다. 

 

(1) 먼저, 일원분산분석을 하려는 모든 숫자형 변수와 집단 변수에 대한 가능한 조합의 MultiIndex 를 생성해줍니다. 

 

# make a multiindex for possible combinations of Xs and Group
num_col = ['x1','x2', 'x3', 'x4']
cat_col =  ['grp']
mult_idx = pd.MultiIndex.from_product([num_col, cat_col],
                                   names=['x', 'grp'])

print(mult_idx)
[Out]
MultiIndex([('x1', 'grp'),
            ('x2', 'grp'),
            ('x3', 'grp'),
            ('x4', 'grp')],
           names=['x', 'grp'])
           

 

 

(2) for loop 순환문(for x, grp in mult_idx:)으로 model = ols('{} ~ {}'.format(x, grp) 의 선형모델의  y, x 부분의 변수 이름을 바꾸어가면서 sm.stats.anova_lm(model, typ=1) 로 일원분산분석을 수행합니다. 이렇게 해서 나온 일원분산분석 결과 테이블을 anova_tables.append(anova_table) 로 순차적으로 append 해나가면 됩니다.  

 

# ANOVA test for multiple combinations of X and Group
import statsmodels.api as sm
from statsmodels.formula.api import ols

anova_tables = []
for x, grp in mult_idx:
    model = ols('{} ~ {}'.format(x, grp), data=df).fit()
    anova_table = sm.stats.anova_lm(model, typ=1)
    anova_tables.append(anova_table)

df_anova_tables = pd.concat(anova_tables, keys=mult_idx, axis=0)

df_anova_tables
[Out]

df	sum_sq	mean_sq	F	PR(>F)
x1	grp	grp	1.0	0.235732	0.235732	0.221365	6.391661e-01
Residual	88.0	93.711314	1.064901	NaN	NaN
x2	grp	grp	1.0	0.448662	0.448662	0.415853	5.206912e-01
Residual	88.0	94.942885	1.078896	NaN	NaN
x3	grp	grp	1.0	6375.876120	6375.876120	259.202952	5.779374e-28
Residual	88.0	2164.624651	24.598007	NaN	NaN
x4	grp	grp	1.0	13760.538009	13760.538009	256.515180	8.145953e-28
Residual	88.0	4720.684932	53.644147	NaN	NaN

anova tables

 

 

만약 특정 변수에 대한 일원분산분석 결과만을 조회하고 싶다면, 아래처럼 DataFrame의 MultiIndex 에 대해 인덱싱을 해오면 됩니다. 가령, 'x3' 에 대한 집단별 평균 차이 여부를 검정한 결과는 아래처럼 인덱싱해오면 됩니다. 

 

## Getting values of 'x3' from ANOVA tables
df_anova_tables.loc[('x3', 'grp', 'grp')]
[Out]

df         1.000000e+00
sum_sq     6.375876e+03
mean_sq    6.375876e+03
F          2.592030e+02
PR(>F)     5.779374e-28
Name: (x3, grp, grp), dtype: float64

 

 

F 통계량과 p-value 에 대해서 조회하고 싶으면 위의 결과에서 DataFrame 의 칼럼 이름으로 선택해오면 됩니다. 

 

# F-statistic
df_anova_tables.loc[('x3', 'grp', 'grp')]['F']
[Out]
259.2029515179077


# P-value
df_anova_tables.loc[('x3', 'grp', 'grp')]['PR(>F)']
[Out]
5.7793742588216585e-28

 

 

 

MultiIndex 를 인덱싱해오는게 좀 불편할 수 도 있는데요, 이럴 경우  df_anova_tables.reset_index() 로  MultiIndex 를 칼럼으로 변환해서 사용할 수도 있습니다. 

# resetting index to columns
df_anova_tables_2 = df_anova_tables.reset_index().dropna()


df_anova_tables_2
[Out]

level_0	level_1	level_2	df	sum_sq	mean_sq	F	PR(>F)
0	x1	grp	grp	1.0	0.235732	0.235732	0.221365	6.391661e-01
2	x2	grp	grp	1.0	0.448662	0.448662	0.415853	5.206912e-01
4	x3	grp	grp	1.0	6375.876120	6375.876120	259.202952	5.779374e-28
6	x4	grp	grp	1.0	13760.538009	13760.538009	256.515180	8.145953e-28

 

 

Greenplum DB에서 PL/Python (또는 PL/R)을 사용하여 여러개의 숫자형 변수에 대해 일원분산분석을 분산병렬처리하는 방법(one-way ANOVA in parallel using PL/Python on Greenplum DB)은 rfriend.tistory.com/640 를 참고하세요. 

 

 

[reference] 

* ANOVA test using Python statsmodels
 
: https://www.statsmodels.org/stable/generated/statsmodels.stats.anova.anova_lm.html

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 2개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 2개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  이론적인 부분에 대한 좀더 자세한 내용은 https://rfriend.tistory.com/131 를 참고하세요. 

 

one-way ANOVA  일원분산분석

 

이번 포스팅에서는 Python의  scipy 모듈의 stats.f_oneway() 메소드를 사용하여 샘플의 크기가 서로 다른 3개 그룹 간 평균에 차이가 존재하는지 여부를 일원분산분석(one-way ANOVA)으로 분석하는 방법을 소개하겠습니다. 

 

분산분석(Analysis Of Variance) 검정의 3가지 가정사항을 고려해서, 샘플 크기가 서로 다른 가상의 3개 그룹의 예제 데이터셋을 만들어보겠습니다. 

 

[ 분산분석  검정의 가정사항 (assumptions of ANOVA test) ]

  (1) 독립성: 각 샘플 데이터는 서로 독립이다. 
  (2) 정규성: 각 샘플 데이터는 정규분포를 따르는 모집단으로 부터 추출되었다. 
  (3) 등분산성: 그룹들의 모집단의 분산은 모두 동일하다. 

 

먼저, 아래의 예제 샘플 데이터셋은 그룹1과 그룹2의 평균은 '0'으로 같고, 그룹3의 평균은 '5'로서 나머지 두 그룹과 다르게 난수를 발생시켜 가상으로 만든 것입니다. 

 

# 3 groups of dataset with different sized samples 
import numpy as np
import pandas as pd
np.random.seed(1004)

data1 = np.random.normal(0, 1, 50)
data2 = np.random.normal(0, 1, 40)
data3 = np.random.normal(5, 1, 30) # different mean

data123 = [data1, data2, data3]


print(data123)
[Out]
[array([ 0.59440307,  0.40260871, -0.80516223,  0.1151257 , -0.75306522,
       -0.7841178 ,  1.46157577,  1.57607553, -0.17131776, -0.91448182,
        0.86013945,  0.35880192,  1.72965706, -0.49764822,  1.7618699 ,
        0.16901308, -1.08523701, -0.01065175,  1.11579838, -1.26497153,
       -1.02072516, -0.71342119,  0.57412224, -0.45455422, -1.15656742,
        1.29721355, -1.3833716 ,  0.3205909 , -0.59086187, -1.43420648,
        0.60998011,  0.51266756,  1.9965168 ,  1.42945668,  1.82880165,
       -1.40997132,  0.49433367,  0.9482873 , -0.35274099, -0.15359935,
       -1.18356064, -0.75440273, -0.85981073,  1.14256322, -2.21331694,
        0.90651805,  2.23629   ,  1.00743665,  1.30584548,  0.46669171]), array([-0.49206651, -0.08727244, -0.34919043, -1.11363541, -1.71982966,
       -0.14033817,  0.90928317, -0.60012686,  1.03906073, -0.03332287,
       -1.03424396,  0.15929405,  0.33053582,  0.02563551, -0.09213904,
       -0.91851177,  0.3099129 , -1.24211638, -0.33113027, -1.64086666,
       -0.27539834, -0.05489003,  1.50604364, -1.37756156, -1.25561652,
        0.16120867, -0.42121705,  0.2341905 , -1.20155195,  1.48131392,
        0.29105321,  0.4022031 , -0.41466037,  1.00502917,  1.45376705,
       -0.07038153,  0.52897801, -2.37895295, -0.75054747,  1.10641762]), array([5.91098191, 4.14583073, 4.71456131, 3.88528941, 5.23020779,
       5.12814125, 3.44610618, 5.36308351, 4.69463298, 7.49521024,
       5.41246681, 3.8724271 , 4.60265531, 4.60082925, 4.9074518 ,
       3.8141367 , 6.4457503 , 4.27553455, 3.63173152, 6.25540542,
       3.77536981, 7.19435668, 6.25339789, 6.43469547, 4.27431061,
       5.16694916, 7.21065725, 5.68274021, 4.81732021, 3.81650656])]

 

 

 

위의 3개 그룹의 커널밀도추정 그래프 (Kernel Density Estimates Plot)를 겹쳐서 그려보면, 아래와 같이 2개 집단은 서로 평균이 비슷하고 나머지 1개 집단은 평균이 확연히 다르다는 것을 직관적으로 알 수 있습니다. 

# Kernel Density Estimate Plot
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['figure.figsize'] = [10, 6]

sns.kdeplot(data1)
sns.kdeplot(data2)
sns.kdeplot(data3)
plt.show()

 

 

상자 그래프 (Box Plot) 으로 3개 집단 간 평균의 위치와 퍼짐 정도를 비교해보면, 역시 아래와 같이 그룹1과 그룹2는 서로 중심위치가 서로 비슷하고 그룹3만 중심위치가 확연히 다름을 알 수 있습니다. 

 

# Boxplot
sns.boxplot(data=data123)
plt.xlabel("Group", fontsize=14)
plt.ylabel("Value", fontsize=14)
plt.show()

 

 

이제 scipy 모듈의 scipy.stats.f_oneway() 메소드를 사용해서 서로 다른 개수의 샘플을 가진 3개 집단에 대해 평균의 차이가 존재하는지 여부를 일원분산분석을 통해 검정을 해보겠습니다. 

 

  - 귀무가설(H0): 3 집단의 평균은 모두 같다. (mu1 = mu2 = m3)

  - 대립가설(H1):  3 집단의 평균은 같지 않다.(적어도 1개 집단의 평균은 같지 않다) (Not H0)

 

F통계량은 매우 큰 값이고 p-value가 매우 작은 값이 나왔으므로 유의수준 5% 하에서 귀무가설을 기각하고 대립가설을 채택합니다. 즉, 3개 집단 간 평균의 차이가 존재한다고 평가할 수 있습니다

# ANOVA with different sized samples using scipy
from scipy import stats

stats.f_oneway(data1, data2, data3)
[Out] 
F_onewayResult(statistic=262.7129127080777, pvalue=5.385523527223916e-44)

 

 

F통계량과 p-value 를 각각 따로 따로 반환받을 수도 있습니다. 

# returning f-statistic and p-value
f_val, p_val = stats.f_oneway(*data123)

print('f-statistic:', f_val)
print('p-vale:', p_val)
[Out]
f-statistic: 262.7129127080777
p-vale: 5.385523527223916e-44

 

 

scipy 모듈의  stats.f_oneway() 메소드를 사용할 때 만약 데이터에 결측값(NAN)이 포함되어 있으면  'NAN'을 반환합니다. 위에서 만들었던 'data1'  numpy array 의 첫번째 값을  np.nan 으로 대체한 후에 scipy.stats.f_oneway() 로 일원분산분석 검정을 해보면  'NAN'(Not A Number)을 반환한 것을 볼 수 있습니다. 

 

# if there is 'NAN', then returns 'NAN'
data1[0] = np.nan
print(data1)
[Out]
array([        nan,  0.40260871, -0.80516223,  0.1151257 , -0.75306522,
       -0.7841178 ,  1.46157577,  1.57607553, -0.17131776, -0.91448182,
        0.86013945,  0.35880192,  1.72965706, -0.49764822,  1.7618699 ,
        0.16901308, -1.08523701, -0.01065175,  1.11579838, -1.26497153,
       -1.02072516, -0.71342119,  0.57412224, -0.45455422, -1.15656742,
        1.29721355, -1.3833716 ,  0.3205909 , -0.59086187, -1.43420648,
        0.60998011,  0.51266756,  1.9965168 ,  1.42945668,  1.82880165,
       -1.40997132,  0.49433367,  0.9482873 , -0.35274099, -0.15359935,
       -1.18356064, -0.75440273, -0.85981073,  1.14256322, -2.21331694,
        0.90651805,  2.23629   ,  1.00743665,  1.30584548,  0.46669171])

# returns 'nan'
stats.f_oneway(data1, data2, data3)
[Out] 
F_onewayResult(statistic=nan, pvalue=nan)

 

 

샘플 데이터에 결측값이 포함되어 있는 경우, 결측값을 먼저 제거해주고 일원분산분석 검정을 실시해주시기 바랍니다. 

 

# get rid of the missing values before applying ANOVA test
stats.f_oneway(data1[~np.isnan(data1)], data2, data3)
[Out]
F_onewayResult(statistic=260.766426640122, pvalue=1.1951277551195217e-43)

 

 

위의  일원분산분석(one-way ANOVA) 는 2개 이상의 그룹 간 평균의 차이가 존재하는지만을 검정할 뿐이며, 집단 간 평균의 차이가 존재한다는 대립가설을 채택하게 된 경우 어느 그룹 간 차이가 존재하는지는 사후검정 다중비교(post-hoc pair-wise multiple comparisons)를 통해서 알 수 있습니다. (rfriend.tistory.com/133)

 

pandas DataFrame 데이터셋에서 여러개의 숫자형 변수에 대해 for loop 순환문을 사용하여 집단 간 평균 차이 여부를 검정하는 방법은 rfriend.tistory.com/639 포스팅을 참고하세요. 

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. 

 

[reference]

* scipy.stats.f_oneway : https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html

 

scipy.stats.f_oneway — SciPy v1.6.3 Reference Guide

G.W. Heiman, “Understanding research methods and statistics: An integrated introduction for psychology”, Houghton, Mifflin and Company, 2001.

docs.scipy.org

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  

 

 

 

[ 일원분산분석 (one-way ANOVA)의 개념 ]

 

 

 

 

분산분석 (ANOVA)은 영국의 통계학자 피셔(R.A. Fisher)가 고안한 분석 기법으로서, 최초에는 농사실험과 관련된 연구에서 사용되었으나 요즘에는 사회과학, 공학, 의학 등 폭넓게 적용되고 있습니다.

 

 

[ 분산분석 (ANOVA)의 분류 ]

 

 

 

3개 이상의 집단에 대한 비교 시 모집단이 정규분포를 따르지 않을 경우에는 비모수 추정과 검정 기법을 이용해야 하며, 이에 대해서는 추후에 별도로 다루도록 하겠습니다.

 

 

분산분석은 측정하고자 하는 값에 영향을 미치는 요인(Factor)의 수에 따라서 구분하는데요, 가령, 작업자(Worker)에 따라서 생산량(Volume of Production)의 차이가 있는가를 비교 분석한다고 했을 때 '작업자(Worker)'가 바로 요인(Factor)이 되겠습니다. 그리고 작업자 3명 '홍길동', '김철희', '최영희'를 요인 수준 (Factor Level) 혹은 처리 (treatment) 라고 합니다.

 

측정값에 영향을 미치는 요인(Factor)이 1 개인 실험은 '일원분산분석' (one-way ANOVA) 라고 하며, 측정값에 영향을 미치는 요인(Factor)이 2 개인 실험은 '이원분산분석' (two-way ANOVA) 라고 부릅니다.

 

 

[ ANOVA Model and Description in R ]

 

n-way ANOVA

Model

Description

one-way ANOVA

y ~ x1

  y is explained by x1 only

two-way ANOVA

y ~ x1 + x2

  y is explained by x1 and x2

two-way ANOVA

y ~ x1 * x2

  y is explained by x1, x2 and the

interaction   
  between them

three-way ANOVA

y ~ x1 + x2 + x3

  y is explained by x1, x2 and x3

 

 

 

이번 포스팅에서는 요인(Factor)이 1개인 일원분산분석(one-way ANOVA)의 이론과 R의 aov() 함수의 사용법에 대해서 설명해보겠습니다.

 

요인(Factor)이 1개이면서 r개의 요인 수준(Factor Level)을 가지고 있고, 각 수준에서 박복 측정수가 n개인 일원분산분석(one-way ANOVA)는 다음과 같은 형태의 자료 형태와 모형을 가지게 됩니다.

 

 

 

[ 일원분산분석 데이터 형태 (dataset for one-way ANOVA) ]

 

 

 

 

[ 일원분산분석 모형 (one-way ANOVA model) ] 

 

one-way ANOVA model

 

분산분석(ANOVA)은 측정된 자료값들의 전체 변동을 비교하고자 하는 요인 수준(Factor Level) 간의 차이에 의해서 발생하는 변동과 그 밖의 요인에 의해서 발생하는 변동으로 나누어서 자료를 분석하는 것이 기본 원리가 되겠습니다.

 

 

[ 분산분석 기본 원리 1 ]

 

 

 

위에서 제시한 식의 양변을 제곱하여 전체 측정값들에 대해서 모두 더하여 식을 정리하면 아래와 같이 됩니다. (자세한 과정은 생략함)

 

 

[ 분산분석 기본 원리 2 ]

 

 

 

 

위의 총제곱합(SST)와 처리제곱합(SSTR), 오차제곱합(SSE) 간의 기본원리를 이용하여 분산분석에서 사용하는 통계량들을 표로 일목요연하게 정리한 것이 아래의 일원분산분석표(one-way ANOVA table) 가 되겠습니다.

 

 

[ 원분산분석표 (one-way ANOVA table) ]

 

 구분

 제곱합

(Sum of Squares)

자유도

(degrees of freedom) 

평균제곱

(Mean Square Error) 

검정통계량 F0

(F statistics)

 처리

(Treatment)

 SSTR

r-1 

MSTR 

MSTR/MSE 

 오차

(Error)

 SSE

 nT - r

 MSE

 

 전체

(Total)

 SST

 nT - 1

 

 

 

 

분산분석표의 제일 오른쪽에 있는 F0 통계량을 이용해서 전체 수준들간의 평균이 같은지 아니면 다른지를 검정합니다.  기본 개념을 설명하자면, F0 통계량은 처리평균제곱 (MSTR)의 크기에 영향을 받으므로 처리평균제곱 (MSTR)이 커지면 오차평균제곱(MSE)은 작아지게 되며, 따라서 F0 통계량은 분자가 커지고 분모가 작아지므로 당연히 커지게 됩니다. 즉, F0 통계량 값이 크다는 것은 수준 평균들간에 평균의 차이가 크다는 의미가 되겠습니다.

 

 

그러면, 요인효과에 대한 검정을 위해서 분산분석에서는 아래와 같은 귀무가설과 대립가설을 사용하며, 검정통계량으로는 F 를 사용하여 기각역 또는 P-value 에 의해서 검정을 하게 됩니다.  

 

 

  [ 가설 ]

   귀무가설 H0 : μ1 = μ2 = ... = μr

   대립가설 H1 : 모든 μi 는 같지 않다.  i = 1, 2, ..., r

 

  [ 기각역 혹은 P-value에 의한 검정 ]

   표본으로부터 계산된 검정통계량의 값 f0가 유의수준(significance level) α 에서

 

   f0 >= Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) < α 이면, H0 기각 (H0 reject)

 

   f0 < Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) >= α 이면, H0 채택 (H0 accept)

 

 

 

 

이론 설명이 무척 길어졌는데요, 이제 드디어, 아래 문제를 R의 aov() 함수를 사용해서 풀어보도록 하겠습니다.

 

 

 


 

 

문제 ) 정유사에서 온도(Factor, one-way)에 따라서 휘발유 생산량에 변화가 있는지 (즉, 영향이 있는지) 알아보기 위하여 온도를 3가지 조건(3 Factor Levels)으로 실험설계를 하여 10번에 걸쳐서 휘발유 생산량을 측정하였습니다. 관찰값이 아래와 같을 때 조사되었을 때 온도의 조건에 따라서 휘발유 생산량에 차이가 있는지 유의수준 α = 10% 로 검정하시오.

 

 

 

 10번 실험에 의한 측정값

 요인 수준

(Factor Level)

 1

2

 10

 온도 조건 1

(Temp condition 1)

 50.5

52.1

51.9

52.4 

50.6

51.4 

51.2 

52.2 

51.5 

50.8 

 온도 조건 2

(Temp condition 2)

 47.5

47.7 

46.6 

47.1 

47.2 

47.8 

45.2 

47.4 

45.0 

47.9 

 온도 조건 3

(Temp condition 3)

 46.0

47.1

45.6

47.1

47.2

46.4

45.9 

47.1

44.9

46.2 

 


 

R에 (1) 위의 관측값들을 벡터로 입력하고,

      (2) boxplot() 함수와 summary()를 이용하여 탐색적 분석을 해본 후에,

      (3) aov() 함수를 이용하여 one-way ANOVA 분석을 하였으며 
         : aov(y ~ group, data = dataset)

      (4) 기본 가정 중에 오차의 등분산성을 검정하기 위해 Bartlett 검정

         : bartlett.test(y ~ group, data = dataset)

을 실시하였습니다.

 

 

이때 조심해야 할 것이 있는데요, dataset을 데이터프레임으로 해서 그룹(group, factor level)에 해당하는 변수는 반드시 요인(Factor)형이어야 하므로, 만약 요인(Factor)형이 아니라면 먼저 transfrom(factor()) 함수를 이용해서 변환을 시켜주고 ANOVA 분석을 수행해야 합니다.

 

 

> ##---------------------------------------------------------- 
> ## One-way ANOVA : aov(), oneway.test 
> ##---------------------------------------------------------- 
>  
> ##--- Are there any daily outcome differences among temperature conditions? 
> # group 1 : temperature condition 1  
> # group 2 : temperature condition 2 
> # group 3 : temperature condition 3 
>  
> # daily outcome by tmep condition (group 1/2/3) 
> y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8) 
> y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9) 
> y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2) 
>  
> y <- c(y1, y2, y3) 
> y  
[1] 50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 47.5 47.7 46.6 47.1 47.2 47.8 45.2 47.4 45.0 
[20] 47.9 46.0 47.1 45.6 47.1 47.2 46.4 45.9 47.1 44.9 46.2 
>  
> n <- rep(10, 3) > n [1] 10 10 10 
>  
> group <- rep(1:3, n) 
> group  [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 
>  
> # combining into data.frame 
> group_df <- data.frame(y, group) 
> group_df       
y group 
1  1 50.5     1 2  52.1     1 3  51.9     1 4  52.4     1 5  50.6     1 6  51.4     1 7  51.2     1 8  52.2     1 9  51.5     1 10 50.8     1 11 47.5     
2 12 47.7     2 13 46.6     2 14 47.1     2 15 47.2     2 16 47.8     2 17 45.2     2 18 47.4     2 19 45.0     2 20 47.9     2 21 46.0     
3 22 47.1     3 23 45.6     3 24 47.1     3 25 47.2     3 26 46.4     3 27 45.9     3 28 47.1     3 29 44.9     3 30 46.2     3 
>  
> sapply(group_df, class)         
y     group  
"numeric" "integer"  
>  
> # transform from 'integer' to 'factor' 
> group_df <- transform(group_df, group = factor(group)) 
> sapply(group_df, class)         
y     group  
"numeric"  "factor"  
>  
>  
> # boxplot 
> attach(group_df) 
The following objects are masked _by_ .GlobalEnv:      group, y  
> boxplot(y ~ group,  
+         main = "Boxplot of Daily Outcome by Temperature condition 1/2/3",  
+         xlab = "Factor Levels : Temperature condition 1/2/3",  
+         ylab = "Daily Outcome") 
>  

 

 

 

>  
> # descriptive statistics by group 
> tapply(y, group, summary) 
$`1`    
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
50.50   50.90   51.45   51.46   52.05   52.40   
$`2`    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
45.00   46.72   47.30   46.94   47.65   47.90   
$`3`    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
44.90   45.92   46.30   46.35   47.10   47.20   
>
>
> detach(group_df) 
>  
> # one-wayANOVA 
> aov(y ~ group, data = group_df) 
Call:    aov(formula = y ~ group, data = group_df)  
Terms:                   
group Residuals Sum of Squares  156.302    19.393 
Deg. of Freedom       2        27  
Residual standard error: 0.8475018 Estimated effects may be unbalanced 

 

> summary(aov(y ~ group, data = group_df))             
Df Sum Sq Mean Sq F value  Pr(>F)     
group        2 156.30   78.15   108.8 1.2e-13 *** Residuals   27  19.39    0.72                     
--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

일원분산분석(one-way ANOVA) 결과 검정통계량 F-value가 108.8으로 나왔으며, P-value 값은 '1.2e-13'으로서 매우 작게 나와 유의수준 10% 에서 귀무가설을 기각하고 대립가설을 채택하게 되어 "온도 조건 1/2/3에 따라서 휘발유 생산량에 차이가 있다"라고 판단할 수 있겠습니다.

 

 오차의 등분산성 가정에 대해 Bartlett 검정 결과 P-value가 0.4368로서 유의수준 10%에서 귀무가설을 채택하여 "오차의 등분산성 가정을 만족한다"고 할 수 있겠습니다.

 

많은 도움 되었기를 바랍니다.

 

질문은 댓글로 남겨주세요.

 

  • 쌍을 이룬 집단 간 평균 다중비교 (multiple comparisons)

Tukey's HSD(honestly significant difference) test 참조

Duncan's LSR(least significant range) test 참고

 

 

  • 대비 (contrast)

샤페 검정법 (scheffe test) 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,