Rozdział 8 Porównanie zmiennych niezależnych

---

8.1 Porównanie średnich

Test t-Studenta / Test t-Welcha. Do porównania dwóch średnich tj. do zweryfikowania hipotezy \(H_0:\mu_1=\mu_2\) najczęsciej proponowany jest test t-Studenta Wymaga on spełnienia dwóch warunków: normalność rozkładu oraz jednorodności wariancji. Statystyka klasycznego testu dla dwóch średnich: \((\bar{x}_1-\bar{x}_2)/\sqrt{d_1+d_2}\) ma rozkład t-Studenta ze stopniami swobody \(df=n_1+n_2-2\). Jeśli wariancje w próbkach nie są równe to zalecane jest stosowanie poprawki Welcha (Derrick i White 2016) która polega na modyfikacji stopni swobody: \[\begin{equation} df_{\mathrm{Welch}}=\frac{(d_1+d_2)^2}{\frac{d_1}{n_1-1}+\frac{d_2}{n_2-1}} \tag{8.1} \end{equation}\] gdzie: \(s_k^2\) to wariancja, \(n_k\) to liczebność próby dla \(k=1,2\) oraz \(d_k=s^2_k/n_k\).

Dzięki funkcji pingouin.ttest jest dostępna klasyczna wersja testu t-Studenta oraz z poprawką Welcha.

import scipy.stats as stats
import pingouin as pg

x = stats.norm.rvs(0, 1, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)

print(pg.ttest(x,y, correction=True)[['T','dof','CI95%','p-val']])

##             T    dof           CI95%     p-val
## T-test -3.448  23.26  [-2.95, -0.74]  0.002164

Anova / Welch-Anova. Klasyczna analiza wariancji - inaczej ANOVA to rozwinięcie testu t-Studenta dla więcej niż dwóch zmiennych niezależnych. Inaczej mówiąc w przypadku porównania średnich z dwóch grup wyniki z obu procedur są tożsame. Funkcja pingouin.anova realizuje jedno lub dwuczynnikową analizę wariancji z interakcją. Natomiast do funkcji pingouin.welch_anova jest zaimplementowana metoda Welcha dla jednej zmiennej grupującej.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

y = np.concatenate((stats.norm.rvs(0, 1, size=20, random_state=2305),
                    stats.norm.rvs(1.5, 2, size=20, random_state=4101),
                    stats.norm.rvs(1.5, 2, size=20, random_state=4026)))
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)

print(pg.welch_anova(dv='y', between='g', data=pd.DataFrame({'y':y,'g':g})))

##   Source  ddof1   ddof2      F     p-unc
## 0      g      2  30.937  6.448  0.004562

Dalsza analiza. Po odrzuceniu hipotezy zerowej w analizie wariancji stosujemy testy do porównań wielokrotnych. Jednym z bardziej popularnych tzw. testów po fakcie dla grup niezależnych jest procedura Tukeya lub seria testów t-Studenta z odpowiednią korektą p-wartości. Są one zaimplementowane odpowiednio do funkcji pingouin.pairwise_tukey oraz pingouin.pairwise_ttests. Natomiast w warunkach heteroskedastyczności można wykonać serię testów t-Welcha z odpowiednią korektą p-wartości lub test Gamesa-Howella. Są one dostępne odpowiednio w funkcji pingouin.pairwise_ttests oraz pingouin.pairwise_gameshowell. Wiele ciekawych rozwiązań np. test Tamhane T2 zostało zaimplementowanych do pakietu scikit-posthocs który bazuje na bibliotece PMCMRplus/PMCMR dla programu R.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

y = np.concatenate((stats.norm.rvs(0, 1, size=20, random_state=2305),
                    stats.norm.rvs(1.5, 2, size=20, random_state=4101),
                    stats.norm.rvs(1.5, 2, size=20, random_state=4026)))
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)

print(pg.pairwise_gameshowell(dv='y', between='g',
                           data=pd.DataFrame({'y':y,'g':g})))

##      A    B  mean(A)  mean(B)   diff  ...      T      df      pval  efsize  eftype
## 0  1.0  2.0    0.155    2.001 -1.846  ... -3.448  23.261  0.002070  -1.069  hedges
## 1  1.0  3.0    0.155    0.857 -0.702  ... -1.539  25.039  0.275781  -0.477  hedges
## 2  2.0  3.0    2.001    0.857  1.143  ...  1.730  36.818  0.196416   0.536  hedges
## 
## [3 rows x 12 columns]

8.2 Porównanie rang

Test Manna-Whitneya. Przy założeniu, że dwa badane rozkłady mają ten sam kształt (takie same wariancje, skośność itp.) można zweryfikować hipotezę zerową o postaci \(H_{0}:\;F(x)=G(y+\Delta)\) w której parametr \(\Delta\) określa przesunięcie dystrybuanty \(G(y)\) względem dystrybuanty \(F(x)\) (Divine i in. 2018). Inaczej mówiąc rozmieszczenie rozkładów \(F(x)\) i \(G(y)\) różni się w zależności od \(\Delta\). Parametr przesunięcia można oszacować za pomocą estymatorora Hodgesa-Lehmanna: \[\begin{equation} \hat{\Delta}=\mbox{mediana}\{x_i-y_j\;:\;i=1,\;\dots n_1\;;\;j=1,\;\dots n_2\} \tag{8.2} \end{equation}\]

Warto zaznaczyć, że parametr \(\Delta\) w funkcji scipy.stats.mannwhitneyu oraz pingouin.mwu ma stałą wartość równą zero. Zatem rozważana hipoteza zerowa ma postać: \[\begin{equation} H_{0}:\;\Delta=0\quad\textrm{vs.}\quad H_{1}:\;\Delta\neq 0 \tag{8.3} \end{equation}\] Równoważnym zapisem może być: \[\begin{equation} H_{0}:\;F(x)=G(y)\quad\textrm{vs.}\quad H_{1}:\;F(x)\neq G(y) \tag{8.4} \end{equation}\]

Statystyka testowa: \[\begin{equation} Z=\frac{|W-\frac{n_1n_2}{2}|-0,5}{\sqrt{\frac{n_1n_2(n_1+n_2+1)}{12}-\frac{n_1n_2\sum_{i=1}^{c}(t^3-t)}{12(n_1+n_2)(n_1+n_2-1)}}} \tag{8.5} \end{equation}\] gdzie: \(c\) to liczba grup pomiarów wiązanych, \(t_i\) to liczba pomiarów wiązanych w \(i\)-tej grupie pomiarów wiązanych.

import scipy.stats as stats
import pingouin as pg
import numpy as np

x = stats.norm.rvs(0, 1, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)

hl = np.median(x[:, None] - y)
df = pg.mwu(x,y)
df['LH-median'] = hl
print(df)

##      U-val     p-val   RBC  CLES  LH-median
## MWU   96.0  0.005115  0.52  0.76  -2.003235

Test Brunera-Munzela. Dobrą alternatywną dla testu sumy rang Wilcoxona w warunakch heteroskedastyczności może być test Brunera-Munzela (Brunner i Munzel 2000) dostępny dzięki funkcji scipy.stats.brunnermunzel. Wersja permutacyjna tego testu (Neubert i Brunner 2007) jest zalecana dla przypadku małolicznych próbek o nierównych liczebnościach. W tym teście hipoteza zerowa dla równości stochastycznej ma postać: \[\begin{equation} H_0:\;p=0,5\quad\textrm{vs.}\quad H_{1}:\;p\neq 0,5 \tag{8.6} \end{equation}\] gdzie \(p\) określa prawdopodobieństwo tego, że obserwacje w grupie pierwszej są zazwyczaj mniejsze niż w grupie drugiej.

Wynika z tego, że prawdopodobieństwo zdarzenia przeciwnego (obserwacje w grupie pierwszej \(x\) są zazwyczaj większe niż w grupie drugiej \(y\)) jest także równe \(0,5\). Zatem w hipotezie zerowej zakładamy, że wartości w obu próbkach mają porównywalne wartości tzn. wartości z pierwszej próbki nie mają tendencji do mniejszych/większych wartości niż w próbce drugiej. Estymację tego prawdopodobieństwa można dokonać w dwojaki sposób: \[\begin{equation} \hat{p}=P(x<y)+0,5\cdot P(x=y)\quad\mbox{lub} \quad \hat{p}=\frac{\bar{r}_2-(n_2+1)\cdot 0,5}{n_1} \tag{8.7} \end{equation}\] gdzie \(\bar{r}_2\) to średnia ranga dla drugiej zmiennej a rangi są liczone na podstawie próbki zbiorczej.

Warto dodać, że na podstawie estymatora prawdopodobieństwa \(\hat{p}\) można obliczyć statystykę testu sumy rang Wilcoxona na podstawie wzoru: \[\begin{equation} W=(1-\hat{p})n_1n_2 \tag{8.8} \end{equation}\]

Statystyka testu Brunera-Munzela: \[\begin{equation} BM=\frac{n_1n_2(\bar{r}_1-\bar{r}_2)}{(n_1+n_2)\sqrt{n_1s_1^2+n_2s_2^2}} \tag{8.9} \end{equation}\] ma rozkład t-Studenta ze stopniami swobody według formuły: \[\begin{equation} df_{\mathrm{Satterthwaite}}=\frac{(d_1+d_2)^2}{\frac{d_1}{n_1-1}+\frac{d_2}{n_2-1}} \tag{8.10} \end{equation}\] gdzie \(d_k=n_k\cdot s^2_k\) to iloczyn liczebności próby \(n_k\) oraz wariancji \(s_k^2\) dla każdej \(k\)-tej grupy.

Wariancja jest zdefiniowana w następujący sposób: \[\begin{equation} s_k^2=\frac{1}{n_k-1}\sum_{i=1}^{n_k}\left(r_{ki}-w_{ki}-\bar{r}_k+\frac{n_k+1}{2}\right)^2 \tag{8.11} \end{equation}\] gdzie: \(\bar{r}_k\) oznacza średnią rangę \(k\)-tej grupy z próbki zbiorczej, \(r_k\) to rangi dla \(k\)-tej grupy z próbki zbiorczej, \(w_k\) to rangi dla \(k\)-tej grupy.

import warnings
warnings.filterwarnings("ignore")
import scipy.stats as stats
import PyNonpar
from PyNonpar import *
    
x = stats.norm.rvs(0, 1, size=20, random_state=2305).tolist()
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101).tolist()

res = PyNonpar.twosample.brunner_munzel_test(x,y)
print("BM= %.4f, df= %.4f, pvalue= %.4f" % (res[1],res[2],res[3]))

## BM= 2.9380, df= 20.6352, pvalue= 0.0080

Test Kruskala-Wallisa. Nieparametrycznym odpowiednikiem analizy wariancji jest test Kruskala-Walisa jako rozszerzenie testu sumy rang Wilcoxona na kilka grup. W tej metodzie zakładamy, że próbki pochodzą z tego samego rozkładu o dowolnym kształcie. Oznacza to, że rozkład w grupach nie musi być normalny ale w dalszym ciągu zakładamy homoskedastyczność wariancji. Dokładny rozkład statystyki Kruskala-Wallisa można przybliżać za pomocą metod permutacyjnych lub takich dystrybuant jak: chi-kwadrat, F-Snedecora oraz beta (Meyer i Seaman 2013).

Statystyka testowa: \[\begin{equation} \chi^2_{KW}=\left(1-\frac{\sum_{i=1}^{c}(t_i^3-t_i)}{n^3-n}\right)^{-1}\left[\frac{12}{n(n+1)}\left(\sum_{j=1}^{k}\frac{R_j^2}{n_j}\right)-3(n+1)\right] \tag{8.12} \end{equation}\] gdzie: \(n\) to liczebność z wszystkich \(k\) grup, \(n_j\) to liczebność w \(j\)-tej grupie, \(R_j\) to suma rang w \(j\)-tej grupie, \(c\) to liczba grup pomiarów wiązanych, \(t_i\) to liczba pomiarów wiązanych w \(i\)-tej grupie pomiarów wiązanych.

Poniżej implementacja wersji permutacyjnej testu Kruskala-Wallisa:

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

y = np.concatenate((stats.expon.rvs(0, 1, size=20, random_state=2305),
                    stats.expon.rvs(0.5, 1, size=20, random_state=4101),
                    stats.expon.rvs(1, 1, size=20, random_state=4026)))
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)

kw = pg.kruskal(dv='y', between='g', data=pd.DataFrame({'y':y,'g':g}))

H = kw["H"][0]
B = 1000
h = [list(pg.kruskal(dv='y', between='g',\
          data=pd.DataFrame({'y':y,'g':np.random.choice(g,size=60)}))["H"])[0]\
          for i in range(B)]
perm = np.greater(h,[H]).mean()
kw["p-perm"] = perm
print(kw)

##         Source  ddof1     H     p-unc  p-perm
## Kruskal      g      2  7.61  0.022254   0.012

ANOVA-rank. Warto zauważyć, że problem heterogeniczności wariancji można uwzględnić za pomocą testu Brunner-Dette-Munk (Brunner, Dette, i Munk 1997) w którym można także testować interakcję w dwuczynnikowej analizie wariancji. Jednak ta metoda nie jest dostępna w pakietach scipy.stats oraz pingouin. Alternatywą może być zastosowanie procedury wykorzystującej rozkład F-Snedecora która polega na porangowaniu danych i zastosowaniu klasycznej metody ANOVA. Innym rozwiązaniem może być wykorzystanie ważonej metody najmniejszych kwadratów lub odpornych błędów standardowych z wykorzystaniem funkcji statsmodels.stats.anova.anova_lm.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

y = np.concatenate((stats.expon.rvs(0, 1, size=20, random_state=2305),
                    stats.expon.rvs(0.5, 1, size=20, random_state=4101),
                    stats.expon.rvs(1, 1, size=20, random_state=4026)))
r = stats.rankdata(y)
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)
d = pd.DataFrame({"y":y,"g":g,"r":r})
print(pg.anova(dv='r', between='g', data=d))

##   Source  ddof1  ddof2      F     p-unc    np2
## 0      g      2     57  4.221  0.019527  0.129

Dalsza analiza. Po odrzuceniu hipotezy zerowej w teście Kruskala-Wallisa można dokonać bardziej szczególowej analizy czyli przeprowadzić porównania wielokrotne. Popularnym rozwiązaniem jest zastosowanie serii testów sumy rang Wilcoxona. Ta metoda jest dostępna dzięki funkcji pingouin.pairwise_ttests z zaznaczeniem opcji parametric=False. Jednak szerszy zestaw testów post hoc dla grup niezależnych znajdziemy w pakiecie scikit-posthocs. Poniżej przykład testu Conovera.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg
from scikit_posthocs import posthoc_conover

y = np.concatenate((stats.expon.rvs(0, 1, size=20, random_state=2305),
                    stats.expon.rvs(0.5, 1, size=20, random_state=4101),
                    stats.expon.rvs(1, 1, size=20, random_state=4026)))
r = stats.rankdata(y)
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)
d = pd.DataFrame({"y":y,"g":g,"r":r})
print(posthoc_conover(d, val_col='y', group_col='g', p_adjust='holm'))

##           1.0       2.0       3.0
## 1.0 -1.000000  0.221096  0.015936
## 2.0  0.221096 -1.000000  0.221096
## 3.0  0.015936  0.221096 -1.000000

8.3 Porównanie wariancji

Test Z-diff / Test Z-ratio. Jeśli chcemy porównać dwie wariancje to rozważamy hipotezy statystyczne o postaci: \[\begin{equation} H_0:\;\sigma_1^2=\sigma_2^2\quad\mbox{vs.}\quad H_1:\;\sigma_1^2\neq\sigma_2^2 \tag{8.13} \end{equation}\] Zauważmy, że powyższą hipotezę statystyczną można sprowadzić do zapisu: \[\begin{equation} H_{0}:\;\sigma^2_1/\sigma^2_2=1\quad\textrm{vs.}\quad H_{1}:\;\sigma^2_1/\sigma^2_2\neq1 \tag{8.14} \end{equation}\] Statystyka testowa: \[\begin{equation} Z_{ratio}=\frac{(s^2_1/s^2_2)-1}{SE_{ratio}} \tag{8.15} \end{equation}\] gdzie: \(SE_{ratio}=\frac{1}{s_2^{2}}\sqrt{SE_1^2+r_0^2\cdot SE_2^2}\) to błąd standardowy ilorazu dwóch wariancji oraz \(SE=\sqrt{s^2/n}\) to błąd standardowy wariancji \(s^2\) dla przekształconej zmiennej \((x_i-\bar{x})^2\), \(r_0^2\) to iloraz wariancji podniesiony do drugiej potęgi.

import scipy.stats as stats
import numpy as np

x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)
  
z1 = (x-np.mean(x))**2
z2 = (y-np.mean(y))**2
ratV = np.var(x,ddof=1)/np.var(y,ddof=1)
SE = np.sqrt(np.var(z1,ddof=1)/len(x)+ratV**2*np.var(z2,ddof=1)/len(y))/np.var(y,ddof=1)
conf = [stats.norm.ppf(i,ratV,SE) for i in [0.025,0.975]]
h0 = 1
p = stats.norm.cdf(h0,ratV,SE)

print("iloraz wariancji:",ratV,", błąd:",SE)
print("95% przedział ufności:",conf)
print("\nH0: rVar = %.0f vs. H1: rVar != %.0f" % (h0,h0))
print("p-wartość:",2*min(p,1-p))

## iloraz wariancji: 0.2555508988591833 , błąd: 0.0975515172257403
## 95% przedział ufności: [0.06435343845949357, 0.446748359258873]
## 
## H0: rVar = 1 vs. H1: rVar != 1
## p-wartość: 2.3314683517128287e-14

Równoważnym zapisem powyższych hipotez statystycznych (8.13) oraz (8.14) będzie zapis: \[\begin{equation} H_{0}:\;\sigma^2_1-\sigma^2_2=0\quad\textrm{vs.}\quad H_{1}:\;\sigma^2_1-\sigma^2_2\neq0 \tag{8.16} \end{equation}\] Statystyka testowa: \[\begin{equation} Z_{diff}=\frac{(s^2_1-s^2_2)-0}{SE_{diff}} \tag{8.17} \end{equation}\]

gdzie: \(SE_{diff}=\sqrt{SE_{1}^2+\rho^2\cdot SE_{2}^2}\) to błąd standardowy różnicy dwóch wariancji oraz \(SE=\sqrt{s^2/n}\) to błąd standardowy wariancji \(s^2\) dla przekształconej zmiennej \((x_i-\bar{x})^2\), \(\rho^2\) to opcjonalny parametr do osłabienia/wzmocnienia udziału drugiej wariancji.

import scipy.stats as stats
import numpy as np

x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)

z1 = (x-np.mean(x))**2
z2 = (y-np.mean(y))**2
difV = np.var(x,ddof=1)-np.var(y,ddof=1)
SE = np.sqrt(np.var(z1,ddof=1)/len(x)+1*np.var(z2,ddof=1)/len(y))
conf = [stats.norm.ppf(i,difV,SE) for i in [0.025,0.975]]
h0 = 0
p = stats.norm.cdf(h0,difV,SE)

print("różnica wariancji:",difV,", błąd:",SE)
print("95% przedział ufności:",conf)
print("\nH0: dVar = %.0f vs. H1: dVar != %.0f" % (h0,h0))
print("p-wartość:",2*min(p,1-p))

## różnica wariancji: -3.8307446353496024 , błąd: 1.267036271883883
## 95% przedział ufności: [-6.314090095347914, -1.347399175351292]
## 
## H0: dVar = 0 vs. H1: dVar != 0
## p-wartość: 0.0024995998590693347

Test Bartletta / Test Levene. Badanie równości wariancji można wykonać również za pomocą testu Fligner-Killen lub testu Levene które w przeciwieństwie do testu Bartletta są mało wrażliwe na odchylenia od rozkładu normalnego w próbkach. Przeważnie są one stosowane do badania równości kilku wariancji ale nic nie stoi na przeszkodzie aby wykorzystać je do porównania dwóch wariancji na podstawie hipotezy (8.13). Dodajmy, że test Levene i Fligner-Killeen mogą występować w trzech wariantach tzn. za parametr lokalizacji można przyjąć średnią, średnią uciętą lub medianę. Taki wybór oferują funkcje scipy.stats.levene oraz scipy.stats.fligner. Natomiast funkcja pingouin.homoscedasticity jako parametr lokalizacji stosuje medianę. Jeśli zmienne mają rozkład normalny to podawany jest wynik testu Bartletta.

import scipy.stats as stats
import numpy as np
import pingouin as pg

x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)

print(stats.levene(x,y))
print(stats.fligner(x,y))
print(stats.bartlett(x,y))
print(pg.homoscedasticity(x,y))

## LeveneResult(statistic=8.994663638939507, pvalue=0.0047575161444811925)
## FlignerResult(statistic=6.954492770953616, pvalue=0.008360900023288296)
## BartlettResult(statistic=8.019535289470467, pvalue=0.004627544281207727)
## (False, 0.005)

8.4 Porównanie rozkładów

Test normalności Andersona-Darlinga. Założenie normalności zmiennych to jedno z głównych założeń w klasycznej statystyce. W związku z tym zostało opracowanych wiele metod porównywania dystrybuanty empirycznej z rozkładem normalnym. Jednym z bardziej popularnych rozwiązań jest test Shapiro-Wilka który wymaga aby liczebność próby nie przekraczała 5000 elementów. Inne metody jak np. test Jarque-Bera, D’Agostino-Pearsona czy Andersona-Darlinga nie mają tego ograniczenia. Dodajmy jeszcze, że wysoka moc testu może być dobrym uzasadnieniem wyboru konkretnej metody (Biecek 2013). Przykładowo test normalności Andersona-Darlinga może być ciekawą alternatywą dla testu Shapiro-Wilka w przypadku wielomodalności lub występowania grubych ogonów (Biecek 2017, str. 244-246).

Statystyka testu Andersona-Darlinga ma postać: \[\begin{equation} AD = -n-\frac{1}{n}\sum_{i=1}^n(2i-1)\big(\ln(z_i)+\ln(1-z_{n+1-i})\big) \tag{8.18} \end{equation}\] gdzie: \(z_i\) to wartości wyznaczone na podstawie dystrybuanty rozkładu normalnego \(\Phi(x_i,\bar{x},s)\) dla posortowanych rosnąco elementów próby \(x_i\).

W przypadku badania normalności o nieznanych parametrach \(\mu\) oraz \(\sigma\) jest stosowana poprawka: \[\begin{equation} A1=AD\left(1+\frac{0,75}{n}+\frac{2,25}{n^2}\right) \tag{8.19} \end{equation}\]

Weryfikację hipotezy zerowej można wykonać w oparciu o otrzymaną p-wartość która jest uzależniona od wartości statystyki testu (8.19).

• jeżeli \(A1 < 0,2\) to:

\[\begin{equation} p-value=1-\exp(-13,436+101,14\,A1-223,73\,A1^2) \tag{8.20} \end{equation}\]

• jeżeli \(0,2\leq A1<0,34\) to:

\[\begin{equation} p-value=1-\exp(-8,318+42,796\,A1-59,938\,A1^2) \tag{8.21} \end{equation}\]

• jeżeli \(0,34\leq A1 < 0,6\) to:

\[\begin{equation} p-value=\exp(0,9177-4,279\,A1-1,38\,A1^2) \tag{8.22} \end{equation}\]

• jeżeli \(A1\geq 0,6\) to:

\[\begin{equation} p-value= \exp(1,2937-5,709\,A1+0,0186\,A1^2) \tag{8.23} \end{equation}\]

import scipy.stats as stats
from statsmodels.stats.diagnostic import normal_ad
  
x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)
ad = normal_ad(x)
print('AD = %.4f, p-value = %.4f' % (ad[0],ad[1]))

## AD = 0.2568, p-value = 0.6848

W stosunkowo prosty sposób można wygenerować wartości krytyczne na podstawie wzoru:

\[\begin{equation} A_{crit}=a\left(1-\frac{b}{n}-\frac{d}{n^2}\right) \tag{8.24} \end{equation}\] gdzie \(n\) to liczebności próby oraz \(a\), \(b\) i \(d\) to parametry które zależą od poziomu istotności \(\alpha\): \[\begin{equation} \begin{array}{c|llllll} \alpha & 0.005 & 0.01 & 0.025 & 0.05 & 0.10 & 0.20\\ \hline\hline a & 1.1578 & 1.0348 & 0.8728 & 0.7514 & 0.6305 & 0.5091\\ b & 1.063 & 1.013 & 0.881 & 0.795 & 0.750 & 0.756\\ d & 1.34 & 0.93 & 0.94 & 0.89 & 0.80 & 0.39 \end{array} \tag{8.25} \end{equation}\] Poniżej przykład wygenerowania różnych wartości krytycznych testu normalności Andersona-Darlinga.

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

def q(alpha=0.05,n=10):
    if alpha == 0.005:\
    return 1.1578*(1-1.063/n-1.34/n**2)
    elif alpha == 0.01:\
    return 1.0348*(1-1.013/n-0.93/n**2)
    elif alpha == 0.025:\
    return 0.8728*(1-0.881/n-0.94/n**2)
    elif alpha==0.05:\
    return 0.7514*(1-0.795/n-0.89/n**2)
    elif alpha == 0.1:\
    return 0.6305*(1-0.750/n-0.80/n**2)
    elif alpha == 0.2:\
    return 0.5091*(1-0.756/n-0.39/n**2)

n = [20,50,100,150,300,900,1500]
q0_01  = [q(alpha=0.01, n=i) for i in n]
q0_025 = [q(alpha=0.025,n=i) for i in n]
q0_05  = [q(alpha=0.05, n=i) for i in n]
q0_1   = [q(alpha=0.1,  n=i) for i in n]
print(pd.DataFrame({'1%':q0_01,'2.5%':q0_025,'5%':q0_05,'10%':q0_1},index=n))

##             1%      2.5%        5%       10%
## 20    0.979981  0.832302  0.719860  0.605595
## 50    1.013450  0.857093  0.739185  0.620841
## 100   1.024221  0.865029  0.745359  0.625721
## 150   1.027769  0.867637  0.747388  0.627325
## 300   1.031295  0.870228  0.749401  0.628918
## 900   1.033634  0.871945  0.750735  0.629974
## 1500  1.034101  0.872287  0.751001  0.630185

Poniżej wygenerujemy w sposób symulacyjny wartości krytyczne dla \(n=20\):

from statsmodels.stats.diagnostic import normal_ad
import numpy as np
import scipy.stats as stats
import pandas as pd
  
res = [normal_ad(stats.norm.rvs(0, 1.5, size=20))[0] for i in range(10000)]
q = np.percentile(res,[99,97.5,95,90])
print(pd.DataFrame({'1%':q[0],'2.5%':q[1],'5%':q[2],'10%':q[3]},index=['20']))

##          1%      2.5%        5%       10%
## 20  0.97439  0.826135  0.711582  0.600183

Test zgodności Andersona-Darlinga. Oprócz rozkładu normalnego dystrybuantę empiryczną można porównywać również z innymi dystrybuantami teoretycznymi. Do badania zgodności z rozkładami ciągłymi można wykorzystać test Kołmogorowa który został zaimplementowany do funkcji scipy.stats.kstest. Alternatywą do tego rozwiązania jest test Andersona-Darlinga dostępny dzięki funkcji scipy.stats.anderson. W tej implementacji zamiast p-wartości są podawane wartości krytyczne \(A_{crit}\) które określają granicę prawostronnego obszaru odrzucenia. Inaczej mówiąc jeśli \(AD>A_{crit}\) to hipotezę zerową o zadanej dystrybuancie należy odrzucić. Dodajmy jeszcze, że wartości krytyczne zależą od liczebności próby \(n\), poziomu istotności \(\alpha\) oraz roważanego rozkładu. W zaimplementowanej funkcji można założyć rozkład np. normalny, wykładniczy, logistyczny, gumbela.

W przypadku rozkładu normalnego wartości krytyczne są obliczane za pomocą wzoru: \[\begin{equation} A_{crit}=k(\alpha)/\left(1 + \frac{4}{n} - \frac{25}{n^2}\right) \tag{8.26} \end{equation}\] gdzie wartość współczynnika \(k(\alpha)\) jest uzależniona od tego czy znane są parametry rozkładu. Poniżej wykaz współczynników dla różnych wariantów. \[\begin{equation} \begin{array}{c|c|lllll} \mbox{wariant} & \alpha & 0.15 & 0.10 & 0.05 & 0.025 & 0.01\\ \hline\hline N(\mu,\sigma) & k(\alpha) & 1.610 & 1.993 & 2.492 & 3.070 & 3.857\\ N(?,?) & k(\alpha) & 0.576 & 0.656 & 0.787 & 0.918 & 1.092 \end{array} \tag{8.27} \end{equation}\]

Wartości krytyczne dla dwóch pozostałych rozkładów np. wykładniczegi oraz logistycznego obliczamy za pomocą wzoru: \[\begin{equation} A_{crit}=k(\alpha)/ \left(1 + \frac{v}{n}\right) \tag{8.27} \end{equation}\] gdzie odpowiednie współczynniki \(k(\alpha)\) dla danej liczebności próby \(n\) są przedstawione poniżej: \[\begin{equation} \begin{array}{c|l|c|llll} \mbox{wariant} & v & \alpha & 0.10 & 0.05 & 0.025 & 0.01\\ \hline\hline Expon & 0.6 & k(\alpha) & 1.065 & 1.325 & 1.587 & 1.934\\ Logist & 0.25 & k(\alpha) & 0.56 & 0.657 & 0.765 & 0.901 \end{array} \tag{8.28} \end{equation}\]

Poniżej przykład jak wygenerować tablicę z wartościami krytycznymi dla rozkładu normalnego, wykładniczego i logistycznego przy założonym \(\alpha=0,05\) oraz różnych liczebności próby \(n\). Dodajmy jeszcze, że funkcja podaje wartości krytyczne dla przypadku gdy parametry rozkładu nie są znane i trzeba je oszacować.

import scipy.stats as stats
import numpy as np
import pandas as pd
  
n = [10,20,30,50,70,100,150,300]
nor = [stats.anderson(stats.norm.rvs(size=i), dist='norm')[1][2] for i in n]
exp = [stats.anderson(stats.norm.rvs(size=i), dist='expon')[1][2] for i in n]
logis = [stats.anderson(stats.norm.rvs(size=i), dist='logistic')[1][2] for i in n]
gumbel = [stats.anderson(stats.norm.rvs(size=i), dist='gumbel')[1][2] for i in n]
print(pd.DataFrame({'nor_0.05':nor,'exp_0.05':exp,
                    'logis_0.05':logis,'gumbel_0.05':gumbel},index=n))

##      nor_0.05  exp_0.05  logis_0.05  gumbel_0.05
## 10      0.684     1.265       0.644        0.712
## 20      0.692     1.302       0.652        0.725
## 30      0.712     1.315       0.655        0.730
## 50      0.736     1.325       0.657        0.736
## 70      0.748     1.330       0.658        0.739
## 100     0.759     1.333       0.658        0.742
## 150     0.767     1.336       0.659        0.745
## 300     0.777     1.338       0.659        0.748

Poniżej przykład wywołania funkcji scipy.stats.anderson w celu zbadania zgodności rozkładu empirycznego z rozkładem normalnym.

import scipy.stats as stats
import pandas as pd

x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)

ad = stats.anderson(x, dist='norm')
print(pd.DataFrame({'20':ad[1]},index=ad[2]/100).T)
print('\nad: %.4f'% (ad[0]))

##     0.150  0.100  0.050  0.025  0.010
## 20  0.506  0.577  0.692  0.807   0.96
## 
## ad: 0.2568

Otrzymana wartość statystyki testu Andersona-Darlinga nie przekracza wartości krytycznej nawet dla \(\alpha=0.15\) więc brak jest podstaw do odrzucenia hipotezy zerowej. Warto dodać, że w tym teście można zweryfikować hipotezę zerową w oparciu o p-wartość wyznaczoną w sposób analityczny (Jäntschi i Bolboacă 2018) lub symulacyjny. Jedna z propozycji (Marsaglia i Marsaglia 2004) została zaimplementowana do pakietu ADGofTest dla środowiska R.

import scipy.stats as stats
import numpy as np

x = stats.norm.rvs(0, 1.5, size=20, random_state=2305)
  
A = stats.anderson(x,dist='norm')
ad = [stats.anderson(np.random.choice(x,size=len(x),replace=True), dist='norm')[0] \
      for i in range(1000)]
print('AD: %.4f, p-value: %.4f' % (A[0], np.mean(np.greater(ad,[A[0]]))))

## AD: 0.2568, p-value: 0.9820

Test zgodności Cressie-Read. Badanie zgodności rozkładu empirycznego z założonym rozkładem teoretycznym (ciągłym lub dyskretnym) o zdefiniowanych parametrach można wykonać za pomocą testu chi-kwadrat lub jego uogólnionej wersji tzn. testu Cressie-Reada. Do tego celu można wykorzystać odpowiednio funkcje scipy.stats.chisquare oraz scipy.stats.power_divergence w których argumentami są wartości empiryczne f_obs=fi oraz teoretyczne f_exp=ei. Jeśli w teście Cressie-Reada ustalimy, że parametr lambda będzie równy "1" lub przypiszemy mu nazwę "pearson" to zostanie wykonany test chi-kwadrat.

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

x = stats.poisson.rvs(1.5, size=80, random_state=2305)
    
def goodfitPois(x): 
    t = pd.Series(x).value_counts(sort=False)
    fi = t.values
    xi = list(t.index)
    pi = [stats.poisson.pmf(i,np.mean(x)) for i in xi]
    pi.append(1-sum(pi))
    ei = np.asarray(pi) * len(x)
    e = ei[-1]+ei[-2]
    ei = ei[:-2]
    ei = list(ei)
    ei.append(e)
    return stats.power_divergence(fi, ei, ddof=1, lambda_=2/3)
  
print(goodfitPois(x))

## Power_divergenceResult(statistic=0.47455822372954887, pvalue=0.9759300084452401)

Test Andersona-Darlinga dla k prób. Test Kołmogorowa-Smirnowa (funkcja scipy.stats.ks_2samp) to częsty wybór do weryfikacji hipotezy zerowej w której zakładamy, że dwie dystrybuanty są takie same. Inaczej mówiąc badamy czy dwie zmienne losowe pochodzą z tego samego ciągłego rozkładu o takich samych parametrach. Gdy porównujemy dwie próbki warto zwrócić uwagę także na test Eppsa-Singletona który jest zaimplementowany do funkcji scipy.stats.epps_singleton_2samp. Ta metoda charakteryzuje się między innymi tym, że ma większą moc niż test Kołmogorowa-Smirnowa oraz może porównywać także rozkłady dyskretne. Alternatytwnym rozwiązaniem jest test Andersona-Darlinga (funkcja scipy.stats.anderson_ksamp) który można stosować dla dwóch lub większej liczby próbek z rozkładu ciągłego. Dodatkowo po odrzuceniu hipotezy zerowej można sprawdzić które zmienne różnią się między sobą za pomocą testów post hoc – porównania wielokrotne.

import scipy.stats as stats
import numpy as np
import pandas as pd
from scikit_posthocs import posthoc_anderson

x1 = stats.expon.rvs(0, 1, size=20, random_state=2305)
x2 = stats.expon.rvs(0.5, 1, size=20, random_state=4101)
x3 = stats.expon.rvs(1, 1, size=20, random_state=4026)
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)
d = pd.DataFrame({"y":np.concatenate((x1,x2,x3)),"g":g})
adk = stats.anderson_ksamp([x1,x2,x3])

print('ad = %.4f, p-wartość = %.4f' % (adk[0],adk[2]),'\n')
print(posthoc_anderson(d,val_col='y',group_col='g'))

## ad = 3.8139, p-wartość = 0.0066 
## 
##           1.0       2.0       3.0
## 1.0 -1.000000  0.103108  0.001627
## 2.0  0.103108 -1.000000  0.121164
## 3.0  0.001627  0.121164 -1.000000

8.5 Moc testu

Test t-Studenta. Standardowy rozkład t-Studenta ma swój ogólniejszy odpowiednik tzn. niecentralny rozkład t-Studenta z dodatkowym parametrem ncp – non-centrality parameter. Dla \(ncp = 0\) niecentralny rozkład t-Studenta jest tożsamy z centralnym rozkładem t-Studenta – takie szczególne przypadki mają także rozkłady chi-kwadrat oraz F-Snedecora. Rozkłady niecentralne są często wykorzystywane do obliczania mocy testów np. funkcja pingouin.power_ttest oblicza moc testu t-Studenta dla dwóch niezależnych prób (test dwustronny) według wzoru: \[\begin{equation} \mbox{moc}=P(T\leq t_{crit},df,ncp) \tag{8.29} \end{equation}\] gdzie: \(t_{crit}\) to kwantyl rzędu \(1-\alpha/2\) z rozkładu t-Studenta o stopniach swobody \(df=2n-2\) oraz \(ncp=|d|\cdot\sqrt{\frac{n}{2}}\) to non-centrality parameter.

Wielkość efektu \(d\) można obliczyć na podstawie wzoru: \[\begin{equation} d=t_{val}\cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}. \tag{8.30} \end{equation}\] gdzie: \(t_{val}\) to statystyka testu t-Studenta dla dwóch niezależnych prób, \(n_1\) oraz \(n_2\) to liczenbość odpowiednio dla pierwszej i drugiej próby.

import scipy.stats as stats
import numpy as np
import pingouin as pg
  
x = stats.norm.rvs(0, 1, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)

tval, n1, n2 = stats.ttest_ind(x,y)[0], len(x), len(y)
d = pg.compute_effsize_from_t(tval, nx=n1, ny=n2, eftype='cohen')
power = pg.power_ttest(d=d, n=len(x), contrast='two-samples')
print("Efekt: %.4f, Moc: %.4f" % (d, power))

## Efekt: -1.0903, Moc: 0.9192

ANOVA. Moc testu dla klasycznej wersji jednoczynnikowej ANOVY można obliczyć za pomocą nie centralnego rozkładu F-Snedecora czyli z dodatkowym parametrem \(ncp\): \[\begin{equation} \mbox{moc}=P(F\geq F_{crit},df_1,\, df_2,\, ncp) \tag{8.31} \end{equation}\] gdzie: \(F_{crit}\) to kwantyl rzędu \(1-\alpha\) z rozkładu F-Snedecora o stopniach swobody \(df1=k-1\), \(df2=n-3\) oraz \(npc=f^2 N\) to non-centrality parameter.

Wielkość efektu \(f\) można obliczyć według formuły: \[\begin{equation} f=\sqrt{\frac{\sum_{i=1}^{k}p_i(\mu_i-\mu)^2}{\sigma^2}}=\sqrt{\frac{SS_{betveen}}{MS_{residuals}\cdot N}} \tag{8.32} \end{equation}\] gdzie: \(p_i=n_i/N\), \(n_i\) to liczba obserwacji w \(i\)-tej grupie, \(N\) to suma wszystkich obserwacji, \(\mu_i\) to średnia w \(i\)-tej grupie, \(\mu\) to ogólna średnia, \(\sigma^2\) to wariancja błędu w obrębie grupy (\(MS_{residuals}\) - mean squares for resuduals).

Metoda zaimplementowana do funkcji pingouin.power_anova bazuje na obliczeniu wielkości efektu \(f\) według wzoru: \[\begin{equation} f=\sqrt{\frac{\eta^2}{1-\eta^2}} \tag{8.33} \end{equation}\] gdzie: \(\eta^2\) to wielkość efektu dla jednoczynnikowej analizy wariancji która jest tożsama z współczynnikiem determinacji \(R^2\) dla regresji liniowej. \[\begin{equation} \eta^2=\frac{df_1\cdot F}{df_1\cdot F+df_2}\quad\mbox{lub} \quad \eta^2= \frac{SS_{between}}{SS_{residuals}+SS_{betveen}} \tag{8.34} \end{equation}\] gdzie: \(SS_{between}\) to suma kwadratów dla czynnika, \(SS_{total}\) to suma kwadratów dla czynnika oraz reszt.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

y = np.concatenate((stats.expon.rvs(0, 1, size=20, random_state=2305),
                    stats.expon.rvs(0.5, 1, size=20, random_state=4101),
                    stats.expon.rvs(1, 1, size=20, random_state=4026)))
g = np.repeat(np.linspace(1,3,3), [20,20,20], axis=0)
d = pd.DataFrame({"y":y,"g":g})
eta2 = pg.anova(dv='y', between='g', data=d)['np2']
f = np.sqrt(eta2/(1-eta2))
power = pg.power_anova(eta=eta2, k=3, n=20)[0]
print("Efekt: %.4f, Moc: %.4f" % (f[0], power))

## Efekt: 0.1534, Moc: 0.1636

Jeśli metoda analityczna do obliczenia mocy wybranego testu nie jest dostępne to wygodnym rozwiązaniem może być symulacja komputerowa.

import numpy as np
import scipy.stats as stats
import pandas as pd
import pingouin as pg

x = stats.norm.rvs(0, 1, size=20, random_state=2305)
y = stats.norm.rvs(1.5, 2, size=20, random_state=4101)
z = np.concatenate((x,y))
g = np.repeat(np.linspace(1,2,2), [20,20], axis=0)
  
def pvalA(x,y):
    nx = len(x)
    ny = len(y)
    z = np.concatenate((np.random.choice(x,nx),np.random.choice(y,ny)))
    g = np.repeat(np.linspace(1,2,2), [nx,ny], axis=0)
    return pg.welch_anova(dv='z', between='g',data=pd.DataFrame({"z":z,"g":g}))['p-unc'][0]

m = [pvalA(x,y) for i in range(1000)]
print("Moc: ", np.less(m,[0.05]).mean(), "dla 1000 symulacji testu Welch-Anova")

## Moc:  0.914 dla 1000 symulacji testu Welch-Anova

Bibliografia

Biecek, P. 2013. „Wybrane testy normalności”. SmarterPoland. http://tofesi.mimuw.edu.pl/~cogito/smarterpoland/samouczki/testyNormalnosci/testyNormalnosci.pdf.

Biecek, P. 2017. Przewodnik po pakiecie R. Oficyna Wydawnicza "GIS". http://www.biecek.pl/R/.

Brunner, Edgar, Holger Dette, i Axel Munk. 1997. „Box-Type Approximations in Nonparametric Factorial Designs”. Journal of the American Statistical Association 92 (440). Taylor & Francis: 1494–1502. https://doi.org/10.1080/01621459.1997.10473671.

Brunner, E., i U. Munzel. 2000. „The Nonparametric Behrens‐Fisher Problem: Asymptotic Theory and a Small‐Sample Approximation”. Biometrical Journal 42 (1): 17–25. https://doi.org/10.1002/(SICI)1521-4036(200001)42:1<17::AID-BIMJ17>3.0.CO;2-U.

Derrick, B., i P. White. 2016. „Why Welch’s test is Type I error robust”. The Quantitative Methods for Psychology 12 (1). TQMP: 30–38. https://www.tqmp.org/RegularArticles/vol12-1/p030/.

Divine, George W., H. James Norton, Anna E. Barón, i Elizabeth Juarez-Colunga. 2018. „The Wilcoxon–Mann–Whitney Procedure Fails as a Test of Medians”. The American Statistician 0 (0). Taylor & Francis: 1–9. https://doi.org/10.1080/00031305.2017.1305291.

Jäntschi, Lorentz, i Sorana D. Bolboacă. 2018. „Computation of Probability Associated with Anderson–Darling Statistic”. Mathematics 6 (6). https://doi.org/10.3390/math6060088.

Marsaglia, John, i George Marsaglia. 2004. „Evaluating the Anderson-Darling Distribution”. Journal of Statistical Software 09 (luty). https://doi.org/10.18637/jss.v009.i02.

Meyer, J. Patrick, i Michael A. Seaman. 2013. „A Comparison of the Exact Kruskal-Wallis Distribution to Asymptotic Approximations for All Sample Sizes up to 105”. The Journal of Experimental Education 81 (2). Routledge: 139–56. https://doi.org/10.1080/00220973.2012.699904.

Neubert, Karin, i Edgar Brunner. 2007. „A studentized permutation test for the non-parametric Behrens–Fisher problem”. Computational Statistics and Data Analysis 51 (10): 5192–5204. https://doi.org/10.1016/j.csda.2006.05.024.