import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
Policy Learning#
https://bookdown.org/stanfordgsbsilab/ml-ci-tutorial/policy-evaluation-i—binary-treatment.html
https://bookdown.org/stanfordgsbsilab/ml-ci-tutorial/policy-learning-i—binary-treatment.html
이 둘을 파이썬으로 바꿔야하지 않을까 생각합니다.
causal ml 책에 있는 내용은 사실상 내용이 너무 적어서 굳이 볼 필요는 없는 것 같습니다. 우리가 파이썬 코드를 만들어야 합니다.
질문이나 의견을 남겨주세요.#
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import matplotlib.patches as mpatches
# Set random seed for reproducibility
np.random.seed(42)
# Generate data
n = 1000
p = 4
X = np.random.uniform(0, 1, (n, p))
W = np.random.binomial(1, 0.5, n) # Independent from X and Y
Y = 0.5 * (X[:, 0] - 0.5) + (X[:, 1] - 0.5) * W + 0.1 * np.random.randn(n)
# Normalize Y for plotting
y_norm = 1 - (Y - Y.min()) / (Y.max() - Y.min())
# First plot: All data points
fig1, ax1 = plt.subplots(1, 1, figsize=(8, 6))
for i in range(n):
if W[i] == 1:
ax1.scatter(X[i, 0], X[i, 1], marker='o', s=100,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
else:
ax1.scatter(X[i, 0], X[i, 1], marker='D', s=80,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
ax1.set_xlabel('X1', fontsize=12)
ax1.set_ylabel('X2', fontsize=12)
ax1.set_title('All Data Points (○: Treated, ◇: Untreated)', fontsize=14)
plt.show()
# Second plot: Separated by treatment
fig2, (ax2, ax3) = plt.subplots(1, 2, figsize=(14, 6))
# Untreated group
untreated_idx = W == 0
ax2.scatter(X[untreated_idx, 0], X[untreated_idx, 1], marker='D', s=80,
c=y_norm[untreated_idx], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
ax2.set_xlabel('X1', fontsize=12)
ax2.set_ylabel('X2', fontsize=12)
ax2.set_title('Untreated', fontsize=14)
# Treated group
treated_idx = W == 1
ax3.scatter(X[treated_idx, 0], X[treated_idx, 1], marker='o', s=100,
c=y_norm[treated_idx], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
ax3.set_xlabel('X1', fontsize=12)
ax3.set_ylabel('X2', fontsize=12)
ax3.set_title('Treated', fontsize=14)
plt.show()
# Third plot: Policy regions
fig3, ax4 = plt.subplots(1, 1, figsize=(8, 6))
# Define colors with transparency
col1 = (0.9960938, 0.7539062, 0.0273438, 0.35) # Yellow-ish
col2 = (0.250980, 0.690196, 0.650980, 0.35) # Teal-ish
# Draw policy regions
rect1 = Rectangle((-0.1, -0.1), 0.6, 1.2, linewidth=0,
edgecolor='none', facecolor=col1, hatch='///')
rect2 = Rectangle((0.5, -0.1), 0.6, 0.6, linewidth=0,
edgecolor='none', facecolor=col1, hatch='///')
rect3 = Rectangle((0.5, 0.5), 0.6, 0.6, linewidth=0,
edgecolor='none', facecolor=col2, hatch='///')
ax4.add_patch(rect1)
ax4.add_patch(rect2)
ax4.add_patch(rect3)
# Plot data points
for i in range(n):
if W[i] == 1:
ax4.scatter(X[i, 0], X[i, 1], marker='o', s=100,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
else:
ax4.scatter(X[i, 0], X[i, 1], marker='D', s=80,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
# Add text labels
ax4.text(0.75, 0.75, 'TREAT (A)', fontsize=16, ha='center', va='center')
ax4.text(0.25, 0.25, 'DO NOT TREAT (A^C)', fontsize=16, ha='left', va='center')
ax4.set_xlabel('X1', fontsize=12)
ax4.set_ylabel('X2', fontsize=12)
ax4.set_xlim(-0.1, 1.1)
ax4.set_ylim(-0.1, 1.1)
ax4.set_title('Policy Regions', fontsize=14)
plt.show()
# Policy Evaluation Methods
print("=" * 60)
print("POLICY EVALUATION RESULTS")
print("=" * 60)
# Method 1: Value of policy A (only valid in randomized setting)
A = (X[:, 0] > 0.5) & (X[:, 1] > 0.5)
value_estimate = np.mean(Y[A & (W == 1)]) * np.mean(A) + \
np.mean(Y[~A & (W == 0)]) * np.mean(~A)
value_stderr = np.sqrt(
np.var(Y[A & (W == 1)]) / np.sum(A & (W == 1)) * np.mean(A)**2 +
np.var(Y[~A & (W == 0)]) / np.sum(~A & (W == 0)) * np.mean(~A)**2
)
print(f"\nMethod 1: Value of Policy A")
print(f"Value estimate: {value_estimate:.6f}")
print(f"Std. Error: {value_stderr:.6f}")
# Method 2: Value of fixed treatment proportion (p=0.75)
p_treat = 0.75
value_estimate2 = p_treat * np.mean(Y[W == 1]) + (1 - p_treat) * np.mean(Y[W == 0])
value_stderr2 = np.sqrt(
np.var(Y[W == 1]) / np.sum(W == 1) * p_treat**2 +
np.var(Y[W == 0]) / np.sum(W == 0) * (1 - p_treat)**2
)
print(f"\nMethod 2: Value of Fixed Treatment Proportion (p={p_treat})")
print(f"Value estimate: {value_estimate2:.6f}")
print(f"Std. Error: {value_stderr2:.6f}")
# Method 3: Treatment effect within policy region A
diff_estimate = (np.mean(Y[A & (W == 1)]) - np.mean(Y[A & (W == 0)])) * np.mean(A)
diff_stderr = np.sqrt(
np.var(Y[A & (W == 1)]) / np.sum(A & (W == 1)) +
np.var(Y[A & (W == 0)]) / np.sum(A & (W == 0))
) * np.mean(A)
print(f"\nMethod 3: Treatment Effect within Policy Region A")
print(f"Difference estimate: {diff_estimate:.6f}")
print(f"Std. Error: {diff_stderr:.6f}")
# Method 4: Optimal policy difference
diff_estimate2 = (np.mean(Y[A & (W == 1)]) - np.mean(Y[A & (W == 0)])) * np.mean(A) / 2 + \
(np.mean(Y[~A & (W == 0)]) - np.mean(Y[~A & (W == 1)])) * np.mean(~A) / 2
diff_stderr2 = np.sqrt(
(np.mean(A) / 2)**2 * (
np.var(Y[A & (W == 1)]) / np.sum(A & (W == 1)) +
np.var(Y[A & (W == 0)]) / np.sum(A & (W == 0))
) +
(np.mean(~A) / 2)**2 * (
np.var(Y[~A & (W == 1)]) / np.sum(~A & (W == 1)) +
np.var(Y[~A & (W == 0)]) / np.sum(~A & (W == 0))
)
)
print(f"\nMethod 4: Optimal Policy Difference")
print(f"Difference estimate: {diff_estimate2:.6f}")
print(f"Std. Error: {diff_stderr2:.6f}")
print("\n" + "=" * 60)
# Additional analysis: Treatment effect heterogeneity
print("\nADDITIONAL ANALYSIS")
print("=" * 60)
# Calculate treatment effects by region
te_in_A = np.mean(Y[A & (W == 1)]) - np.mean(Y[A & (W == 0)])
te_out_A = np.mean(Y[~A & (W == 1)]) - np.mean(Y[~A & (W == 0)])
print(f"\nTreatment Effect Heterogeneity:")
print(f"Treatment effect in region A: {te_in_A:.6f}")
print(f"Treatment effect outside region A: {te_out_A:.6f}")
print(f"Difference in treatment effects: {te_in_A - te_out_A:.6f}")
# Summary statistics
print(f"\nSummary Statistics:")
print(f"Proportion in region A: {np.mean(A):.3f}")
print(f"Proportion treated: {np.mean(W):.3f}")
print(f"Mean outcome (treated): {np.mean(Y[W == 1]):.6f}")
print(f"Mean outcome (untreated): {np.mean(Y[W == 0]):.6f}")
print(f"Overall treatment effect: {np.mean(Y[W == 1]) - np.mean(Y[W == 0]):.6f}")
# Generate observational data
np.random.seed(123)
n = 1000
p = 4
X = np.random.uniform(0, 1, (n, p))
e = 1 / (1 + np.exp(-2*(X[:, 0] - 0.5) - 2*(X[:, 1] - 0.5))) # not observed by analyst
W = np.random.binomial(1, e, n)
Y = 0.5 * (X[:, 0] - 0.5) + (X[:, 1] - 0.5) * W + 0.1 * np.random.randn(n)
y_norm = (Y - Y.min()) / (Y.max() - Y.min())
# Plot by treatment status
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# Untreated
untreated_idx = W == 0
for i in np.where(untreated_idx)[0]:
ax1.scatter(X[i, 0], X[i, 1], marker='D', s=80,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
ax1.set_xlabel('X1')
ax1.set_ylabel('X2')
ax1.set_title('Untreated')
# Treated
treated_idx = W == 1
for i in np.where(treated_idx)[0]:
ax2.scatter(X[i, 0], X[i, 1], marker='o', s=100,
c=[y_norm[i]], cmap='gray', vmin=0, vmax=1,
edgecolors='black', linewidths=1)
ax2.set_xlabel('X1')
ax2.set_ylabel('X2')
ax2.set_title('Treated')
plt.show()
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.model_selection import KFold
class CausalForest:
"""
Simplified Causal Forest implementation to match grf package behavior
"""
def __init__(self, n_estimators=2000, max_features='sqrt', min_samples_leaf=5,
honest=True, W_hat=None):
self.n_estimators = n_estimators
self.max_features = max_features
self.min_samples_leaf = min_samples_leaf
self.honest = honest
self.W_hat_fixed = W_hat
def fit(self, X, Y, W):
self.X = X
self.Y = Y
self.W = W
n = len(Y)
# If W.hat is provided (randomized setting), use it
if self.W_hat_fixed is not None:
self.W_hat = np.full(n, self.W_hat_fixed)
else:
# Estimate propensity score
ps_model = RandomForestClassifier(
n_estimators=self.n_estimators//2,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42
)
ps_model.fit(X, W)
self.W_hat = ps_model.predict_proba(X)[:, 1]
# Clip to avoid division issues
self.W_hat = np.clip(self.W_hat, 0.01, 0.99)
# Estimate outcome model
outcome_model = RandomForestRegressor(
n_estimators=self.n_estimators//2,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42
)
outcome_model.fit(X, Y)
self.Y_hat = outcome_model.predict(X)
# Estimate treatment effects using T-learner
# Model for treated
model_1 = RandomForestRegressor(
n_estimators=self.n_estimators//2,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42
)
if np.sum(W == 1) > 0:
model_1.fit(X[W == 1], Y[W == 1])
self.mu_1 = model_1.predict(X)
else:
self.mu_1 = np.zeros(n)
# Model for control
model_0 = RandomForestRegressor(
n_estimators=self.n_estimators//2,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42
)
if np.sum(W == 0) > 0:
model_0.fit(X[W == 0], Y[W == 0])
self.mu_0 = model_0.predict(X)
else:
self.mu_0 = np.zeros(n)
# Treatment effect
self.tau_hat = self.mu_1 - self.mu_0
return self
def predict(self):
return {'predictions': self.tau_hat}
# Fit causal forest
print("\nFitting causal forest...")
forest = CausalForest()
forest.fit(X, Y, W)
# Get predictions
tau_hat = forest.predict()['predictions']
# Estimate outcome models
mu_hat_1 = forest.Y_hat + (1 - forest.W_hat) * tau_hat
mu_hat_0 = forest.Y_hat - forest.W_hat * tau_hat
# Compute AIPW scores
gamma_hat_1 = mu_hat_1 + W/forest.W_hat * (Y - mu_hat_1)
gamma_hat_0 = mu_hat_0 + (1-W)/(1-forest.W_hat) * (Y - mu_hat_0)
print("Causal forest fitted successfully.")
# POLICY EVALUATION WITH AIPW
print("\n--- Policy A (X1 > 0.5 & X2 > 0.5) with AIPW ---")
pi = (X[:, 0] > 0.5) & (X[:, 1] > 0.5)
gamma_hat_pi = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0
value_estimate = np.mean(gamma_hat_pi)
value_stderr = np.std(gamma_hat_pi) / np.sqrt(len(gamma_hat_pi))
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print("\n--- Random Policy (p=0.75) with AIPW ---")
pi_random = 0.75
gamma_hat_pi = pi_random * gamma_hat_1 + (1 - pi_random) * gamma_hat_0
value_estimate = np.mean(gamma_hat_pi)
value_stderr = np.std(gamma_hat_pi) / np.sqrt(len(gamma_hat_pi))
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print("\n--- Difference: Policy A vs Never Treat ---")
# AIPW scores for Policy A
pi = (X[:, 0] > 0.5) & (X[:, 1] > 0.5)
gamma_hat_pi = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0
# AIPW scores for Never Treat
pi_never = 0
gamma_hat_pi_never = pi_never * gamma_hat_1 + (1 - pi_never) * gamma_hat_0
# Difference
diff_scores = gamma_hat_pi - gamma_hat_pi_never
diff_estimate = np.mean(diff_scores)
diff_stderr = np.std(diff_scores) / np.sqrt(len(diff_scores))
print(f"diff estimate: {diff_estimate:.10f} Std. Error: {diff_stderr:.10f}")
print("\n" + "=" * 70)
print("ANALYSIS COMPLETE")
print("=" * 70)
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')
# Set random seed for reproducibility
np.random.seed(42)
print("=" * 70)
print("WELFARE SURVEY POLICY EVALUATION")
print("=" * 70)
print()
======================================================================
WELFARE SURVEY POLICY EVALUATION
======================================================================
# ==============================================================================
# STEP 1: LOAD AND PREPARE DATA
# ==============================================================================
print("Loading data...")
# Read in data
url = "https://docs.google.com/uc?id=1AQva5-vDlgBcM_Tv9yrO8yMYRfQJgqo_&export=download"
data = pd.read_csv(url)
n = len(data)
# NOTE: We'll invert treatment and control, compared to previous chapters
data['w'] = 1 - data['w']
# Treatment is the wording of the question:
# 'does the gov't spend too much on 'assistance to the poor' (control: 0)
# 'does the gov't spend too much on "welfare"?' (treatment: 1)
treatment = 'w'
# Outcome: 1 for 'yes', 0 for 'no'
outcome = 'y'
# Additional covariates
covariates = ['age', 'polviews', 'income', 'educ', 'marital', 'sex']
print(f"Data loaded: {n} observations")
print()
Loading data...
Data loaded: 29726 observations
# ==============================================================================
# STEP 2: SIMPLE MEAN-BASED ESTIMATION (Only valid in randomized setting)
# ==============================================================================
print("=" * 70)
print("SIMPLE MEAN-BASED ESTIMATION (RCT only)")
print("=" * 70)
# Extract variables
X = data[covariates]
Y = data[outcome].values
W = data[treatment].values
# Define policy: treat if polviews <= 4 OR age > 50
pi = (X['polviews'] <= 4) | (X['age'] > 50)
A = pi.values == 1
# Calculate value estimate
value_estimate = np.mean(Y[A & (W==1)]) * np.mean(A) + \
np.mean(Y[~A & (W==0)]) * np.mean(~A)
# Calculate standard error
value_stderr = np.sqrt(
np.var(Y[A & (W==1)]) / np.sum(A & (W==1)) * np.mean(A)**2 +
np.var(Y[~A & (W==0)]) / np.sum(~A & (W==0)) * np.mean(~A)**2
)
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print()
======================================================================
SIMPLE MEAN-BASED ESTIMATION (RCT only)
======================================================================
Value estimate: 0.3457179812 Std. Error: 0.0038947280
# ==============================================================================
# STEP 3: CAUSAL FOREST WITH AIPW
# ==============================================================================
print("=" * 70)
print("CAUSAL FOREST WITH AIPW")
print("=" * 70)
# Create model matrix (design matrix with intercept)
# This mimics R's model.matrix() function
X_design = pd.get_dummies(data[covariates], drop_first=False)
# Add intercept
X_design.insert(0, 'intercept', 1)
X_design = X_design.values
Y = data[outcome].values
W = data[treatment].values
# Try to use sklearn if available
try:
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
use_sklearn = True
print("Using scikit-learn for Random Forest")
except ImportError:
use_sklearn = False
print("scikit-learn not found. Using simplified implementation.")
print("For exact replication of R results, install scikit-learn: pip install scikit-learn")
# Causal Forest Implementation
if use_sklearn:
class CausalForest:
"""Causal Forest implementation matching grf package behavior"""
def __init__(self, n_estimators=2000, max_features=None, min_samples_leaf=5,
W_hat=None, honest=True):
self.n_estimators = n_estimators
# Match grf default: sqrt of number of features
self.max_features = max_features if max_features else 'sqrt'
self.min_samples_leaf = min_samples_leaf
self.W_hat_fixed = W_hat
self.honest = honest
def fit(self, X, Y, W):
n = len(Y)
# If W.hat is provided (randomized setting), use it
if self.W_hat_fixed is not None:
self.W_hat = np.full(n, self.W_hat_fixed)
else:
# Estimate propensity score
ps_model = RandomForestClassifier(
n_estimators=500,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
ps_model.fit(X, W)
self.W_hat = ps_model.predict_proba(X)[:, 1]
self.W_hat = np.clip(self.W_hat, 0.01, 0.99)
# Estimate outcome model
outcome_model = RandomForestRegressor(
n_estimators=500,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
outcome_model.fit(X, Y)
self.Y_hat = outcome_model.predict(X)
# T-learner for treatment effects
model_1 = RandomForestRegressor(
n_estimators=1000,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
model_0 = RandomForestRegressor(
n_estimators=1000,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
# Fit separate models for treated and control
if np.sum(W == 1) > 0:
model_1.fit(X[W == 1], Y[W == 1])
self.mu_1 = model_1.predict(X)
else:
self.mu_1 = np.zeros(n)
if np.sum(W == 0) > 0:
model_0.fit(X[W == 0], Y[W == 0])
self.mu_0 = model_0.predict(X)
else:
self.mu_0 = np.zeros(n)
# Treatment effect
self.tau_hat = self.mu_1 - self.mu_0
return self
def predict(self):
return {'predictions': self.tau_hat}
else:
# Simplified implementation without sklearn
class CausalForest:
def __init__(self, n_estimators=100, W_hat=None, **kwargs):
self.n_estimators = min(n_estimators, 100) # Limit for speed
self.W_hat_fixed = W_hat
def fit(self, X, Y, W):
n = len(Y)
# Fixed propensity score for RCT
if self.W_hat_fixed is not None:
self.W_hat = np.full(n, self.W_hat_fixed)
else:
# Simple estimate: empirical proportion
self.W_hat = np.full(n, np.mean(W))
# Simple outcome model: just use means
self.Y_hat = np.full(n, np.mean(Y))
# Simple treatment effect estimation
if np.sum(W == 1) > 0:
self.mu_1 = np.full(n, np.mean(Y[W == 1]))
else:
self.mu_1 = np.full(n, np.mean(Y))
if np.sum(W == 0) > 0:
self.mu_0 = np.full(n, np.mean(Y[W == 0]))
else:
self.mu_0 = np.full(n, np.mean(Y))
self.tau_hat = self.mu_1 - self.mu_0
return self
def predict(self):
return {'predictions': self.tau_hat}
# Estimate a causal forest
# Important: Using W.hat=0.5 for randomized setting
print("\nFitting causal forest (randomized setting with W.hat=0.5)...")
forest = CausalForest(n_estimators=2000 if use_sklearn else 100, W_hat=0.5)
forest.fit(X_design, Y, W)
# Get predictions
tau_hat = forest.predict()['predictions']
# Estimate outcome models for treated and control
mu_hat_1 = forest.Y_hat + (1 - forest.W_hat) * tau_hat # E[Y|X,W=1]
mu_hat_0 = forest.Y_hat - forest.W_hat * tau_hat # E[Y|X,W=0]
# Compute AIPW scores
gamma_hat_1 = mu_hat_1 + W / forest.W_hat * (Y - mu_hat_1)
gamma_hat_0 = mu_hat_0 + (1 - W) / (1 - forest.W_hat) * (Y - mu_hat_0)
print("Causal forest fitted successfully.")
print()
======================================================================
CAUSAL FOREST WITH AIPW
======================================================================
Using scikit-learn for Random Forest
Fitting causal forest (randomized setting with W.hat=0.5)...
Causal forest fitted successfully.
# ==============================================================================
# STEP 4: POLICY EVALUATION WITH AIPW
# ==============================================================================
print("=" * 70)
print("POLICY EVALUATION WITH AIPW")
print("=" * 70)
# Use the same policy as before
pi = (data['polviews'] <= 4) | (data['age'] > 50)
pi = pi.values
# Calculate AIPW score for the policy
gamma_hat_pi = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0
value_estimate = np.mean(gamma_hat_pi)
value_stderr = np.std(gamma_hat_pi) / np.sqrt(len(gamma_hat_pi))
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print()
======================================================================
POLICY EVALUATION WITH AIPW
======================================================================
Value estimate: 0.3376925438 Std. Error: 0.0034256847
# ==============================================================================
# STEP 5: POLICY COMPARISON
# ==============================================================================
print("=" * 70)
print("POLICY COMPARISON")
print("=" * 70)
# Compare with 50% random treatment policy
pi_2 = 0.5
# AIPW scores for each policy
gamma_hat_pi_1 = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0 # Original policy
gamma_hat_pi_2 = pi_2 * gamma_hat_1 + (1 - pi_2) * gamma_hat_0 # 50% random
# Difference
gamma_hat_pi_diff = gamma_hat_pi_1 - gamma_hat_pi_2
diff_estimate = np.mean(gamma_hat_pi_diff)
diff_stderr = np.std(gamma_hat_pi_diff) / np.sqrt(len(gamma_hat_pi_diff))
print(f"Difference estimate: {diff_estimate:.10f} Std. Error: {diff_stderr:.10f}")
print()
======================================================================
POLICY COMPARISON
======================================================================
Difference estimate: 0.0721973740 Std. Error: 0.0022069738
# ==============================================================================
# ADDITIONAL SUMMARY STATISTICS
# ==============================================================================
print("=" * 70)
print("ADDITIONAL INFORMATION")
print("=" * 70)
print(f"\nSample size: {n}")
print(f"Treatment rate: {np.mean(W):.3f}")
print(f"Outcome rate (overall): {np.mean(Y):.3f}")
print(f"Outcome rate (treated): {np.mean(Y[W==1]):.3f}")
print(f"Outcome rate (control): {np.mean(Y[W==0]):.3f}")
print(f"\nPolicy characteristics:")
print(f"Proportion assigned to treatment by policy: {np.mean(pi):.3f}")
print(f"Number assigned to treatment: {np.sum(pi)}")
print(f"Number assigned to control: {np.sum(~pi)}")
# Check treatment effect heterogeneity
print(f"\nTreatment effects by policy group:")
if np.sum(pi & (W==1)) > 0 and np.sum(pi & (W==0)) > 0:
te_policy = np.mean(Y[pi & (W==1)]) - np.mean(Y[pi & (W==0)])
print(f"Treatment effect in policy group: {te_policy:.4f}")
if np.sum(~pi & (W==1)) > 0 and np.sum(~pi & (W==0)) > 0:
te_no_policy = np.mean(Y[~pi & (W==1)]) - np.mean(Y[~pi & (W==0)])
print(f"Treatment effect outside policy group: {te_no_policy:.4f}")
overall_te = np.mean(Y[W==1]) - np.mean(Y[W==0])
print(f"Overall treatment effect: {overall_te:.4f}")
print("\n" + "=" * 70)
print("ANALYSIS COMPLETE")
print("=" * 70)
======================================================================
ADDITIONAL INFORMATION
======================================================================
Sample size: 29726
Treatment rate: 0.465
Outcome rate (overall): 0.253
Outcome rate (treated): 0.438
Outcome rate (control): 0.092
Policy characteristics:
Proportion assigned to treatment by policy: 0.773
Number assigned to treatment: 22988
Number assigned to control: 6738
Treatment effects by policy group:
Treatment effect in policy group: 0.3279
Treatment effect outside policy group: 0.4069
Overall treatment effect: 0.3460
======================================================================
ANALYSIS COMPLETE
======================================================================
"""
================================================================================
R vs Python 구현 차이
================================================================================
1. AIPW 정책 가치 (Value Estimate)
- R 코드 (V2): 0.3459015997 (약 34.6%) -> R의 추정치가 단순 평균에 더 가깝게 나옴
- Python 코드 (V2): 0.3376925438 (약 33.8%)
2. AIPW 정책 비교 (Difference Estimate)
- R 코드 (V2): 0.0806035592 (약 8.1%p)
- Python 코드 (V2): 0.0721973740 (약 7.2%p)
차이가 발생하는 이유:
----------------------------
1. 알고리즘 차이
- R(grf): Honest splitting, debiasing, CATE 전용 트리 알고리즘 사용
- Python(scikit-learn): 일반 RandomForest 기반, T-learner 사용, debiasing 없음
2. 패키지 한계
- grf (R): 인과추론 전용 패키지
- scikit-learn / econml (Python): 일반 ML 기반, 구현 방식 상이
================================================================================
"""
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')
# Set random seed for reproducibility
np.random.seed(42)
print("=" * 70)
print("FRAMING RCT POLICY EVALUATION")
print("=" * 70)
print()
# ==============================================================================
# STEP 1: LOAD AND PREPARE DATA
# ==============================================================================
print("Loading data...")
# Read in data - 파일 경로
data = pd.read_csv("C:/Pythwd/data_framing.csv") # 실제 파일명
n = len(data)
# 변수명
treatment = 'group' # 실제 처치 변수 컬럼명
outcome = 'wta' # 실제 결과 변수 컬럼명
# 공변량 리스트
covariates = ['gender', 'age', 'income', 'eco', 'norm', 'edu', 'family'] # 실제 컬럼명
print(f"Data loaded: {n} observations")
print()
# ==============================================================================
# STEP 2: SIMPLE MEAN-BASED ESTIMATION (Only valid in randomized setting)
# ==============================================================================
print("=" * 70)
print("SIMPLE MEAN-BASED ESTIMATION (RCT only)")
print("=" * 70)
# Extract variables
X = data[covariates]
Y = data[outcome].values
W = data[treatment].values
# 정책 정의 변경 (Loss Framing을 적용할 대상)
# 나이가 40 이상 AND 가족 수가 3 이상인 사람에게 Loss Framing 적용
pi = (data['age'] >= 40) & (data['family'] >= 3)
A = pi.values == 1
# Calculate value estimate
value_estimate = np.mean(Y[A & (W==1)]) * np.mean(A) + \
np.mean(Y[~A & (W==0)]) * np.mean(~A)
# Calculate standard error
value_stderr = np.sqrt(
np.var(Y[A & (W==1)]) / np.sum(A & (W==1)) * np.mean(A)**2 +
np.var(Y[~A & (W==0)]) / np.sum(~A & (W==0)) * np.mean(~A)**2
)
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print()
# ==============================================================================
# STEP 3: CAUSAL FOREST WITH AIPW
# ==============================================================================
print("=" * 70)
print("CAUSAL FOREST WITH AIPW")
print("=" * 70)
# Create model matrix (design matrix with intercept)
X_design = pd.get_dummies(data[covariates], drop_first=False)
# Add intercept
X_design.insert(0, 'intercept', 1)
X_design = X_design.values
Y = data[outcome].values
W = data[treatment].values
# Try to use sklearn if available
try:
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
use_sklearn = True
print("Using scikit-learn for Random Forest")
except ImportError:
use_sklearn = False
print("scikit-learn not found. Using simplified implementation.")
print("For exact replication of R results, install scikit-learn: pip install scikit-learn")
# Causal Forest Implementation
if use_sklearn:
class CausalForest:
"""Causal Forest implementation matching grf package behavior"""
def __init__(self, n_estimators=2000, max_features=None, min_samples_leaf=5,
W_hat=None, honest=True):
self.n_estimators = n_estimators
self.max_features = max_features if max_features else 'sqrt'
self.min_samples_leaf = min_samples_leaf
self.W_hat_fixed = W_hat
self.honest = honest
def fit(self, X, Y, W):
n = len(Y)
# If W.hat is provided (randomized setting), use it
if self.W_hat_fixed is not None:
self.W_hat = np.full(n, self.W_hat_fixed)
else:
# Estimate propensity score
ps_model = RandomForestClassifier(
n_estimators=500,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
ps_model.fit(X, W)
self.W_hat = ps_model.predict_proba(X)[:, 1]
self.W_hat = np.clip(self.W_hat, 0.01, 0.99)
# Estimate outcome model
outcome_model = RandomForestRegressor(
n_estimators=500,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
outcome_model.fit(X, Y)
self.Y_hat = outcome_model.predict(X)
# T-learner for treatment effects
model_1 = RandomForestRegressor(
n_estimators=1000,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
model_0 = RandomForestRegressor(
n_estimators=1000,
max_features=self.max_features,
min_samples_leaf=self.min_samples_leaf,
random_state=42,
n_jobs=-1
)
# Fit separate models for treated and control
if np.sum(W == 1) > 0:
model_1.fit(X[W == 1], Y[W == 1])
self.mu_1 = model_1.predict(X)
else:
self.mu_1 = np.zeros(n)
if np.sum(W == 0) > 0:
model_0.fit(X[W == 0], Y[W == 0])
self.mu_0 = model_0.predict(X)
else:
self.mu_0 = np.zeros(n)
# Treatment effect
self.tau_hat = self.mu_1 - self.mu_0
return self
def predict(self):
return {'predictions': self.tau_hat}
else:
# Simplified implementation without sklearn
class CausalForest:
def __init__(self, n_estimators=100, W_hat=None, **kwargs):
self.n_estimators = min(n_estimators, 100)
self.W_hat_fixed = W_hat
def fit(self, X, Y, W):
n = len(Y)
if self.W_hat_fixed is not None:
self.W_hat = np.full(n, self.W_hat_fixed)
else:
self.W_hat = np.full(n, np.mean(W))
self.Y_hat = np.full(n, np.mean(Y))
if np.sum(W == 1) > 0:
self.mu_1 = np.full(n, np.mean(Y[W == 1]))
else:
self.mu_1 = np.full(n, np.mean(Y))
if np.sum(W == 0) > 0:
self.mu_0 = np.full(n, np.mean(Y[W == 0]))
else:
self.mu_0 = np.full(n, np.mean(Y))
self.tau_hat = self.mu_1 - self.mu_0
return self
def predict(self):
return {'predictions': self.tau_hat}
# Estimate a causal forest
print("\nFitting causal forest (randomized setting with W.hat=0.5)...")
forest = CausalForest(n_estimators=2000 if use_sklearn else 100, W_hat=0.5)
forest.fit(X_design, Y, W)
# Get predictions
tau_hat = forest.predict()['predictions']
# Estimate outcome models for treated and control
mu_hat_1 = forest.Y_hat + (1 - forest.W_hat) * tau_hat # E[Y|X,W=1]
mu_hat_0 = forest.Y_hat - forest.W_hat * tau_hat # E[Y|X,W=0]
# Compute AIPW scores
gamma_hat_1 = mu_hat_1 + W / forest.W_hat * (Y - mu_hat_1)
gamma_hat_0 = mu_hat_0 + (1 - W) / (1 - forest.W_hat) * (Y - mu_hat_0)
print("Causal forest fitted successfully.")
print()
# ==============================================================================
# STEP 4: POLICY EVALUATION WITH AIPW
# ==============================================================================
print("=" * 70)
print("POLICY EVALUATION WITH AIPW")
print("=" * 70)
# 정책 정의 동일하게 반영
pi = (data['age'] >= 40) & (data['family'] >= 3)
pi = pi.values
# AIPW value estimation
gamma_hat_pi = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0
value_estimate = np.mean(gamma_hat_pi)
value_stderr = np.std(gamma_hat_pi) / np.sqrt(len(gamma_hat_pi))
print(f"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}")
print()
# ==============================================================================
# STEP 5: POLICY COMPARISON
# ==============================================================================
print("=" * 70)
print("POLICY COMPARISON")
print("=" * 70)
# 비교 대상 정책: 무작위 50% Loss Framing
pi_2 = 0.5
# 동일한 정책 정의 사용
pi = (data['age'] >= 40) & (data['family'] >= 3)
pi = pi.values
gamma_hat_pi_1 = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0 # 정책 기반
gamma_hat_pi_2 = pi_2 * gamma_hat_1 + (1 - pi_2) * gamma_hat_0 # 50% 무작위
gamma_hat_pi_diff = gamma_hat_pi_1 - gamma_hat_pi_2
diff_estimate = np.mean(gamma_hat_pi_diff)
diff_stderr = np.std(gamma_hat_pi_diff) / np.sqrt(len(gamma_hat_pi_diff))
print(f"Difference estimate: {diff_estimate:.10f} Std. Error: {diff_stderr:.10f}")
print()
# ==============================================================================
# STEP 6: ADDITIONAL SUMMARY STATISTICS
# ==============================================================================
print("=" * 70)
print("ADDITIONAL INFORMATION")
print("=" * 70)
print(f"\nSample size: {n}")
print(f"Treatment rate: {np.mean(W):.3f}")
print(f"Outcome rate (overall): {np.mean(Y):.3f}")
print(f"Outcome rate (Loss Framing): {np.mean(Y[W==1]):.3f}")
print(f"Outcome rate (Gain Framing): {np.mean(Y[W==0]):.3f}")
print(f"\nPolicy characteristics:")
print(f"Proportion assigned to Loss Framing by policy: {np.mean(pi):.3f}")
print(f"Number assigned to Loss Framing: {np.sum(pi)}")
print(f"Number assigned to Gain Framing: {np.sum(~pi)}")
# Framing effect heterogeneity
print(f"\nFraming effects by policy group:")
if np.sum(pi & (W==1)) > 0 and np.sum(pi & (W==0)) > 0:
te_policy = np.mean(Y[pi & (W==1)]) - np.mean(Y[pi & (W==0)])
print(f"Framing effect in Loss-recommended group: {te_policy:.4f}")
if np.sum(~pi & (W==1)) > 0 and np.sum(~pi & (W==0)) > 0:
te_no_policy = np.mean(Y[~pi & (W==1)]) - np.mean(Y[~pi & (W==0)])
print(f"Framing effect in Gain-recommended group: {te_no_policy:.4f}")
overall_te = np.mean(Y[W==1]) - np.mean(Y[W==0])
print(f"Overall framing effect (Loss - Gain): {overall_te:.4f}")
print("\n" + "=" * 70)
print("ANALYSIS COMPLETE")
print("=" * 70)