statsmodels 기초

smf.ols()는 formula API입니다. "y ~ x"처럼 식을 문자열로 적으면 회귀식을 바로 만들 수 있어서, 처음 배울 때 가장 이해하기 쉽습니다.

작은 예시 데이터부터 보면 아래와 같습니다.

import pandas as pd
import statsmodels.formula.api as smf

df = pd.DataFrame(
    {
        "x": [1, 2, 3, 4],
        "y": [3, 5, 7, 9],
    }
)

result = smf.ols("y ~ x", data=df).fit()

print(result.params)
print(round(result.predict(pd.DataFrame({"x": [5]}))[0], 1))

실행 결과:

Intercept    1.0
x            2.0
dtype: float64
11.0

절편 1.0, 기울기 2.0이라는 뜻입니다. 그래서 x=5를 넣으면 1 + 2*5 = 11에 가까운 예측이 나옵니다.

외부 데이터에 적용하면 아래처럼 실제 회귀를 볼 수 있습니다.

import pandas as pd
import statsmodels.formula.api as smf

data_path = "/data/body.csv"

df = pd.read_csv(data_path)

result = smf.ols('Q("체중 : kg") ~ Q("신장 : cm")', data=df).fit()

print(round(result.params['Q("신장 : cm")'], 3))
print(round(result.rsquared, 3))
print(round(result.predict(pd.DataFrame({"신장 : cm": [170]}))[0], 1))

실행 결과:

1.036
0.533
68.9

신장 : cm 계수 1.036은 키가 1cm 커질 때 예측 몸무게가 약 1.036kg 커지는 방향으로 학습되었다는 뜻입니다. rsquared는 이 모델이 몸무게 변동을 얼마나 설명하는지 보여 주고, 0.533이면 절반 정도를 설명한다고 읽을 수 있습니다.

formula API에서는 절편이 기본으로 자동 포함됩니다. 그래서 Q("체중 : kg") ~ Q("신장 : cm")처럼 써도 절편과 기울기를 함께 추정합니다.

sm.OLS()는 class API입니다. formula 문자열 대신 X와 y를 따로 넣어야 합니다. 이때 절편을 자동으로 넣어 주지 않으므로, sm.add_constant()로 상수항 열을 직접 추가해야 합니다.

import pandas as pd
import statsmodels.api as sm

data_path = "/data/body.csv"

df = pd.read_csv(data_path)

X = sm.add_constant(df[["신장 : cm", "측정나이", "악력D : kg"]])
y = df["체중 : kg"]

result = sm.OLS(y, X).fit()

new_data = pd.DataFrame({"신장 : cm": [170], "측정나이": [30], "악력D : kg": [35]})
new_data = sm.add_constant(new_data, has_constant="add")

print(X.columns.tolist())
print(result.params.round(3))
print(round(result.predict(new_data)[0], 1))

실행 결과:

['const', '신장 : cm', '측정나이', '악력D : kg']
const    -72.043
신장 : cm      0.720
측정나이         0.097
악력D : kg     0.394
dtype: float64
67.0

첫 번째 줄의 const가 절편용 열입니다. smf.ols()에서는 자동이었지만, sm.OLS()에서는 직접 넣어야 한다는 점이 가장 큰 차이입니다.

계수는 각각 아래처럼 읽으면 됩니다.

• 신장 : cm: 키가 1 커질 때 예측 몸무게가 얼마나 변하는지
• 측정나이: 나이가 1살 늘 때 예측 몸무게가 얼마나 변하는지
• 악력D : kg: 악력이 1 커질 때 예측 몸무게가 얼마나 변하는지

여기서는 여러 설명변수가 함께 들어간 다중회귀 예시라, 계수 하나를 읽을 때도 “다른 변수들을 함께 고려한 상태에서의 변화량”으로 읽는 것이 맞습니다.

statsmodels를 배우는 핵심은 결과 객체를 읽는 것입니다. 초보자가 가장 먼저 익혀야 할 속성은 아래 여섯 가지입니다.

• params: 계수
• fittedvalues: 모델이 맞춘 예측값
• resid: 실제값 - 예측값인 잔차
• rsquared: 설명력
• pvalues: 각 계수가 통계적으로 유의한지 볼 때 참고하는 값
• conf_int(): 계수의 신뢰구간

import pandas as pd
import statsmodels.api as sm

data_path = "/data/body.csv"

df = pd.read_csv(data_path)

X = sm.add_constant(df[["신장 : cm", "측정나이", "악력D : kg"]])
y = df["체중 : kg"]

result = sm.OLS(y, X).fit()

print(round(result.rsquared, 3))
print(round(result.resid.iloc[0], 1))
print((result.pvalues < 0.05).to_dict())
print(result.conf_int().round(3))

실행 결과:

0.603
-5.0
{'const': True, '신장 : cm': True, '측정나이': True, '악력D : kg': True}
                0       1
const    -75.579 -68.506
신장 : cm    0.697   0.743
측정나이       0.087   0.106
악력D : kg   0.376   0.412

잔차 -5.0은 첫 번째 행에서 실제 몸무게가 모델 예측보다 약 5kg 작았다는 뜻입니다. pvalues < 0.05가 True라는 것은 이 예시에서는 각 계수가 통계적으로 유의하다고 볼 수 있다는 뜻입니다.

conf_int()는 계수의 추정 범위를 보여 줍니다. 예를 들어 신장 : cm의 신뢰구간이 0.697 ~ 0.743이라면, 이 구간 안쪽에서 계수가 형성되었을 가능성이 높다고 읽을 수 있습니다.

summary()도 아주 중요하지만, 표가 길기 때문에 처음에는 위 속성들을 직접 확인하는 방식이 더 이해하기 쉽습니다.

smf.logit()은 formula API 방식의 로지스틱 회귀입니다. 결과를 보면 각 계수와 함께, 새 데이터가 1일 확률을 바로 확인할 수 있습니다.

작은 예시 데이터부터 보겠습니다.

import pandas as pd
import statsmodels.formula.api as smf

df = pd.DataFrame(
    {
        "study_hours": [1, 2, 3, 4, 5, 6, 7, 8],
        "passed": [0, 0, 0, 1, 0, 1, 1, 1],
    }
)

result = smf.logit("passed ~ study_hours", data=df).fit(disp=False)

print(result.params.round(3))
print(round(result.predict(pd.DataFrame({"study_hours": [6]}))[0], 3))

실행 결과:

Intercept     -5.770
study_hours    1.282
dtype: float64
0.873

study_hours 계수가 양수이므로, 공부 시간이 늘수록 합격 확률이 올라가는 방향으로 모델이 학습되었다고 볼 수 있습니다. study_hours=6일 때 예측 확률이 0.873이라는 것은, 합격할 확률이 약 87.3%라는 뜻입니다.

외부 데이터로 보면 아래와 같습니다.

import pandas as pd
import statsmodels.formula.api as smf

data_path = "/data/body.csv"

df = pd.read_csv(data_path)

df["is_male"] = (df["측정회원성별"] == "M").astype(int)

result = smf.logit('is_male ~ Q("신장 : cm") + Q("악력D : kg")', data=df).fit(disp=False)

new_data = pd.DataFrame({"신장 : cm": [175], "악력D : kg": [40]})

print(round(result.params['Q("신장 : cm")'], 3))
print(round(result.predict(new_data)[0], 3))

실행 결과:

0.218
0.992

여기서 0.992는 새 데이터가 is_male=1, 즉 남성일 확률이 약 99.2%라는 뜻입니다. statsmodels의 로지스틱 회귀는 이런 식으로 클래스 확률을 해석하는 데 매우 좋습니다.

로지스틱 회귀의 계수는 바로 확률이 아닙니다. 계수에 np.exp()를 적용하면 오즈비(odds ratio)로 바꿔서 볼 수 있습니다. 오즈비는 “설명변수가 1 늘 때 odds가 몇 배가 되는가”를 보여 줍니다.

먼저 확률과 오즈의 관계를 아주 짧게 보면 아래와 같습니다.

prob = 0.8
odds = prob / (1 - prob)

print(round(odds, 1))

실행 결과:

4.0

확률이 0.8이면 오즈는 4.0입니다. 즉 일어날 쪽이 일어나지 않을 쪽보다 4배라는 뜻입니다.

외부 데이터에서는 이렇게 계수를 오즈비로 바꿉니다.

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

data_path = "/data/body.csv"

df = pd.read_csv(data_path)

df["is_male"] = (df["측정회원성별"] == "M").astype(int)

result = smf.logit('is_male ~ Q("신장 : cm") + Q("악력D : kg")', data=df).fit(disp=False)

print(round(np.exp(result.params['Q("악력D : kg")']), 3))

실행 결과:

1.463

악력D : kg의 오즈비가 1.463이라는 것은, 악력이 1 증가할 때 is_male=1 쪽의 오즈가 약 1.463배가 되는 방향이라는 뜻입니다. 오즈비가 1보다 크면 양의 방향, 1보다 작으면 음의 방향으로 이해하면 됩니다.

포켓몬 데이터에서도 같은 방식으로 볼 수 있습니다.

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

data_path = "/data/pokemon.csv"

df = pd.read_csv(data_path)
df["is_legendary"] = df["Legendary"].astype(int)

result = smf.logit("is_legendary ~ Total + Speed", data=df).fit(disp=False)

print(round(np.exp(result.params["Total"]), 3))

실행 결과:

1.031

이 값은 Total이 1 증가할 때 전설 포켓몬일 오즈가 약 1.031배가 되는 방향으로 읽을 수 있다는 뜻입니다.

sm.Logit()도 sm.OLS()처럼 class API입니다. 그래서 입력값 X와 정답 y를 따로 만들고, sm.add_constant()로 상수항을 추가해야 합니다.

import pandas as pd
import statsmodels.api as sm

data_path = "/data/pokemon.csv"

df = pd.read_csv(data_path)
df["is_legendary"] = df["Legendary"].astype(int)

X = sm.add_constant(df[["Total", "Speed"]])
y = df["is_legendary"]

result = sm.Logit(y, X).fit(disp=False)

new_data = pd.DataFrame({"Total": [600], "Speed": [110]})
new_data = sm.add_constant(new_data, has_constant="add")

print(round(result.params["Total"], 3))
print(round(result.predict(new_data)[0], 3))
print(round(result.pvalues["Speed"], 3))

실행 결과:

0.03
0.388
0.064

Total 계수가 양수라서 총합 능력치가 클수록 전설일 확률이 커지는 방향으로 학습되었다고 볼 수 있습니다. 반면 Speed의 pvalue가 0.064라서, 이 예시에서는 0.05 기준으로 아주 강하게 유의하다고 보기는 어렵습니다.

이런 식으로 statsmodels는 단순 예측뿐 아니라 “어떤 변수가 더 의미 있어 보이는가”를 읽는 데 강합니다.

로지스틱 회귀는 기본적으로 확률을 예측합니다. 하지만 분류 문제에서는 마지막에 그 확률을 0과 1, 혹은 클래스 이름으로 나누어야 할 때가 많습니다. 가장 흔한 기준은 0.5입니다.

import numpy as np

probs = np.array([0.2, 0.51, 0.8, 0.49])
labels = (probs >= 0.5).astype(int)

print(labels.tolist())

실행 결과:

[0, 1, 1, 0]

0.5 이상이면 1, 작으면 0으로 바꾼 결과입니다. 이 기준은 가장 흔하지만, 항상 정답은 아닙니다. 문제에 따라 0.3이나 0.7 같은 다른 기준을 쓰기도 합니다.

외부 데이터 예측 확률을 같은 방식으로 나누면 아래처럼 읽을 수 있습니다.

body_prob = 0.992
pokemon_prob = 0.388

print(int(body_prob >= 0.5))
print(int(pokemon_prob >= 0.5))

실행 결과:

1
0

첫 번째 값은 남성으로 분류되고, 두 번째 값은 전설 포켓몬이 아닌 쪽으로 분류됩니다.