Section 3.7 — Inventory of statistical tests¶

This notebook contains the code examples from Section 3.7 Inventory of statistical tests from the No Bullshit Guide to Statistics.

Notebook setup¶

In [1]:

# load Python modules
import os
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# load Python modules import os import numpy as np import pandas as pd import seaborn as sns import matplotlib.pyplot as plt

In [2]:

# Plot helper functions
from ministats import plot_pdf

# Plot helper functions from ministats import plot_pdf

In [3]:

# Figures setup
plt.clf()  # needed otherwise `sns.set_theme` doesn't work
from plot_helpers import RCPARAMS
RCPARAMS.update({'figure.figsize': (10, 3)})   # good for screen
# RCPARAMS.update({'figure.figsize': (5, 1.6)})  # good for print
sns.set_theme(
    context="paper",
    style="whitegrid",
    palette="colorblind",
    rc=RCPARAMS,
)

# Useful colors
snspal = sns.color_palette()
blue, orange, purple = snspal[0], snspal[1], snspal[4]

# High-resolution please
%config InlineBackend.figure_format = 'retina'

# Where to store figures
DESTDIR = "figures/stats/inventory"

# Figures setup plt.clf() # needed otherwise `sns.set_theme` doesn't work from plot_helpers import RCPARAMS RCPARAMS.update({'figure.figsize': (10, 3)}) # good for screen # RCPARAMS.update({'figure.figsize': (5, 1.6)}) # good for print sns.set_theme( context="paper", style="whitegrid", palette="colorblind", rc=RCPARAMS, ) # Useful colors snspal = sns.color_palette() blue, orange, purple = snspal[0], snspal[1], snspal[4] # High-resolution please %config InlineBackend.figure_format = 'retina' # Where to store figures DESTDIR = "figures/stats/inventory"

<Figure size 640x480 with 0 Axes>

In [4]:

# set random seed for repeatability
np.random.seed(42)

# set random seed for repeatability np.random.seed(42)

In [5]:

#######################################################

#######################################################

Definitions¶

In [ ]:

Assumptions¶

In [ ]:

NHST procedure¶

In [ ]:

Categorization of statistical test recipes¶

In [ ]:

Z-Tests¶

One-sample $z$-test¶

See the examples/one_sample_z-test.ipynb notebook.

In [ ]:

Proportion tests¶

One-sample $z$-test for proportions¶

Binomial test¶

Two-sample $z$-test for proportions¶

In [6]:

from statsmodels.stats.proportion import proportions_ztest

from statsmodels.stats.proportion import proportions_ztest

In [ ]:

T-tests¶

One sample $t$-test¶

See the examples/one_sample_t-test.ipynb notebook.

Welch's two-sample $t$-test¶

(explain pooled variance as special case "Two sample t-test", but inferior)

Paired $t$-test¶

In [ ]:

Chi-square tests¶

Chi-square test for goodness of fit¶

Example: are digits of $\pi$ random?¶

In [7]:

pidigits = [99959,  99757, 100026, 100230, 100230, 100359,  99548,  99800, 99985, 100106]
# obtained using   np.bincount(list(str(sympy.N(sympy.pi, 1_000_000)).replace('.','')))

os = pidigits           # observed
es = [1_000_000/10]*10  # expected (uniform)

from scipy.stats import chisquare
chisquare(f_obs=os, f_exp=es)

pidigits = [99959, 99757, 100026, 100230, 100230, 100359, 99548, 99800, 99985, 100106] # obtained using np.bincount(list(str(sympy.N(sympy.pi, 1_000_000)).replace('.',''))) os = pidigits # observed es = [1_000_000/10]*10 # expected (uniform) from scipy.stats import chisquare chisquare(f_obs=os, f_exp=es)

Out[7]:

Power_divergenceResult(statistic=5.51852, pvalue=0.7869706202650393)

See original blog post for useful historical context about this https://probabilityandstats.wordpress.com/2017/03/14/are-digits-of-pi-random/

In [ ]:

Chi-square test of independence¶

Chi-square test for homogeneity¶

Chi-square test for the population variance¶

In [ ]:

Analysis of variance (ANOVA) tests¶

One-way analysis of variance (ANOVA)¶

Two-way ANOVA¶

In [ ]:

Nonparametric tests¶

Use when assumptions for other tests not valid

Sign test for the population median¶

via https://vitalflux.com/sign-test-hypothesis-python-examples/

In [8]:

from scipy.stats import binomtest

n_pos = 6
n_neg = 9
n_min = min(n_pos, n_neg)
n_tot = n_pos + n_neg

# Calculate p-value (two-tailed) using the binomial test
binomtest(k=n_min, n=n_tot, p=0.5, alternative='two-sided')

from scipy.stats import binomtest n_pos = 6 n_neg = 9 n_min = min(n_pos, n_neg) n_tot = n_pos + n_neg # Calculate p-value (two-tailed) using the binomial test binomtest(k=n_min, n=n_tot, p=0.5, alternative='two-sided')

Out[8]:

BinomTestResult(k=6, n=15, alternative='two-sided', statistic=0.4, pvalue=0.6072387695312499)

In [9]:

n_max = max(n_pos, n_neg)
binomtest(k=n_max, n=n_tot, p=0.5, alternative='two-sided')

n_max = max(n_pos, n_neg) binomtest(k=n_max, n=n_tot, p=0.5, alternative='two-sided')

Out[9]:

BinomTestResult(k=9, n=15, alternative='two-sided', statistic=0.6, pvalue=0.6072387695312499)

In [ ]:

In [ ]:

One-sample Wilcoxon signed-rank test¶

In [ ]:

Mann-Whitney U-test¶

example via https://www.reneshbedre.com/blog/mann-whitney-u-test.html

In [10]:

dfw = pd.read_csv("https://reneshbedre.github.io/assets/posts/mann_whitney/genotype.csv")
dfw.shape
# dfw

dfw = pd.read_csv("https://reneshbedre.github.io/assets/posts/mann_whitney/genotype.csv") dfw.shape # dfw

Out[10]:

(23, 2)

In [11]:

from scipy.stats import mannwhitneyu

mannwhitneyu(x=dfw["A"], y=dfw["B"], alternative="two-sided")

from scipy.stats import mannwhitneyu mannwhitneyu(x=dfw["A"], y=dfw["B"], alternative="two-sided")

Out[11]:

MannwhitneyuResult(statistic=489.5, pvalue=7.004695394561307e-07)

In [ ]:

Kruskal-Wallis analysis of variance by ranks¶

In [ ]:

Resampling methods¶

Simulation tests¶

In [12]:

from ministats.hypothesis_tests import simulation_test

%psource simulation_test

from ministats.hypothesis_tests import simulation_test %psource simulation_test

def simulation_test(sample, rvH0, estfunc, alt="two-sided"):
    """
    Compute the p-value of the observed estimate `estfunc(sample)` under H0
    described by the random variable `rvH0`.
    """
    # 1. Compute the observed value of `estfunc`
    obsest = estfunc(sample)
    n = len(sample)

    # 2. Get sampling distribution of `estfunc` under H0
    sampl_dist_H0 = gen_sampling_dist(rvH0, estfunc, n)

    # 3. Compute the p-value
    tails = tailvalues(sampl_dist_H0, obsest, alt=alt)
    pvalue = len(tails) / len(sampl_dist_H0)
    return pvalue

In [ ]:

Two-sample permutation test¶

In [13]:

from ministats.hypothesis_tests import permutation_test

%psource permutation_test

from ministats.hypothesis_tests import permutation_test %psource permutation_test

def permutation_test(xsample, ysample, estfunc, P=10000):
    """
    Compute the p-value of the observed estimate `estfunc(xsample,ysample)`
    under the null hypothesis where the group membership is randomized.
    """
    # 1. Compute the observed value of `estfunc`
    obsest = estfunc(xsample, ysample)

    # 2. Get sampling dist. of `estfunc` under H0
    pestimates = []
    for i in range(0, P):
        rsx, rsy = resample_under_H0(xsample, ysample)
        pestimate = estfunc(rsx, rsy)
        pestimates.append(pestimate)

    # 3. Compute the p-value
    tails = tailvalues(pestimates, obsest)
    pvalue = len(tails) / len(pestimates)
    return pvalue

In [ ]:

Permutation ANOVA¶

In [14]:

from ministats import permutation_anova

%psource permutation_anova

from ministats import permutation_anova %psource permutation_anova

def permutation_anova(samples, P=10000, alt="greater"):
    """
    Compute the p-value of the observed F-statistic for `samples` list
    under the null hypothesis where the group membership is randomized.
    """
    ns = [len(sample) for sample in samples]

    # 1. Comptue the observed F-statistic
    obsfstat, _ = f_oneway(*samples)

    # 2. Get sampling dist. of F-statistic under H0
    pfstats = []
    for i in range(0, P):
        values = np.concatenate(samples)
        pvalues = np.random.permutation(values)
        psamples = []
        nstart = 0
        for nstep in ns:
            psample = pvalues[nstart:nstart+nstep]
            psamples.append(psample)
            nstart = nstart + nstep
        pfstat, _ = f_oneway(*psamples)
        pfstats.append(pfstat)

    # 3. Compute the p-value
    tails = tailvalues(pfstats, obsfstat, alt=alt)
    pvalue = len(tails) / len(pfstats)
    return pvalue

In [15]:

# test on three samples
from scipy.stats import norm

# Random samples
np.random.seed(43)
sample1 = norm(loc=0).rvs(size=30)
sample2 = norm(loc=0).rvs(size=30)
sample3 = norm(loc=0.7).rvs(size=30)

np.random.seed(45)
permutation_anova([sample1, sample2, sample3])

# test on three samples from scipy.stats import norm # Random samples np.random.seed(43) sample1 = norm(loc=0).rvs(size=30) sample2 = norm(loc=0).rvs(size=30) sample3 = norm(loc=0.7).rvs(size=30) np.random.seed(45) permutation_anova([sample1, sample2, sample3])

Out[15]:

0.0297

In [16]:

# compare with analytical formula
from scipy.stats import f_oneway

f_oneway(sample1, sample2, sample3)

# compare with analytical formula from scipy.stats import f_oneway f_oneway(sample1, sample2, sample3)

Out[16]:

F_onewayResult(statistic=3.6808227678358856, pvalue=0.029206733498721497)

In [ ]:

Equivalence tests¶

See examples/two_sample_equivalence_test.ipynb for an example.

In [ ]:

Distribution checks¶

Kolmogorov-Smirnov test¶

Shapiro-Wilk normality test¶

Discussion¶

In [ ]:

In [ ]:

Discussion¶

In [ ]:

Exercises¶

In [ ]:

Links¶

In [ ]: