# Aspirin RCT (toy simulation) — minimal, vectorized

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import chi2_contingency

np.random.seed(42)
n = 10_000

# Simulate participants
df = pd.DataFrame({
    "age": np.random.normal(55, 12, n),
    "risk_factor": np.random.uniform(0, 1, n),
    "treatment": np.random.choice(["aspirin", "placebo"], size=n)
})

# Baseline risk and treatment effect (aspirin RR = 0.8)
base = 0.02 + 0.03 * (df["age"] - 40) / 40 + 0.05 * df["risk_factor"]
rr   = np.where(df["treatment"].eq("aspirin"), 0.8, 1.0)
p    = np.clip(base * rr, 0, 1)

# Cancer outcome
df["cancer"] = (np.random.rand(n) < p).astype(int)

# Summary by arm
summary = df.groupby("treatment")["cancer"].agg(cases="sum", total="count", rate="mean")
print("\n=== RCT RESULTS ===")
print(summary)

# Effect size (relative risk)
asp, plc = summary.loc["aspirin","rate"], summary.loc["placebo","rate"]
RR  = asp / plc
RRR = 1 - RR
print(f"\nRelative Risk (RR): {RR:.3f}")
print(f"Risk Reduction (1−RR): {RRR*100:.1f}%")

# Minimal visualization
plt.figure(figsize=(5,4))
plt.bar(["Aspirin","Placebo"], [asp, plc], alpha=0.8)
plt.ylabel("Cancer rate")
plt.title("Cancer rate by treatment")
for i, v in enumerate([asp, plc]):
    plt.text(i, v + 0.001, f"{v:.3f}", ha="center")
plt.ylim(0, max(asp, plc)*1.2)
plt.tight_layout()
plt.show()

=== RCT RESULTS ===
           cases  total      rate
treatment                        
aspirin      243   5087  0.047769
placebo      257   4913  0.052310

Relative Risk (RR): 0.913
Risk Reduction (1−RR): 8.7%

fig, ax = plt.subplots(figsize=(8, 1.2))
left = 0.0
for wi in w:
    ax.barh(0, wi, left=left, height=1, alpha=0.3)
    ax.text(left + wi/2, 0, f"{int(wi)}", ha='center', va='center')
    left += wi
ax.set_xlim(0, T); ax.set_yticks([]); ax.set_xlabel("minutes (total 40)")
ax.set_title("One long gap (37) + three short gaps (1,1,1)")
plt.tight_layout(); plt.show()

# Inspection paradox in a 40-min window with 4 buses: [37, 1, 1, 1]
import numpy as np
import matplotlib.pyplot as plt

w = np.array([37, 1, 1, 1], dtype=float)
T = w.sum()

# --- Naive unweighted observed samples (equal-probability buses) ---
naive_wait = (w.mean()) / 2

# --- Student's biased observation (length-biased arrivals) ---
p_len = w / T                          # P(arrive in interval i)
student_wait = (w/2 @ p_len)           # E[wait] under length bias

# ---Importance / "inverse-weighted" check ---
# p_x: equal-probability buses; p_y: length-biased over intervals
p_x = np.full(len(w), 1/len(w))
p_y = (p_x * w) / (p_x @ w)            # ∝ w; normalizes to length-bias

print("Purdue's collected expected wait (mean-interval/2):        {:.2f} min".format(naive_wait))
print("Students' observed expected wait (length-biased):          {:.2f} min".format(student_wait))

# --- Monte Carlo (simulated) sanity check (uniform arrival in [0, T]) ---
rng = np.random.default_rng(0)
N = 100_000
arrivals = rng.uniform(0, T, size=N)
edges = np.concatenate([[0], np.cumsum(w)])
idx = np.searchsorted(edges, arrivals, side='right') - 1
mc_wait = edges[idx+1] - arrivals
print("Simulated (Monte Carlo) E[wait] (uniform in [0,T]):        {:.2f} min".format(mc_wait.mean()))

Purdue's collected expected wait (mean-interval/2):        5.00 min
Students' observed expected wait (length-biased):          17.15 min
Simulated (Monte Carlo) E[wait] (uniform in [0,T]):        17.17 min

# --- Importance weighting to recover the unweighted (equal-bus) average from MC samples ---
# We sampled arrivals ~ uniform on [0,T] ⇒ interval index 'idx' is drawn from p_y(i) ∝ w_i
# Weight each sample by p_x(i)/p_y(i) ∝ 1/w_i to estimate the equal-bus expectation.
imp_raw = 1.0 / w[idx]                     # ∝ p_x(idx)/p_y(idx)
imp_w   = imp_raw / imp_raw.sum()          # self-normalized weights
naive_from_MC = np.sum(imp_w * mc_wait)    # ≈ naive_wait

print("Unweighted bus-average via Simulation + Inverse weighted estimation:          {:.2f} min".format(naive_from_MC))

Unweighted bus-average via Simulation + Inverse weighted estimation:          5.04 min

Age	Total Conditional Individuals (approx.)	Not Vax %	Fully Vax %	Severe Cases Not Vax per 100k	Severe Cases Fully Vax per 100k	Score function Conditional Efficacy
[12,15]	650,000	62.1%	29.9%	0.3	0.0	100.0%
[16,19]	600,000	21.9%	73.5%	1.6	0.0	100.0%
[20,29]	1,200,000	20.5%	76.2%	1.5	0.0	100.0%
[30,39]	1,050,000	16.2%	80.9%	6.2	0.2	96.8%
[40,49]	900,000	13.2%	84.4%	16.5	1.0	94.2%
[50,59]	750,000	10.0%	88.0%	40.2	2.9	93.2%
[60,69]	550,000	8.8%	89.8%	76.6	8.7	89.8%
[70,79]	350,000	4.2%	94.6%	190.1	19.8	90.6%
[80,89]	120,000	5.6%	92.6%	252.3	47.9	84.0%
90+	50,000	6.1%	90.5%	510.9	38.6	93.0%

Data Mining & Machine Learning¶

CS37300¶

Purdue University¶

Learning Outcomes¶

Primer on hypotheses¶

The Scientific Method¶

What is a hypothesis?¶

Types of hypotheses¶

Broad categories¶

Specific characteristics¶

From claims to testable hypotheses¶

From claims to testable hypotheses¶

From claims to testable hypotheses¶

From claims to testable hypotheses¶

Using Data to Test Hypotheses¶

Example:¶

Directional Relational Hypothesis¶

Causal Claims without Experiments are Possible¶

Is Aspirin effective in reducing cancer risk?¶

Q3: Traditional Gender Roles¶

The Strong Effect of Data Biases¶

Observational Data Biases¶

Example of data bias: Inspection Paradox¶

Data bias explanation¶

Fixing data biases: Inverse weighted estimator¶

Summary¶

References¶