# numerical calculation & data framesimport numpy as npimport pandas as pd# visualizationimport matplotlib.pyplot as pltimport seaborn as snsimport seaborn.objects as so# statisticsimport statsmodels.api as sm# pandas optionspd.set_option('mode.copy_on_write', True) # pandas 2.0pd.options.display.float_format ='{:.2f}'.format# pd.reset_option('display.float_format')pd.options.display.max_rows =7# max number of rows to display# NumPy optionsnp.set_printoptions(precision =2, suppress=True) # suppress scientific notation# For high resolution displayimport matplotlib_inlinematplotlib_inline.backend_inline.set_matplotlib_formats("retina")
Generalized Linear Models
Ordinary least squares (OLS)
Mean function: \(E(Y | X = x_i) = \beta_0 + \beta_1 x_i\)
1999년 7월 4일, 미네소타 북동부의 바운더리 워터스 카누 지역 야생지대에 시속 90마일이 넘는 폭풍이 몰아쳐 숲에 심각한 피해를 입혔습니다. Rich 등(2007)은 이 지역에 대한 매우 광범위한 지상 조사를 통해 이 폭풍의 영향을 연구하여 3600여 그루의 나무의 생사 여부를 확인했습니다.
d: 나무의 지름(cm), y: 나무가 죽은 경우 1, 그렇지 않으면 0 spp: 나무의 종, s: 특정 지역의 폭풍의 심각도
blowbs_bn["log2d"] = np.log2(blowbs_bn["d"])( so.Plot(blowbs_bn, x='log2d', y='p') .add(so.Dot(), pointsize='m') .add(so.Dots(alpha=.3, color=".6"), so.Jitter(y=.1), x=blowbs.log2d, y=blowbs.y) .add(so.Line(color="#87bc45"), so.PolyFit(1)) .limit(y=(-0.2, 1.2)) .layout(size=(8, 5)) .label(title='Emperical Probability', y ='Proportion of Died', x='Log of Diameter (cm)'))
위에서 살펴본 OLS의 문제들 즉,
잔차에 패턴이 보인다는 것은 충분히 좋은 모형이 아니라는 것을 의미하고,
예측값이 확률을 의미하지 못할 수 있음.
잔차의 분산이 x값에 따라 달라져 모집단에 대한 추론을 어렵게 함.
이런 문제들을 해결하고 예측값이 분명한 “확률”의 의미를 품도록 여러 방식이 제시되는데 주로 사용되는 것이 logistic regression임.
Important
Binary outcome을 예측하는 모형은 binary 값을 예측하는 것이 아니고, 확률 값을 예측하는 것임.
이후에 이를 이용해 binary outcome을 예측.
예를 들어, 두께가 5cm인 (특정 종의) 나무가 태풍에 쓰러질 확률(true probability)을 파악하고자 함.
이 때, 관측값은 5cm인 나무 중 쓰러진 나무의 “비율”이고, 이 관측치들로부터 true probability를 추정하고자 함.
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: 50if x == np.inf else x)blowbs_bn
파라미터 \(b_0, b_1\)의 추정은 잔차들의 제곱의 합을 최소로 하는 OLS 방식은 부적절하며, 대신에 Maximum Likelihood Estimation을 사용함.
아이디어는 관측치가 전체적으로 관찰될 likelihood가 최대가 되도록 \(b_0, b_1\)을 선택하는 것임
이를 위해서 적절한 확률모형을 결합시켜야 함.
선택하는 확률 모형은 Bernoulli 분포임; 평균 \(E(Y|X) = p\)이고, \(logit(E(Y|X)) = b_0 + b_1 \cdot log_{2}(d)\)
발생 비율값으로 변환해 Binomial distribution (이항분포)로 전개하는 방식도 있음; binomial logit model
관찰값은 Bernoulli 분포로부터 발생했다고 가정함으로써, 실제 관찰값들이 관찰될 확률/가능도를 기준으로 파라미터를 추정할 수 있음.
blowbs.sort_values("log2d")
d s y spp log2d label
2102 5.00 0.45 0 black spruce 2.32 alive
724 5.00 0.18 0 black spruce 2.32 alive
723 5.00 0.18 1 black spruce 2.32 died
... ... ... .. ... ... ...
1784 29.00 0.38 1 black spruce 4.86 died
1079 29.00 0.25 0 black spruce 4.86 alive
3455 32.00 0.80 1 black spruce 5.00 died
[659 rows x 6 columns]
각 likelihood는 관측치들이 모두 독립적으로 발생했다고 가정했을 때의 확률값이고,
모든 데이터가 관찰될 likelihood = \(p_1^7 (1-p_1)^{83} \cdot p_2^7 (1-p_2)^{85} \cdot p_3^{18} \cdot (1-p_3)^{73} \cdots\)
높은 precision은 사용자가 검색 결과에서 불필요한 정보를 적게 얻고, 대부분 유용한 정보를 얻는다는 것을 의미
하지만, 관련은 있지만 놓친 정보는 많을 수 있음.
검색 결과의 정확도를 평가
검사를 통해 암으로 진단받은 사람이 실제로 암을 가지고 있을 확률
높은 precision은 암 진단을 받은 사람에게 잘못된 암 진단을 내리지 않는 것을 의미
Recall (재현율): \(y=1\)일 때, \(\hat{y}=1\)일 확률
실제 스팸 이메일 중에서 시스템이 스팸으로 올바르게 분류한 비율
높은 recall은 대부분의 스팸 이메일이 정확히 스팸으로 분류된다는 것을 의미
하지만, 정상 이메일이 스팸함으로 분류될 수 있음.
정보검색에서 관련 문서들 중에서 실제로 시스템이 반환한 문서의 비율
높은 recall은 사용자가 찾고자 하는 모든 관련 정보를 검색 결과에서 얻을 수 있다는 것을 의미
하지만, 관련 없는 정보도 많이 포함될 수 있음.
검색 시스템의 탐색 능력을 평가
실제로 암을 가진 사람 중에서 검사를 통해 암으로 진단받은 사람의 비율
높은 recall은 암 환자를 놓치지 않고 모두 찾아낸다는 것을 의미
참 양성, true positive rate (TPR), 또는 sensitivity라고도 함
precision & recall trade-off: 이 둘 사이에는 종종 상충 관계가 있음. Precision을 높이기 위해 예측 확률에 대한 더 엄격한 기준을 사용하면 recall이 낮아질 수 있고, 반대로 recall을 높이기 위해 더 느슨한 기준을 사용하면 precision이 낮아질 수 있음.
정보 검색의 맥락: precision을 높여 사용자가 불필요한 정보를 받지 않게 해주며, recall을 높여 필요한 정보를 놓치지 않도록 함.
스팸 필터링의 맥락: precision을 높여 자주 스팸 메일함을 확인하지 않아도 되도록 하며, recall을 높여 대부분의 스팸 이메일이 차단되어 정상 메일함에서 보이지 않도록 쾌적한 환경을 제공할 수 있음.
Sensitivity & specificity: 의학 분야에서 진단 테스트의 정확성을 평가하는데에서 유래
Precision & recall이 주로 양성(positive) 클래스에 초점을 두는 반면, sensitivity & specificity는 양성(positive)와 음성(negative) 클래스 모두에 초점을 둠.
특히, 이상치 탐지라든가 희귀 질병 진단 등 불균형 데이터셋에서는 precision & recall이 더 유용할 수 있음.
Sensitivity (민감도): \(y=1\)일 때, \(\hat{y}=1\)일 확률
Recall, true positive rate (TPR)
높은 민감도는 질병이 있는 사람을 놓치지 않고 모두 찾아내는 것(참 양성)을 의미
Specificity (특이도): \(y=0\)일 때, \(\hat{y}=0\)일 확률
True negative rate (TNR)
실제로 질병이 없는 사람을 얼마나 잘 (질병이 없다고) 식별하는지(참 음성)를 나타냄
높은 특이도는 질병이 없는 사람에게 질병이 있다고 잘못 진단하지 않는 것을 의미함.
희귀 질병을 진단하는 경우, 정상인(negative)을 올바로 식별하는 것은 매우 쉽기 때문에 질병이 없는 사람에 대한 판별력이 과대추정될 수 있음; precision이 더 유용할 수 있음.
1 - FPR (False Positive Rate; false alarm)
0, 1을 바꿨을 때의 즉, 0을 기준으로 했을 때의 recall
특이도가 높은 테스트라면 질병이 없다는 판정은 신뢰할 만함. 즉, 0(음성)을 기준으로한 sensitivity가 높은 것이 됨.
거짓 음성 비율: False negative rate (FNR); \(y=1\)일 때, \(\hat{y}=0\)일 확률: 1 - sensitivity
- 정상이라고 진단된 환자가 실제로 암인 확률
거짓 양성 비율: False positive rate (FPR) 또는 False alarm; \(y=0\)일 때, \(\hat{y}=1\)일 확률: 1 - specificity
- 정상인 사람이 검사를 통해 암으로 진단받을 확률
예측된 확률에 대한 임계치를 조정함에 따라 옳은 예측과 틀린 예측의 비율이 어떻게 달라지는지 살펴봄으로써 임계치를 설정하는데 도움을 줌
잘못된 예측에 대한 비용이 다르다면, 특정 임계치를 선택하는 것을 고려해야 함
code
from sklearn.metrics import roc_curvefpr, tpr, thresholds = roc_curve(test_pred.y, test_pred.pred_prob)roc = pd.DataFrame( {"thresholds": thresholds.round(2),"False Pos": fpr,"sensitivity(TPR)": tpr,"specificity(TNR)": 1- fpr })roc.sort_values("thresholds").head(10)## ROC curve# 위에서 얻은 roc 데이터을 사용하여 그리거나# RocCurveDisplay를 이용from sklearn.metrics import RocCurveDisplayroc_display = RocCurveDisplay(fpr=fpr, tpr=tpr).plot()# scikit-learn의 visualization 문서 참조# https://scikit-learn.org/stable/visualizations.html
code
from sklearn.metrics import precision_recall_curveprecision, recall, thresholds = precision_recall_curve(test_pred.y, test_pred.pred_prob)pr = pd.DataFrame( {"thresholds": np.append(thresholds.round(2), 1),"precision": precision,"recall": recall,"f1-score": 2* precision * recall / (precision + recall) }).query("recall > 0.1")pr.sort_values("thresholds").head(10)# precision-recall curve# 위에서 얻은 pr 데이터을 사용하여 그리거나# PrecisionRecallDisplay를 이용from sklearn.metrics import PrecisionRecallDisplaypr_display = PrecisionRecallDisplay(precision=precision, recall=recall).plot()# scikit-learn의 visualization 문서 참조# https://scikit-learn.org/stable/visualizations.html
위에서 언급한 클래스 불균형의 예들: 이상치 탐지, 희귀 질병 진단 등에서는 precision & recall이 더 유용할 수 있음.
Classifier의 전체적 성능에 대한 지표
ROC curve:
AUC: Area Under the Curve = Concordance Index
각 specificity값에 대한 sensitivity의 합; 모형(classifier) 대한 전반적 평가
0.5: random guess, 1: perfect prediction
Concordance Index(c-index): 모든 서로 다른 클래스의 Y쌍, 예를 들어 \((0_i, 1_j)\)에 대해서 해당하는 예측된 확률의 크기가 \(p_i < p_j\) 인 비율, 즉 순서가 맞는(concordance) 비율
위 blowdown 데이터셋 대한 모형의 AUC, concordance index 계산
Code
from sklearn.metrics import roc_auc_score, aucroc_auc = roc_auc_score(test_pred.y, test_pred.pred_prob)print(f"AUC: {roc_auc:.2f}") # 또는 auc(fpr, tpr)
AUC: 0.81
Code
def concordance(y, pred_prob): concord =0 total =0 test_pred = pd.DataFrame({"y": y, "pred_prob": pred_prob})for idx, case in test_pred.iterrows(): other_cases = test_pred[test_pred.index != idx]# 같은 케이스 제외 df = other_cases[case["y"] != other_cases["y"]] concord += ( (case["y"] - df["y"] >0) == (case["pred_prob"] - df["pred_prob"] >0) ).sum() total +=len(df)return concord / totalprint(f"Concordance index: {concordance(test_pred.y, test_pred.pred_prob):.2f}")
Concordance index: 0.81
Precision-recall:
Average precision: the area under the precision-recall curve
Code
from sklearn.metrics import average_precision_scoreaverage_precision = average_precision_score(test_pred.y, test_pred.pred_prob)print(f"Average Precision: {average_precision:.2f}")
Average Precision: 0.68
The classification report in scikit-learn Thershold: 0.4인 경우,
code
from sklearn.metrics import classification_report, recall_scorethreshold =.4pred_class = test_pred["pred_prob"] >= thresholdprint(classification_report(test_pred["y"], pred_class))# specificity: 양성을 0으로 보고 recall 계산specificity = recall_score(test_pred["y"], pred_class, pos_label=0)print(f"specificity (recall for negative): {specificity:.2f}")
왼편의 logit 함수는 증가함수이므로, threshold를 정하면 선형 결정 경계(decision boundary)가 나타남: hyperplane
비선형 함수로 모형을 세우면, decision boundary가 비선형으로 나타남.
예를 들어, \(b_0 + b_1 \cdot X + b_2 \cdot X^2\) 등
밑은 두 예측변수 \(X_1= log2d\) (diameter)와 \(X_2=s\) (severity)로 나무가 쓰러질지 여부를 선형함수로 만들었을 때의 decision boundary (threshold=0.5)
즉, Log odds(logit): \(\displaystyle log\left(\frac{\hat{p}}{1 - \hat{p}}\right) = b_0 + b_1 \cdot X_1 + b_2 \cdot X_2\)
밑은 두 예측변수 \(X_1= log2d\) (diameter)와 \(X_2=s\) (severity)에 대한 모든 2, 3차항을 추가했을 때 decision boundary; 10개의 파라미터
예를 들어, \(X_1^2, ~ X_1 \cdot X_2, ~X_1^3, ~X_1^2 \cdot X_2, ~X_2^3\)
Multiclass (K > 2)의 경우
Multinomial logistic regression
클래스 \(C_k\)에 속할 확률을 예측하는 모형: Base 클래스를 기준으로 나머지 클래스에 대한 확률을 예측
Binary 경우에서의 확률 \(\displaystyle P(Y=1|X=x) = \frac{1}{1+e^{-(b_0 + b_1 \cdot x)}}=\frac{e^{b_0 + b_1 \cdot x}}{1 + e^{b_0 + b_1 \cdot x}}\)을 다음과 같이 확장할 수 있음.
# split the data into training and test setsX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)lr = LogisticRegression(max_iter=1000)lr.fit(X_train, y_train);
