교호작용 항¶

import pandas as pd

# Advertising: 광고비(TV, Radio, Newspaper)와 매출(Sales), 200개 시장
advertising = pd.read_csv("https://www.statlearning.com/s/Advertising.csv",
                          index_col=0)
advertising = advertising.rename(columns={"radio": "Radio",
                                          "newspaper": "Newspaper",
                                          "sales": "Sales"})

# Credit: 신용카드 잔액(Balance)과 소득(Income), 학생 여부(Student), 400명
credit = pd.read_csv("https://raw.githubusercontent.com/vincentarelbundock/"
                     "Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/ISLR/Credit.csv", index_col=0)

print(advertising[["TV", "Radio", "Sales"]].describe().round(2).to_string())
print()
print(credit.groupby("Student")[["Income", "Balance"]].mean().round(2).to_string())

           TV   Radio   Sales
count  200.00  200.00  200.00
mean   147.04   23.26   14.02
std     85.85   14.85    5.22
min      0.70    0.00    1.60
25%     74.38    9.98   10.38
50%    149.75   22.90   12.90
75%    218.82   36.52   17.40
max    296.40   49.60   27.00

         Income  Balance
Student                 
No        44.99   480.37
Yes       47.29   876.82

print("Credit 의 학생 수:", credit['Student'].value_counts().to_dict())
tab = credit.groupby("Student")[["Income", "Balance"]].agg(['count', 'mean', 'std'])
print(tab.round(2).to_string())
print(f"Income 의 범위: {credit['Income'].min():.2f} ~ {credit['Income'].max():.2f}")
print(f"Balance 가 0 인 사람: {(credit['Balance'] == 0).sum()}명")
print(f"TV x Radio 의 범위: {(advertising['TV'] * advertising['Radio']).min():.2f} ~ "
      f"{(advertising['TV'] * advertising['Radio']).max():.2f}")
print(f"요약표의 최댓값으로 셈하면 296.40 x 49.60 = {296.40 * 49.60:.2f}")

Credit 의 학생 수: {'No': 360, 'Yes': 40}
        Income               Balance                
         count   mean    std   count    mean     std
Student                                             
No         360  44.99  34.91     360  480.37  439.41
Yes         40  47.29  38.55      40  876.82  490.00
Income 의 범위: 10.35 ~ 186.63
Balance 가 0 인 사람: 90명
TV x Radio 의 범위: 0.00 ~ 13540.41
요약표의 최댓값으로 셈하면 296.40 x 49.60 = 14701.44

import statsmodels.formula.api as smf

# 수식 문법. R 의 것을 그대로 따른다
# 별표는 주효과와 교호작용을 한꺼번에 넣는다는 뜻이다
model = smf.ols('Sales ~ TV * Radio', data=advertising).fit()
# summary()는 실행 날짜와 시각을 함께 찍으므로 계수 표만 인쇄한다.
print(model.summary().tables[1])
print(f"R^2 = {model.rsquared:.4f}")

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      6.7502      0.248     27.233      0.000       6.261       7.239
TV             0.0191      0.002     12.699      0.000       0.016       0.022
Radio          0.0289      0.009      3.241      0.001       0.011       0.046
TV:Radio       0.0011   5.24e-05     20.727      0.000       0.001       0.001
==============================================================================
R^2 = 0.9678

b = model.params
print("계수:", {k: round(v, 6) for k, v in b.items()})
for r in [0, advertising['Radio'].mean(), 30, advertising['Radio'].max()]:
    print(f"  Radio={r:6.2f}: dSales/dTV = {b['TV'] + b['TV:Radio'] * r:.4f}")
for t_ in [0, advertising['TV'].mean(), advertising['TV'].max()]:
    print(f"  TV={t_:7.2f}: dSales/dRadio = {b['Radio'] + b['TV:Radio'] * t_:.4f}")
print(f"dSales/dTV = 0 이 되는 Radio = {-b['TV'] / b['TV:Radio']:.2f}"
      f"   (Radio 의 범위 {advertising['Radio'].min()} ~ {advertising['Radio'].max()})")

additive = smf.ols('Sales ~ TV + Radio', data=advertising).fit()
F = (additive.ssr - model.ssr) / (model.ssr / model.df_resid)
t = model.tvalues['TV:Radio']
print(f"가법모형 R^2 = {additive.rsquared:.4f},  교호작용 R^2 = {model.rsquared:.4f}")
print(f"중첩 F = {F:.4f},  t^2 = {t ** 2:.4f},  p = {model.pvalues['TV:Radio']:.3e}")

계수: {'Intercept': 6.75022, 'TV': 0.019101, 'Radio': 0.02886, 'TV:Radio': 0.001086}
  Radio=  0.00: dSales/dTV = 0.0191
  Radio= 23.26: dSales/dTV = 0.0444
  Radio= 30.00: dSales/dTV = 0.0517
  Radio= 49.60: dSales/dTV = 0.0730
  TV=   0.00: dSales/dRadio = 0.0289
  TV= 147.04: dSales/dRadio = 0.1886
  TV= 296.40: dSales/dRadio = 0.3509
dSales/dTV = 0 이 되는 Radio = -17.58   (Radio 의 범위 0.0 ~ 49.6)
가법모형 R^2 = 0.8972,  교호작용 R^2 = 0.9678
중첩 F = 429.5905,  t^2 = 429.5905,  p = 2.758e-51

# 범주형 변수가 든 예
# statsmodels 는 범주형을 알아서 가변수로 바꾼다
model = smf.ols('Balance ~ Income + C(Student) + Income:C(Student)',
                data=credit).fit()
print(model.summary().tables[1])
print(f"R^2 = {model.rsquared:.4f}")

============================================================================================
                               coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------------------
Intercept                  200.6232     33.698      5.953      0.000     134.373     266.873
C(Student)[T.Yes]          476.6758    104.351      4.568      0.000     271.524     681.827
Income                       6.2182      0.592     10.502      0.000       5.054       7.382
Income:C(Student)[T.Yes]    -1.9992      1.731     -1.155      0.249      -5.403       1.404
============================================================================================
R^2 = 0.2799

import numpy as np

pr = model.params
print("비학생:", f"{pr['Intercept']:.4f} + {pr['Income']:.4f} * Income")
print("학생  :", f"{pr['Intercept'] + pr['C(Student)[T.Yes]']:.4f} + "
                 f"{pr['Income'] + pr['Income:C(Student)[T.Yes]']:.4f} * Income")
print(f"두 직선이 만나는 Income = {pr['C(Student)[T.Yes]'] / -pr['Income:C(Student)[T.Yes]']:.2f}"
      f"   (자료의 최댓값 {credit['Income'].max():.2f})")

# 집단별로 따로 적합한다.
no = smf.ols('Balance ~ Income', data=credit[credit['Student'] == 'No']).fit()
yes = smf.ols('Balance ~ Income', data=credit[credit['Student'] == 'Yes']).fit()
print(f"따로 적합 비학생: 절편 {no.params['Intercept']:.4f}  기울기 {no.params['Income']:.4f}  (n={int(no.nobs)})")
print(f"따로 적합 학생  : 절편 {yes.params['Intercept']:.4f}  기울기 {yes.params['Income']:.4f}  (n={int(yes.nobs)})")
print(f"계수의 최대 차이 = "
      f"{max(abs(no.params['Income'] - pr['Income']), abs(yes.params['Income'] - (pr['Income'] + pr['Income:C(Student)[T.Yes]']))):.2e}")
print(f"잔차표준오차: 묶음 {np.sqrt(model.mse_resid):.2f},  "
      f"비학생 {np.sqrt(no.mse_resid):.2f},  학생 {np.sqrt(yes.mse_resid):.2f}")
print(f"학생 기울기의 SE: 묶은 모형 쪽은 따로 적합의 {yes.bse['Income']:.4f} 과 다르다")
print(f"교호작용 t = {model.tvalues['Income:C(Student)[T.Yes]']:.4f},  "
      f"p = {model.pvalues['Income:C(Student)[T.Yes]']:.4f}")
print(f"교호작용 없는 모형 R^2 = {smf.ols('Balance ~ Income + C(Student)', data=credit).fit().rsquared:.4f},  "
      f"넣은 모형 {model.rsquared:.4f}")

비학생: 200.6232 + 6.2182 * Income
학생  : 677.2990 + 4.2190 * Income
두 직선이 만나는 Income = 238.44   (자료의 최댓값 186.63)
따로 적합 비학생: 절편 200.6232  기울기 6.2182  (n=360)
따로 적합 학생  : 절편 677.2990  기울기 4.2190  (n=40)
계수의 최대 차이 = 1.78e-14
잔차표준오차: 묶음 391.62,  비학생 382.59,  학생 468.27
학생 기울기의 SE: 묶은 모형 쪽은 따로 적합의 1.9452 과 다르다
교호작용 t = -1.1547,  p = 0.2489
교호작용 없는 모형 R^2 = 0.2775,  넣은 모형 0.2799