R for Data Science by Wickham & Grolemund
# numerical calculation & data frames
import numpy as np
import pandas as pd
# visualization
import matplotlib.pyplot as plt
import seaborn as sns
import seaborn.objects as so
from sbcustom import *
# statistics
import statsmodels.api as sm# pandas options
pd.set_option("mode.copy_on_write", True)
pd.options.display.precision = 2
pd.options.display.float_format = '{:.2f}'.format # pd.reset_option('display.float_format')
pd.options.display.max_rows = 8
# Numpy options
np.set_printoptions(precision = 2, suppress=True)Mean function: \(E(Y | X = x_i) = \beta_0 + \beta_1 x_i\)
Variance function: \(Var(Y | X = x_i) = \sigma^2\)
Distribution: \((Y | X = x_i) \sim N(\beta_0 + \beta_1 x_i, \sigma^2)\)
Source: p.107, Applied Regression Analysis and Generalized Linear Models (3e), by John Fox
Source: p.482, Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences (3e), by Cohen, J., Cohen, P., West, S. G., & Aiken, L. S.
Response 변수가 binary인 경우,
예제: Blowdown
p.274, Applied Linear Regression (4e) by Sanford Weisberg
Data: blowdown.csv
On July 4, 1999, a storm with winds exceeding 90 miles per hour hit the Boundary Waters Canoe Area Wilderness in northeastern Minnesota, causing serious damage to the forest. Rich et al. (2007) studied the effects of this storm using a very extensive ground survey of the area, determining status, either alive or dead, of more than 3600 trees.
blowdown = pd.read_csv('data/blowdown.csv')
blowdown d s y spp
0 9.00 0.02 0 balsam fir
1 14.00 0.02 0 balsam fir
2 18.00 0.02 0 balsam fir
3 23.00 0.02 0 balsam fir
... ... ... .. ...
3662 20.00 0.98 1 jackpine
3663 19.00 0.98 1 jackpine
3664 37.00 0.98 1 black ash
3665 48.00 0.98 1 black ash
[3666 rows x 4 columns]# black spruce에 제한하고, diameter d를 log2로 변환
blowbs = blowdown.query('spp == "black spruce"')
blowbs["d"] = np.floor(blowbs["d"])
blowbs["log2d"] = np.log2(blowbs["d"])blowbs d s y spp log2d
17 9.00 0.02 0 black spruce 3.17
24 11.00 0.03 0 black spruce 3.46
25 9.00 0.03 0 black spruce 3.17
50 9.00 0.03 0 black spruce 3.17
... ... ... .. ... ...
3636 5.00 0.94 0 black spruce 2.32
3646 9.00 0.94 1 black spruce 3.17
3647 17.00 0.94 1 black spruce 4.09
3661 8.00 0.98 1 black spruce 3.00
[659 rows x 5 columns]기존의 ordinary least squares (OLS) 방법으로 fit을 하면, (linear probability model)
(
so.Plot(blowbs, x='d', y='y')
.add(so.Dots(alpha=.3), so.Jitter(y=.1))
.add(so.Line(), so.PolyFit(5))
.add(so.Line(), so.PolyFit(1))
.limit(y=(-0.2, 1.2))
.layout(size=(8, 5))
)
from statsmodels.formula.api import ols
mod_ols = ols('y ~ d', data = blowbs).fit()
mod_ols.paramsIntercept -0.16
d 0.05
dtype: float64blowbs_plot = blowbs.assign(resid = mod_ols.resid)
(
so.Plot(blowbs_plot, x='d', y='resid')
.add(so.Dots(), so.Jitter(y=.1))
.label(title='Residual Plot')
)
Distributions for d
(
so.Plot(blowbs_plot.query('d < 15'), x='resid')
.add(so.Bar(), so.Hist(stat="proportion"))
.facet("d", wrap=5)
.layout(size=(8, 4))
.label(title='d = {}'.format)
)
이제 y를 직접 예측하기보다, 확률을 예측하는 방식을 취하면,
m)와died)를 고려하면,p= died/m)을 계산할 수 있음. 이를 emperical probability라고 함.p는 binary response (0, 1)의 conditional mean(평균)인데,blowbs_bn = blowbs.groupby("d")["y"].agg([("died", "sum"), ("m", "count"), ("p", "mean")]).reset_index()
# blowbs_bn = blowbs_bn.assign(
# p = (blowbs_bn.died + 0.1)/ (blowbs_bn.m + 0.2), # 확률이 0이거나 1이 되는 것을 피하기 위해 작은 값을 추가
# )
blowbs_bn d died m p
0 5.00 7 90 0.08
1 6.00 7 92 0.08
2 7.00 18 91 0.20
3 8.00 13 70 0.19
.. ... ... .. ...
21 27.00 1 1 1.00
22 28.00 2 2 1.00
23 29.00 1 2 0.50
24 32.00 1 1 1.00
[25 rows x 4 columns](
so.Plot(blowbs_bn, x='d', y='p')
.add(so.Dots(color="orangered"), pointsize='m')
.add(so.Dots(alpha=.3), so.Jitter(y=.1), x=blowbs.d, y=blowbs.y)
.add(so.Line(), so.PolyFit(5))
.add(so.Line(), so.PolyFit(1))
.limit(y=(-0.2, 1.2))
.layout(size=(8, 5))
.label(title='Emperical Probability')
)
우선 d를 log2 변환해서 살펴보면,
blowbs_bn["log2d"] = np.log2(blowbs_bn["d"])
(
so.Plot(blowbs_bn, x='log2d', y='p')
.add(so.Dots(color="orangered"), pointsize='m')
.add(so.Dots(alpha=.3), so.Jitter(y=.1), x=blowbs.log2d, y=blowbs.y)
.add(so.Line(), so.PolyFit(5))
.add(so.Line(), so.PolyFit(1))
.limit(y=(-0.2, 1.2))
.layout(size=(8, 5))
.label(title='Emperical Probability')
)
만약, 이 proportion p (= died/m)를 log2d로 예측하는 모형을 만들어, 그 잔차를 살펴보면,
mod_prob = ols('p ~ log2d', data = blowbs_bn).fit()
mod_prob_plot = blowbs_bn.assign(resid = mod_prob.resid)
(
so.Plot(mod_prob_plot, x='log2d', y='resid')
.add(so.Dot())
.add(so.Line(color=".3", linestyle=":"), so.Agg(lambda x: 0))
.label(title='Residual Plot')
)
위에서 살펴본 OLS의 문제들 즉,
이런 문제들을 해결하고 예측값이 분명한 “확률”의 의미를 품도록 여러 방식이 제시되는데 주로 사용되는 것이 logistic regression임.
Binary outcome을 예측하는 모형은 binary 값을 예측하는 것이 아니고, 확률 값을 예측하는 것임.
이후에 이를 이용해 binary outcome을 예측.
Odds의 정의: 실패할 확률 대비 성공할 확률의 비율
\(\displaystyle odds = \frac{p}{1-p}\)
예를 들어, 5cm 두께의 나무는 90그루 중 7그루가 죽었으므로 83그루는 살았음.
즉, 죽음:생존 = 7:83 \(\approx\) 1:12 이고 odds = 7/83 = 0.084; 생존할 가능성 대비 죽을 가능성이 8.4%임.
확률로 표현하면, \(odds = \frac{\frac{7}{90}}{1 - \frac{7}{90}} = \frac{7}{90-7} = \frac{7}{83}\)
확률과 odds, logit(log odds)의 관계
blowbs_bn = blowbs_bn.assign(odds = lambda x: x.p / (1 - x.p))
# p = 1인 경우 odds가 무한대가 되므로 편의상 inf 값을 50으로 대체
blowbs_bn["odds"] = blowbs_bn["odds"].apply(lambda x: 50 if x == np.inf else x)
blowbs_bn d died m p log2d odds
0 5.00 7 90 0.08 2.32 0.08
1 6.00 7 92 0.08 2.58 0.08
2 7.00 18 91 0.20 2.81 0.25
3 8.00 13 70 0.19 3.00 0.23
.. ... ... .. ... ... ...
21 27.00 1 1 1.00 4.75 50.00
22 28.00 2 2 1.00 4.81 50.00
23 29.00 1 2 0.50 4.86 1.00
24 32.00 1 1 1.00 5.00 50.00
[25 rows x 6 columns]이 odds를 선형모형으로 나무두께로 예측하려고 하는 것인데, 예를 들어,
\(\widehat{odds} = b_0 + b_1 \cdot log_{2}(d)\)
odds 값의 범위는 (0, \(\infty\))이므로, 선형모형으로 fit하는 것이 적절하지 않음.
한편, odds를 log 변환하면, 그 범위는 (\(-\infty\), \(\infty\))가 되어 선형모형으로 fit하는 것이 적절해짐.
\(\displaystyle logit(p):= log(odds) = log\left(\frac{p}{1-p}\right)\)
이 때, 예측모형은 \(\displaystyle log\left(\frac{\hat{p}}{1-\hat{p}}\right) = b_0 + b_1 \cdot log_{2}(d)\)
즉, logit은 예측변수 \(x\) 와 선형적으로 연결되는데 반해, 확률 \(p\) 는 예측변수 \(x\) 와 비선형적 관계를 맺음
log odds값을 구해보면,
blowbs_bn = blowbs_bn.assign(log_odds = lambda x: np.log(x.odds))
blowbs_bn d died m p log2d odds log_odds
0 5.00 7 90 0.08 2.32 0.08 -2.47
1 6.00 7 92 0.08 2.58 0.08 -2.50
2 7.00 18 91 0.20 2.81 0.25 -1.40
3 8.00 13 70 0.19 3.00 0.23 -1.48
.. ... ... .. ... ... ... ...
21 27.00 1 1 1.00 4.75 50.00 3.91
22 28.00 2 2 1.00 4.81 50.00 3.91
23 29.00 1 2 0.50 4.86 1.00 0.00
24 32.00 1 1 1.00 5.00 50.00 3.91
[25 rows x 7 columns]정리하면, 아래와 같은 관측값들과 emperical probability/log odds을 함께 살펴보면,
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
sns.scatterplot(x=blowbs_bn.log2d, y=blowbs_bn.p, size=blowbs_bn.m, sizes=(20, 200), ax=ax)
def jitter(values, j):
return values + np.random.normal(0, j, values.shape)
sns.scatterplot(x=blowbs.log2d, y=jitter(blowbs.y, 0.02), alpha=.3, ax=ax)
for i, row in blowbs_bn.iterrows():
ax.annotate(f"{row.died:n}/{row.m:n}", xy=(row.log2d, row.p), xytext=(row.log2d, row.p-0.05), size=9)
ax.set_xticks(blowbs_bn.log2d.unique())
ax.set_xticklabels(blowbs_bn.log2d.unique().round(2))
ax.tick_params(axis='x', rotation=45)
ax.set_title("Emperical Probability")
ax.legend_.remove()
Logit 값을 추가해서 그리면,
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
sns.scatterplot(x=blowbs_bn.log2d, y=blowbs_bn.p, size=blowbs_bn.m, sizes=(20, 200), ax=ax, legend=False)
sns.scatterplot(x=blowbs_bn.log2d, y=blowbs_bn.log_odds, size=blowbs_bn.m, sizes=(20, 200), color="red", ax=ax)
def jitter(values, j):
return values + np.random.normal(0, j, values.shape)
sns.scatterplot(x=blowbs.log2d, y=jitter(blowbs.y, 0.02), alpha=.3, ax=ax)
for i, row in blowbs_bn.iterrows():
ax.annotate(f"logit{row.died:n}/{row.m:n}", xy=(row.log2d, row.log_odds), xytext=(row.log2d, row.log_odds-0.4), size=9)
ax.set_xticks(blowbs_bn.log2d.unique())
ax.set_xticklabels(blowbs_bn.log2d.unique().round(2))
ax.tick_params(axis='x', rotation=45)
ax.set_ylabel("P & Logit")
ax.set_title("Emperical Probability and Logit")
ax.legend_.remove()
# # polyfit 5
# x = np.linspace(blowbs_bn.log2d.min(), blowbs_bn.log2d.max(), 100)
# y = np.polyval(np.polyfit(blowbs_bn.log2d, blowbs_bn.log_odds, 5), x)
# sns.lineplot(x=x, y=y, ax=ax, color=".6")
# # polyfit 1
# y = np.polyval(np.polyfit(blowbs_bn.log2d, blowbs_bn.log_odds, 1), x)
# sns.lineplot(x=x, y=y, ax=ax, color=".6")
위의 logit값을 선형모형으로 예측하는 모형: \(\displaystyle log\left(\frac{\hat{p}}{1-\hat{p}}\right) = b_0 + b_1 \cdot log_{2}(d)\)
파라미터 \(b_0, b_1\)의 추정은 잔차들의 제곱의 합을 최소로 하는 OLS 방식은 부적절하며, 대신에 Maximum Likelihood Estimation을 사용함.
참고: 이항분포 (Y: 사건의 횟수)
blowbs_bn.head(3) d died m p log2d odds log_odds
0 5.00 7 90 0.08 2.32 0.08 -2.47
1 6.00 7 92 0.08 2.58 0.08 -2.50
2 7.00 18 91 0.20 2.81 0.25 -1.40
각 likelihood는 관측치들이 모두 독립적으로 발생했다고 가정했을 때의 확률값이고,
모든 데이터가 관찰될 likelihood:
Lik = \(_{90}C_7 \cdot p_1^7 \cdot (1-p_1)^{83} \cdot _{92}C_7 \cdot p_2^7 \cdot (1-p_2)^{85} \cdot _{91}C_{18} \cdot p_3^{18} \cdot (1-p_3)^{73} \cdots\)
이 때, 모형을 예측변수 x의 1차 다항함수로 fit한다면, \(\displaystyle log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 \cdot x\) 인데,
변형하면 \(\displaystyle p = \frac{1}{1+e^{-(\beta_0 + \beta_1 \cdot x)}}\)
Lik \(\propto\) \(\displaystyle \left(\frac{1}{1+e^{-(\beta_0 + \beta_1 \cdot 2.32)}}\right)^7 \cdot \left(1 - \frac{1}{1+e^{-(\beta_0 + \beta_1 \cdot 2.32)}}\right)^{83} \cdot \left(\frac{1}{1+e^{-(\beta_0 + \beta_1 \cdot 2.58)}}\right)^7 \cdot \left(\frac{1}{1+e^{-(\beta_0 + \beta_1 \cdot 2.58)}}\right)^{85} \cdots\)
이 Likelihood를 최대가 되도록 \(\beta_0, \beta_1\)의 추정치를 찾는 것이 Maximum Likelihood Estimation임.
from statsmodels.formula.api import glm, logit
# Binomial response
mod_bn = glm('died + I(m-died) ~ log2d', data=blowbs_bn, family=sm.families.Binomial()).fit()
# Binary response: 동일한 분석
mod = logit('y ~ log2d', data=blowbs).fit() # 아래와 동일하나 다른 methods를 제공
mod = glm('y ~ log2d', data=blowbs, family=sm.families.Binomial()).fit()
print(mod.summary())Optimization terminated successfully.
Current function value: 0.499165
Iterations 6
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 659
Model: GLM Df Residuals: 657
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -328.95
Date: Mon, 05 Jun 2023 Deviance: 657.90
Time: 10:32:39 Pearson chi2: 679.
No. Iterations: 5 Pseudo R-squ. (CS): 0.2599
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -7.8162 0.628 -12.437 0.000 -9.048 -6.584
log2d 2.2408 0.190 11.773 0.000 1.868 2.614
==============================================================================Model deviance:
\(D_k = -2[log(likelihood_k) - log(likelihood_{perfect})]\)
# deviance for an intercept only model
mod.null_deviance856.2073760911842위 fitted model의 예측값들: fitted line
fig, ax = plt.subplots(1, 1, figsize=(9, 5))
sns.scatterplot(x=blowbs_bn.log2d, y=blowbs_bn.p, size=blowbs_bn.m, sizes=(20, 200), ax=ax, legend=False)
def jitter(values, j):
return values + np.random.normal(0, j, values.shape)
sns.scatterplot(x=blowbs.log2d, y=jitter(blowbs.y, 0.02), alpha=.3, ax=ax)
for i, row in blowbs_bn.iterrows():
ax.annotate(f"{row.died:n}/{row.m:n}", xy=(row.log2d, row.p), xytext=(row.log2d, row.p-0.05), size=9)
ax.set_xticks(blowbs_bn.log2d.unique())
ax.set_xticklabels(blowbs_bn.log2d.unique().round(2))
# x-axis with 45 degree rotation
ax.tick_params(axis='x', rotation=45)
# fitted line
sns.lineplot(x=blowbs.log2d, y=mod.fittedvalues, ax=ax, color=".6")
plt.show()
Logistic regression에서는 세 가지 타입의 예측값들이 있음.
blowbs_bn_pred = blowbs_bn.assign(
pred_prob = mod_bn.fittedvalues,
pred_odds = lambda x: x.pred_prob / (1 - x.pred_prob),
pred_logit = lambda x: np.log(x.pred_odds)
)
blowbs_bn_pred.iloc[:, 1:] died m p log2d odds log_odds pred_prob pred_odds pred_logit
0 7 90 0.08 2.32 0.08 -2.47 0.07 0.07 -2.61
1 7 92 0.08 2.58 0.08 -2.50 0.12 0.13 -2.02
2 18 91 0.20 2.81 0.25 -1.40 0.18 0.22 -1.53
3 13 70 0.19 3.00 0.23 -1.48 0.25 0.33 -1.09
.. ... .. ... ... ... ... ... ... ...
21 1 1 1.00 4.75 50.00 3.91 0.94 17.09 2.84
22 2 2 1.00 4.81 50.00 3.91 0.95 19.22 2.96
23 1 2 0.50 4.86 1.00 0.00 0.96 21.53 3.07
24 1 1 1.00 5.00 50.00 3.91 0.97 29.59 3.39
[25 rows x 9 columns]예를 들어, 두께(log2d)가 3인 나무의 경우,
Odds의 비율 (odds ratio, OR)을 통해 해석
\(\displaystyle odds: \frac{\hat{p}}{1 - \hat{p}} = e^{b_0 + b_1 \cdot x}=e^{b_0}\cdot e^{b_1 \cdot x}\) 로부터
\(x\)가 1 증가할 때 odds의 비율: \(\displaystyle odds ~ratio: \frac{\frac{\hat{p_2}}{1 - \hat{p_2}}}{\frac{\hat{p_1}}{1 - \hat{p_1}}} = e^{b_1 \cdot (x+1) - b_1 \cdot x} = e^{b_1}\)
즉, \(x\)가 1 증가하면, odds가 “몇 배”로 증가하는지를 나타냄.
따라서, odds ratio가 1보다 크면 (\(b_1\)이 양수) \(x\)가 1 증가할 때, event의 odds가 커지며,
odds ratio가 1보다 작으면 (\(b_1\)이 음수) event의 odds가 줄어듬.
위의 경우, \(\displaystyle \widehat{odds} = e^{-7.82 + 2.24 \cdot log_2(d)}\) 이므로 odds ratio = \(e^{2.24} = 9.4\)
\(b_0\): d = 1일 때의 odds: \(\displaystyle e^{-7.82 + 2.34 \cdot 0} = e^{-7.82} = 0.004\), 즉 두께가 1cm 일 때 태풍에 나무가 쓰러지지 않을 가능성 대비 쓰러질 가능성은 0.004임.
이제 이 모형이 좋은 모형인지 살펴보기 위해 residual, 잔차를 살펴볼 수 있는가?
Binomial version
Binary version
Deviance를 이용해 OLS에서의 \(R^2\)와 비슷한 개념을 구성
“Coefficient of discrimination” (Tjur, 2009): average \(\hat{p}\) when \(y=1\) - average \(\hat{p}\) when \(y=0\)
Binomial response에 대해 residual들을 살펴보면,
glm('died + I(m-died) ~ log2d', data=blowbs_bn, family=sm.families.Binomial())
blowbs_bn_resid = blowbs_bn.assign(
prob_pred = mod_bn.fittedvalues,
died_pred = lambda x: x.prob_pred * x.m,
count_resid = lambda x: x.died - x.died_pred,
pearson_resid = mod_bn.resid_pearson,
deviance_resid = mod_bn.resid_deviance,
)
blowbs_bn_resid d died m p log2d odds log_odds prob_pred died_pred
0 5.00 7 90 0.08 2.32 0.08 -2.47 0.07 6.15 \
1 6.00 7 92 0.08 2.58 0.08 -2.50 0.12 10.74
2 7.00 18 91 0.20 2.81 0.25 -1.40 0.18 16.26
3 8.00 13 70 0.19 3.00 0.23 -1.48 0.25 17.56
.. ... ... .. ... ... ... ... ... ...
21 27.00 1 1 1.00 4.75 50.00 3.91 0.94 0.94
22 28.00 2 2 1.00 4.81 50.00 3.91 0.95 1.90
23 29.00 1 2 0.50 4.86 1.00 0.00 0.96 1.91
24 32.00 1 1 1.00 5.00 50.00 3.91 0.97 0.97
count_resid pearson_resid deviance_resid
0 0.85 0.36 0.35
1 -3.74 -1.21 -1.29
2 1.74 0.48 0.47
3 -4.56 -1.26 -1.30
.. ... ... ...
21 0.06 0.24 0.34
22 0.10 0.32 0.45
23 -0.91 -3.13 -1.88
24 0.03 0.18 0.26
[25 rows x 12 columns](
so.Plot(blowbs_bn_resid, x='log2d', y='pearson_resid')
.add(so.Dot())
.add(so.Dot(color="orangered", alpha=.5), y = 'deviance_resid')
.label(y="Pearson / Deviance Residuals")
)
(
so.Plot(blowbs_bn_resid, x='prob_pred', y='p',)
.add(so.Dot())
.add(so.Line(color=".6"), y = 'prob_pred') # y = x line
.label(x="Predicted Probability", y="Emperical Probability")
)
Binary response에 대해 residual들을 살펴보면,
glm('y ~ log2d', data=blowbs, family=sm.families.Binomial())
blowbs_resid = blowbs.assign(
prob_pred = mod.fittedvalues,
prob_resid = mod.resid_response,
pearson_resid = mod.resid_pearson,
deviance_resid = mod.resid_deviance
)
blowbs_resid.sort_values("d") d s y spp log2d prob_pred prob_resid pearson_resid
2102 5.00 0.45 0 black spruce 2.32 0.07 -0.07 -0.27 \
724 5.00 0.18 0 black spruce 2.32 0.07 -0.07 -0.27
723 5.00 0.18 1 black spruce 2.32 0.07 0.93 3.69
1687 5.00 0.37 0 black spruce 2.32 0.07 -0.07 -0.27
... ... ... .. ... ... ... ... ...
2634 28.00 0.57 1 black spruce 4.81 0.95 0.05 0.23
1784 29.00 0.38 1 black spruce 4.86 0.96 0.04 0.22
1079 29.00 0.25 0 black spruce 4.86 0.96 -0.96 -4.64
3455 32.00 0.80 1 black spruce 5.00 0.97 0.03 0.18
deviance_resid
2102 -0.38
724 -0.38
723 2.32
1687 -0.38
... ...
2634 0.32
1784 0.30
1079 -2.50
3455 0.26
[659 rows x 9 columns](
so.Plot(blowbs_resid, x='log2d')
.pair(y=["prob_resid", "pearson_resid", "deviance_resid"])
.add(so.Dots(alpha=.3), so.Jitter(y=.5))
.add(so.Line(), so.PolyFit(5))
.layout(size=(5, 8))
)
Emperical probability vs. predicted probability
blowbs_resid["y2"] = blowbs_resid.y.map({0: "survived(y=0)", 1: "died(y=1)"})
(
so.Plot(blowbs_resid, x='prob_pred')
.add(so.Bars(), so.Hist())
.facet("y2")
.label(x="Predicted Probability", y="Count")
)
blowbs_resid["prob_cat"] = (
pd.cut(
blowbs_resid.prob_pred,
bins=np.arange(0, 11) * 0.1,
include_lowest=True,
labels=False, # labels=False: 0, 1, 2, ... 9,
)
* 0.1 + 0.05 # 각 percentile의 중앙값을 만듦
)
blowbs_resid[["prob_pred", "prob_cat"]]
# prob_pred prob_cat
# 17 0.33 0.35
# 24 0.48 0.45
# 25 0.33 0.35
# 50 0.33 0.35
# ... ... ...
# 3636 0.07 0.05
# 3646 0.33 0.35
# 3647 0.79 0.75
# 3661 0.25 0.25
# 각 percentile에서 실제 event(y=1)의 관찰 비율을 계산
prob_plot = blowbs_resid.groupby("prob_cat")["y"].mean().reset_index()
(
so.Plot(prob_plot, x='prob_cat', y='y')
.add(so.Dot())
.add(so.Line(color=".6"), y = 'prob_cat') # y = x line
.label(x="Predicted Probability", y="Observed Probability")
)Calibration plot
예측된 확률을 기반으로 binary outcome으로 예측하여, 모형의 예측력을 평가
| Predicted | ||
|---|---|---|
| Actual | survied | died |
| survied | 377 | 49 |
| died | 107 | 126 |
| Predicted | ||
|---|---|---|
| Actual | survied | died |
| survied | True Negative | False Positive |
| died | False Negative | True Positive |
mod_lg = logit('y ~ log2d', data=blowbs).fit() # logit has .pred_table() methodOptimization terminated successfully.
Current function value: 0.499165
Iterations 6# accuracy rate with threshold 0.3
cm = mod_lg.pred_table(0.3) # confusion matrix
display(cm)
(cm[0, 0] + cm[1, 1]) / cm.sum()array([[298., 128.],
[ 45., 188.]])0.7374810318664643# accuracy rate with threshold 1: no information
cm = mod_lg.pred_table(1) # confusion matrix
display(cm)
(cm[0, 0] + cm[1, 1]) / cm.sum()array([[426., 0.],
[233., 0.]])0.6464339908952959옳은 예측의 비율
틀린 예측의 비율
Receiver operating characteristic (ROC) curve
예측된 확률에 대한 임계치를 조정함에 따라 옳은 예측과 틀린 예측의 비율이 어떻게 달라지는지 살펴봄으로써 임계치를 설정하는데 도움을 줌
from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(blowbs_resid.y, blowbs_resid.prob_pred)
roc = pd.DataFrame(
{
"thresholds": thresholds.round(2),
"False Pos": fpr,
"sensitivity(True Pos)": tpr,
"specificity(True Neg)": 1 - fpr,
}
)
roc.sort_values("thresholds").head(10) thresholds False Pos sensitivity(True Pos) specificity(True Neg)
25 0.07 1.00 1.00 0.00
24 0.12 0.81 0.97 0.19
23 0.18 0.61 0.94 0.39
22 0.25 0.43 0.86 0.57
.. ... ... ... ...
19 0.48 0.15 0.64 0.85
18 0.55 0.12 0.54 0.88
17 0.62 0.09 0.48 0.91
16 0.67 0.08 0.38 0.92
[10 rows x 4 columns]
AUC: Area Under the Curve = concordance index
from sklearn.metrics import roc_auc_score, auc
roc_auc = roc_auc_score(blowbs_resid.y, blowbs_resid.prob_pred)
roc_auc # auc(fpr, tpr)0.8161306897177054Classifier로서 전반적인 모형의 성능 vs. 특정 임계치에서의 모형의 성능 vs. 확률모형
Binary response에 대해 그들의 conditional 평균인 확률 \(p(x_i)\)를 추정
Distribution: \((Y_i | X = x_i) \sim Bin(m_i, p(x_i))\) 의 분포가 binomial distribution을 따름. 따라서 \(Y_i\)의 평균과 분산은 \(m_i\)와 \(p(x_i)\)에 따라 결정되어 변함.
Linear predictors: 확률 \(p(x_i)\)는 link function (위에서는 logit)를 통해 예측변수들의 선형 결합으로 결정됨.
\(\displaystyle log\left(\frac{\hat{p}}{1-\hat{p}}\right) = b_0 + b_1 \cdot X_1 + b_2 \cdot X_2\)
Probit link function으로 쓸 수도 있음: \(\displaystyle \Phi^{-1}(\hat{p}) = b_0 + b_1 \cdot X_1 + b_2 \cdot X_2\)
위의 아이디어는 generalized linear model (GLM)로 확장되는데,
Source: p.274, Applied Linear Regression (4e) by Sanford Weisberg
위에서 다룬 blowdown 데이터셋에서 black spruce라는 종에 대해서만 살펴봤는데, 이번에는 black spruce와 aspen 두 종의 나무를 모두 포함한 데이터로 비슷한 모형을 세워 분석해보세요.
spp를 추가하고,d)에 따라 태풍에 쓰러질(y) 확률이 두 종에서 큰 차이를 보이는지 살펴보고,d) 쓰러질 odds ratio(OR)가 얼마나 증가하는지 파라미터를 통해 구해보고,blow2 = blowdown.query('spp in ["black spruce", "aspen"]')
blow2["d"] = np.floor(blow2["d"])
blow2["log2d"] = np.log2(blow2["d"])
blow2 d s y spp log2d
17 9.00 0.02 0 black spruce 3.17
24 11.00 0.03 0 black spruce 3.46
25 9.00 0.03 0 black spruce 3.17
50 9.00 0.03 0 black spruce 3.17
... ... ... .. ... ...
3640 29.00 0.94 1 aspen 4.86
3646 9.00 0.94 1 black spruce 3.17
3647 17.00 0.94 1 black spruce 4.09
3661 8.00 0.98 1 black spruce 3.00
[1095 rows x 5 columns]blow2_bn = blow2.groupby(["log2d", "spp"])["y"].agg([("died", "sum"), ("m", "count"), ("p", "mean")]).reset_index()
blow2_bn log2d spp died m p
0 2.32 aspen 0 2 0.00
1 2.32 black spruce 7 90 0.08
2 2.58 aspen 1 2 0.50
3 2.58 black spruce 7 92 0.08
.. ... ... ... .. ...
66 5.61 aspen 1 1 1.00
67 5.67 aspen 0 1 0.00
68 5.75 aspen 1 1 1.00
69 5.81 aspen 0 1 0.00
[70 rows x 5 columns]Source: p.172, Beyond Multiple Linear Regression, by Paul Roback, Julie Legler.
다음에서는 체중감량을 하고자하는 학생들의 의도에 영향을 줄 수 있는 3가지 요소들을 탐색합니다.
sex)bmipct)media)A sample of 500 teens from data collected in 2009 through the U.S. Youth Risk Behavior Surveillance System (YRBSS) [Centers for Disease Control and Prevention, 2009]. The YRBSS is an annual national school-based survey conducted by the Centers for Disease Control and Prevention (CDC) and state, territorial, and local education and health agencies and tribal governments.
Data: diet.csv
lose_wt_01
Q66. Which of the following are you trying to do about your weight?
A. Lose weight (1)
B. Gain weight (0)
C. Stay the same weight (0)
D. I am not trying to do anything about my weight (0)
media
Q81. On an average school day, how many hours do you watch TV?
A. I do not watch TV on an average school day (0)
B. Less than 1 hour per day (0.5)
C. 1 hour per day (1)
D. 2 hours per day (2)
E. 3 hours per day (3)
F. 4 hours per day (4)
G. 5 or more hours per day (5)
bmipct
The percentile for a given BMI for members of the same sex.
diet = pd.read_csv("data/diet.csv")
diet lose_wt lose_wt_01 sex media bmipct
0 Lose weight 1 Male 3.00 98
1 No weight loss 0 Female 1.00 41
2 No weight loss 0 Male 3.00 6
3 No weight loss 0 Male 3.00 41
.. ... ... ... ... ...
441 No weight loss 0 Female 0.50 43
442 No weight loss 0 Male 3.00 40
443 Lose weight 1 Female 1.00 39
444 No weight loss 0 Male 2.00 34
[445 rows x 5 columns]lose_wt_01을 response로 하는 logistic regression 모형을 세워 분석하고,
bmipct를 10개의 percentile 구간으로 나누어 (즉, 10%씩 10개 구간) emperical logit값을 구해 그려볼 것pd.qcut을 이용
media)을 discrete하게 보고 각 시간 구간에서 emperical logit을 구해 살펴볼 것
다음 모형에 대해
lose_wt_01 ~ sex + bmipct + media