{
"cells": [
{
"cell_type": "markdown",
"id": "md-intro",
"metadata": {},
"source": [
"# Case Study: Cox Frailty Model for Kidney Stone Recurrence\n",
"\n",
"This notebook walks through a realistic multi-center clinical trial analysis using `interlace.coxme()` — a Cox proportional hazards model with Gaussian frailty.\n",
"\n",
"## Why Cox frailty over standard Cox PH?\n",
"\n",
"Standard Cox PH assumes all patients share the same baseline hazard, ignoring clustering. When patients are nested in centers (hospitals, clinics), unobserved center-level differences inflate the apparent treatment effect and produce overconfident standard errors. A **frailty model** adds a random intercept per center on the log-hazard scale, capturing this heterogeneity explicitly.\n",
"\n",
"| Feature | Standard Cox | Cox Frailty (`coxme`) |\n",
"|---|---|---|\n",
"| Handles clustering | No | Yes |\n",
"| Center-level BLUPs | No | Yes |\n",
"| Correct SE under clustering | No | Yes |\n",
"| Variance components | No | Yes |\n",
"\n",
"## Dataset: Kidney Stone Recurrence Trial\n",
"\n",
"- **300 patients** across **20 clinical centers** (15 per center)\n",
"- **Outcome**: time-to-recurrence in months, right-censored\n",
"- **Predictors**: `treatment` (0=standard, 1=new), `age_c` (age mean-centred at 58, divided by 10)\n",
"- **Random effect**: `center` (frailty, Gaussian on log-hazard scale)\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "code-imports-data",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-05T08:16:24.847047Z",
"iopub.status.busy": "2026-05-05T08:16:24.846954Z",
"iopub.status.idle": "2026-05-05T08:16:25.960135Z",
"shell.execute_reply": "2026-05-05T08:16:25.959721Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Observations : 300\n",
"Centers : 20\n",
"Event rate : 57.0%\n",
"\n",
"Event rate by treatment:\n",
" events rate\n",
"trt_label \n",
"New treatment 72 0.464516\n",
"Standard 99 0.682759\n"
]
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import interlace\n",
"\n",
"rng = np.random.default_rng(99)\n",
"\n",
"TRUE_LAMBDA0 = 0.04\n",
"TRUE_BETA_TRT = np.log(0.65) # -0.431\n",
"TRUE_BETA_AGE = np.log(1.25) # 0.223 (per 10 years)\n",
"TRUE_SD_CENTER = 0.5\n",
"\n",
"n_centers = 20\n",
"n_per_center = 15\n",
"\n",
"u_center = rng.normal(0, TRUE_SD_CENTER, n_centers)\n",
"\n",
"rows = []\n",
"for j in range(n_centers):\n",
" for i in range(n_per_center):\n",
" trt = int(rng.binomial(1, 0.5))\n",
" age_raw = rng.normal(58, 12)\n",
" age_c = (age_raw - 58) / 10\n",
"\n",
" h = TRUE_LAMBDA0 * np.exp(TRUE_BETA_TRT * trt + TRUE_BETA_AGE * age_c + u_center[j])\n",
" surv_time = rng.exponential(1 / h)\n",
" cens_time = rng.uniform(12, 36)\n",
" obs_time = min(surv_time, cens_time)\n",
" event = int(surv_time <= cens_time)\n",
"\n",
" rows.append({\n",
" 'patient': f'P{j*n_per_center+i+1:03d}',\n",
" 'center': f'C{j+1:02d}',\n",
" 'treatment': trt,\n",
" 'age': round(age_raw, 1),\n",
" 'age_c': round(age_c, 3),\n",
" 'time': round(obs_time, 2),\n",
" 'event': event,\n",
" })\n",
"\n",
"df = pd.DataFrame(rows)\n",
"df['trt_label'] = df['treatment'].map({0: 'Standard', 1: 'New treatment'})\n",
"\n",
"print(f\"Observations : {len(df)}\")\n",
"print(f\"Centers : {df['center'].nunique()}\")\n",
"print(f\"Event rate : {df['event'].mean():.1%}\")\n",
"print()\n",
"print(\"Event rate by treatment:\")\n",
"print(df.groupby('trt_label')['event'].agg(['sum', 'mean']).rename(columns={'sum': 'events', 'mean': 'rate'}))\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "code-head",
"metadata": {
"execution": {
"iopub.execute_input": "2026-05-05T08:16:25.962460Z",
"iopub.status.busy": "2026-05-05T08:16:25.962224Z",
"iopub.status.idle": "2026-05-05T08:16:25.968550Z",
"shell.execute_reply": "2026-05-05T08:16:25.968180Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" patient | \n",
" center | \n",
" treatment | \n",
" age | \n",
" time | \n",
" event | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" P001 | \n",
" C01 | \n",
" 0 | \n",
" 53.8 | \n",
" 33.36 | \n",
" 0 | \n",
"
\n",
" \n",
" | 1 | \n",
" P002 | \n",
" C01 | \n",
" 0 | \n",
" 39.8 | \n",
" 14.98 | \n",
" 1 | \n",
"
\n",
" \n",
" | 2 | \n",
" P003 | \n",
" C01 | \n",
" 1 | \n",
" 57.0 | \n",
" 8.43 | \n",
" 1 | \n",
"
\n",
" \n",
" | 3 | \n",
" P004 | \n",
" C01 | \n",
" 0 | \n",
" 50.3 | \n",
" 1.64 | \n",
" 1 | \n",
"
\n",
" \n",
" | 4 | \n",
" P005 | \n",
" C01 | \n",
" 1 | \n",
" 68.1 | \n",
" 13.38 | \n",
" 1 | \n",
"
\n",
" \n",
" | 5 | \n",
" P006 | \n",
" C01 | \n",
" 0 | \n",
" 75.7 | \n",
" 25.74 | \n",
" 0 | \n",
"
\n",
" \n",
" | 6 | \n",
" P007 | \n",
" C01 | \n",
" 1 | \n",
" 78.4 | \n",
" 30.41 | \n",
" 1 | \n",
"
\n",
" \n",
" | 7 | \n",
" P008 | \n",
" C01 | \n",
" 1 | \n",
" 54.1 | \n",
" 2.87 | \n",
" 1 | \n",
"
\n",
" \n",
"
\n",
"