Isolation forest을 이용한 이상탐지

1. 이상 탐지(Anomaly detection)

이상 탐지(anomaly detection)란 자료에서 다른 패턴을 보이는 자료를 찾는 것을 말합니다. 이런 데이터를 이상값(anomaly)라 하며 이상 탐지는 사기 탐지, 침입 탐지, 안전 관리를 포함한 다양한 분야에 널리 활용된다.

2. 이상의 종류

2.1. Point anomaly

데이터셋 내에 하나의 데이터가 나머지에 대해 이상하다고 판단되는 경우, 흔히 아웃라이어(Outlier)라고 부른다.

2.2. Collective anomaly

데이터셋 내에 여러 데이터 포인트가 이상하다고 판단되는 경우

2.3. Contextual anomaly

전체적인 데이터셋의 맥락을 고려했을때 이상하다고 판단되는 경우

3. Isolation forest

기계학습으로 이상을 탐지하는 다양한 알고리즘이 존재하고 문제 마다 적합한 알고리즘을 선택하는 것이 중요하다. 여기에서는 밀도기반으로 이상 탐지를 하는 Isolation forest의 예제를 배운다.

Isolation forest는 기본적으로 데이터셋을 의사결정나무(Decision Tree) 형태로 표현해 정상값을 분리하기 위해서는 의사결정나무를 깊숙하게 타고 내려가야 하고, 반대로 이상값은 의사결정나무 상단부에서 분리할 수 있다는 것을 이용한다.

이 특성을 사용해 의사결정나무를 몇 회 타고 내려가야 분리되는가를 기준으로 정상과 이상을 분리한다.

Isolation forest의 장점

  • 군집기반 이상탐지 알고리즘에 비해 계산량이 매우 적다
  • 강건한(Robust)한 모델을 만들 수 있다

4. 예제 살펴보기

4.1. 가장 단순한 예제

먼저 사용할 라이브러리를 불러온다.

In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%matplotlib inline
import seaborn as sns

sns.set_style("darkgrid")
from sklearn.ensemble import IsolationForest

간단한 데이터를 생성한다.

In [2]:
rng = np.random.RandomState(42)
# Generating training data
X_train = 0.2 * rng.randn(1000, 2)
X_train = np.r_[X_train + 3, X_train]
X_train = pd.DataFrame(X_train, columns=["x1", "x2"])
# Generating new, 'normal' observation
X_test = 0.2 * rng.randn(200, 2)
X_test = np.r_[X_test + 3, X_test]
X_test = pd.DataFrame(X_test, columns=["x1", "x2"])
# Generating outliers
X_outliers = rng.uniform(low=-1, high=5, size=(50, 2))
X_outliers = pd.DataFrame(X_outliers, columns=["x1", "x2"])

생성한 데이터셋을 시각화해 살펴본다.

In [3]:
plt.rcParams["figure.figsize"] = [10, 10]
p1 = plt.scatter(
    X_train.x1,
    X_train.x2,
    c="white",
    s=20 * 4,
    edgecolor="k",
    label="training observations",
)
# p2 = plt.scatter(X_test.x1, X_test.x2, c='green', s=20*4, edgecolor='k', label='new regular obs.')
p3 = plt.scatter(
    X_outliers.x1,
    X_outliers.x2,
    c="red",
    s=20 * 4,
    edgecolor="k",
    label="new abnormal obs.",
)

plt.legend()
Out[3]:
<matplotlib.legend.Legend at 0xd059518>
No description has been provided for this image

scikit-learnIsolationForest 함수를 이용해 학습 모델을 만든다.

In [4]:
clf = IsolationForest(max_samples=100, contamination=0.1, random_state=42)
clf.fit(X_train)
y_pred_train = clf.predict(X_train)
y_pred_test = clf.predict(X_test)
y_pred_outliers = clf.predict(X_outliers)

학습된 모델을 가지고 X_outliers 데이터셋의 이상 탐지를 해보고 시각화를 통해 나타낸다.

In [5]:
X_outliers = X_outliers.assign(y=y_pred_outliers)
p1 = plt.scatter(
    X_train.x1,
    X_train.x2,
    c="white",
    s=20 * 4,
    edgecolor="k",
    label="training observations",
)
p2 = plt.scatter(
    X_outliers.loc[X_outliers.y == -1, ["x1"]],
    X_outliers.loc[X_outliers.y == -1, ["x2"]],
    c="red",
    s=20 * 4,
    edgecolor="k",
    label="detected outliers",
)
p3 = plt.scatter(
    X_outliers.loc[X_outliers.y == 1, ["x1"]],
    X_outliers.loc[X_outliers.y == 1, ["x2"]],
    c="green",
    s=20 * 4,
    edgecolor="k",
    label="detected regular obs",
)
plt.legend()
Out[5]:
<matplotlib.legend.Legend at 0xd2db630>
No description has been provided for this image

학습 모델의 정확도를 수치로 나타내본다.

In [6]:
print("테스트 데이터셋에서 정확도:", list(y_pred_test).count(1) / y_pred_test.shape[0])
print(
    "이상치 데이터셋에서 정확도:",
    list(y_pred_outliers).count(-1) / y_pred_outliers.shape[0],
)
테스트 데이터셋에서 정확도: 0.905
이상치 데이터셋에서 정확도: 0.96

4.2. 더 복잡한 예제

좀 더 실용적인 예제로 직업에 따른 급여를 조사한 데이터셋을 가지고 이상 탐지를 해보자. 먼저 필요한 라이브러리와 데이터셋(jobtitle_pay.csv)을 불러온다.

In [7]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest

%matplotlib inline

salary = pd.read_csv("../data/jobtitle_pay.csv")
salary.tail()
df = salary.copy()

생성한 데이터프레임의 모양을 살펴보자.

In [8]:
df.head()
Out[8]:
JobTitle TotalPay
0 Transit Operator 148711.28
1 Transit Operator 146680.39
2 Transit Operator 145024.94
3 Transit Operator 144704.29
4 Transit Operator 144557.78

데이터셋의 첫번째 열에는 직업의 이름(JobTitle)이 두번째 열에는 연봉(TotalPay)가 달러화 기준으로 나와 있다. 이제 데이터셋에 존재하는 직업의 종류와 각각 몇명의 데이터가 있는지 출력한다.

In [9]:
df["JobTitle"].value_counts()
Out[9]:
Transit Operator                9424
Special Nurse                   4389
Registered Nurse                3736
Public Svc Aide-Public Works    2518
Police Officer 3                2421
Custodian                       2418
Firefighter                     2359
Recreation Leader               1971
Patient Care Assistant          1945
Name: JobTitle, dtype: int64

위 결과를 통해 총 9개의 직업이 많개는 9424명의 데이터에서 적게는 1945명의 데이터를 포함하고 있다는 것을 알 수 있다.

각 직업별로 박스플롯을 그려 평균 연봉이 어느정도 되는지 시각화해보자.

In [10]:
plt.figure(figsize=(10, 8))
plt.xticks(rotation=45, horizontalalignment="right")
sns.boxplot(x="JobTitle", y="TotalPay", data=df)
Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0xd118908>
No description has been provided for this image

위 시각화를 통해 미국에서는 소방관(Firefighter)의 연봉이 높다는 것과 연봉의 분포가 일정하지 않다는 것을 알 수 있다.

기계학습을 위해 위의 데이터를 LabelEncoder를 통해 전처리 해줘야한다.

In [11]:
from sklearn.preprocessing import LabelEncoder

job_encode = LabelEncoder()
salary["JobTitle"] = job_encode.fit_transform(salary["JobTitle"])
salary.tail()
Out[11]:
JobTitle TotalPay
31176 8 15.35
31177 6 6.00
31178 0 0.00
31179 0 0.00
31180 0 0.00

기계 학습모델을 만들어 학습을 수행한다.

In [12]:
model = IsolationForest(
    n_estimators=100, max_samples="auto", n_jobs=-1, max_features=2, contamination=0.01
)
model.fit(salary.to_numpy())
Out[12]:
IsolationForest(behaviour='deprecated', bootstrap=False, contamination=0.01,
                max_features=2, max_samples='auto', n_estimators=100, n_jobs=-1,
                random_state=None, verbose=0, warm_start=False)

데이터프레임에 평가 점수(score)와 이상(anomaly) 판단 여부에 대한 값을 추가한다.

In [13]:
score = model.decision_function(salary.to_numpy())
anomaly = model.predict(salary.to_numpy())
df["scores"] = score
df["anomaly"] = anomaly
anomaly_data = df.loc[df["anomaly"] == -1]  # 이상값은 -1으로 나타낸다.
anomaly_data
Out[13]:
JobTitle TotalPay scores anomaly
2388 Firefighter 290076.13 -0.063016 -1
2389 Firefighter 267951.41 -0.061552 -1
2390 Firefighter 259740.36 -0.061552 -1
2391 Firefighter 245399.19 -0.055934 -1
2392 Firefighter 246369.02 -0.059329 -1
... ... ... ... ...
31108 Custodian 192.04 -0.024968 -1
31155 Custodian 54.39 -0.024505 -1
31178 Custodian 0.00 -0.024505 -1
31179 Custodian 0.00 -0.024505 -1
31180 Custodian 0.00 -0.024505 -1

309 rows × 4 columns

총 31180개 데이터 중 312개의 데이터가 이상치로 판별된다. 이제 직업별로 이상치가 얼마나 되는지 살펴보자.

In [14]:
anomaly_data["JobTitle"].value_counts()
Out[14]:
Firefighter         183
Custodian            71
Police Officer 3     51
Special Nurse         3
Registered Nurse      1
Name: JobTitle, dtype: int64

총 9개의 직업에서 5개 직업에 이상치가 발견되며, 소방관에서 가장 많이 분포한다. 이제 시각화를 통해 이상치가 어떻게 분포하는지 살펴보자.

In [15]:
plt.figure(figsize=(10, 8))
plt.xticks(rotation=45, horizontalalignment="right")
sns.stripplot(x="JobTitle", y="TotalPay", data=anomaly_data, jitter=True)
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0xd2ae278>
No description has been provided for this image

시각화를 해보면 소방관과 후견인(Custodian)에서 연봉이 유난히 낮은 사람들을 볼 수 있다. 후견인의 경우는 생업을 위한 직업이라고는 볼 수 없기 때문에 그런것으로 보여진다.

In [21]:
anomaly_data[
    (anomaly_data["JobTitle"] == "Custodian") & (anomaly_data["TotalPay"] <= 1000)
].sort_values(by=["TotalPay"], axis=0)
Out[21]:
JobTitle TotalPay scores anomaly
31180 Custodian 0.00 -0.024505 -1
31178 Custodian 0.00 -0.024505 -1
11784 Custodian 0.00 -0.024505 -1
19727 Custodian 0.00 -0.024505 -1
19754 Custodian 0.00 -0.024505 -1
31179 Custodian 0.00 -0.024505 -1
21319 Custodian 41.95 -0.024505 -1
31155 Custodian 54.39 -0.024505 -1
31106 Custodian 135.58 -0.024968 -1
21226 Custodian 152.85 -0.024968 -1
31093 Custodian 165.48 -0.024968 -1
21268 Custodian 167.17 -0.024968 -1
31108 Custodian 192.04 -0.024968 -1
11670 Custodian 380.36 -0.021584 -1
31042 Custodian 427.50 -0.021584 -1
21204 Custodian 462.26 -0.021124 -1
30978 Custodian 597.20 -0.020204 -1
21158 Custodian 644.77 -0.020204 -1
21146 Custodian 685.91 -0.020204 -1
21128 Custodian 755.10 -0.020204 -1
30856 Custodian 840.26 -0.020204 -1
21091 Custodian 841.94 -0.020204 -1
11445 Custodian 845.72 -0.020204 -1
30930 Custodian 916.07 -0.019354 -1

위 결과는 후견인들중 연봉이 1000달러 미만인 사람의 데이터를 정렬한 것이다. 연봉 1000달러도 사실 너무 낮지만 연봉을 받지 않는 사람들도 존재하고 있다. 이들의 데이터는 상식적으로 직업별 연봉 데이터셋에 불필요한 이상치라고 볼 수 있다.

5. 마치며,

이 글에서 이상이 무엇인지, 그리고 Isolation forest 알고리즘을 사용해 감지하는 방법을 배웠다. 박스 플롯을 사용해 데이터셋을 살펴보고 마지막으로 Isolation forest 알고리즘을 구현했다.

추가적으로 더 알아보고 싶다면 다음 링크를 확인하라.

통계 분석 주석넣기

0. statannot 소개

statannot은 파이썬 시각화 라이브러리 seaborn로 그린 박스플롯에 통계분석에 대한 주석을 자동으로 달아주는 파이썬 패키지입니다. 자세한 것은 공식 홈페이지를 참고하세요.

막대그래프에 주석을 추가하는 것도 가능하지만 아직 완전하지 않습니다.

1. statannot 특징

  • scipy.stats 메소드를 사용한 통계분석:
    • Mann-Whitney
    • t-test (independent and paired)
    • Welch's t-test
    • Levene test
    • Wilcoxon test
    • Kruskal-Wallis test
  • 플롯의 내부 또는 외부에 주석 넣기
  • 통계분석의 주석을 별표(*), p-value 등으로 사용자가 지정할 수 있음
  • 선택적으로 사용자 정의 p- value를 입력값으로 제공(이 경우 통계 분석은 수행되지 않음).

2. statannot 의존성

    Python >= 3.5
    numpy >= 1.12.1
    seaborn >= 0.8.1
    matplotlib >= 2.2.2
    pandas >= 0.23.0
    scipy >= 1.1.0

3. statannot 설치

PyPI를 통해 설치:

pip install statannot

개발중인 최신 버전은 다음 명령어로 설치:

pip install git+https://github.com/webermarcolivier/statannot.git

4. 사용 예시

먼저 필요한 라이브러리를 불러옵니다.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from statannot import add_stat_annotation

sns.set(style="ticks")

4.1 Boxplot non-hue

4.1.1 다중 비교 교정

기본적으로 본페로니 교정이 적용됩니다.

In [2]:
df = sns.load_dataset("tips")
x = "day"
y = "total_bill"
order = ["Sun", "Thur", "Fri", "Sat"]
ax = sns.boxplot(data=df, x=x, y=y, order=order)
ax, test_results = add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    order=order,
    box_pairs=[("Thur", "Fri"), ("Thur", "Sat"), ("Fri", "Sun")],
    test="Mann-Whitney",
    text_format="star",
    loc="outside",
    verbose=2,
)
sns.despine()  # 필요없는 axis border 제거
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

Thur v.s. Fri: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.000e+00 U_stat=6.305e+02
Thur v.s. Sat: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.407e-01 U_stat=2.180e+03
Sun v.s. Fri: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=8.041e-02 U_stat=9.605e+02
No description has been provided for this image

4.1.2 통계 분석 결과

add_stat_annotation 함수는 튜플 ax, test_results 를 반환합니다. 여기서 test_results는 원래 데이터와 통계 분석 결과 (p-value 등)를 모두 포함하는StatResult 리스트 객체입니다.

In [3]:
for res in test_results:
    print(res)

print("\nStatResult attributes:", test_results[0].__dict__.keys())
Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.000e+00 U_stat=6.305e+02
Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.407e-01 U_stat=2.180e+03
Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=8.041e-02 U_stat=9.605e+02

StatResult attributes: dict_keys(['test_str', 'test_short_name', 'stat_str', 'stat', 'pval', 'box1', 'box2'])

4.1.3. 다중 비교 교정이 없는 경우

In [4]:
x = "day"
y = "total_bill"
order = ["Sun", "Thur", "Fri", "Sat"]
ax = sns.boxplot(data=df, x=x, y=y, order=order)
test_results = add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    order=order,
    box_pairs=[("Thur", "Fri"), ("Thur", "Sat"), ("Fri", "Sun")],
    test="Mann-Whitney",
    comparisons_correction=None,
    text_format="star",
    loc="outside",
    verbose=2,
)
sns.despine()  # 필요없는 axis border 제거
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

Thur v.s. Fri: Mann-Whitney-Wilcoxon test two-sided, P_val=6.477e-01 U_stat=6.305e+02
Thur v.s. Sat: Mann-Whitney-Wilcoxon test two-sided, P_val=4.690e-02 U_stat=2.180e+03
Sun v.s. Fri: Mann-Whitney-Wilcoxon test two-sided, P_val=2.680e-02 U_stat=9.605e+02
No description has been provided for this image

4.1.4. 주석의 위치

통계 분석 주석은 플롯 안쪽(loc = 'inside') 또는 바깥쪽(loc ='outside')에 위치할 수 있습니다. 아래 예시는 안쪽 있는 예시입니다.

In [5]:
x = "day"
y = "total_bill"
order = ["Sun", "Thur", "Fri", "Sat"]
ax = sns.boxplot(data=df, x=x, y=y, order=order)
add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    order=order,
    box_pairs=[("Sun", "Thur"), ("Sun", "Sat"), ("Fri", "Sun")],
    perform_stat_test=False,
    pvalues=[0.1, 0.1, 0.001],
    test=None,
    text_format="star",
    loc="inside",
    verbose=2,
)
sns.despine()  # 필요없는 axis border 제거
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

Sun v.s. Thur: Custom statistical test, P_val:1.000e-01
Sun v.s. Fri: Custom statistical test, P_val:1.000e-03
Sun v.s. Sat: Custom statistical test, P_val:1.000e-01
No description has been provided for this image

4.2. Boxplot with hue

ymax 위치가 다른 상자 플롯을 만듭니다.

In [6]:
df = sns.load_dataset("diamonds")
df = df[df["color"].map(lambda x: x in "EIJ")]
# Modifying data to yield unequal boxes in the hue value
df.loc[df["cut"] == "Ideal", "price"] = df.loc[df["cut"] == "Ideal", "price"].map(
    lambda x: min(x, 5000)
)
df.loc[df["cut"] == "Premium", "price"] = df.loc[df["cut"] == "Premium", "price"].map(
    lambda x: min(x, 7500)
)
df.loc[df["cut"] == "Good", "price"] = df.loc[df["cut"] == "Good", "price"].map(
    lambda x: min(x, 15000)
)
df.loc[df["cut"] == "Very Good", "price"] = df.loc[
    df["cut"] == "Very Good", "price"
].map(lambda x: min(x, 3000))
df.head()
Out[6]:
carat cut color clarity depth table price x y z
0 0.23 Ideal E SI2 61.5 55.0 326 3.95 3.98 2.43
1 0.21 Premium E SI1 59.8 61.0 326 3.89 3.84 2.31
2 0.23 Good E VS1 56.9 65.0 327 4.05 4.07 2.31
3 0.29 Premium I VS2 62.4 58.0 334 4.20 4.23 2.63
4 0.31 Good J SI2 63.3 58.0 335 4.34 4.35 2.75
In [7]:
x = "color"
y = "price"
hue = "cut"
hue_order = ["Ideal", "Premium", "Good", "Very Good", "Fair"]
box_pairs = [
    (("E", "Ideal"), ("E", "Very Good")),
    (("E", "Ideal"), ("E", "Premium")),
    (("E", "Ideal"), ("E", "Good")),
    (("I", "Ideal"), ("I", "Premium")),
    (("I", "Ideal"), ("I", "Good")),
    (("J", "Ideal"), ("J", "Premium")),
    (("J", "Ideal"), ("J", "Good")),
    (("E", "Good"), ("I", "Ideal")),
    (("I", "Premium"), ("J", "Ideal")),
]
ax = sns.boxplot(
    data=df,
    x=x,
    y=y,
    hue=hue,
)
add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    hue=hue,
    box_pairs=box_pairs,
    test="Mann-Whitney",
    loc="inside",
    verbose=2,
)
plt.legend(loc="upper left", bbox_to_anchor=(1.03, 1))
sns.despine()  # 필요없는 axis border 제거
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

E_Ideal v.s. E_Premium: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.404e-30 U_stat=3.756e+06
I_Ideal v.s. I_Premium: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=4.627e-60 U_stat=1.009e+06
J_Ideal v.s. J_Premium: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=3.616e-36 U_stat=2.337e+05
E_Ideal v.s. E_Good: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=4.681e-18 U_stat=1.480e+06
I_Ideal v.s. I_Good: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=4.507e-12 U_stat=4.359e+05
J_Ideal v.s. J_Good: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=9.056e-04 U_stat=1.174e+05
E_Ideal v.s. E_Very Good: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.562e-01 U_stat=4.850e+06
E_Good v.s. I_Ideal: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=1.000e+00 U_stat=9.882e+05
I_Premium v.s. J_Ideal: Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction, P_val=4.643e-26 U_stat=8.084e+05
No description has been provided for this image

4.3. 버킷 카테고리가 포함 된 박스 플롯

In [8]:
df = sns.load_dataset("tips")
df["tip_bucket"] = pd.cut(df["tip"], 3)
df.head()
Out[8]:
total_bill tip sex smoker day time size tip_bucket
0 16.99 1.01 Female No Sun Dinner 2 (0.991, 4.0]
1 10.34 1.66 Male No Sun Dinner 3 (0.991, 4.0]
2 21.01 3.50 Male No Sun Dinner 3 (0.991, 4.0]
3 23.68 3.31 Male No Sun Dinner 2 (0.991, 4.0]
4 24.59 3.61 Female No Sun Dinner 4 (0.991, 4.0]
In [9]:
# In this case we just have to pass the list of categories objects to the add_stat_annotation function.
tip_bucket_list = df["tip_bucket"].unique()
tip_bucket_list
Out[9]:
[(0.991, 4.0], (4.0, 7.0], (7.0, 10.0]]
Categories (3, interval[float64]): [(0.991, 4.0] < (4.0, 7.0] < (7.0, 10.0]]
In [10]:
x = "day"
y = "total_bill"
hue = "tip_bucket"
data = df
ax = sns.boxplot(data=df, x=x, y=y, hue=hue)
add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    hue=hue,
    box_pairs=[(("Sat", tip_bucket_list[2]), ("Fri", tip_bucket_list[0]))],
    test="t-test_ind",
    loc="inside",
    verbose=2,
)
plt.legend(loc="upper left", bbox_to_anchor=(1.03, 1))
sns.despine()  # 필요없는 axis border 제거
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

Fri_(0.991, 4.0] v.s. Sat_(7.0, 10.0]: t-test independent samples with Bonferroni correction, P_val=6.176e-07 stat=-7.490e+00
No description has been provided for this image

4.4. y 오프셋 조절하기

In [11]:
df = sns.load_dataset("tips")
x = "day"
y = "total_bill"
hue = "smoker"
ax = sns.boxplot(data=df, x=x, y=y, hue=hue)
add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    hue=hue,
    box_pairs=[
        (("Thur", "No"), ("Fri", "No")),
        (("Sat", "Yes"), ("Sat", "No")),
        (("Sun", "No"), ("Thur", "Yes")),
    ],
    test="t-test_ind",
    text_format="full",
    loc="inside",
    comparisons_correction=None,
    line_offset_to_box=0.2,
    line_offset=0.1,
    line_height=0.05,
    text_offset=8,
    verbose=2,
)
plt.legend(loc="upper left", bbox_to_anchor=(1.03, 1))
sns.despine()  # 필요없는 axis border 제거
Sat_Yes v.s. Sat_No: t-test independent samples, P_val=4.304e-01 stat=7.922e-01
Thur_No v.s. Fri_No: t-test independent samples, P_val=7.425e-01 stat=-3.305e-01
Thur_Yes v.s. Sun_No: t-test independent samples, P_val=5.623e-01 stat=-5.822e-01
No description has been provided for this image

4.5. 사용자 정의 p-value를 입력으로 사용하기

In [12]:
df = sns.load_dataset("iris")
x = "species"
y = "sepal_length"

box_pairs = [
    ("setosa", "versicolor"),
    ("setosa", "virginica"),
    ("versicolor", "virginica"),
]

from scipy.stats import bartlett

test_short_name = "Bartlett"
pvalues = []
for pair in box_pairs:
    data1 = df.groupby(x)[y].get_group(pair[0])
    data2 = df.groupby(x)[y].get_group(pair[1])
    stat, p = bartlett(data1, data2)
    print(
        "Performing Bartlett statistical test for equal variances on pair:",
        pair,
        "stat={:.2e} p-value={:.2e}".format(stat, p),
    )
    pvalues.append(p)
print("pvalues:", pvalues)
Performing Bartlett statistical test for equal variances on pair: ('setosa', 'versicolor') stat=6.89e+00 p-value=8.66e-03
Performing Bartlett statistical test for equal variances on pair: ('setosa', 'virginica') stat=1.60e+01 p-value=6.38e-05
Performing Bartlett statistical test for equal variances on pair: ('versicolor', 'virginica') stat=2.09e+00 p-value=1.48e-01
pvalues: [0.008659557933879902, 6.378941946712463e-05, 0.14778816016231236]
In [13]:
ax = sns.boxplot(data=df, x=x, y=y)
sns.despine()  # 필요없는 axis border 제거
test_results = add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    box_pairs=box_pairs,
    perform_stat_test=False,
    pvalues=pvalues,
    test_short_name=test_short_name,
    text_format="star",
    verbose=2,
)
p-value annotation legend:
ns: 5.00e-02 < p <= 1.00e+00
*: 1.00e-02 < p <= 5.00e-02
**: 1.00e-03 < p <= 1.00e-02
***: 1.00e-04 < p <= 1.00e-03
****: p <= 1.00e-04

setosa v.s. versicolor: Custom statistical test, P_val:8.660e-03
versicolor v.s. virginica: Custom statistical test, P_val:1.478e-01
setosa v.s. virginica: Custom statistical test, P_val:6.379e-05
No description has been provided for this image

4.6. 사용자 정의 주석문구 넣기

In [14]:
df = sns.load_dataset("tips")
x = "day"
y = "total_bill"
order = ["Sun", "Thur", "Fri", "Sat"]
ax = sns.boxplot(data=df, x=x, y=y, order=order)
sns.despine()  # 필요없는 axis border 제거
test_results = add_stat_annotation(
    ax,
    data=df,
    x=x,
    y=y,
    order=order,
    box_pairs=[("Thur", "Fri"), ("Thur", "Sat"), ("Fri", "Sun")],
    text_annot_custom=["first pair", "second pair", "third pair"],
    perform_stat_test=False,
    pvalues=[0, 0, 0],
    loc="outside",
    verbose=0,
)
No description has been provided for this image

5. 마치며

기존에는 번거롭게 플랏에 통계분석에 대한 주석을 수작업으로 넣었지만, statannot 패키지를 사용하면 간편하게 주석처리를 해줍니다. 게다가 기본 설정값도 충실해서 앞으로는 statannot을 계속 사용하게 될 것 같습니다.