""" 04-12 Capstone — Business Statistics Project Every chapter, on one data set. Run from the project folder: python python/analysis.py """ import numpy as np import pandas as pd from scipy import stats import statsmodels.api as sm from statsmodels.formula.api import ols from statsmodels.stats.multicomp import pairwise_tukeyhsd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt pd.set_option("display.width", 130) pd.set_option("display.max_columns", 30) df = pd.read_csv("data/customers.csv") n = len(df) spend = df["annual_spend"] def hr(title): print("\n" + "=" * 74) print(title) print("=" * 74) # ════════════════════════════════════════════════════════════════════════════ hr("PART 1 — ORGANIZE & DESCRIBE") # ════════════════════════════════════════════════════════════════════════════ print("\n-- data quality --") print(f"rows {n}") print(f"duplicate ids {df['customer_id'].duplicated().sum()}") print(f"missing values {df.isna().sum().sum()}") print(f"annual_spend range {spend.min():.2f} to {spend.max():.2f}") print(f"visits range {df['visits_per_year'].min()} to " f"{df['visits_per_year'].max()}") print("\ncategory levels:") for c in ("region", "channel", "satisfied", "returned_item"): print(f" {c:<15} {', '.join(sorted(df[c].unique()))}") print("\n-- annual spend --") q1, med, q3 = spend.quantile([0.25, 0.5, 0.75]) sd = spend.std(ddof=1) desc = pd.Series({ "n": n, "mean": spend.mean(), "median": med, "sd": sd, "cv_pct": 100 * sd / spend.mean(), "min": spend.min(), "Q1": q1, "Q3": q3, "max": spend.max(), "IQR": q3 - q1, "skew_SK": 3 * (spend.mean() - med) / sd, }) print(desc.round(3).to_string()) sk = desc["skew_SK"] shape = ("roughly SYMMETRIC" if abs(sk) < 0.5 else "RIGHT-SKEWED" if sk > 0 else "LEFT-SKEWED") print(f"\nshape: {shape} -> report the " + ("MEAN and SD" if shape == "roughly SYMMETRIC" else "MEDIAN and IQR")) iqr = q3 - q1 lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr print(f"outlier fences: {lo:.2f} to {hi:.2f} -> " f"{((spend < lo) | (spend > hi)).sum()} outlier(s)") print("\n-- by region --") print(df.groupby("region")["annual_spend"] .agg(n="count", mean="mean", sd=lambda x: x.std(ddof=1), median="median") .round(2).to_string()) print("\n-- channel x satisfaction --") print(pd.crosstab(df["channel"], df["satisfied"], margins=True).to_string()) # ════════════════════════════════════════════════════════════════════════════ hr("PART 2 — ESTIMATE") # ════════════════════════════════════════════════════════════════════════════ ci_lo, ci_hi = stats.t.interval(0.95, n - 1, spend.mean(), stats.sem(spend)) print(f"\nMean annual spend ${spend.mean():.2f}") print(f"95% CI ${ci_lo:.2f} to ${ci_hi:.2f}") print(f"margin of error ${(ci_hi - ci_lo) / 2:.2f}") print('\n=> "We are 95% confident that the mean annual spend of all customers') print(f' lies between ${ci_lo:.2f} and ${ci_hi:.2f}."') x = int((df["satisfied"] == "Yes").sum()) ph = x / n se_p = np.sqrt(ph * (1 - ph) / n) cp_lo, cp_hi = ph - 1.959964 * se_p, ph + 1.959964 * se_p print(f"\nProportion satisfied {ph:.4f} ({x} of {n})") print(f"condition check n*p = {n*ph:.0f}, n*(1-p) = {n*(1-ph):.0f} (both >= 5)") print(f"95% CI {cp_lo:.4f} to {cp_hi:.4f} " f"({100*cp_lo:.1f}% to {100*cp_hi:.1f}%)") need = int(np.ceil(ph * (1 - ph) * (1.959964 / 0.02) ** 2)) print(f"\nFor a +/-2% margin we would need n = {need} (currently {n})") # ════════════════════════════════════════════════════════════════════════════ hr("PART 3 — TEST") # ════════════════════════════════════════════════════════════════════════════ # ---- TEST 1 --------------------------------------------------------------- print("\n--- TEST 1: does mean spend exceed the $250 target? ---") print("H0: mu <= 250 H1: mu > 250 right-tailed, alpha = 0.05\n") t1 = stats.ttest_1samp(spend, 250, alternative="greater") d1 = (spend.mean() - 250) / sd print(f"t({n-1}) = {t1.statistic:.4f}, p = {t1.pvalue:.6f}, Cohen's d = {d1:.4f}") print(f"two-sided 95% CI: ${ci_lo:.2f} to ${ci_hi:.2f}") print(f"=> {'REJECT H0' if t1.pvalue <= 0.05 else 'FAIL TO REJECT H0'}. " f"Mean spend {'EXCEEDS' if t1.pvalue <= 0.05 else 'is not shown to exceed'} " "the $250 target.") # ---- TEST 2 --------------------------------------------------------------- print("\n--- TEST 2: does mean spend differ across regions? ---") print("H0: mu_N = mu_S = mu_E = mu_W H1: at least one differs\n") groups = [g["annual_spend"].to_numpy() for _, g in df.groupby("region")] region_names = list(df.groupby("region").groups.keys()) vars_ = [g.var(ddof=1) for g in groups] print(f"variance ratio (max/min) = {max(vars_)/min(vars_):.3f} -> " + ("pooled ANOVA is appropriate" if max(vars_) / min(vars_) < 4 else "unequal; prefer Welch")) print(f"Levene test p = {stats.levene(*groups).pvalue:.4f}") aov = ols("annual_spend ~ C(region)", data=df).fit() atab = sm.stats.anova_lm(aov, typ=2) print(atab.round(4).to_string()) Fp = atab["PR(>F)"].iloc[0] eta2 = atab["sum_sq"].iloc[0] / atab["sum_sq"].sum() print(f"\neta squared = {eta2:.4f} -> {100*eta2:.1f}% of spending variation " "explained by region") print(f"Shapiro-Wilk on residuals: p = {stats.shapiro(aov.resid).pvalue:.4f}") if Fp <= 0.05: print("\nF IS significant -> post-hoc comparisons:") tukey = pairwise_tukeyhsd(df["annual_spend"], df["region"], alpha=0.05) print(tukey) else: tukey = None print("\nF is NOT significant -> no post-hoc test.") # ---- TEST 3 --------------------------------------------------------------- print("\n--- TEST 3: are channel and satisfaction associated? ---") print("H0: channel and satisfaction are INDEPENDENT H1: associated\n") tab = pd.crosstab(df["channel"], df["satisfied"]) chi2, pchi, dof, expected = stats.chi2_contingency(tab) expected = pd.DataFrame(expected, index=tab.index, columns=tab.columns) print("observed:") print(pd.crosstab(df["channel"], df["satisfied"], margins=True).to_string()) print("\nexpected:") print(expected.round(2).to_string()) print(f"\nall expected >= 5? {(expected.values >= 5).all()}") print("\ncell contributions:") print((((tab - expected) ** 2) / expected).round(4).to_string()) print(f"\nchi-square = {chi2:.4f}, df = {dof}, p = {pchi:.4f}") V = np.sqrt(chi2 / (n * (min(tab.shape) - 1))) print(f"Cramer's V = {V:.4f}") print("\nsatisfaction rate by channel:") print((100 * tab.div(tab.sum(axis=1), axis=0)).round(1).to_string()) print("=> " + ("REJECT H0. Channel and satisfaction ARE associated." if pchi <= 0.05 else "FAIL TO REJECT H0. No evidence of an association.")) # ════════════════════════════════════════════════════════════════════════════ hr("PART 4 — MODEL") # ════════════════════════════════════════════════════════════════════════════ pr = stats.pearsonr(df["visits_per_year"], spend) print("\ncorrelation of spend with visits:") print(f" r = {pr.statistic:.4f}, r^2 = {pr.statistic**2:.4f}, p = {pr.pvalue:.3e}") print("\n--- simple regression: annual_spend ~ visits_per_year ---") s1 = ols("annual_spend ~ visits_per_year", data=df).fit() print(pd.DataFrame({"estimate": s1.params, "std_err": s1.bse, "t": s1.tvalues, "p_value": s1.pvalues}).round(4).to_string()) print(f"R^2 = {s1.rsquared:.4f} s_e = {np.sqrt(s1.mse_resid):.2f}") c1 = s1.conf_int().loc["visits_per_year"] print(f"each extra visit is worth ${s1.params['visits_per_year']:.2f} " f"(95% CI: ${c1[0]:.2f} to ${c1[1]:.2f})") print("\n--- multiple regression: + channel + region ---") s2 = ols("annual_spend ~ visits_per_year + C(channel) + C(region)", data=df).fit() print(pd.DataFrame({"estimate": s2.params, "std_err": s2.bse, "t": s2.tvalues, "p_value": s2.pvalues}).round(4).to_string()) print(f"\nR^2 = {s2.rsquared:.4f} Adjusted R^2 = {s2.rsquared_adj:.4f} " f"s_e = {np.sqrt(s2.mse_resid):.2f}") print(f"reference levels: channel = {sorted(df['channel'].unique())[0]}, " f"region = {sorted(df['region'].unique())[0]}") print("\nmodel comparison (adjusted R^2 and AIC, never raw R^2):") print(pd.DataFrame({ "model": ["visits only", "+ channel + region"], "adj_r2": [round(s1.rsquared_adj, 4), round(s2.rsquared_adj, 4)], "AIC": [round(s1.aic, 2), round(s2.aic, 2)], }).to_string(index=False)) print("\nnested F-test:") print(sm.stats.anova_lm(s1, s2).round(4).to_string()) print("\ndiagnostics on the full model:") print(f" Shapiro-Wilk on residuals p = {stats.shapiro(s2.resid).pvalue:.4f}") print(f" cor(|resid|, fitted) " f"{np.corrcoef(np.abs(s2.resid), s2.fittedvalues)[0,1]:+.4f} " "(near 0 = equal variance)") print(f" Durbin-Watson {sm.stats.durbin_watson(s2.resid):.4f}") cooks = s2.get_influence().cooks_distance[0] print(f" influential points {(cooks > 4/n).sum()}") new = pd.DataFrame({"visits_per_year": [12], "channel": [sorted(df["channel"].unique())[0]], "region": [sorted(df["region"].unique())[0]]}) pf = s2.get_prediction(new).summary_frame(alpha=0.05) print(f"\nprediction for a 12-visit, {new['channel'][0]}, {new['region'][0]} customer:") print(f" point estimate ${pf['mean'].iloc[0]:,.2f}") print(f" 95% CI (mean) ${pf['mean_ci_lower'].iloc[0]:,.2f} to " f"${pf['mean_ci_upper'].iloc[0]:,.2f}") print(f" 95% PI (one person) ${pf['obs_ci_lower'].iloc[0]:,.2f} to " f"${pf['obs_ci_upper'].iloc[0]:,.2f}") print("Use the PREDICTION interval when the question is about ONE customer.") # ════════════════════════════════════════════════════════════════════════════ hr("PART 5 — EXECUTIVE SUMMARY") # ════════════════════════════════════════════════════════════════════════════ print(f"\n1. The average customer spends ${spend.mean():.0f} a year " f"(95% CI: ${ci_lo:.0f}-${ci_hi:.0f}),") print(f" with {100*ph:.0f}% satisfied (95% CI: {100*cp_lo:.0f}%-{100*cp_hi:.0f}%).") print(f"\n2. Mean spend " f"{'EXCEEDS' if t1.pvalue <= 0.05 else 'does NOT clearly exceed'} " f"the $250 target (one-sample t, p = {t1.pvalue:.4f}).") print(f"\n3. Regional spending " f"{'DIFFERS significantly' if Fp <= 0.05 else 'shows no significant difference'} " f"(ANOVA p = {Fp:.4f}, eta-sq = {eta2:.3f}).") print(f"\n4. Channel and satisfaction " f"{'ARE associated' if pchi <= 0.05 else 'are NOT associated'} " f"(chi-square p = {pchi:.4f}, V = {V:.3f}).") print(f"\n5. Visits, channel and region explain {100*s2.rsquared_adj:.0f}% of " "spending variation;") print(f" a typical prediction is within about ${np.sqrt(s2.mse_resid):.0f}.") print("\nLIMITATION: this is OBSERVATIONAL data. Every finding is an ASSOCIATION.") print("Only a randomized experiment would license a causal claim.") # ── PLOTS ─────────────────────────────────────────────────────────────────── fig, axes = plt.subplots(2, 3, figsize=(15, 9)) ax = axes.ravel() ax[0].hist(spend, bins=25, color="#8A5FBF", edgecolor="white") ax[0].axvline(spend.mean(), color="#0FA3A3", lw=2) ax[0].axvline(med, color="#0B7A7A", lw=2, ls="--") ax[0].set(title="Annual spend", xlabel="$") ax[1].boxplot(groups, tick_labels=region_names) ax[1].set(title="Spend by region", ylabel="$") ax[1].tick_params(axis="x", rotation=30) ch_groups = [g["annual_spend"].to_numpy() for _, g in df.groupby("channel")] ch_names = list(df.groupby("channel").groups.keys()) ax[2].boxplot(ch_groups, tick_labels=ch_names) ax[2].set(title="Spend by channel", ylabel="$") ax[2].tick_params(axis="x", rotation=30) sat_rate = 100 * tab["Yes"] / tab.sum(axis=1) ax[3].bar(sat_rate.index, sat_rate.values, color="#0B7A7A") ax[3].set(title="Satisfaction rate by channel", ylabel="% satisfied", ylim=(0, 100)) ax[3].tick_params(axis="x", rotation=30) ax[4].scatter(df["visits_per_year"], spend, color="#5B2A86") xs = np.linspace(df["visits_per_year"].min(), df["visits_per_year"].max(), 50) ax[4].plot(xs, s1.params["Intercept"] + s1.params["visits_per_year"] * xs, color="#0FA3A3", lw=2) ax[4].set(title="Spend vs visits", xlabel="visits per year", ylabel="$") ax[5].scatter(s2.fittedvalues, s2.resid, color="#0B7A7A") ax[5].axhline(0, ls="--", color="grey") ax[5].set(title="Full model: residuals vs fitted", xlabel="fitted", ylabel="residual") plt.tight_layout() plt.savefig("capstone_plots.png", dpi=110) print("\nWrote capstone_plots.png")