회귀 진단 (주택)¶

import pandas as pd

# "Practical Statistics for Data Scientists" 저장소의 자료. 탭으로 구분되어 있다.
url = ("https://raw.githubusercontent.com/gedeck/"
       "practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/house_sales.csv")
house = pd.read_csv(url, sep='\t')
house_98105 = house.loc[house['ZipCode'] == 98105, :]

print(f"전체 {len(house)}건 중 98105 지역 {len(house_98105)}건")
print(house_98105[['AdjSalePrice', 'SqFtTotLiving', 'Bedrooms']].describe().round(1).to_string())

전체 22687건 중 98105 지역 313건
       AdjSalePrice  SqFtTotLiving  Bedrooms
count         313.0          313.0     313.0
mean       756277.3         2069.9       3.4
std        391463.3          905.4       1.1
min        119748.0          490.0       1.0
25%        521101.0         1410.0       3.0
50%        630041.0         1850.0       3.0
75%        851954.0         2570.0       4.0
max       3013254.0         5570.0       9.0

import numpy as np
from scipy.stats import skew

price = house_98105['AdjSalePrice']
print(f"평균 {price.mean():,.0f},  중앙값 {price.median():,.0f},  "
      f"차 {price.mean() - price.median():,.0f}")
print(f"변동계수 = sd/mean = {price.std() / price.mean():.4f}")
print(f"최대/최소 = {price.max() / price.min():.1f},  최대/중앙값 = {price.max() / price.median():.2f}")
print(f"\n왜도:  가격 {skew(price):+.4f},  로그가격 {skew(np.log(price)):+.4f}")
print(f"sd(log 가격) = {np.log(price).std():.4f}   (변동계수의 근사값)")

평균 756,277,  중앙값 630,041,  차 126,236
변동계수 = sd/mean = 0.5176
최대/최소 = 25.2,  최대/중앙값 = 4.78

왜도:  가격 +2.2581,  로그가격 +0.6970
sd(log 가격) = 0.4199   (변동계수의 근사값)

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import OLSInfluence

# assign(const=1) 이 절편 열을 만든다. statsmodels 의 OLS 는 절편을
# 자동으로 넣지 않는다.
predictors = ['SqFtTotLiving', 'SqFtLot', 'Bathrooms', 'Bedrooms', 'BldgGrade']
X = house_98105[predictors].assign(const=1)
y = house_98105['AdjSalePrice']

results = sm.OLS(y, X).fit()
influence = OLSInfluence(results)

# 세 진단량을 함께 본다. 스튜던트화 잔차는 그 점이 얼마나 벗어났는지,
# 지렛값은 설명변수 쪽에서 얼마나 외따로 있는지, Cook 거리는 그 둘을
# 합쳐 결론을 얼마나 흔드는지를 잰다.
studentized_resids = influence.resid_studentized_internal
hat_values = influence.hat_matrix_diag
cooks_dist, _ = influence.cooks_distance

print(results.summary().tables[1])
print(f"n = {len(y)}, 문턱 4/n = {4/len(y):.4f}")
print(f"Cook 거리 최댓값 = {cooks_dist.max():.4f} (관측 {cooks_dist.argmax()})")
print(f"문턱을 넘는 관측값 수 = {(cooks_dist > 4/len(y)).sum()}")

=================================================================================
                    coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
SqFtTotLiving   209.6023     24.408      8.587      0.000     161.574     257.631
SqFtLot          38.9333      5.330      7.305      0.000      28.445      49.421
Bathrooms      2282.2641      2e+04      0.114      0.909    -3.7e+04    4.16e+04
Bedrooms      -2.632e+04   1.29e+04     -2.043      0.042   -5.17e+04    -973.867
BldgGrade        1.3e+05   1.52e+04      8.533      0.000       1e+05     1.6e+05
const         -7.725e+05   9.83e+04     -7.861      0.000   -9.66e+05   -5.79e+05
=================================================================================
n = 313, 문턱 4/n = 0.0128
Cook 거리 최댓값 = 0.5608 (관측 152)
문턱을 넘는 관측값 수 = 20

import numpy as np

r = np.asarray(studentized_resids)
h = np.asarray(hat_values)
D = np.asarray(cooks_dist)
e = np.asarray(results.resid)
n, k = len(y), X.shape[1]
s2 = results.mse_resid

print(f"sum h_ii = {h.sum():.10f}   (p = {k} 여야 한다)")
print(f"평균 h = {h.mean():.6f} = p/n,   문턱 2p/n = {2 * k / n:.6f}")
print(f"지렛값이 문턱을 넘는 관측 = {(h > 2 * k / n).sum()}건")

# 두 가지 Cook 공식이 모두 같은 값을 주는가
D1 = r ** 2 / k * h / (1 - h)
D2 = e ** 2 * h / (k * s2 * (1 - h) ** 2)
print(f"\n|D1 - cooks| 최대 = {np.abs(D1 - D).max():.2e}")
print(f"|D2 - cooks| 최대 = {np.abs(D2 - D).max():.2e}")

# 가장 영향이 큰 집을 뜯어본다
i = int(D.argmax())
print(f"\n위치 {i}:  D = {D[i]:.4f},  r = {r[i]:.4f},  h = {h[i]:.4f}"
      f"  (평균 h 의 {h[i] / h.mean():.1f}배)")
print(f"  실제 가격 {y.iloc[i]:,.0f},  적합값 {results.fittedvalues.iloc[i]:,.0f},"
      f"  잔차 {e[i]:,.0f}")
print(f"  두 인자:  r^2/p = {r[i] ** 2 / k:.4f},"
      f"  h/(1-h) = {h[i] / (1 - h[i]):.4f},  곱 = {r[i] ** 2 / k * h[i] / (1 - h[i]):.4f}")

sum h_ii = 6.0000000000   (p = 6 여야 한다)
평균 h = 0.019169 = p/n,   문턱 2p/n = 0.038339
지렛값이 문턱을 넘는 관측 = 19건

|D1 - cooks| 최대 = 5.55e-17
|D2 - cooks| 최대 = 1.11e-16

위치 152:  D = 0.5608,  r = 4.9420,  h = 0.1211  (평균 h 의 6.3배)
  실제 가격 3,013,254,  적합값 2,186,265,  잔차 826,989
  두 인자:  r^2/p = 4.0705,  h/(1-h) = 0.1378,  곱 = 0.5608

# 문턱을 넘는 관측값을 빼고 다시 적합해 계수가 얼마나 달라지는지 본다.
# 크게 달라진다면 결론이 몇 채의 집에 기대고 있다는 뜻이다.
threshold_cooks = 4 / len(y)
mask_keep = cooks_dist < threshold_cooks

X_filtered = X[mask_keep]
y_filtered = y[mask_keep]
results_filtered = sm.OLS(y_filtered, X_filtered).fit()

# 계수를 견준다
comparison = pd.DataFrame({
    'Original': results.params,
    'Filtered': results_filtered.params,
})
print(comparison.round(2).to_string())
print(f"\n제거된 관측값: {(~mask_keep).sum()}건")
print(f"R^2: {results.rsquared:.4f} -> {results_filtered.rsquared:.4f}")

                Original   Filtered
SqFtTotLiving     209.60     201.83
SqFtLot            38.93      42.59
Bathrooms        2282.26    1664.95
Bedrooms       -26320.27  -23734.89
BldgGrade      130000.10  111175.45
const         -772549.86 -644099.89

제거된 관측값: 20건
R^2: 0.7954 -> 0.8415

import numpy as np

h = np.asarray(influence.hat_matrix_diag)
e = np.asarray(results.resid)
D = np.asarray(cooks_dist)
n, k = len(y), X.shape[1]
s2 = results.mse_resid
i = int(D.argmax())

# 한 점만 빼는 공식 (셔먼-모리슨)
XtXinv = np.linalg.inv(X.values.T @ X.values)
beta_formula = results.params.values - XtXinv @ X.values[i] * e[i] / (1 - h[i])

# 실제로 빼고 다시 적합
keep_one = np.ones(n, dtype=bool)
keep_one[i] = False
beta_refit = sm.OLS(y[keep_one], X[keep_one]).fit().params.values

print(f"{'':>14}{'공식':>16}{'실제 재적합':>16}")
for name, a, b in zip(X.columns, beta_formula, beta_refit):
    print(f"{name:>14}{a:>16.4f}{b:>16.4f}")
print(f"\n최대 상대오차 = {np.abs((beta_formula - beta_refit) / beta_refit).max():.2e}")

# Cook 거리가 정말 그 변화의 크기인가
d = results.params.values - beta_formula
print(f"\n이차형식으로 계산한 D = {d @ (X.values.T @ X.values) @ d / (k * s2):.12f}")
print(f"cooks_distance        = {D[i]:.12f}")

# 20건을 뺐을 때 R^2 이 오른 까닭
print(f"\nSSE: {results.ssr:.4g} -> {results_filtered.ssr:.4g}"
      f"  ({results_filtered.ssr / results.ssr:.1%})")
print(f"SST: {results.centered_tss:.4g} -> {results_filtered.centered_tss:.4g}"
      f"  ({results_filtered.centered_tss / results.centered_tss:.1%})")

                            공식          실제 재적합
 SqFtTotLiving        217.5874        217.5874
       SqFtLot         31.5277         31.5277
     Bathrooms       9031.4566       9031.4566
      Bedrooms     -28072.1958     -28072.1958
     BldgGrade     120524.2661     120524.2661
         const    -693113.5065    -693113.5065

최대 상대오차 = 1.77e-14

이차형식으로 계산한 D = 0.560788525781
cooks_distance        = 0.560788525781

SSE: 9.781e+12 -> 4.069e+12  (41.6%)
SST: 4.781e+13 -> 2.568e+13  (53.7%)

from statsmodels.stats.diagnostic import het_breuschpagan

# 집값 자료는 비싼 집일수록 오차도 커지는 것이 보통이라, 등분산 검정이
# 기각되는 일이 흔하다. 그럴 때는 로그 변환이나 로버스트 표준오차로 간다.
bp_stat, bp_pval, _, _ = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_pval:.4f}")

Breusch-Pagan p-value: 0.0000

import numpy as np
import statsmodels.api as sm
from scipy.stats import chi2
from statsmodels.stats.diagnostic import het_breuschpagan

e = np.asarray(results.resid)
n = len(y)

aux = sm.OLS(e ** 2, X).fit()
lm = n * aux.rsquared
print(f"보조회귀 R^2 = {aux.rsquared:.6f}")
print(f"LM = n R^2   = {lm:.6f}")
print(f"statsmodels  = {het_breuschpagan(results.resid, X)[0]:.6f}")
print(f"p = chi2(5).sf(LM) = {chi2.sf(lm, 5):.4e}   (0.0000 이 아니다)")

# 로그 변환이 얼마나 고치는가
res_log = sm.OLS(np.log(y), X).fit()
lm_log, p_log, _, _ = het_breuschpagan(res_log.resid, X)
print(f"\n로그 모형:  LM = {lm_log:.4f},  p = {p_log:.4f}")
print(f"로그 모형의 R^2 = {res_log.rsquared:.4f}  (원 모형 {results.rsquared:.4f})")

보조회귀 R^2 = 0.199511
LM = n R^2   = 62.446794
statsmodels  = 62.446794
p = chi2(5).sf(LM) = 3.7899e-12   (0.0000 이 아니다)

로그 모형:  LM = 12.6167,  p = 0.0272
로그 모형의 R^2 = 0.7651  (원 모형 0.7954)

ext_resids = influence.resid_studentized_external
int_resids = influence.resid_studentized_internal
discrepancy = np.abs(ext_resids - int_resids)
top_idx = np.argsort(discrepancy)[-5:]

resid_sq = results.resid ** 2
aux_model = sm.OLS(resid_sq, X).fit()
bp_stat_manual = len(y) * aux_model.rsquared
from scipy import stats
bp_pval_manual = 1 - stats.chi2.cdf(bp_stat_manual, X.shape[1] - 1)