Hypothesis Testing Examples¶

import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.proportion import proportions_ztest
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import matplotlib.pyplot as plt
import seaborn as sns
pd.set_option('display.float_format', lambda x: f"{x:,.6f}")

from statsmodels.stats.proportion import proportion_confint

# Load dataset (raw CSV via GitHub)
url_ab = "https://raw.githubusercontent.com/beery4010/Analyze-AB-Test-Results/master/ab_data.csv"
ab = pd.read_csv(url_ab)
ab

# Align groups and landing pages as done in the Udacity project
ab = ab[((ab['group'] == 'control') & (ab['landing_page'] == 'old_page')) |
        ((ab['group'] == 'treatment') & (ab['landing_page'] == 'new_page'))]

# Compute conversions and sample sizes
conv_control = ab.loc[ab['group']=='control','converted'].sum()
n_control    = ab.loc[ab['group']=='control','converted'].count()
conv_treat   = ab.loc[ab['group']=='treatment','converted'].sum()
n_treat      = ab.loc[ab['group']=='treatment','converted'].count()
print('In our control group, we had', n_control,' out of whom ', conv_control, 'converted, making a probability of p=', round(conv_control/n_control, 4))
print("In our treatment group we had", n_treat, " out of which ", conv_treat, 'converted., making a probability of p=', round(conv_treat/n_treat, 4))

In our control group, we had 145274  out of whom  17489 converted, making a probability of p= 0.1204
In our treatment group we had 145311  out of which  17264 converted., making a probability of p= 0.1188

# Two‑sided two‑proportion z‑test
z_stat, p_val = proportions_ztest([conv_treat, conv_control], [n_treat, n_control], alternative='two-sided')
print(f"Two‑proportion z‑test: z = {z_stat:.4f}, p = {p_val:.6f}")

Two‑proportion z‑test: z = -1.3116, p = 0.189653

# Wald 95% CI for each proportion
ci_treat = proportion_confint(conv_treat, n_treat, method='normal')
ci_ctrl  = proportion_confint(conv_control, n_control, method='normal')

# Difference in proportions & CI (Wald)
p_hat_t = conv_treat / n_treat
p_hat_c = conv_control / n_control
diff = p_hat_t - p_hat_c
se = np.sqrt(p_hat_t*(1-p_hat_t)/n_treat + p_hat_c*(1-p_hat_c)/n_control)
ci_diff = (diff - 1.96*se, diff + 1.96*se)

print(f"Counts (treat/control): {conv_treat}/{n_treat} vs {conv_control}/{n_control}")
print(f"p_treat = {p_hat_t:.5f} 95% CI {ci_treat}")
print(f"p_ctrl  = {p_hat_c:.5f} 95% CI {ci_ctrl}")
print(f"Diff (treat - control) = {diff:.5f} 95% CI {ci_diff}")

Counts (treat/control): 17264/145311 vs 17489/145274
p_treat = 0.11881 95% CI (0.11714362162601945, 0.12047087417952866)
p_ctrl  = 0.12039 95% CI (0.11871294722381814, 0.12205966177710426)
Diff (treat - control) = -0.00158 95% CI (-0.003938713688889012, 0.0007806004935147211)

def plot_z_hypothesis(
    data_list,
    pop_mean=0.0,
    pop_sd=1.0,
    alternative='two.sided',   # 'two.sided' | 'greater' | 'less'
    type_test = 'mean', #'mean' | 'prob'
    alpha=0.05,
    label='Sampling distribution',
    title='z-test (sampling distribution of x̄)',
    x_label = "Sampling distribution of x̄",
    figsize=(8, 5)
):
    """
    Visualize z-test decision regions for the sampling distribution of the sample mean.

    Parameters
    ----------
    data_list : array-like
        Sample values (used only to compute x̄ and n).
    pop_mean : float
        Hypothesized population mean (μ under H0).
    pop_sd : float
        Known population standard deviation (σ).
    alternative : str
        'two.sided', 'greater', or 'less' (same semantics as your R code).
    alpha : float
        Significance level for critical regions.
    label : str
        Label for the distribution curve.
    title : str
        Plot title.
    figsize : tuple
        Figure size.

    Returns
    -------
    fig, ax : matplotlib Figure and Axes
    """
    x = np.asarray(data_list)
    n = len(x)
    xbar = np.mean(x)

    # Standard error of the mean
    se = pop_sd / np.sqrt(n)

    # Range for plotting ±4 SE around μ (as in your R function)
    grid = np.linspace(pop_mean - 4 * se, pop_mean + 4 * se, 4001)
    pdf = stats.norm.pdf(grid, loc=pop_mean, scale=se)

    # Compute critical cutoffs under H0 for the chosen alternative
    if alternative == 'two.sided':
        # symmetric cutoffs: (alpha/2) and (1 - alpha/2)
        lower_cut = stats.norm.ppf(alpha / 2, loc=pop_mean, scale=se)
        upper_cut = stats.norm.ppf(1 - alpha / 2, loc=pop_mean, scale=se)

        # Retain between cutoffs; reject outside
        retain_mask = (grid >= lower_cut) & (grid <= upper_cut)
        reject_mask = ~retain_mask

    elif alternative == 'greater':
        # reject on right tail
        cutoff = stats.norm.ppf(1 - alpha, loc=pop_mean, scale=se)
        retain_mask = (grid <= cutoff)
        reject_mask = (grid > cutoff)
        lower_cut, upper_cut = None, cutoff

    elif alternative == 'less':
        # reject on left tail
        cutoff = stats.norm.ppf(alpha, loc=pop_mean, scale=se)
        retain_mask = (grid >= cutoff)
        reject_mask = (grid < cutoff)
        lower_cut, upper_cut = cutoff, None
    else:
        raise ValueError("alternative must be one of {'two.sided','greater','less'}")

    # Build a DataFrame like the R pipeline for convenience (optional)
    df = pd.DataFrame({'x': grid, 'pdf': pdf, 'retain': retain_mask})

    # Plot
    fig, ax = plt.subplots(figsize=figsize)
    ax.plot(df['x'], df['pdf'], color='black', lw=1.2, label=label)

    # Shade retain region
    ax.fill_between(df['x'], 0, df['pdf'], where=df['retain'], color='#69b3a2', alpha=0.4, label='Retain H₀')

    # Shade reject region(s)
    ax.fill_between(df['x'], 0, df['pdf'], where=~df['retain'], color='#e76f51', alpha=0.4, label='Reject H₀')

    # Critical lines
    if alternative == 'two.sided':
        ax.axvline(lower_cut, color='#e76f51', ls='--', lw=1)
        ax.axvline(upper_cut, color='#e76f51', ls='--', lw=1)
        ax.text(lower_cut, ax.get_ylim()[1]*0.3, f"Lower crit\n{lower_cut:.2f}", ha='right', va='top', fontsize=9)
        ax.text(upper_cut, ax.get_ylim()[1]*0.3, f"Upper crit\n{upper_cut:.2f}", ha='left', va='top', fontsize=9)
    elif alternative == 'greater':
        ax.axvline(upper_cut, color='#e76f51', ls='--', lw=1)
        ax.text(upper_cut, ax.get_ylim()[1]*0.9, f"Crit\n{upper_cut:.2f}", ha='left', va='top', fontsize=9)
    elif alternative == 'less':
        ax.axvline(lower_cut, color='#e76f51', ls='--', lw=1)
        ax.text(lower_cut, ax.get_ylim()[1]*0.9, f"Crit\n{lower_cut:.2f}", ha='right', va='top', fontsize=9)

    # x̄ line and annotation
    if(type_test == 'prob'):
        ax.axvline(xbar, color='#264653', lw=2, ls='-', label='z')
        ax.annotate(f"z = {xbar:.2f}",
                xy=(xbar, stats.norm.pdf(xbar, loc=pop_mean, scale=se)),
                xytext=(xbar, ax.get_ylim()[1]*0.6),
                arrowprops=dict(arrowstyle='->', color='#264653'),
                ha='center', color='#264653')
    else:
        ax.axvline(xbar, color='#264653', lw=2, ls='-', label='Sample mean (x̄)')
        ax.annotate(f"x̄ = {xbar:.2f}",
                xy=(xbar, stats.norm.pdf(xbar, loc=pop_mean, scale=se)),
                xytext=(xbar, ax.get_ylim()[1]*0.6),
                arrowprops=dict(arrowstyle='->', color='#264653'),
                ha='center', color='#264653')

    ax.set_title(title)
    ax.set_xlabel(x_label)
    ax.set_ylabel("Density")
    ax.legend(loc='upper right', frameon=False)
    ax.grid(alpha=0.15)
    plt.show();
plot_z_hypothesis([z_stat], type_test = 'prob', title='Two proportion z-test', label = 'Standard Normal', x_label = '')

# Load ToothGrowth
url_tg = "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/ToothGrowth.csv"
tg = pd.read_csv(url_tg)
# Drop the Rdatasets index column
for c in tg.columns:
    if 'Unnamed' in c:
        tg = tg.drop(columns=[c])
tg.head()

# Groups
oj = tg.loc[(tg['supp']=='OJ') & (tg.dose == 0.5),'len']
vc = tg.loc[(tg['supp']=='VC') & (tg.dose == 0.5),'len']

# Welch's t-test
t_stat, p_val = stats.ttest_ind(oj, vc, equal_var=False)

# Cohen's d (using pooled SD with group sizes)
def cohens_d(x, y):
    nx, ny = len(x), len(y)
    sx2, sy2 = np.var(x, ddof=1), np.var(y, ddof=1)
    sp2 = ((nx-1)*sx2 + (ny-1)*sy2) / (nx+ny-2)
    d = (np.mean(x) - np.mean(y)) / np.sqrt(sp2)
    return d

d = cohens_d(oj, vc)

print(f"Welch t‑test: t = {t_stat:.4f}, p = {p_val:.6f}")
print(f"Mean's: OJ = {np.mean(oj):.3f}, VC = {np.mean(vc):.3f}")
print(f"Cohen's d = {d:.3f}")

Welch t‑test: t = 3.1697, p = 0.006359
Mean's: OJ = 13.230, VC = 7.980
Cohen's d = 1.418

def plot_bivariate_histograms(dataset, con_col, cat_col, title='', x_label = ''):
    bi_con_cat = dataset.groupby([cat_col])[con_col].plot.hist(alpha = 0.5)
    plt.xlabel(con_col)
    plt.legend(dataset.groupby([cat_col])[con_col].count().axes[0].tolist())
    plt.title(title)
    plt.xlabel(x_label)
    plt.show();
plot_bivariate_histograms(tg[tg.dose==0.5], 'len', 'supp', title='Joint histogram', x_label = 'Length')

# Load Supermarket Sales
df = pd.read_csv("https://raw.githubusercontent.com/selva86/datasets/master/supermarket_sales.csv")

# Tidy column names
df.columns = [c.strip().replace(' ', '_').lower() for c in df.columns]

df

# One-way ANOVA
model = ols('gross_income ~ C(branch)', data=df).fit()
anova_tbl = sm.stats.anova_lm(model, typ=2)
print(anova_tbl)

                  sum_sq         df        F   PR(>F)
C(branch)     242.602644   2.000000 0.884583 0.413210
Residual  136,716.894906 997.000000      NaN      NaN

# Levene test for homogeneity of variances
from scipy.stats import levene
branch_a = df[df['branch']=='A']['gross_income']
branch_b = df[df['branch']=='B']['gross_income']
branch_c = df[df['branch']=='C']['gross_income']
print("Levene p-value:", levene(branch_a, branch_b, branch_c).pvalue)

Levene p-value: 0.08946425577002974

# Tukey HSD post-hoc
tukey = pairwise_tukeyhsd(endog=df['gross_income'], groups=df['branch'], alpha=0.05)
print(tukey.summary())

Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     A      B    0.358 0.9171 -1.7627 2.4788  False
     A      C   1.1784 0.3954 -0.9489 3.3057  False
     B      C   0.8203 0.6405 -1.3195 2.9602  False
---------------------------------------------------

bi_variate_boxplot = sns.boxplot(x="branch", y="gross_income", data=df)
bi_variate_boxplot.set(title = 'Box Chart of gross income across branches');

import pandas as pd
from scipy.stats import chi2_contingency
import numpy as np

# Load dataset
url = "https://raw.githubusercontent.com/gokcengiz/Shopping-data-analysis/main/customer_shopping_data.csv"
df = pd.read_csv(url)
df

# Build contingency table: Gender vs Payment Method
ct = pd.crosstab(df['category'], df['gender'])
chi2, p, dof, expected = chi2_contingency(ct)
print("Contingency Table:\n", ct)
print(f"Chi-square = {chi2:.4f}, p-value = {p:.6f}, dof = {dof}")

Contingency Table:
 gender           Female   Male
category                      
Books              2906   2075
Clothing          20652  13835
Cosmetics          9070   6027
Food & Beverage    8804   5972
Shoes              5967   4067
Souvenir           3017   1982
Technology         2981   2015
Toys               6085   4002
Chi-square = 7.5679, p-value = 0.372234, dof = 7

# Compute Cramér's V
n = ct.values.sum()
phi2 = chi2 / n
r, c = ct.shape
cramers_v = np.sqrt(phi2 / min(r-1, c-1))

print(f"Cramér's V = {cramers_v:.3f}")

Cramér's V = 0.009

from statsmodels.graphics.mosaicplot import mosaic
def plot_mosaics(data, x_col, y_col, title='', colors_list =[]):
    dict_of_tuples = {}

    # create the clean set of percentages to print
    for x_col_ in data[x_col].unique():
        for y_col_ in data[y_col].unique():
            n = len(data[(data[x_col]==x_col_)&(data[y_col]==y_col_)][x_col])
            d = len(data[(data[x_col]==x_col_)][x_col])
            len_ = len(data[x_col])
            if((d==0) or (n/d<=0.04)):
                # if the percentage within a class is less than 4%, do not print the percentage
                dict_of_tuples[(str(x_col_), str(y_col_))] = ''
            elif(n/len_<=0.02):
                # If its a tiny class with less than 2% of the total data, do not print
                dict_of_tuples[(str(x_col_), str(y_col_))] = ''
            else:
                dict_of_tuples[(str(x_col_), str(y_col_))] = str(int(n/d*100))+"%"

    dict_of_colors = dict_of_tuples.copy()
    if(len(colors_list)>0):
        # create a clean set of colors
        for i, x_col_ in enumerate(data[x_col].unique()):
            for y_col_ in data[y_col].unique():
                dict_of_colors[(str(x_col_), str(y_col_))] = {'color':colors_list[i], 'alpha':0.8}

    # Plot the mosaic plot
    labelizer = lambda k: dict_of_tuples[k]    
    fig, ax = plt.subplots(figsize=(8,6))
    if(len(colors_list)>0):
        mosaic(data.sort_values([x_col, y_col]), [x_col, y_col], 
               statistic = False, axes_label = True, label_rotation = [90, 0],
               labelizer=labelizer, properties=dict_of_colors, gap=0.008, ax=ax)
    else:
        mosaic(data.sort_values([x_col, y_col]), [x_col, y_col], 
               statistic = False, axes_label = True, label_rotation = [90, 0],
               labelizer=labelizer, gap=0.008, ax=ax)
    if(title==''):
        plt.title(str(y_col) + ' percentages across ' + str(x_col))
    else:
        plt.title(title)
    plt.show();
plot_mosaics(df, 'category', 'gender')