Causal inference-এর পরিচয় — DAG, matching আর difference-in-differences

AI-এর ভিত্তি কোর্সের পর্ব ৩৭-এ দেখেছেন, correlation মানেই কারণ নয়, আর একটি তৃতীয় চলক (confounder) দুটি জিনিসকে একসঙ্গে নাড়াতে পারে। এই পর্ব সেখান থেকে এগোয়: কারণের প্রভাব ঠিক কী সংখ্যা, কোন অনুমানে তা ডেটা থেকে মাপা যায়, আর মাপার পদ্ধতিগুলো কীভাবে চলে। চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python-এ NumPy, pandas, statsmodels, scikit-learn আর matplotlib দিয়ে)। প্রথম তিনটি একই namespace ভাগ করে; প্রথমবার scikit-learn নামতে প্রায় ১০ সেকেন্ড লাগে।

১ · সম্পর্ক থেকে কারণ — potential outcomes, ATE আর ATT

এই ভাবনার আনুষ্ঠানিক রূপ potential outcomes কাঠামো। Jerzy Neyman ১৯২৩ সালে কৃষি-পরীক্ষার প্রসঙ্গে এর প্রথম রূপ লেখেন, আর Donald Rubin ১৯৭৪ সাল থেকে একে পর্যবেক্ষণের ডেটায় ছড়িয়ে দেন। প্রতিটি গ্রাহক i-এর দুটি সম্ভাব্য ফল: Y_i(১) মানে coupon পেলে পরের মাসের খরচ, Y_i(০) মানে না পেলে। T_i = ১ হলে coupon পেয়েছেন। ব্যক্তিগত প্রভাব Y_i(১) − Y_i(০), কিন্তু দেখা যায় শুধু একটি ফল। Paul Holland ১৯৮৬ সালে একে বলেন causal inference-এর মৌলিক সমস্যা।

তাই আমরা ব্যক্তি ছেড়ে গড়ের দিকে যাই:

ATE = E[Y(১) − Y(০)]  ·  ATT = E[Y(১) − Y(০) | T = ১]

E[·] মানে গড় (প্রত্যাশিত মান)। ATE সবার গড় প্রভাব: “সবাইকে coupon দিলে বনাম কাউকে না দিলে” পার্থক্য। ATT শুধু যাঁরা আসলে coupon পেয়েছেন তাঁদের গড় প্রভাব: “যাঁদের দিয়েছি, তাঁদের কাজে লেগেছে কি?” প্রশ্ন আলাদা, তাই সংখ্যাও আলাদা হতে পারে। আর সরল তুলনাটি একটি অভেদ (identity) মেনে ভেঙে যায়:

E[Y | T=১] − E[Y | T=০] = ATT + (E[Y(০) | T=১] − E[Y(০) | T=০])

দ্বিতীয় অংশ selection bias: coupon না পেলেও যাঁরা পেয়েছেন আর যাঁরা পাননি, তাঁদের খরচ কতটা আলাদা হতো। অনুগত গ্রাহক বেশি coupon পেলে এই অংশ বড় ধনাত্মক। নিচের stepper চারজন কাল্পনিক গ্রাহকে পুরো হিসাবটি দেখায়।

এই কাঠামোর একটি নীরব শর্ত SUTVA: একজনের coupon অন্যজনের ফল বদলায় না, আর coupon-এর একটিই রূপ। বন্ধুকে coupon শেয়ার করলে শর্তটি ভাঙে; এমন network প্রভাব পর্ব ৬২ (A/B test-এর পরিকল্পনা — power, sample size আর peeking-এর ফাঁদ)-এর বিষয়।

২ · Randomisation কেন কাজ করে, আর বাজারঘরের coupon-এ confounding

মুদ্রা ছুড়ে coupon দিলে T আর (Y(০), Y(১)) পরস্পর-নিরপেক্ষ হয়: অনুগত আর অননুগত, দুই দলেই সমান হারে পড়েন। তখন E[Y(০) | T=১] = E[Y(০) | T=০], selection bias শূন্য, আর সরল পার্থক্যই ATE-র নিরপেক্ষ অনুমান। Ronald Fisher ১৯২০–৩০-এর দশকে কৃষি-পরীক্ষায় এই যুক্তি দাঁড় করান (The Design of Experiments, ১৯৩৫)। এ কারণেই randomised controlled trial বা A/B test কারণ-প্রশ্নের সোনার মান; পরিকল্পনা আর বিশ্লেষণ পর্ব ৬২ আর পর্ব ৬৩ (A/B test-এর বিশ্লেষণ — t-test, chi-square, CUPED আর ভুল ব্যাখ্যা)-এ।

কিন্তু বাজারঘর গত মাসে coupon দিয়েছে marketing দলের পছন্দে, এলোমেলোভাবে নয়। নিচের seed দেওয়া ডেটায় অর্ডার বেশি (অনুগত) হলে coupon পাওয়ার সম্ভাবনা বেশি, আর একই আনুগত্য খরচও বাড়ায়: অর্ডার সংখ্যা তাই একটি confounder চলক। Simulation দুটি সম্ভাব্য ফল (y0, y1) রাখে, তাই আসল ATE জানা; বাস্তবে হাতে থাকে শুধু spend.

def coupon_data(n=2000, seed=64):
    r = np.random.default_rng(seed)
    orders = r.poisson(4, n)                       # গত ছয় মাসের অর্ডার: আনুগত্যের মাপ
    tenure = r.integers(1, 49, n)                  # অ্যাকাউন্টের বয়স, মাস
    dhaka = (r.random(n) < 0.4).astype(int)        # ঢাকা শহরের ঠিকানা
    logit = -2.2 + 0.5 * orders + 0.015 * tenure - 0.4 * dhaka
    coupon = (r.random(n) < 1 / (1 + np.exp(-logit))).astype(int)   # অনুগতরা বেশি পান
    y0 = 250 + 70 * orders + 4 * orders ** 2 + 3 * tenure + 80 * dhaka + r.normal(0, 120, n)
    y1 = y0 + 40 + 15 * orders                     # coupon-এর আসল প্রভাব: অনুগতদের বেশি
    spend = np.where(coupon == 1, y1, y0).round()  # দেখা যায় শুধু একটি
    return pd.DataFrame({"orders": orders, "tenure": tenure, "dhaka": dhaka,
                         "coupon": coupon, "spend": spend, "y0": y0, "y1": y1})

নিচের ব্লক তিনটি অনুমান পাশাপাশি রাখে: সরল পার্থক্য, regression adjustment (OLS-এ coupon-এর সঙ্গে confounder-গুলো রেখে coupon-এর সহগ), আর একই মানুষদের মধ্যে মুদ্রা ছুড়ে coupon দিলে কী হতো।

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from bazar_causal import coupon_data

df = coupon_data()
true_ate = (df.y1 - df.y0).mean()                   # শুধু simulation জানে
true_att = (df.y1 - df.y0)[df.coupon == 1].mean()
print(df.groupby("coupon")[["orders", "tenure", "dhaka"]].mean().round(2))

naive = df.spend[df.coupon == 1].mean() - df.spend[df.coupon == 0].mean()
reg = smf.ols("spend ~ coupon + orders + tenure + dhaka", data=df).fit()
reg2 = smf.ols("spend ~ coupon + orders + I(orders**2) + tenure + dhaka", data=df).fit()
print(f"সরল পার্থক্য        {naive:7.1f}")
print(f"regression (রৈখিক)  {reg.params['coupon']:7.1f}  SE {reg.bse['coupon']:.1f}")
print(f"regression (+orders²) {reg2.params['coupon']:5.1f}  SE {reg2.bse['coupon']:.1f}")
print(f"আসল ATE {true_ate:.1f} | আসল ATT {true_att:.1f}")

rng = np.random.default_rng(1)                      # একই মানুষ, কিন্তু coupon মুদ্রা ছুড়ে
coin = rng.random(len(df)) < 0.5
spend_r = np.where(coin, df.y1, df.y0)
print(f"randomised: সরল পার্থক্য {spend_r[coin].mean() - spend_r[~coin].mean():.1f}, "
      f"গড় orders {df.orders[coin].mean():.2f} বনাম {df.orders[~coin].mean():.2f}")

        orders  tenure  dhaka
coupon                       
0         3.16   24.22   0.48
1         4.69   25.78   0.36
সরল পার্থক্য          267.6
regression (রৈখিক)     98.7  SE 6.0
regression (+orders²) 102.4  SE 5.8
আসল ATE 98.8 | আসল ATT 110.4
randomised: সরল পার্থক্য 88.8, গড় orders 3.93 বনাম 3.92

Regression adjustment দুটি অনুমানে দাঁড়ায়। (১) Unconfoundedness: মাপা confounder X ধরে রাখলে coupon পাওয়া আর সম্ভাব্য ফল পরস্পর-নিরপেক্ষ, অর্থাৎ X-এর বাইরে আর কোনো লুকানো কারণ নেই। (২) Overlap: X-এর প্রতিটি মানে coupon পাওয়া আর না পাওয়া দুই ধরনের মানুষই আছেন। প্রথমটি ডেটা দিয়ে প্রমাণ করা যায় না; দ্বিতীয়টি যাচাই করা যায় (অংশ ৪)। কোন চলক X-এ রাখবেন, তার উত্তর দেয় পরের অংশের ছবি।

৩ · DAG — confounder, mediator, collider আর backdoor criterion

Causal DAG (directed acyclic graph) হলো চলকগুলোর একটি ছবি, যেখানে A → B মানে “A সরাসরি B-কে প্রভাবিত করে”, আর কোনো বৃত্ত নেই। AI-এর ভিত্তি কোর্সের পর্ব ৫২-এর Bayesian network-ও DAG; তফাত হলো, causal DAG-এর তীর কারণের দাবি, শুধু সম্ভাবনার গঠন নয়। Judea Pearl ১৯৯০-এর দশকে এই ভাষা দাঁড় করান (Biometrika-য় ১৯৯৫ সালের প্রবন্ধ আর Causality বই, ২০০০)। তিন চলকের তিনটি মৌলিক গঠন:

এখন প্রশ্ন: কোন চলকগুলো ধরলে T → Y-এর মোট প্রভাব পাওয়া যায়? Pearl-এর backdoor criterion: একটি সেট S বৈধ, যদি (ক) S-এর কোনো চলক T-এর বংশধর (descendant, T থেকে তীর ধরে পৌঁছানো যায় এমন) না হয়, আর (খ) S, T-তে ঢোকা তীর দিয়ে শুরু হওয়া (“পেছনের দরজা”) প্রতিটি পথ বন্ধ করে। পথ বন্ধ হয়, যদি পথের কোনো fork বা chain-এর মাঝের চলক S-এ থাকে, অথবা কোনো collider (বা তার বংশধর) S-এর বাইরে থাকে।

শেষ ধাপের collider ফাঁদটি সবচেয়ে অবাক করা, কারণ “বেশি চলক ধরলে ক্ষতি কী?” ভাবনাটি স্বাভাবিক। Joseph Berkson ১৯৪৬ সালে হাসপাতালের রোগীদের মধ্যে দুটি রোগের ভুয়া সম্পর্ক দেখিয়ে এর একটি রূপ বর্ণনা করেন, যা এখন collider bias বা Berkson's paradox নামে পরিচিত। নিচের simulation-এ coupon এলোমেলো, delivery দেরি এলোমেলো, আর দুটি পরস্পর-নিরপেক্ষ; কিন্তু দুটিই review লেখার সম্ভাবনা বাড়ায়।

import numpy as np
import pandas as pd

r = np.random.default_rng(64)
n = 20000
coupon = r.random(n) < 0.5                   # মুদ্রা ছুড়ে coupon
late = r.random(n) < 0.3                     # delivery দেরি: coupon-এর সঙ্গে সম্পর্কহীন
p_review = 0.05 + 0.25 * coupon + 0.40 * late   # দুটিই review লেখায় ঠেলে দেয়
review = r.random(n) < p_review              # collider: coupon → review ← late
d = pd.DataFrame({"coupon": coupon, "late": late, "review": review})

def late_rate(s):                            # coupon পাওয়া আর না-পাওয়াদের দেরির হার
    return s.groupby("coupon").late.mean().round(3).to_dict()
print("সবাই:            ", late_rate(d), " corr", round(np.corrcoef(coupon, late)[0, 1], 3))
rv = d[d.review]
print("শুধু review লেখা:", late_rate(rv), " corr",
      round(np.corrcoef(rv.coupon, rv.late)[0, 1], 3), f" n={len(rv)}")
nr = d[~d.review]
print("review লেখেননি:  ", late_rate(nr), " corr", round(np.corrcoef(nr.coupon, nr.late)[0, 1], 3))

সবাই:             {False: 0.303, True: 0.303}  corr -0.001
শুধু review লেখা: {False: 0.812, True: 0.496}  corr -0.289  n=5883
review লেখেননি:   {False: 0.201, True: 0.162}  corr -0.049

৪ · Propensity score — matching, IPW আর balance যাচাই

Regression-এ ফলের (খরচের) একটি model ঠিকঠাক লিখতে হয়। বিকল্প পথ: ফলের model না লিখে coupon পাওয়ার model লেখা। Paul Rosenbaum আর Donald Rubin ১৯৮৩ সালে (Biometrika) দেখান, propensity score e(x) = P(T = ১ | X = x) একটি “balancing score”: একই e(x)-এর মানুষদের মধ্যে X-এর বণ্টন দুই দলে সমান। তাই অনেক চলকের বদলে একটি সংখ্যায় তুলনা করা যায়। e(x) অনুমান করা হয় AI-এর ভিত্তি কোর্সের পর্ব ৫৯-এর logistic regression দিয়ে।

দুটি ব্যবহার:

• Matching (ATT): প্রতিটি coupon-পাওয়া গ্রাহকের জন্য সবচেয়ে কাছের e(x)-এর একজন না-পাওয়া গ্রাহক খোঁজা (১-nearest neighbour, একই control একাধিকবার চলতে পারে)। জোড়ার পার্থক্যের গড় ATT. খোঁজা হয় k-d tree দিয়ে (algorithm="kd_tree"; ব্রাউজারের Python-এ default পথটি একটি অপ্রাসঙ্গিক warning ছাপে)।
• IPW (ATE): coupon-পাওয়াদের ওজন ১/e আর না-পাওয়াদের ১/(১ − e), যেখানে e তাঁর propensity score. যাঁর coupon পাওয়ার কথা কম ছিল অথচ পেয়েছেন, তিনি বিরল, তাই ভারী: তিনি একা অনেক “অদেখা” মানুষের প্রতিনিধি। ধারণাটি জরিপ-পরিসংখ্যানের Horvitz–Thompson (১৯৫২) থেকে। এখানে দলভেদে ওজনসহ গড় নেওয়া হয় (Hájek রূপ), আর e-কে ০.০২ থেকে ০.৯৮-এ কেটে চরম ওজন ঠেকানো হয়।

কাজ শেষে balance যাচাই বাধ্যতামূলক। প্রতিটি চলকের standardised mean difference: SMD = (x̄_treated − x̄_control) ÷ √((s²_treated + s²_control) ÷ ২), যেখানে x̄ গড় আর s² variance. |SMD| < ০.১ একটি প্রচলিত হাতে-কলমে সীমা (Austin-এর লেখায় জনপ্রিয়), কোনো পরীক্ষার নিয়ম নয়। ছবিতে এগুলো পাশাপাশি দেখানোর ছবিকে বলে love plot.

import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors
X = ["orders", "tenure", "dhaka"]

def estimate(d):                                   # PS → matching-এ ATT, IPW-তে ATE
    t, y = d.coupon.to_numpy() == 1, d.spend.to_numpy()
    e = LogisticRegression(C=1e10, max_iter=1000).fit(d[X], t).predict_proba(d[X])[:, 1]
    e = np.clip(e, 0.02, 0.98)                     # চরম ওজন ঠেকাতে PS কাটা
    nn = NearestNeighbors(n_neighbors=1, algorithm="kd_tree").fit(e[~t, None])
    m = np.flatnonzero(~t)[nn.kneighbors(e[t, None])[1][:, 0]]   # প্রতিটি treated-এর জোড়া
    w = np.where(t, 1 / e, 1 / (1 - e))
    ate = np.average(y[t], weights=w[t]) - np.average(y[~t], weights=w[~t])
    return y[t].mean() - y[m].mean(), ate, t, m, w

def smd(a, b, wa=None, wb=None):                   # standardised mean difference
    return (np.average(a, weights=wa) - np.average(b, weights=wb)) / np.sqrt((a.var() + b.var()) / 2)

att, ate, t, m, w = estimate(df)
bal = pd.DataFrame({c: [smd(df[c][t], df[c][~t]), smd(df[c][t], df[c].to_numpy()[m]),
                        smd(df[c][t], df[c][~t], w[t], w[~t])] for c in X},
                   index=["before", "matched", "IPW"]).T.round(3)
print(bal, "\nসবচেয়ে বড় ওজন", w.max().round(1), "| জোড়ায় আলাদা control", len(set(m)))
boot, rng = [], np.random.default_rng(0)
for _ in range(200):                               # bootstrap: PS-সহ পুরো পদ্ধতি আবার
    boot.append(estimate(df.iloc[rng.integers(0, len(df), len(df))])[:2])
lo, hi = np.percentile(boot, [2.5, 97.5], axis=0)
print(f"matching ATT {att:.1f} [{lo[0]:.1f}, {hi[0]:.1f}] | আসল ATT {true_att:.1f}")
print(f"IPW ATE      {ate:.1f} [{lo[1]:.1f}, {hi[1]:.1f}] | আসল ATE {true_ate:.1f}")
fig, ax = plt.subplots(figsize=(6, 2.8))
for k, mk in zip(bal.columns, "osD"):
    ax.plot(bal[k].abs(), bal.index, mk, label=k)
ax.axvline(0.1, ls="--", c="gray"); ax.legend(); fig.subplots_adjust(left=0.18, bottom=0.2)
ax.set(xlabel="|standardised mean difference|", title="Covariate balance (love plot)")
plt.show()

        before  matched    IPW
orders   0.839    0.011  0.013
tenure   0.113   -0.022 -0.012
dhaka   -0.242   -0.002  0.002 
সবচেয়ে বড় ওজন 16.6 | জোড়ায় আলাদা control 346
matching ATT 115.0 [95.9, 136.9] | আসল ATT 110.4
IPW ATE      106.5 [90.8, 121.0] | আসল ATE 98.8

Overlap যাচাই: দুই দলের e-এর histogram পাশাপাশি আঁকুন। কোনো অঞ্চলে শুধু এক দলের মানুষ থাকলে (যেমন e > ০.৯৫-এ প্রায় কেউ coupon ছাড়া নেই), সেখানে তুলনার মতো কেউ নেই; তখন e কাটা (clipping) বা ওই মানুষদের বাদ দেওয়া (trimming, যেমন e ০.১ থেকে ০.৯-এর বাইরে বাদ) প্রচলিত পথ। কিন্তু trimming করলে estimand বদলায়: তখন উত্তর শুধু “তুলনীয়” মানুষদের জন্য, আর রিপোর্টে সেটি স্পষ্ট লিখতে হয়। Overlap দুর্বল হলে regression-ও নিরাপদ নয়; সে তখন ডেটা নেই এমন জায়গায় সরলরেখা টেনে অনুমান করে।

৫ · Difference-in-differences — delivery fee বাড়ানোর প্রভাব

সপ্তাহ ৯ থেকে বাজারঘর ১৬টি জেলার ৬টিতে delivery fee বাড়িয়েছে। প্রশ্ন: সাপ্তাহিক অর্ডার (প্রতি ১,০০০ গ্রাহকে) কতটা কমল? দুটি সরল তুলনাই ভুল। শুধু আগে-পরে তুলনায় সময়ের প্রভাব (মৌসুম, ঈদ, সাধারণ প্রবৃদ্ধি) মিশে যায়। শুধু fee-বাড়ানো বনাম না-বাড়ানো জেলার তুলনায় জেলার নিজস্ব পার্থক্য মিশে যায় (fee বাড়ানো হয়েছে যেখানে চাহিদা আগে থেকেই বেশি)।

Difference-in-differences (DiD) দুটি পার্থক্য একসঙ্গে বাদ দেয়:

DiD = (ȳ_{treated, পরে} − ȳ_{treated, আগে}) − (ȳ_{control, পরে} − ȳ_{control, আগে})

ȳ মানে সেই ঘরের গড় অর্ডার। Control জেলার পরিবর্তনটি “fee না বাড়ালে treated জেলায় কী হতো”-র প্রতিনিধি। এর মূল অনুমান parallel trends: fee না বাড়ালে দুই দলের গড় একই হারে বদলাত। মাত্রা আলাদা হতে পারে, ধারা নয়। ক্রিকেটের ভাষায়: দুই দলের রান রেট একই হারে বাড়ছিল; ১০ ওভারের পর এক দলের ব্যাটার বদলালে, রান রেটের ফাঁক কতটা বদলাল সেটিই প্রভাব।

একই হিসাব regression-এ: y = β₀ + β₁·treated + β₂·post + β₃·(treated × post) + ε, যেখানে treated আর post ০/১ চলক, আর ε বাকি এলোমেলোতা। Interaction-এর সহগ β₃ ঠিক DiD. Regression-এর সুবিধা SE: একই জেলার ১৬ সপ্তাহের অর্ডার পরস্পর-সম্পর্কিত, তাই সাধারণ SE অতি-আশাবাদী। Bertrand, Duflo আর Mullainathan (QJE, ২০০৪) দেখান, এমন ডেটায় সাধারণ SE ব্যবহার করলে ভুয়া “significant” ফল অনেক বেড়ে যায়। সমাধান হলো জেলা ধরে cluster-robust SE নেওয়া। নিচের ব্লক দুটি জগৎ চালায়: একটিতে parallel trends সত্যি, অন্যটিতে fee-বাড়ানো জেলার চাহিদা আগে থেকেই প্রতি সপ্তাহে ০.৬ করে নামছিল। আসল প্রভাব দুটিতেই −৫. শেষে একটি placebo test: শুধু সপ্তাহ ১–৮ নিয়ে ভান করা হয় fee বেড়েছে সপ্তাহ ৫-এ।

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from bazar_causal import fee_panel

def did(p, name):                        # 2×2 হাতে, তারপর OLS-এর interaction
    t = p.groupby(["treated", "post"]).orders.mean().unstack()
    by_hand = (t.loc[1, 1] - t.loc[1, 0]) - (t.loc[0, 1] - t.loc[0, 0])
    fit = smf.ols("orders ~ treated * post", data=p).fit(
        cov_type="cluster", cov_kwds={"groups": p.district})
    lo, hi = fit.conf_int().loc["treated:post"]
    print(f"{name:13s} হাতে {by_hand:6.2f} | OLS {fit.params['treated:post']:6.2f}, "
          f"cluster SE {fit.bse['treated:post']:.2f}, 95% CI [{lo:.2f}, {hi:.2f}]")
    return t.round(2)

good = fee_panel()                       # parallel trends সত্যি; আসল প্রভাব -5
bad = fee_panel(trend_gap=-0.6)          # treated জেলার চাহিদা আগে থেকেই নামছিল
print(did(good, "সমান্তরাল"))
did(bad, "ধারা ভাঙা")
for name, p in [("placebo ঠিক", good), ("placebo ভাঙা", bad)]:   # শুধু সপ্তাহ 1–8, ভান 5-এ
    pre = p[p.week <= 8]
    did(pre.assign(post=(pre.week >= 5).astype(int)), name)
fig, axs = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
for ax, p, ttl in zip(axs, [good, bad], ["Parallel pre-trends", "Diverging pre-trends"]):
    wk = p.groupby(["week", "treated"]).orders.mean().unstack()
    ax.plot(wk.index, wk[0], "o-", ms=3, label="control districts")
    ax.plot(wk.index, wk[1], "s-", ms=3, label="fee raised (week 9)")
    ax.axvline(8.5, ls="--", c="gray"); ax.set(title=ttl, xlabel="week")
axs[0].set_ylabel("orders per 1,000 customers"); axs[0].legend(fontsize=8)
fig.subplots_adjust(bottom=0.17, wspace=0.08)
plt.show()

সমান্তরাল     হাতে  -5.35 | OLS  -5.35, cluster SE 0.50, 95% CI [-6.33, -4.36]
post         0      1
treated              
0        58.27  60.20
1        71.03  67.61
ধারা ভাঙা     হাতে -10.15 | OLS -10.15, cluster SE 0.50, 95% CI [-11.13, -9.16]
placebo ঠিক   হাতে   0.71 | OLS   0.71, cluster SE 0.51, 95% CI [-0.29, 1.72]
placebo ভাঙা  হাতে  -1.69 | OLS  -1.69, cluster SE 0.51, 95% CI [-2.69, -0.68]

DiD-এর একটি বিখ্যাত প্রয়োগ David Card আর Alan Krueger-এর ১৯৯৪ সালের ন্যূনতম মজুরির গবেষণা: New Jersey-তে মজুরি বাড়ার আগে-পরে পাশের Pennsylvania-র fast-food দোকানের সঙ্গে তুলনা। এর আগের এক উদাহরণ হিসেবে প্রায়ই John Snow-এর ১৮৫০-এর দশকের লন্ডনের কলেরা-বিশ্লেষণের কথা বলা হয়, যেখানে দুটি পানি-সরবরাহ কোম্পানির এলাকার মৃত্যুহার তুলনা করা হয়েছিল।

৬ · আরও দুটি নকশা — regression discontinuity আর instrumental variable

এই দুটি পদ্ধতির পুরো চর্চা এখানে নয়; ধারণা আর শর্তগুলো জানা থাকলে কখন বিশেষজ্ঞের সাহায্য নেবেন বুঝবেন।

Card, Angrist আর Imbens ২০২১ সালে অর্থনীতির নোবেল-স্মারক পুরস্কার পান, মূলত এমন “প্রাকৃতিক পরীক্ষা” থেকে কারণ খোঁজার পদ্ধতি নিয়ে কাজের জন্য। মূল কথা একই: কোথাও একটি এলোমেলো-ঘেঁষা কারণে treatment বদলেছে, এমন জায়গা খোঁজা।

৭ · লুকানো confounder, DoWhy/EconML আর পর্যবেক্ষণের সীমা

অংশ ২–৪-এর সব পদ্ধতি দাঁড়িয়ে আছে “সব confounder মাপা হয়েছে” অনুমানের ওপর। বাস্তবে আয়, পরিবারের আকার, প্রতিযোগীর অফার প্রায়ই অজানা। Sensitivity analysis প্রশ্ন বদলায়: “লুকানো confounder কত শক্তিশালী হলে আমার সিদ্ধান্ত উল্টে যেত?” এর একটি জনপ্রিয় রূপ VanderWeele আর Ding-এর (২০১৭) E-value: risk ratio RR > ১ হলে E = RR + √(RR × (RR − ১))। যেমন RR = ২ হলে E = ২ + √২ ≈ ৩.৪১: লুকানো confounder-কে treatment আর ফল দুটির সঙ্গেই অন্তত ৩.৪১ গুণ risk ratio-য় জড়িত থাকতে হবে, তবেই পুরো প্রভাব ব্যাখ্যা হতো। এটি একটি নির্দেশক, প্রমাণ নয়; ধারণাটির পূর্বসূরি Cornfield ও সহকর্মীদের ১৯৫৯ সালের ধূমপান-ক্যান্সার যুক্তি।

DoWhy আর EconML দিয়ে সত্যিকারের run

DoWhy লাইব্রেরি কাজটিকে চার ধাপে বাঁধে: DAG দিয়ে model, identify (backdoor সেট খোঁজা), estimate, আর refute (ভুয়া treatment দিয়ে ফল ভেঙে পড়ে কি না দেখা)। শুরু Microsoft Research-এ, এখন PyWhy নামের open-source উদ্যোগের অংশ। EconML মানুষভেদে প্রভাব (CATE) অনুমানের পদ্ধতি দেয়, যেমন double machine learning (Chernozhukov ও সহকর্মী, ২০১৮)। দুটির কোনোটিই ব্রাউজারের Python-এ নেই, তাই ডেস্কটপে একটি আলাদা uv project-এ একই bazar_causal.py দিয়ে সত্যিকারের run (একটি যন্ত্র, একটি run; লেখার সময়ে dowhy ০.১৪-এর ঘোষিত Python সীমা ৩.১৪-এর নিচে, তাই Python ৩.১৩):

import sys, warnings, logging
import numpy as np, networkx as nx
import dowhy, econml
from dowhy import CausalModel
from econml.dml import LinearDML
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
from bazar_causal import coupon_data

warnings.filterwarnings("ignore"); logging.getLogger("dowhy").setLevel(logging.ERROR)
print(f"Python {sys.version.split()[0]} | dowhy {dowhy.__version__} | econml {econml.__version__}")
df = coupon_data()
data = df[["orders", "tenure", "dhaka", "coupon", "spend"]].astype({"coupon": bool})
g = nx.DiGraph([(c, "coupon") for c in ("orders", "tenure", "dhaka")] +
               [(c, "spend") for c in ("orders", "tenure", "dhaka")] + [("coupon", "spend")])
model = CausalModel(data=data, treatment="coupon", outcome="spend", graph=g)
est = model.identify_effect()
print("backdoor set:", est.get_backdoor_variables())
for m, unit in [("backdoor.linear_regression", "ate"),
                ("backdoor.propensity_score_weighting", "ate"),
                ("backdoor.propensity_score_matching", "att")]:
    e = model.estimate_effect(est, method_name=m, target_units=unit)
    print(f"{m:38s} {unit}: {e.value:7.1f}")
lin = model.estimate_effect(est, method_name="backdoor.linear_regression")
ref = model.refute_estimate(est, lin, method_name="placebo_treatment_refuter",
                            placebo_type="permute", num_simulations=50, random_seed=0)
print(f"placebo refuter: new effect {ref.new_effect:.1f}, p-value {ref.refutation_result['p_value']:.2f}")
dml = LinearDML(model_y=GradientBoostingRegressor(random_state=0),
                model_t=GradientBoostingClassifier(random_state=0),
                discrete_treatment=True, random_state=0)
dml.fit(df.spend, df.coupon, X=df[["orders"]], W=df[["tenure", "dhaka"]])
print(f"LinearDML: effect = {dml.intercept_:.1f} + {dml.coef_[0]:.1f} x orders  (truth 40 + 15 x orders)")
lo, hi = dml.effect_interval(np.array([[2], [8]]), alpha=0.05)
print("effect at orders=2, 8:", dml.effect(np.array([[2], [8]])).round(1), "95% CI", lo.round(1), hi.round(1))

~/bazarghar-causal$ uv init --python 3.13 --vcs none
…
~/bazarghar-causal$ uv add dowhy econml
…
~/bazarghar-causal$ uv run python causal_check.py
Python 3.13.5 | dowhy 0.14 | econml 0.17.0
backdoor set: ['tenure', 'dhaka', 'orders']
backdoor.linear_regression             ate:    98.7
backdoor.propensity_score_weighting    ate:   106.7
backdoor.propensity_score_matching     att:   123.2
placebo refuter: new effect -0.3, p-value 0.48
LinearDML: effect = 61.2 + 10.5 x orders  (truth 40 + 15 x orders)
effect at orders=2, 8: [ 82.3 145.3] 95% CI [ 65.2 117. ] [ 99.3 173.7]

DoWhy DAG থেকে নিজেই backdoor সেট (তিনটি confounder) বের করেছে; linear regression পাতার ব্লকের সঙ্গে হুবহু (৯৮.৭)। IPW ১০৬.৭, পাতার ১০৬.৫-এর কাছে (ভেতরের propensity model-এর সেটিং আলাদা হতে পারে)। Matching-এর ATT ১২৩.২, পাতার ১১৫.০ থেকে বেশ দূরে: DoWhy-র matching-এর ভেতরের খুঁটিনাটি (কী দিয়ে দূরত্ব মাপা হয়) আলাদা, আর matching অনুমান এমন খুঁটিনাটিতে নড়ে। Placebo refuter coupon কলাম এলোমেলো করে ৫০ বার চালায়; প্রভাব −০.৩, প্রত্যাশিতভাবে প্রায় শূন্য। এই পরীক্ষা লুকানো confounder ধরে না, শুধু কোডের বড় ভুল ধরে। LinearDML-এর রেখা ৬১.২ + ১০.৫ × orders, আসল ৪০ + ১৫ × orders; orders = ২ আর ৮-এ আসল মান ৭০ আর ১৬০, দুটিই ৯৫% পরিসরের ভেতরে, তবে ঢাল কম ধরা পড়েছে। লেখার সময়ের version প্রথম লাইনে; API বদলায়, তাই নিজের version-এর docs দেখুন।

পর্যবেক্ষণের ডেটার সৎ সীমা

Robert LaLonde ১৯৮৬ সালে একটি চাকরি-প্রশিক্ষণ কর্মসূচির randomised পরীক্ষার ফলকে মানদণ্ড ধরে দেখান, সাধারণ পর্যবেক্ষণভিত্তিক পদ্ধতিগুলো প্রায়ই সেই ফল থেকে অনেক দূরে, এমনকি উল্টো দিকে যায়। পরে propensity score দিয়ে আরও কাছে পৌঁছানোর দাবি এসেছে, আর সেই দাবির সংবেদনশীলতা নিয়ে তর্কও হয়েছে। চিকিৎসায় menopause-এর পরে hormone therapy নিয়ে পর্যবেক্ষণের গবেষণা হৃদরোগে সুরক্ষার ইঙ্গিত দিয়েছিল; ২০০২ সালে Women's Health Initiative-এর randomised trial সেই সুরক্ষা খুঁজে পায়নি বলে জানায়, আর ব্যাখ্যা নিয়ে আলোচনা এখনো চলে। শিক্ষা: পদ্ধতির নাম যত জটিলই হোক, শক্তি আসে অনুমান থেকে। সিদ্ধান্তের রিপোর্টে লিখুন estimand (ATE না ATT), DAG, কোন অনুমান যাচাই হয়েছে (balance, overlap, pre-trend, placebo) আর কোনটি হয়নি (unconfoundedness, exclusion), আর sensitivity. সম্ভব হলে একটি ছোট A/B test-ই সবচেয়ে সৎ উত্তর।

৮ · AI-এর যুগে — causal দাবি যাচাই

“Coupon-এর প্রভাব বের করে দাও” বললে AI assistant দ্রুত কোড দেয়, কিন্তু ভুলগুলো প্রায়ই চুপচাপ থাকে:

AI-এর দেওয়া causal বিশ্লেষণ হাতে পেলে — যাচাইয়ের চেকলিস্ট

• Estimand (ATE, ATT বা অন্য) প্রশ্ন থেকে বেছে লিখেছি
• DAG এঁকেছি; treatment-এর পরের চলক adjustment-এ রাখিনি
• Overlap আর SMD balance দেখেছি, শুধু শেষ সংখ্যা নয়
• DiD-এ pre-trend ছবি, placebo আর cluster SE দেখেছি
• লাইব্রেরির version-এর docs মিলিয়েছি, আর জানা উত্তরের একটি simulation-এ কোড যাচাই করেছি
• যাচাই-অযোগ্য অনুমান আর sensitivity রিপোর্টে লিখেছি; আসল গ্রাহকের ডেটা AI-তে দিইনি