선형 회귀 | 프리슈-워-로벨 정리 | 더미 변수 : 실무로 통하는 인과추론 학습기(2)

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from matplotlib import pyplot as plt
import seaborn as sns


data = pd.read_csv('student_grades.csv')

formula = 'Grade ~ Sleep_Hours'
model = smf.ols(formula, data=data).fit()
model.summary().tables[1]

formula = 'Grade ~ Sleep_Hours+Tutoring_Hours+Study_Hours+Family_Support'
model = smf.ols(formula, data=data).fit()
model.summary().tables[1]

# 회귀 모델을 생성하고 데이터에 적합시킵니다.
debiasing_model = smf.ols('Sleep_Hours ~ Family_Support + Study_Hours + Tutoring_Hours', data = data).fit()

debiasing_model.resid + data['Sleep_Hours'].mean()

# 회귀 모델의 잔차(residuals)와 평균을 더하여 sleep_hours_res 열을 생성합니다.
data_deb = data.assign(
    sleep_hours_res = (debiasing_model.resid + data['Sleep_Hours'].mean()))

# Sleep_Hours에 대한 편향제거된 회귀모델
model_w_deb_data = smf.ols('Grade ~ sleep_hours_res', data = data_deb).fit()

print(model_w_deb_data.summary().tables[1])

# Grade에 대한 잡음제거
denoising_model = smf.ols(
    'Grade ~ Family_Support + Study_Hours + Tutoring_Hours', data = data).fit()

data_denoise = data_deb.assign(
    grade_res=denoising_model.resid + data_deb["Grade"].mean())

#표준오차 Family_Support + Tutoring_Hours + Study_Hours 세 변수 활용
model_se = smf.ols(
    'Grade ~ Family_Support + Tutoring_Hours + Study_Hours', data = data).fit()

print("SE regression:", model_se.bse["Family_Support"])

>>> SE regression: 0.15239839952386264

model_family_support_aux = smf.ols(
    'Family_Support ~ Tutoring_Hours + Study_Hours', data = data
).fit()

# subtract the degrees of freedom - 4 model parameters - from N.
se_formula = (np.std(model_se.resid)
              /(np.std(model_family_support_aux.resid)*np.sqrt(len(data)-4)))

print("SE formula:   ", se_formula)

>>> SE formula:    0.15239839952386258

# Grade, sleep_hours_res 활용한 최종모델
model_w_orthogonal = smf.ols('grade_res ~ sleep_hours_res', data = data_denoise).fit()

print(model_w_orthogonal.summary().tables[1])

print(model_w_orthogonal.summary())

# ATE는 T의 계수로부터 추정
ATE = model_w_orthogonal.params['sleep_hours_res']
print('ATE:', ATE)

# Study_Hours 값을 10개의 구간으로 나눕니다
data['Study_Hours_Bin'] = pd.cut(data['Study_Hours'], bins=10, labels=False)
data

data['Study_Hours_Bin'].value_counts()

plt.figure(figsize=(15,6))
sns.histplot(data=data,
             x="Study_Hours",
             hue="Study_Hours_Bin")
plt.title("Conditional random experiment")
plt.show()

data_dummies = (data
                     .join(pd.get_dummies(data["Study_Hours_Bin"],
                                          prefix="sb",
                                          drop_first=True)))
data_dummies.head()

model = smf.ols(
    "Grade ~ Sleep_Hours + sb_1+sb_2+sb_3+sb_4+sb_5+sb_6+sb_7+sb_8+sb_9",
    data=data_dummies
).fit()

model.summary().tables[1]

model = smf.ols("Grade ~ Sleep_Hours + C(Study_Hours_Bin)",
                data=data).fit()

model.summary().tables[1]

def regress(df, t, y):
    return round(smf.ols(f'{y} ~ {t}', data = df).fit().params[t], 6)

effect_by_group = data_dummies.groupby('Study_Hours_Bin').apply(regress, y = 'Grade', t = 'Sleep_Hours')
effect_by_group

>>>
Study_Hours_Bin
0    9.015105
1    6.002520
2    0.747068
3   -1.999126
4   -0.464798
5   -0.066598
6   -0.019569
7   -0.000000
8   -0.000000
9    0.000000
dtype: float64