Skip to content

policy aipw v1 #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

policy aipw v1 #23

Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
eb7f133
policy aipw v1
AIgoGracia Sep 21, 2025
1a6025b
policy aipw v2
AIgoGracia Sep 27, 2025
8f7fb9a
policy aipw v3
AIgoGracia Oct 4, 2025
9f69e18
case study v1
AIgoGracia Oct 12, 2025
17ec9ed
case study v2
AIgoGracia Oct 18, 2025
c4fa2ef
case study v3
AIgoGracia Oct 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
384 changes: 384 additions & 0 deletions book/cate_and_policy/policy_learning.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -965,6 +965,390 @@
"print(\"ANALYSIS COMPLETE\")\n",
"print(\"=\" * 70)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"================================================================================\n",
"R vs Python 구현 차이\n",
"================================================================================\n",
"\n",
"1. AIPW 정책 가치 (Value Estimate)\n",
"- R 코드 (V2): 0.3459015997 (약 34.6%) -> R의 추정치가 단순 평균에 더 가깝게 나옴\n",
"- Python 코드 (V2): 0.3376925438 (약 33.8%)\n",
"\n",
"2. AIPW 정책 비교 (Difference Estimate)\n",
"- R 코드 (V2): 0.0806035592 (약 8.1%p)\n",
"- Python 코드 (V2): 0.0721973740 (약 7.2%p)\n",
"\n",
"차이가 발생하는 이유:\n",
"----------------------------\n",
"1. 알고리즘 차이\n",
" - R(grf): Honest splitting, debiasing, CATE 전용 트리 알고리즘 사용\n",
" - Python(scikit-learn): 일반 RandomForest 기반, T-learner 사용, debiasing 없음\n",
"\n",
"2. 패키지 한계\n",
" - grf (R): 인과추론 전용 패키지\n",
" - scikit-learn / econml (Python): 일반 ML 기반, 구현 방식 상이\n",
"\n",
"================================================================================\n",
"\"\"\""
Comment on lines +975 to +999
Copy link
Contributor

@Funbucket Funbucket Oct 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

초반에 연구 주제와 데이터셋에 대한 설명이 함께 들어가면 전체 흐름을 이해하는 데 더 도움이 될 것 같습니다!

]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"# Set random seed for reproducibility\n",
"np.random.seed(42)\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"FRAMING RCT POLICY EVALUATION\")\n",
"print(\"=\" * 70)\n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 1: LOAD AND PREPARE DATA\n",
"# ==============================================================================\n",
"\n",
"print(\"Loading data...\")\n",
"# Read in data - 파일 경로\n",
"data = pd.read_csv(\"C:/Pythwd/data_framing.csv\") # 실제 파일명\n",
"n = len(data)\n",
"\n",
"# 변수명\n",
"treatment = 'group' # 실제 처치 변수 컬럼명\n",
"outcome = 'wta' # 실제 결과 변수 컬럼명\n",
"\n",
"# 공변량 리스트\n",
"covariates = ['gender', 'age', 'income', 'eco', 'norm', 'edu', 'family'] # 실제 컬럼명\n",
"\n",
"print(f\"Data loaded: {n} observations\")\n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 2: SIMPLE MEAN-BASED ESTIMATION (Only valid in randomized setting)\n",
"# ==============================================================================\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"SIMPLE MEAN-BASED ESTIMATION (RCT only)\")\n",
"print(\"=\" * 70)\n",
"\n",
"# Extract variables\n",
"X = data[covariates]\n",
"Y = data[outcome].values\n",
"W = data[treatment].values\n",
"\n",
"# 정책 정의 변경 (Loss Framing을 적용할 대상)\n",
"# 나이가 40 이상 AND 가족 수가 3 이상인 사람에게 Loss Framing 적용\n",
"pi = (data['age'] >= 40) & (data['family'] >= 3)\n",
"A = pi.values == 1\n",
"\n",
"# Calculate value estimate\n",
"value_estimate = np.mean(Y[A & (W==1)]) * np.mean(A) + \\\n",
" np.mean(Y[~A & (W==0)]) * np.mean(~A)\n",
"\n",
"# Calculate standard error\n",
"value_stderr = np.sqrt(\n",
" np.var(Y[A & (W==1)]) / np.sum(A & (W==1)) * np.mean(A)**2 + \n",
" np.var(Y[~A & (W==0)]) / np.sum(~A & (W==0)) * np.mean(~A)**2\n",
")\n",
"\n",
"print(f\"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}\")\n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 3: CAUSAL FOREST WITH AIPW\n",
"# ==============================================================================\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"CAUSAL FOREST WITH AIPW\")\n",
"print(\"=\" * 70)\n",
"\n",
"# Create model matrix (design matrix with intercept)\n",
"X_design = pd.get_dummies(data[covariates], drop_first=False)\n",
"# Add intercept\n",
"X_design.insert(0, 'intercept', 1)\n",
"X_design = X_design.values\n",
"\n",
"Y = data[outcome].values\n",
"W = data[treatment].values\n",
"\n",
"# Try to use sklearn if available\n",
"try:\n",
" from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
" use_sklearn = True\n",
" print(\"Using scikit-learn for Random Forest\")\n",
"except ImportError:\n",
" use_sklearn = False\n",
" print(\"scikit-learn not found. Using simplified implementation.\")\n",
" print(\"For exact replication of R results, install scikit-learn: pip install scikit-learn\")\n",
"\n",
"# Causal Forest Implementation\n",
"if use_sklearn:\n",
" class CausalForest:\n",
" \"\"\"Causal Forest implementation matching grf package behavior\"\"\"\n",
" \n",
" def __init__(self, n_estimators=2000, max_features=None, min_samples_leaf=5, \n",
" W_hat=None, honest=True):\n",
" self.n_estimators = n_estimators\n",
" self.max_features = max_features if max_features else 'sqrt'\n",
" self.min_samples_leaf = min_samples_leaf\n",
" self.W_hat_fixed = W_hat\n",
" self.honest = honest\n",
" \n",
" def fit(self, X, Y, W):\n",
" n = len(Y)\n",
" \n",
" # If W.hat is provided (randomized setting), use it\n",
" if self.W_hat_fixed is not None:\n",
" self.W_hat = np.full(n, self.W_hat_fixed)\n",
" else:\n",
" # Estimate propensity score\n",
" ps_model = RandomForestClassifier(\n",
" n_estimators=500,\n",
" max_features=self.max_features,\n",
" min_samples_leaf=self.min_samples_leaf,\n",
" random_state=42,\n",
" n_jobs=-1\n",
" )\n",
" ps_model.fit(X, W)\n",
" self.W_hat = ps_model.predict_proba(X)[:, 1]\n",
" self.W_hat = np.clip(self.W_hat, 0.01, 0.99)\n",
" \n",
" # Estimate outcome model\n",
" outcome_model = RandomForestRegressor(\n",
" n_estimators=500,\n",
" max_features=self.max_features,\n",
" min_samples_leaf=self.min_samples_leaf,\n",
" random_state=42,\n",
" n_jobs=-1\n",
" )\n",
" outcome_model.fit(X, Y)\n",
" self.Y_hat = outcome_model.predict(X)\n",
" \n",
" # T-learner for treatment effects\n",
" model_1 = RandomForestRegressor(\n",
" n_estimators=1000,\n",
" max_features=self.max_features,\n",
" min_samples_leaf=self.min_samples_leaf,\n",
" random_state=42,\n",
" n_jobs=-1\n",
" )\n",
" model_0 = RandomForestRegressor(\n",
" n_estimators=1000,\n",
" max_features=self.max_features,\n",
" min_samples_leaf=self.min_samples_leaf,\n",
" random_state=42,\n",
" n_jobs=-1\n",
" )\n",
" \n",
" # Fit separate models for treated and control\n",
" if np.sum(W == 1) > 0:\n",
" model_1.fit(X[W == 1], Y[W == 1])\n",
" self.mu_1 = model_1.predict(X)\n",
" else:\n",
" self.mu_1 = np.zeros(n)\n",
" \n",
" if np.sum(W == 0) > 0:\n",
" model_0.fit(X[W == 0], Y[W == 0])\n",
" self.mu_0 = model_0.predict(X)\n",
" else:\n",
" self.mu_0 = np.zeros(n)\n",
" \n",
" # Treatment effect\n",
" self.tau_hat = self.mu_1 - self.mu_0\n",
" \n",
" return self\n",
" \n",
" def predict(self):\n",
" return {'predictions': self.tau_hat}\n",
"else:\n",
" # Simplified implementation without sklearn\n",
" class CausalForest:\n",
" def __init__(self, n_estimators=100, W_hat=None, **kwargs):\n",
" self.n_estimators = min(n_estimators, 100)\n",
" self.W_hat_fixed = W_hat\n",
" \n",
" def fit(self, X, Y, W):\n",
" n = len(Y)\n",
" \n",
" if self.W_hat_fixed is not None:\n",
" self.W_hat = np.full(n, self.W_hat_fixed)\n",
" else:\n",
" self.W_hat = np.full(n, np.mean(W))\n",
" \n",
" self.Y_hat = np.full(n, np.mean(Y))\n",
" \n",
" if np.sum(W == 1) > 0:\n",
" self.mu_1 = np.full(n, np.mean(Y[W == 1]))\n",
" else:\n",
" self.mu_1 = np.full(n, np.mean(Y))\n",
" \n",
" if np.sum(W == 0) > 0:\n",
" self.mu_0 = np.full(n, np.mean(Y[W == 0]))\n",
" else:\n",
" self.mu_0 = np.full(n, np.mean(Y))\n",
" \n",
" self.tau_hat = self.mu_1 - self.mu_0\n",
" \n",
" return self\n",
" \n",
" def predict(self):\n",
" return {'predictions': self.tau_hat}\n",
"\n",
"# Estimate a causal forest\n",
"print(\"\\nFitting causal forest (randomized setting with W.hat=0.5)...\")\n",
"forest = CausalForest(n_estimators=2000 if use_sklearn else 100, W_hat=0.5)\n",
"forest.fit(X_design, Y, W)\n",
"\n",
"# Get predictions\n",
"tau_hat = forest.predict()['predictions']\n",
"\n",
"# Estimate outcome models for treated and control\n",
"mu_hat_1 = forest.Y_hat + (1 - forest.W_hat) * tau_hat # E[Y|X,W=1]\n",
"mu_hat_0 = forest.Y_hat - forest.W_hat * tau_hat # E[Y|X,W=0]\n",
"\n",
"# Compute AIPW scores\n",
"gamma_hat_1 = mu_hat_1 + W / forest.W_hat * (Y - mu_hat_1)\n",
"gamma_hat_0 = mu_hat_0 + (1 - W) / (1 - forest.W_hat) * (Y - mu_hat_0)\n",
"\n",
"print(\"Causal forest fitted successfully.\")\n",
"print()"
Comment on lines +1092 to +1249
Copy link
Contributor

@Funbucket Funbucket Oct 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

이번에 CausalForest를 직접 구현해주셔서 구조를 이해하는 데 큰 도움이 됐습니다! 다만 실무적으로는 econml의 CausalForestDML 같은 기존 라이브러리를 활용하면 코드가 훨씬 간결하고 유지보수에도 용이할 것 같아요.

]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 4: POLICY EVALUATION WITH AIPW\n",
"# ==============================================================================\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"POLICY EVALUATION WITH AIPW\")\n",
"print(\"=\" * 70)\n",
"\n",
"# 정책 정의 동일하게 반영\n",
"pi = (data['age'] >= 40) & (data['family'] >= 3)\n",
"pi = pi.values\n",
"\n",
"# AIPW value estimation\n",
"gamma_hat_pi = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0\n",
"value_estimate = np.mean(gamma_hat_pi)\n",
"value_stderr = np.std(gamma_hat_pi) / np.sqrt(len(gamma_hat_pi))\n",
"\n",
"print(f\"Value estimate: {value_estimate:.10f} Std. Error: {value_stderr:.10f}\")\n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 5: POLICY COMPARISON\n",
"# ==============================================================================\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"POLICY COMPARISON\")\n",
"print(\"=\" * 70)\n",
"\n",
"# 비교 대상 정책: 무작위 50% Loss Framing\n",
"pi_2 = 0.5\n",
"\n",
"# 동일한 정책 정의 사용\n",
"pi = (data['age'] >= 40) & (data['family'] >= 3)\n",
"pi = pi.values\n",
"\n",
"gamma_hat_pi_1 = pi * gamma_hat_1 + (1 - pi) * gamma_hat_0 # 정책 기반\n",
"gamma_hat_pi_2 = pi_2 * gamma_hat_1 + (1 - pi_2) * gamma_hat_0 # 50% 무작위\n",
"\n",
"gamma_hat_pi_diff = gamma_hat_pi_1 - gamma_hat_pi_2\n",
"diff_estimate = np.mean(gamma_hat_pi_diff)\n",
"diff_stderr = np.std(gamma_hat_pi_diff) / np.sqrt(len(gamma_hat_pi_diff))\n",
"\n",
"print(f\"Difference estimate: {diff_estimate:.10f} Std. Error: {diff_stderr:.10f}\")\n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# STEP 6: ADDITIONAL SUMMARY STATISTICS\n",
"# ==============================================================================\n",
"\n",
"print(\"=\" * 70)\n",
"print(\"ADDITIONAL INFORMATION\")\n",
"print(\"=\" * 70)\n",
"\n",
"print(f\"\\nSample size: {n}\")\n",
"print(f\"Treatment rate: {np.mean(W):.3f}\")\n",
"print(f\"Outcome rate (overall): {np.mean(Y):.3f}\")\n",
"print(f\"Outcome rate (Loss Framing): {np.mean(Y[W==1]):.3f}\")\n",
"print(f\"Outcome rate (Gain Framing): {np.mean(Y[W==0]):.3f}\")\n",
"\n",
"print(f\"\\nPolicy characteristics:\")\n",
"print(f\"Proportion assigned to Loss Framing by policy: {np.mean(pi):.3f}\")\n",
"print(f\"Number assigned to Loss Framing: {np.sum(pi)}\")\n",
"print(f\"Number assigned to Gain Framing: {np.sum(~pi)}\")\n",
"\n",
"# Framing effect heterogeneity\n",
"print(f\"\\nFraming effects by policy group:\")\n",
"if np.sum(pi & (W==1)) > 0 and np.sum(pi & (W==0)) > 0:\n",
" te_policy = np.mean(Y[pi & (W==1)]) - np.mean(Y[pi & (W==0)])\n",
" print(f\"Framing effect in Loss-recommended group: {te_policy:.4f}\")\n",
"if np.sum(~pi & (W==1)) > 0 and np.sum(~pi & (W==0)) > 0:\n",
" te_no_policy = np.mean(Y[~pi & (W==1)]) - np.mean(Y[~pi & (W==0)])\n",
" print(f\"Framing effect in Gain-recommended group: {te_no_policy:.4f}\")\n",
"\n",
"overall_te = np.mean(Y[W==1]) - np.mean(Y[W==0])\n",
"print(f\"Overall framing effect (Loss - Gain): {overall_te:.4f}\")\n",
"\n",
"print(\"\\n\" + \"=\" * 70)\n",
"print(\"ANALYSIS COMPLETE\")\n",
"print(\"=\" * 70)"
]
}
],
"metadata": {
Expand Down