Survival analysis — ক্রেতা কখন চলে যাবে, Kaplan-Meier আর Cox

এতদিন label ছিল একটি সংখ্যা বা একটি শ্রেণি। এবার প্রশ্ন “কখন?”: একজন ক্রেতা কত দিন পর কেনা বন্ধ করবেন। পর্ব ০২-এ churn-কে সাধারণত একটি হ্যাঁ/না সমস্যা হিসেবে দেখেছেন। এই পর্বে দেখবেন, সময়ের প্রশ্নকে হ্যাঁ/না-তে নামালে কী হারায়, আর পরিসংখ্যানের একটি পুরোনো শাখা, survival analysis, কীভাবে তা সামলায়। ভিত্তি তিনটি: পর্ব ১৮-এর maximum likelihood, পর্ব ২০-এর exposure-এর ধারণা, আর AI-এর ভিত্তি কোর্সের পর্ব ৩৬-এর distribution. চারটি চালানো যায় এমন ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, pandas, SciPy, statsmodels, scikit-learn আর matplotlib); প্রথমবার package নামতে কয়েক সেকেন্ড, scikit-learn-এ প্রায় দশ সেকেন্ড লাগে।

১ · সময় পর্যন্ত ঘটনা আর censoring — অসম্পূর্ণ গল্পও তথ্য

Time-to-event ডেটায় প্রতিটি সারির label দুটি সংখ্যা: একটি সময় T আর একটি ঘটনা-সূচক δ (1 = ঘটনা দেখা গেছে, 0 = দেখা যায়নি)। বাজারঘরে ঘটনা “ক্রেতা চলে গেছেন” (ধরুন ৯০ দিন কোনো অর্ডার নেই, আর চলে যাওয়ার দিন শেষ অর্ডারের দিন), আর সময় মাপা হয় যোগদানের দিন থেকে। ডেটা কাটার দিন ধরা হলো ৩১ মার্চ ২০২৬. যে ক্রেতা তখনো সক্রিয়, তাঁর সময় হলো যোগদান থেকে কাটার দিন পর্যন্ত, আর δ = 0. একে বলে right censoring: আসল সময় এই সংখ্যার ডানে, কোথাও। এখানে censoring আসে শুধু ডেটা কাটার দিন থেকে, তাই একে administrative censoring বলে।

একটি গুরুত্বপূর্ণ অনুমান সবকিছুর নিচে আছে: non-informative censoring. অর্থাৎ censoring-এর কারণ ক্রেতার চলে যাওয়ার ঝুঁকির সঙ্গে সম্পর্কহীন। ডেটা কাটার দিন ক্যালেন্ডারের ঘটনা, ক্রেতার আচরণের নয়, তাই এখানে অনুমানটি যুক্তিসংগত। কিন্তু যদি দোকান নিজেই “ঝামেলার” ক্রেতাদের অ্যাকাউন্ট বন্ধ করে আর তাঁদের censored লেখে, তাহলে অনুমান ভাঙে, আর এই পর্বের কোনো পদ্ধতিই তা ঠিক করতে পারে না।

২ · সাধারণ পথে কেন ভুল — গড়, classifier আর censoring

প্রথমে একটি ডেটাসেট। নিচের module ৮০০ জন কাল্পনিক বাজারঘর ক্রেতা বানায়: যোগদানের দিন (২০২৪–২০২৫), Plus সদস্য কি না (plan), প্রথম মাসে ডেলিভারি নিয়ে অভিযোগের সংখ্যা (complaints), আর কুপন দিয়ে এসেছেন কি না (promo). আসল চলে যাওয়ার সময় আসে একটি Weibull model থেকে (অংশ ৩-এ ব্যাখ্যা): Plus সদস্যের ঝুঁকি ০.৬ গুণ, প্রতিটি অভিযোগে ১.৩৫ গুণ, আর কুপনের ক্রেতাদের ঝুঁকির আকার আলাদা। true_time কলামটি শুধু simulation-এ জানা “সত্য”, যা বাস্তবে কখনো পাওয়া যায় না; এটি রাখা হয়েছে পদ্ধতিগুলোর ভুল মাপার জন্য।

# বাজারঘর (কাল্পনিক দোকান): ক্রেতা কবে চলে যায়, seeded ডেটা
import numpy as np
import pandas as pd

CUTOFF = pd.Timestamp("2026-03-31")      # ডেটা কাটার দিন

def customers(n=800, seed=55):
    rng = np.random.default_rng(seed)
    start = pd.Timestamp("2024-01-01")
    signup = start + pd.to_timedelta(rng.integers(0, 731, n), unit="D")
    plan = rng.binomial(1, 0.4, n)           # 1 = Plus সদস্য
    complaints = rng.poisson(0.8, n)         # প্রথম মাসে ডেলিভারি-অভিযোগ
    promo = rng.binomial(1, 0.4, n)          # 1 = কুপন দিয়ে আসা
    # আসল model (Weibull): কুপনের ক্রেতার আকার আলাদা
    k = np.where(promo == 1, 0.6, 1.3)
    lam = np.where(promo == 1, 700.0, 520.0)
    risk = np.exp(np.log(0.6) * plan + np.log(1.35) * complaints)
    true_t = lam * (-np.log(rng.uniform(size=n)) / risk) ** (1 / k)
    follow = (CUTOFF - signup).days.to_numpy()   # কাটার দিন পর্যন্ত কত দিন
    event = (true_t <= follow).astype(int)       # 1 = চলে গেছে, 0 = censored
    time = np.where(event == 1, np.ceil(true_t), follow).astype(float)
    return pd.DataFrame({"signup": signup, "plan": plan,
                         "complaints": complaints, "promo": promo,
                         "time": time, "event": event,
                         "true_time": np.ceil(true_t)})

এখন চারটি সাধারণ কিন্তু ভুল পথ। (১) শুধু যাঁরা চলে গেছেন, তাঁদের গড় সময়। (২) censored ক্রেতাদেরও “চলে গেছেন” ধরে সবার গড়। (৩) যোগদানের বছর ধরে “চলে গেছেন কি না”-র হার। (৪) সেই হ্যাঁ/না label-এ একটি logistic regression (AI-এর ভিত্তি কোর্সের পর্ব ৫৯), যেখানে একটি feature হলো ক্রেতাকে কত দিন দেখা হয়েছে। প্রতিটির পাশে simulation-এর সত্য।

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from bazar_surv import customers, CUTOFF

df = customers()
pd.set_option("display.width", 70)
print(df[["signup", "plan", "complaints", "promo", "time", "event"]].head(4))
print("censored:", f"{1 - df.event.mean():.1%}")
# সত্য (শুধু simulation-এ জানা)
print("truth P(churn within 365 d):", f"{(df.true_time <= 365).mean():.3f}")
print("truth median lifetime:", df.true_time.median(), "days")
# ভুল ১ আর ২: গড় সময়
print("mean time, churned only:", round(df.time[df.event == 1].mean()), "days")
print("mean time, all as events:", round(df.time.mean()), "days")

# ভুল ৩: 'চলে গেছে কি না' — যোগদানের বছর ধরে
df["year"] = df.signup.dt.year
naive = df.groupby("year").event.mean()
truth = df.groupby("year").true_time.apply(lambda t: (t <= 365).mean())
print(pd.DataFrame({"churned (naive)": naive, "true 1-yr churn": truth})
      .round(3))

# ভুল ৪: churned/not classifier — কত দিন দেখা হয়েছে, সেটাই শেখে
seen = (CUTOFF - df.signup).dt.days / 100          # পর্যবেক্ষণ, 100 দিনে
X = np.c_[df.plan, df.complaints, df.promo, seen]
clf = LogisticRegression(C=1e10).fit(X, df.event)
print("odds ratio per 100 days observed:", round(np.exp(clf.coef_[0][3]), 2))
same = [[0, 0, 0, 1.0], [0, 0, 0, 7.0]]            # একই ক্রেতা, 100 বনাম 700 দিন
print("P(churn) same customer:", clf.predict_proba(same)[:, 1].round(3))

      signup  plan  complaints  promo   time  event
0 2025-12-05     0           0      1  116.0      0
1 2025-08-31     1           1      0   93.0      1
2 2025-06-15     1           0      0  289.0      0
3 2025-09-28     0           0      1  184.0      0
censored: 43.9%
truth P(churn within 365 d): 0.491
truth median lifetime: 376.0 days
mean time, churned only: 212 days
mean time, all as events: 285 days
      churned (naive)  true 1-yr churn
year
2024            0.711            0.492
2025            0.423            0.490
odds ratio per 100 days observed: 1.44
P(churn) same customer: [0.258 0.757]

৩ · Survival function আর hazard — দুই চোখে একই গল্প

দুটি function পুরো বণ্টনের বর্ণনা দেয়। Survival function S(t) = P(T > t): একজন ক্রেতা t দিনের পরেও থাকবেন, তার সম্ভাবনা। S(0) = ১, আর সময়ের সঙ্গে S কখনো বাড়ে না। দ্বিতীয়টি hazard function:

h(t) = lim_Δt→০ P(t ≤ T < t + Δt | T ≥ t) / Δt = f(t) / S(t)

অর্থাৎ “আজ পর্যন্ত যিনি আছেন, তিনি ঠিক এখন চলে যাওয়ার হার”। এখানে f(t) হলো T-এর density. ছোট Δt-তে h(t)·Δt ≈ পরের Δt দিনে চলে যাওয়ার শর্তসাপেক্ষ সম্ভাবনা। hazard একটি হার, সম্ভাবনা নয়; তাই ১-এর বেশিও হতে পারে (একক: প্রতি দিন)। hazard যোগ করলে পাওয়া যায় cumulative hazard H(t) = ∫_০^t h(u) du, আর তিনটির বাঁধন একটি সূত্রে:

S(t) = exp(−H(t))

সবচেয়ে সরল model exponential: hazard স্থির, h(t) = ১/λ, তাই S(t) = e^−t/λ. এর মানে “স্মৃতিহীন”: দুই বছর থাকা ক্রেতা আর নতুন ক্রেতার আগামী মাসে চলে যাওয়ার ঝুঁকি সমান। বাস্তবে এটি কমই সত্য। Weibull একটি আকার-parameter k যোগ করে:

S(t) = exp(−(t/λ)^k),  h(t) = (k/λ)(t/λ)^k−১

λ হলো scale (দিনে), k হলো shape. k = ১ হলে exponential. k < ১ হলে hazard সময়ের সঙ্গে কমে: নতুন ক্রেতা বেশি ঝুঁকিতে, যাঁরা টিকে যান তাঁরা অনুগত। কুপনের ক্রেতাদের জন্য simulation-এ k = ০.৬, ঠিক এই ছাঁদ: কুপন ফুরালে অনেকে চলে যান। k > ১ হলে hazard বাড়ে, যেমন অভ্যাস ধীরে ক্ষয়ে যাওয়া; সাধারণ ক্রেতাদের জন্য k = ১.৩. Median জীবনকাল হলো S(t) = ০.৫-এর সমাধান: λ(ln ২)^১/k. আর ক্রেতার বৈশিষ্ট্য x ঝুঁকিকে e^xβ গুণ করলে S_x(t) = S₀(t)^exβ, যা অংশ ৫-এ আবার আসবে।

৪ · Kaplan–Meier — হাতে, Greenwood আর statsmodels

১৯৫৮ সালে Edward Kaplan আর Paul Meier যে Kaplan–Meier estimator প্রকাশ করেন, তা উপরের ধারণার সরাসরি রূপ। ঘটনার দিনগুলো t₁ < t₂ < … সাজান। প্রতিটি t_j-এ n_j = তখনো ঝুঁকিতে থাকা মানুষ (risk set: যাঁদের সময় ≥ t_j), আর d_j = সেদিন চলে যাওয়া মানুষ। তাহলে

Ŝ(t) = ∏_{tj ≤ t} (১ − d_j/n_j)

প্রতিটি গুণক “ওই দিনে যারা ছিল, তাদের কত ভাগ টিকল”। Censored মানুষ যতদিন দেখা গেছে ততদিন n-এ গোনা হন, তারপর নিঃশব্দে সরে যান। অনিশ্চয়তার জন্য Major Greenwood ১৯২৬ সালে যে সূত্র দেন, তা আজও আদর্শ: Greenwood formula

Var[Ŝ(t)] ≈ Ŝ(t)² · ∑_{tj ≤ t} d_j / (n_j(n_j − d_j))

নিচের stepper ছয়জন ক্রেতার ছোট একটি দলে পুরো হিসাব করে, যার দুজন censored. সময় মাসে: ২, ৪+, ৫, ৫, ৭+, ৯ (“+” মানে censored).

এবার ৮০০ জনে। দ্বিতীয় ব্লক Kaplan–Meier আর Greenwood নিজে লেখে, তারপর statsmodels-এর SurvfuncRight-এর সঙ্গে মেলায়, তিনটি দিনে S আর ৯৫% পরিসর ছাপায়, Basic আর Plus দলের ধাপ-আকারের বক্ররেখা আঁকে, আর দুই দলকে log-rank test-এ (অংশ ৫-এ ব্যাখ্যা) তুলনা করে।

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.duration.survfunc import SurvfuncRight, survdiff
from bazar_surv import customers

def km(time, event):                 # Kaplan-Meier + Greenwood, নিজে লেখা
    t = np.unique(time[event == 1])
    n = np.array([(time >= u).sum() for u in t])           # ঝুঁকিতে
    d = np.array([((time == u) & (event == 1)).sum() for u in t])
    S = np.cumprod(1 - d / n)
    se = S * np.sqrt(np.cumsum(d / (n * (n - d))))         # Greenwood
    return t, S, se

df = customers()
T, E = df.time.to_numpy(), df.event.to_numpy()
t, S, se = km(T, E)
sf = SurvfuncRight(T, E)
print("max |S - statsmodels|  =", f"{np.max(np.abs(S - sf.surv_prob)):.1e}")
print("max |se - statsmodels| =", f"{np.max(np.abs(se - sf.surv_prob_se)):.1e}")
for day in [90, 365, 730]:
    i = np.searchsorted(t, day, side="right") - 1
    print(f"S({day}) = {S[i]:.3f}, 95% CI [{S[i] - 1.96*se[i]:.3f},"
          f" {S[i] + 1.96*se[i]:.3f}]")
print("median (days):", sf.quantile(0.5))
for p, name in [(0, "Basic"), (1, "Plus")]:
    g = df[df.plan == p]
    tg, Sg, _ = km(g.time.to_numpy(), g.event.to_numpy())
    print(name, "n =", len(g), "median =", SurvfuncRight(g.time, g.event).quantile(0.5))
    plt.step(np.r_[0, tg], np.r_[1, Sg], where="post", label=name)
chi2, pval = survdiff(T, E, df.plan)
print(f"log-rank: chi2 = {chi2:.1f}, p = {pval:.1e}")
plt.xlabel("days since sign-up"); plt.ylabel("S(t): still a customer")
plt.legend(); plt.show()

max |S - statsmodels|  = 4.4e-16
max |se - statsmodels| = 2.1e-17
S(90) = 0.831, 95% CI [0.805, 0.857]
S(365) = 0.506, 95% CI [0.468, 0.543]
S(730) = 0.221, 95% CI [0.170, 0.273]
median (days): 371.0
Basic n = 468 median = 276.0
Plus n = 332 median = 497.0
log-rank: chi2 = 25.7, p = 4.0e-07

৫ · দুটি দল তুলনা — log-rank test

দুটি বক্ররেখা আলাদা দেখালেই পার্থক্য সত্য নয়; এলোমেলো ওঠানামাও হতে পারে। Log-rank test প্রশ্নটি এভাবে করে: null hypothesis “দুই দলের hazard সব সময়ে সমান” সত্য হলে, প্রতিটি ঘটনার দিনে d_j টি ঘটনা দুই দলের মধ্যে তাদের ঝুঁকিতে থাকা মানুষের অনুপাতে ভাগ হওয়ার কথা। দল ১-এর প্রত্যাশিত ঘটনা E_1j = d_j · n_1j/n_j. সব দিনে দেখা O₁ আর প্রত্যাশিত E₁ যোগ করে, (O₁ − E₁)²/V-কে chi-square বণ্টনের (১ degree of freedom) সঙ্গে তুলনা করা হয়, যেখানে V হলো hypergeometric variance-এর যোগফল।

ছোট উদাহরণ: এক দিনে ১০০ জন ঝুঁকিতে, Basic ৬০ আর Plus ৪০, আর সেদিন ৫ জন চলে গেলেন। Null-এ Basic-এর প্রত্যাশা ৫ × ৬০/১০০ = ৩. বাস্তবে ৪ জন Basic হলে O − E = +১. এমন শত শত দিনের যোগফল বড় হলে পার্থক্য এলোমেলো নয়। বাজারঘরে χ² = ২৫.৭, p-value ৪ × ১০^−৭: Basic আর Plus-এর churn-এর ছাঁদ সত্যিই আলাদা। Nathan Mantel ১৯৬৬ সালে এই ধরনের test প্রস্তাব করেন, আর Richard Peto ও Julian Peto ১৯৭২ সালে “log-rank” নাম দেন।

৬ · Cox proportional hazards — hazard ratio আর concordance

অনেকগুলো বৈশিষ্ট্য একসঙ্গে দেখতে চাইলে দল ভাগ করে চলে না। David Cox ১৯৭২ সালে Journal of the Royal Statistical Society, Series B-তে যে model প্রস্তাব করেন, তা আজ survival analysis-এর সবচেয়ে ব্যবহৃত regression: Cox proportional hazards model

h(t | x) = h₀(t) · exp(x₁β₁ + … + x_pβ_p)

h₀(t) হলো baseline hazard, সব x = ০ হলে hazard; তার আকার সম্পূর্ণ অজানা রাখা হয়। বৈশিষ্ট্যগুলো শুধু তাকে একটি স্থির গুণক দিয়ে ওঠায়-নামায়। তাই দুজন ক্রেতার hazard-এর অনুপাত সময়ের ওপর নির্ভর করে না: এটিই “proportional hazards”. e^β হলো hazard ratio (HR): x এক একক বাড়লে hazard যত গুণ হয়। একটি অংশ অজানা আর অন্যটি parametric বলে Cox model-কে বলা হয় semi-parametric.

h₀ না জেনে β শেখা যায় কীভাবে? Cox-এর কৌশল partial likelihood. প্রতিটি ঘটনার দিনে প্রশ্ন করুন: “ঝুঁকিতে থাকা সবার মধ্যে ঠিক এই ক্রেতাই চলে গেলেন, এর সম্ভাবনা কত?” উত্তর তাঁর hazard ভাগ সবার hazard-এর যোগফল, আর h₀(t) উপরে-নিচে কেটে যায়:

L(β) = ∏_{ঘটনা i} e^x_iβ / ∑_{j ∈ R(ti)} e^x_jβ

R(t_i) হলো সেই দিনের risk set. লক্ষ করুন, ঘটনার সময়গুলোর শুধু ক্রম কাজে লাগে, মান নয়। Censored ক্রেতারা কোনো গুণক যোগ করেন না, কিন্তু যতদিন ছিলেন ততদিন হরে থাকেন। তারপর পর্ব ১৮-এর মতো log L-কে Newton পদ্ধতিতে (পর্ব ১৫) সর্বোচ্চ করা হয়। একই দিনে একাধিক ঘটনা (tie) হলে আনুমানিক রূপ লাগে: Breslow-এর সরল রূপ, বা Efron-এর একটু নিখুঁত রূপ। এখন HR পড়তে শিখুন।

Cox model-এর মূল্যায়নে accuracy বা RMSE চলে না, কারণ censored সারির সঠিক উত্তর অজানা। প্রধান মাপ concordance index (Harrell's C), যা Frank Harrell ও সহকর্মীরা ১৯৮২ সালে JAMA-তে ব্যবহার করেন। ধারণাটি ক্রমের: দুজন ক্রেতার জোড়া নিন, যাঁদের মধ্যে কে আগে চলে গেছেন তা নিশ্চিত জানা যায় (তুলনাযোগ্য জোড়া: প্রথমজনের ঘটনা দেখা গেছে, আর দ্বিতীয়জন তখনো ছিলেন)। Model যদি আগে চলে যাওয়া জনকে বেশি ঝুঁকি দেয়, জোড়াটি concordant. C = concordant জোড়া / তুলনাযোগ্য জোড়া; ঝুঁকি সমান হলে আধা নম্বর। C = ০.৫ মানে মুদ্রা ছোড়া, ১ মানে নিখুঁত ক্রম। এটি AI-এর ভিত্তি কোর্সের পর্ব ৬০-এর ROC AUC-এর censoring-সহ আত্মীয়।

তৃতীয় ব্লক statsmodels-এর PHReg দিয়ে তিনটি বৈশিষ্ট্যে Cox model চালায় (default ties="breslow", এখানে স্পষ্ট লেখা), HR আর ৯৫% পরিসর ছাপায়, Harrell's C নিজে গোনে, আর শেষে একটি প্রশ্ন করে: কুপনের ক্রেতাদের HR কি সময়ের সঙ্গে স্থির? তার জন্য দুটি আলাদা fit: প্রথম ১৮০ দিন (১৮০-তে সবাইকে censor করে), আর ১৮০ দিনের পর যাঁরা ছিলেন তাঁদের নিয়ে। দ্বিতীয়টিতে entry বলে দেয় ক্রেতা ১৮০ দিন থেকে ঝুঁকির হিসাবে ঢুকেছেন (left truncation).

import numpy as np
import pandas as pd
from statsmodels.duration.hazard_regression import PHReg
from bazar_surv import customers

df = customers()
cols = ["plan", "complaints", "promo"]
X, T, E = df[cols].to_numpy(), df.time.to_numpy(), df.event.to_numpy()
res = PHReg(T, X, status=E, ties="breslow").fit()
ci = np.exp(res.conf_int())
print(pd.DataFrame({"HR": np.exp(res.params), "lo95": ci[:, 0], "hi95": ci[:, 1],
                    "p": res.pvalues}, index=cols).round(3))
print("truth: plan 0.600, complaints 1.350")

def c_index(time, event, risk):      # Harrell's C, নিজে লেখা
    num = den = 0.0
    for i in np.flatnonzero(event):  # i চলে গেছে; j তখনও ছিল
        later = time > time[i]
        den += later.sum()
        num += (risk[i] > risk[later]).sum() + 0.5 * (risk[i] == risk[later]).sum()
    return num / den

print(f"C-index, Cox xb:          {c_index(T, E, X @ res.params):.3f}")
print(f"C-index, complaints only: {c_index(T, E, X[:, 1]):.3f}")
print(f"C-index, random scores:   {c_index(T, E, np.random.default_rng(1).random(len(T))):.3f}")

# promo-র HR কি সময়ের সঙ্গে স্থির? 180 দিনের আগে আর পরে আলাদা fit
early = PHReg(np.minimum(T, 180), X, status=E * (T <= 180), ties="breslow").fit()
m = T > 180                          # 180 দিন পর্যন্ত যারা টিকে ছিল
late = PHReg(T[m], X[m], status=E[m], entry=np.full(m.sum(), 180.0),
             ties="breslow").fit()
print("promo HR, day 0-180:", round(np.exp(early.params[2]), 2),
      "| after day 180:", round(np.exp(late.params[2]), 2))

              HR   lo95   hi95      p
plan        0.60  0.493  0.730  0.000
complaints  1.37  1.242  1.511  0.000
promo       1.18  0.977  1.425  0.086
truth: plan 0.600, complaints 1.350
C-index, Cox xb:          0.619
C-index, complaints only: 0.575
C-index, random scores:   0.522
promo HR, day 0-180: 1.97 | after day 180: 0.63

৭ · অনুমান যাচাই আর parametric model — log(−log S), Weibull আর stratification

Proportional hazards অনুমান ভাঙলে HR একটি সময়-গড় মাত্র, আর তার ব্যাখ্যা ভুল দিকে যায়। ধরার তিনটি পথ। প্রথম, log(−log S) plot: PH সত্য হলে S₁ = S₀^HR, তাই log(−log S₁(t)) = log(−log S₀(t)) + log HR: দুই দলের Kaplan–Meier থেকে আঁকা বক্ররেখা সমান্তরাল হওয়ার কথা। Weibull-এ log(−log S) = k·log t − k·log λ, অর্থাৎ log t-এর বিপরীতে সরলরেখা, ঢাল k. ঢাল আলাদা মানে hazard-এর আকার আলাদা, PH ভাঙা। দ্বিতীয়, Schoenfeld residual: প্রতিটি ঘটনায় “যিনি চলে গেলেন তাঁর x” বিয়োগ “risk set-এর ঝুঁকি-ওজনের গড় x”. PH সত্য হলে এই Schoenfeld residual সময়ের সঙ্গে কোনো ধারা দেখায় না; ধারা থাকলে β সময়ের সঙ্গে বদলাচ্ছে। Patricia Grambsch আর Terry Therneau ১৯৯৪ সালে এর ওপর একটি test দেন। তৃতীয়, সময় ভাগ করে fit: তৃতীয় ব্লকের শেষ লাইনের মতো।

প্রতিকার কয়েকটি। (১) Stratification: কুপন আর সাধারণ ক্রেতার জন্য আলাদা baseline hazard h_0,s(t), কিন্তু অন্য বৈশিষ্ট্যের β একই। দাম: কুপনের নিজের কোনো HR পাওয়া যায় না, শুধু অন্যদের HR “কুপন সামলে”। (২) সময়ভেদে প্রভাব: x × g(t) পদ যোগ, বা সময় ভাগ করে আলাদা HR (যেমন ০–১৮০ আর ১৮০+). (৩) Parametric model, যেখানে hazard-এর আকার নিজেই শেখা হয়।

শেষটি দেখা যাক। Weibull-এর likelihood-এ censoring বসে স্বাভাবিকভাবে: ঘটনা দেখা গেলে density f(t) = h(t)S(t), censored হলে শুধু S(t) (“অন্তত t পর্যন্ত টিকেছেন”)। log নিলে

log L(k, λ) = ∑_i [ δ_i log h(t_i) − H(t_i) ],  H(t) = (t/λ)^k

এটি পর্ব ১৮-এর maximum likelihood, শুধু censored সারিতে density-র বদলে survival. পর্ব ২০-এর Poisson-এর সঙ্গে মিলটাও লক্ষ করুন: exponential-এ (k = ১) MLE λ̂ = ∑t_i / ∑δ_i, “মোট দেখার সময় ভাগ ঘটনার সংখ্যা”, ঠিক exposure-সহ হারের মতো। চতুর্থ ব্লক scipy-এর minimize দিয়ে দুই দলে Weibull শেখে (k আর λ লগ-স্কেলে, যাতে ধনাত্মক থাকে), median মেলায় Kaplan–Meier-এর সঙ্গে, censored সারি ফেলে দিলে কী হয় দেখায়, log(−log S) আঁকে, আর শেষে কুপন দিয়ে stratified Cox চালায়।

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from statsmodels.duration.survfunc import SurvfuncRight
from statsmodels.duration.hazard_regression import PHReg
from bazar_surv import customers

def weibull_nll(p, t, d):            # censoring-সহ: d log h(t) - H(t)
    k, lam = np.exp(p)               # লগ-স্কেলে, যাতে দুটিই ধনাত্মক
    z = t / lam
    return -np.sum(d * (np.log(k / lam) + (k - 1) * np.log(z)) - z**k)

df = customers()
print("group    k_hat  lam_hat  median_W  median_KM  median_dropC")
for g, name in [(0, "organic"), (1, "promo")]:
    s = df[df.promo == g]
    t, d = s.time.to_numpy(), s.event.to_numpy()
    k, lam = np.exp(minimize(weibull_nll, [0.0, 6.0], args=(t, d)).x)
    kd, ld = np.exp(minimize(weibull_nll, [0.0, 6.0],
                             args=(t[d == 1], d[d == 1])).x)
    sf = SurvfuncRight(t, d)
    print(f"{name:8} {k:5.2f} {lam:8.0f} {lam * np.log(2) ** (1 / k):9.0f}"
          f" {sf.quantile(0.5):10.0f} {ld * np.log(2) ** (1 / kd):11.0f}")
    keep = (sf.surv_prob > 0) & (sf.surv_prob < 1)
    plt.plot(np.log(sf.surv_times[keep]), np.log(-np.log(sf.surv_prob[keep])),
             drawstyle="steps-post", label=name)
plt.xlabel("log(days)"); plt.ylabel("log(-log S(t))"); plt.legend(); plt.show()

X = df[["plan", "complaints"]]
st = PHReg(df.time, X, status=df.event, strata=df.promo, ties="breslow").fit()
print("stratified by promo, HR:", np.exp(st.params).round(3))

group    k_hat  lam_hat  median_W  median_KM  median_dropC
organic   1.27      505       379        392         209
promo     0.65      513       292        282          94
stratified by promo, HR: [0.6   1.369]

৮ · ML survival model — random survival forest, lifelines আর scikit-survival

Cox model-এ x-এর প্রভাব log-hazard-এ রৈখিক আর সময়ে স্থির। Interaction বা বাঁকা সম্পর্ক থাকলে ML পদ্ধতি সাহায্য করতে পারে। Random survival forest (Hemant Ishwaran ও সহকর্মীরা, Annals of Applied Statistics, ২০০৮) পর্ব ৩৩-এর random forest-এর মতোই, শুধু গাছ ভাগ হয় log-rank পরিসংখ্যান দিয়ে আর প্রতিটি পাতায় থাকে একটি cumulative hazard-এর অনুমান। Gradient-boosted Cox পর্ব ৩৫-এর boosting-কে partial likelihood-এর loss দিয়ে চালায়; XGBoost-এ এর জন্য survival:cox objective আছে (পর্ব ৩৬)। Neural network-ভিত্তিক model-ও আছে (যেমন DeepSurv); সেগুলো deep learning-এর কোর্সে আসবে।

এগুলো পাতার Python-এ নেই, তাই ডেস্কটপে একটি আলাদা প্রকল্পে চালানো হলো। uv দিয়ে নতুন প্রকল্প (AI-এর ভিত্তি কোর্সের পর্ব ৭ দেখুন), একই bazar_surv.py. Python দুটি library-তে একই ডেটা: lifelines-এ Cox আর Grambsch–Therneau ধাঁচের proportional hazards test, তারপর scikit-survival-এ ৭০/৩০ ভাগে Cox বনাম random survival forest-এর test C-index.

~/bazarghar-survival$ uv init --bare --vcs none --author-from none
~/bazarghar-survival$ uv add numpy scipy pandas statsmodels matplotlib scikit-learn lifelines scikit-survival
 …
 + lifelines==0.30.3
 + pandas==2.3.3
 + scikit-learn==1.9.1
 + scikit-survival==0.28.0
 + scipy==1.18.1
 + statsmodels==0.15.0
 …

from lifelines import CoxPHFitter
from lifelines.statistics import proportional_hazard_test
from sksurv.util import Surv
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.ensemble import RandomSurvivalForest
from sklearn.model_selection import train_test_split
from bazar_surv import customers

df = customers()[["plan", "complaints", "promo", "time", "event"]]
cph = CoxPHFitter().fit(df, "time", "event")
print(cph.summary[["exp(coef)", "p"]].round(3))
pht = proportional_hazard_test(cph, df, time_transform="rank")
print(pht.summary[["test_statistic", "p"]].round(4))

tr, te = train_test_split(df, test_size=0.3, random_state=55)
y_tr = Surv.from_dataframe("event", "time", tr)
y_te = Surv.from_dataframe("event", "time", te)
X_tr, X_te = tr.iloc[:, :3], te.iloc[:, :3]
cox = CoxPHSurvivalAnalysis().fit(X_tr, y_tr)
rsf = RandomSurvivalForest(n_estimators=300, min_samples_leaf=15,
                           random_state=55).fit(X_tr, y_tr)
print("test C-index  Cox:", round(cox.score(X_te, y_te), 3),
      " RSF:", round(rsf.score(X_te, y_te), 3))

সত্যিকারের run (ডেস্কটপ, Python 3.14.4; লেখার সময়ে lifelines 0.30.3, scikit-survival 0.28.0, scikit-learn 1.9.1, pandas 2.3.3, NumPy 2.5.4):

~/bazarghar-survival$ uv run python desktop_survival.py
            exp(coef)      p
covariate
plan            0.600  0.000
complaints      1.371  0.000
promo           1.180  0.085
            test_statistic       p
complaints          3.2966  0.0694
plan                0.1039  0.7472
promo              47.1825  0.0000
test C-index  Cox: 0.634  RSF: 0.639

পড়া: lifelines-এর HR statsmodels-এর সঙ্গে মেলে (অভিযোগে ১.৩৭১ বনাম ১.৩৬৯-এর ছোট পার্থক্য tie সামলানোর ভিন্ন default থেকে: lifelines Efron, এখানে statsmodels Breslow). PH test কুপনে statistic ৪৭.২, p প্রায় শূন্য: অংশ ৭-এর ছবির সংখ্যাগত রূপ। Plus-এ কোনো লক্ষণ নেই (p = ০.৭৫); অভিযোগে p = ০.০৭, সীমার কাছে কিন্তু আসল model-এ এটি proportional, তাই এটি এলোমেলো ওঠানামা হিসেবেই পড়া উচিত। Forest বনাম Cox: test-এ ০.৬৩৯ বনাম ০.৬৩৪, একটি ভাগে একটি run. এই পার্থক্য সম্ভবত noise-এর মধ্যে; পর্ব ২২-এর মতো কয়েকটি ভাগ বা bootstrap ছাড়া “forest ভালো” বলা যায় না। ডেটা যেখানে প্রায় Cox-এর মতো, সেখানে forest-এর বাড়তি সুযোগ কম।

৯ · AI-এর যুগে — churn আর survival কোড যাচাই

“বাজারঘরের churn predict করো” বললে AI assistant প্রায়ই সোজা একটি classifier দেয়: label = “চলে গেছে কি না”, accuracy ৮০%, শেষ। এই পর্বের প্রথম ব্লক দেখিয়েছে, সেই label-এ কী মিশে থাকে। নিচের ভুলগুলো প্রায়ই চোখে পড়ে।

AI-এর দেওয়া survival কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• ৫–৬ জনের একটি খেলনা ডেটায় Kaplan–Meier হাতে করে কোডের সঙ্গে মিলিয়েছি
• event কলামে 1 মানে “ঘটনা ঘটেছে”, আর censored-রা বাদ পড়েনি, তা নিশ্চিত করেছি
• HR-এর ব্যাখ্যা hazard-এর ভাষায় লিখেছি, আর কারণের দাবি করিনি
• Proportional hazards যাচাই করেছি, আর ভাঙলে stratification বা সময়ভেদে প্রভাব ভেবেছি
• C-index আলাদা test সেটে আর একটি সরল baseline-এর সঙ্গে, কয়েকটি seed বা ভাগে তুলনা করেছি
• কোনো আসল ক্রেতার নাম, ফোন নম্বর বা অর্ডারের ইতিহাস AI-তে পেস্ট করিনি