Bootstrap-t-student/bootstrap-t.ipynb

Projekt - Test t studenta

• Marcin Kostrzewski
• Krystian Wasilewski
• Mateusz Tylka

Test t studenta

Metoda statystyczna służącą do porównania dwóch średnich między sobą gdy znamy liczbę badanych próbek, średnią arytmetyczną oraz wartość odchylenia standardowego lub wariancji. Jest to jeden z mniej skomplikowanych i bardzo często wykorzystywanych testów statystycznych używanych do weryfikacji hipotez. Dzięki niemu możemy dowiedzieć się czy dwie różne średnie są różne niechcący (w wyniku przypadku) czy są różne istotnie statystycznie (np. z uwagi na naszą manipulację eksperymentalna). Wyróżniamy 3 wersję testu t:

1. test t Studenta dla jednej próby
2.  test t Studenta dla prób niezależnych
3.  test t Studenta dla prób zależnych

Wszystkie rodzaje testów są testami parametrycznymi, a co za tym idzie nasze mierzone zmienne ilościowe powinny mieć rozkład normalny.

import numpy as np
import pandas as pd
from math import sqrt
from scipy import stats
from scipy.stats import sem
from scipy.stats import t
import matplotlib.pyplot as plt
from statistics import mean, stdev
from scipy.stats import ttest_ind, ttest_1samp, ttest_rel

dataset = pd.read_csv('experiment_data.csv')

def calculate_p(t_stat, df):
    """Funkcja oblicza wartość *p* na podstawie statystyki testowej i stopni swobody"""
    return (1.0 - t.cdf(abs(t_stat), df)) * 2.0

def calculate_cv(df, alpha=0.05):
    """Funkcja oblicza wartość krytyczną (critical value)"""
    return t.ppf(1.0 - alpha, df)

def t_test(sample_1, sample_2=None, df_fn=df_single, t_stat_fn=t_stat_single, population_mean=None, alpha=0.05):
    """
    Funkcja przeprowadza test T-studenta dla dwóch zmiennych.
    liczba kolumn wynosi 1, test jest przeprowadzany dla jednej zmiennej.
    @param df_fn - funkcja obliczająca stopnie swobody
    @param t_stat_fn - funkcja obliczająca statystykę T
    """
    t_stat_list = get_t_stats(sample_1, sample_2, t_stat_fn, population_mean=population_mean)
    t_stat_sum = sum(t_stat_list)

    data_size = sample_1.shape[0]

    t_stat = t_stat_sum / data_size
    # TODO: dolna i górna opcja dają inne wyniki z jakiegoś powodu (???)
    t_stat = mean(t_stat_list)

    if sample_2 is None:
        df = df_fn(sample_1)
    else:
        df = df_fn(sample_1, sample_2)
    cv = calculate_cv(df, alpha)
    p = calculate_p(t_stat, df)
    return t_stat, df, cv, p, t_stat_list

def get_t_stats(sample_1, sample_2=None, t_stat_fn=t_stat_single, population_mean=None):
    """Funkcja oblicza listę statystyk testowych dla każdej próbki bootstrapowej wybranej na podstawie danych sample_1 i sample_2"""
    t_stat_list = []

    # One sample test
    if t_stat_fn==t_stat_single:
        if not population_mean:
            raise Exception("population_mean not provided")
        for bootstrap in generate_bootstraps(sample_1):
            stat = t_stat_fn(bootstrap, population_mean)
            t_stat_list.append(stat)
        return t_stat_list

    # Two sample test
    for bootstrap_1, bootstrap_2 in zip(generate_bootstraps(sample_1), generate_bootstraps(sample_2)):
        stat = t_stat_fn(bootstrap_1, bootstrap_2)
        t_stat_list.append(stat)
    return t_stat_list

Test Shapiro Wilka

Wszystkie rodzaje testów są testami parametrycznymi, a co za tym idzie nasze mierzone zmienne ilościowe powinny mieć rozkład normalny.

ALPHA = 0.05
female_heights = dataset['Female height'].to_numpy()
shapiro_test = stats.shapiro(female_heights)

if shapiro_test.pvalue >  ALPHA:
    print("Female height: Dane mają rozkład normalny.")
else:
    print("Female height: Dane nie mają rozkładu normalnego.")

male_heights = dataset['Male height'].to_numpy()
shapiro_test = stats.shapiro(male_heights)

if shapiro_test.pvalue >  ALPHA:
    print("Male height: Dane mają rozkład normalny.")
else:
    print("Male height: Dane nie mają rozkładu normalnego.")

weights_before = dataset['Weight before'].to_numpy()
shapiro_test = stats.shapiro(weights_before)

if shapiro_test.pvalue >  ALPHA:
    print("Weight before: Dane mają rozkład normalny.")
else:
    print("Weight before: Dane nie mają rozkładu normalnego.")

weights_after = dataset['Weight after'].to_numpy()
shapiro_test = stats.shapiro(weights_after)

if shapiro_test.pvalue >  ALPHA:
    print("Weight after: Dane mają rozkład normalny.")
else:
    print("Weight after: Dane nie mają rozkładu normalnego.")

Female height: Dane mają rozkład normalny.
Male height: Dane mają rozkład normalny.
Weight before: Dane mają rozkład normalny.
Weight after: Dane mają rozkład normalny.

Testowanie hipotez metodą bootstrap

Bootstrap – metoda szacowania (estymacji) wyników poprzez wielokrotne losowanie ze zwracaniem z próby. Polega ona na utworzeniu nowego rozkładu wyników, na podstawie posiadanych danych, poprzez wielokrotne losowanie wartości z posiadanej próby. Metoda ze zwracaniem polega na tym, że po wylosowaniu danej wartości, “wraca” ona z powrotem do zbioru.

Metoda bootstrapowa znajduje zastosowanie w sytuacji, w której nie znamy rozkładu z populacji z której pochodzi próbka lub w przypadku rozkładów małych lub asymetrycznych. W takim wypadku, dzięki tej metodzie, wyniki testów parametrycznych i analiz opartych o modele liniowe są bardziej precyzyjne. Zazwyczaj losuje się wiele próbek, np. 2000 czy 5000.

def generate_bootstraps(data, n_bootstraps=100):
    data_size = data.shape[0]
    for _ in range(n_bootstraps):
        indices =  np.random.choice(len(data), size=data_size)
        yield data.iloc[indices, :]

Test t studenta dla jednej próby

Test t Studenta dla jednej próby wykorzystujemy gdy chcemy porównać średnią “teoretyczną” ze średnią, którą faktycznie możemy zaobserwować w naszej bazie danych. Średnia teoretyczna to średnia pochodząca z innych badań lub po prostu bez większych uzasadnień pochodząca z naszej głowy.

Wyobraźmy sobie, że mamy dane z takimi zmiennymi jak wzrost pewnej grupy ludzi. Dzięki testowi t Studenta dla jednej próby możemy dowiedzieć się np. czy wzrost naszego młodszego brata wynoszący 155cm odbiega znacząco od średniej wzrostu tej grupy. Hipoteza zerowa w takim badaniu wyglądałaby następująco H0: Badana próba została wylosowana z populacji, w której wzrost osób wynosi średnio 155cm. Z kolei hipoteza alternatywna będzie brzmiała H1: Badana próba nie została wylosowana z populacji gdzie średni wzrost wynosi 155cm

def t_stat_single(sample, population_mean):
    """Funkcja oblicza wartość statystyki testowej dla jednej próbki"""
    if sample.empty:
        raise Exception("Empty sample")
    sample = sample['Height'].values.tolist()
    sample_size = len(sample)
    return (mean(sample) - population_mean) / (stdev(sample) / sqrt(sample_size))

def df_single(sample_1):
    """Funkcja oblicza stopnie swobody dla jednej próbki"""
    # TODO: I have no clue what to return from here
    return len(sample_1)

def bootstrap_one_sample(sample, population_mean):
    return t_test(
        sample_1=sample,
        df_fn=df_single,
        t_stat_fn=t_stat_single,
        population_mean=population_mean
    )

Sprawdzenie czy osoba o wzroście 165cm pasuje do populacji (nie jest odmieńcem)

#TODO: poprawić kod aby można było podawać kolumny

t_stat, df, cv, p, _ = bootstrap_one_sample(dataset, 165)
pretty_print_full_stats(t_stat, df, cv, p)

t: 6.854929920812628, df: 500, cv: 1.6479068539295045, p: 2.1091128843409024e-11

TODO: Wniosek

Test t studenta dla prób niezależnych

Test t Studenta dla prób niezależnych jest najczęściej stosowaną metodą statystyczną w celu porównania średnich z dwóch niezależnych od siebie grup. Wykorzystujemy go gdy chcemy porównać dwie grupy pod względem jakiejś zmiennej ilościowej. Na przykład gdy chcemy porównać średni wzrost kobiet i mężczyzn w danej grupie.

Zazwyczaj dwie średnie z różnych od siebie grup będą się różnić. Test t Studenta powie nam jednak czy owe różnice są istotne statystycznie – czy nie są przypadkowe. Hipoteza zerowa takiego testu będzie brzmiała H0: Średni wzrost w grupie mężczyzn jest taki sam jak średni w grupie kobiet. Hipoteza alternatywna z kolei H1: Kobiety będą różnić się od mężczyzn pod wzrostu. Jeśli wynik testu t Studenta będzie istotny na poziomie p < 0,05 możemy odrzucić hipotezę zerową na rzecz hipotezy alternatywnej.

def t_stat_ind(sample_1, sample_2):
    """Funkcja oblicza wartość statystyki testowej dla dwóch próbek niezależnych"""
    if sample_1.empty or sample_2.empty:
        raise Exception("Empty sample")
    sample_1 = sample_1[0].values.tolist()
    sample_2 = sample_2[0].values.tolist()
    sed = sqrt(sem(sample_1)**2 + sem(sample_2)**2)
    return (mean(sample_1) - mean(sample_2)) / sed

def df_ind(sample_1, sample_2):
    """Funkcja oblicza stopnie swobody dla dwóch próbek niezależnych"""
    return len(sample_1) + len(sample_2) - 2

def bootstrap_independent(sample_1, sample_2):
    return t_test(
        sample_1=sample_1,
        sample_2=sample_2,
        df_fn=df_ind,
        t_stat_fn=t_stat_ind
    )

#TODO: Wyciągnąć wysokości kobiet i mężczyzn oraz poprawić kod aby można było podawać kolumny
t_stat, df, cv, p, _ = bootstrap_independent(dataset, dataset)
pretty_print_full_stats(t_stat, df, cv, p)

TODO: Wniosek

Test t studenta dla prób zależnych

W odróżnieniu od testu t – Studenta dla prób niezależnych, gdzie porównujemy dwie grupy, ten rodzaj testu stosujemy gdy poddajemy analizie tą samą pojedynczą grupę, ale dwukrotnie w czasie. Na przykład gdy chcemy porównać średnie wagi grupy osób przed dietą oraz po diecie, aby sprawdzić czy dieta spowodowała istotne zmiany statystyczne.

Hipoteza zerowa takiego testu będzie brzmiała H0: Średnia waga osób po diecie jest taka sama jak przed dietą. Hipoteza alternatywna z kolei H1: Dieta znacząco wpłynęła na średnią wagę danej grupy.

def t_stat_dep(sample_1, sample_2, mu=0):
    """Funkcja oblicza wartość statystyki testowej dla dwóch próbek zależnych"""
    if sample_1.empty or sample_2.empty:
        raise Exception("Empty sample")
    sample_1 = sample_1[0].values.tolist()
    sample_2 = sample_2[0].values.tolist()
    differences = [x_1 - x_2 for x_1, x_2 in zip(sample_1, sample_2)]
    sample_size = len(sample_1)
    return (mean(differences) - mu) / (stdev(differences) / sqrt(sample_size))

def df_dep(sample_1, sample_2):
    """Funkcja oblicza stopnie swobody dla dwóch próbek zależnych"""
    l1, l2 = len(sample_1), len(sample_2)
    if l1 != l2:
        raise Exception("Samples aren't of equal length")
    return l1

def bootstrap_dependent(sample_1, sample_2):
    return t_test(
        sample_1=sample_1,
        sample_2=sample_2,
        df_fn=df_dep,
        t_stat_fn=t_stat_dep
    )

TODO: Wyciągnąć wagi przed dietą i po oraz poprawić kod aby można było podawać kolumny

t_stat, df, cv, p, _ = bootstrap_dependent(dataset, dataset) pretty_print_full_stats(t_stat, df, cv, p)

TODO: Wniosek

Wykresy

def draw_distribution(stats):
    """
    Funkcja rysuje rozkład statystyki testowej
    @param stats: lista statystyk testowych
    """
    plt.hist(stats)
    plt.xlabel('Test statistic value')
    plt.ylabel('Frequency')
    plt.show()

Testy

# Testy dla samych statystyk testowych
def pretty_print_stats(t_stat_selfmade, t_stat_lib, suffix):
    print(f'Statystyka testowa dla {suffix}:')
    print(t_stat_selfmade, '- z naszej funkcji')
    print(t_stat_lib, '- z gotowej biblioteki')
    print()
    
dummy = pd.DataFrame([1, 2, 3, 4, 5])
dummy2 = pd.DataFrame([4, 5, 6, 7, 8])
dummy3 = pd.DataFrame([1, 3 , 3, 4, 6])

t_stat_selfmade = t_stat_single(dummy, 2)
t_stat_lib, _ = ttest_1samp(dummy, 2)
pretty_print_stats(t_stat_selfmade, t_stat_lib, 'jednej próby')

t_stat_selfmade = t_stat_ind(dummy, dummy2)
t_stat_lib, _ = ttest_ind(dummy, dummy2)
pretty_print_stats(t_stat_selfmade, t_stat_lib, 'dwóch prób niezależnych')

t_stat_selfmade = t_stat_dep(dummy, dummy3)
t_stat_lib, _ = ttest_rel(dummy, dummy3)
pretty_print_stats(t_stat_selfmade, t_stat_lib, 'dwóch prób zależnych')

Statystyka testowa dla jednej próby:
1.414213562373095 - z naszej funkcji
[1.41421356] - z gotowej biblioteki

Statystyka testowa dla dwóch prób niezależnych:
-3.0 - z naszej funkcji
[-3.] - z gotowej biblioteki

Statystyka testowa dla dwóch prób zależnych:
-1.6329931618554525 - z naszej funkcji
[-1.63299316] - z gotowej biblioteki

# Testy z bootstrappowaniem

def pretty_print_full_stats(t_stat, df, cv, p):
    print(f't: {t_stat}, df: {df}, cv: {cv}, p: {p}\n')

print(type(dummy))

print('Statystyki dla jednej próby:')
t_stat, df, cv, p, _ = bootstrap_one_sample(dummy, 2)
pretty_print_full_stats(t_stat, df, cv, p)

print('Statystyki dla dwóch prób zależnych:')
t_stat, df, cv, p, _ = bootstrap_dependent(dummy2, dummy3)
pretty_print_full_stats(t_stat, df, cv, p)

print('Statystyki dla dwóch prób niezależnych:')
t_stat, df, cv, p, _ = bootstrap_independent(dummy2, dummy3)
pretty_print_full_stats(t_stat, df, cv, p)

<class 'pandas.core.frame.DataFrame'>
Statystyki dla jednej próby:
t: 1.8073147056683616, df: 5, cv: 2.015048372669157, p: 0.13052275003443325

Statystyki dla dwóch prób zależnych:
t: 3.0790273716290404, df: 5, cv: 2.015048372669157, p: 0.027500015466573435

Statystyki dla dwóch prób niezależnych:
t: 2.8109511013364576, df: 8, cv: 1.8595480375228421, p: 0.02280961069987497