Mathematik für Biologiestudierende II¶

Sommersemester 2024

07.05.2024

© 2024 Prof. Dr. Rüdiger W. Braun

In [1]:
import numpy as np
import pandas as pd
from scipy import stats
import seaborn as sns
sns.set_theme()

Heteroskedastizität¶

  • Die ANOVA vergleicht die Varianzen innerhalb der einzelnen Gruppen mit der Varianz im gesamten Datensatz, um die Unterschiede zwischen den Gruppen zu untersuchen
  • À priori geht das erstmal nur, wenn die Varianzen innerhalb der Gruppen gleich sind
  • Ein Datensatz ist heteroskedastisch, wenn die verschiedenen Gruppen unterschiedliche Varianz haben

Der Levene-Test¶

Der Levene-Test testet auf Gleichheit der Varianzen

Beispiel: Meerschweinchenzähne

  • Drei Gruppen von Meerschweinchen, je nach täglicher Gabe an Vitamin C
    • kleine Dosis
    • mittlerer Dosis
    • große Dosis
  • Nach 42 Tagen wird die Zahnlänge bestimmt

Quelle: The Statistics of Bioassay

In [2]:
small_dose = np.array([
    4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
    15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
])

medium_dose = np.array([
    16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
    19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
])

large_dose = np.array([
    23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
    25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
])

Test auf Heteroskedastizität

In [3]:
stats.levene(small_dose, medium_dose, large_dose)
Out[3]:
LeveneResult(statistic=0.6457341109631506, pvalue=0.5280694573759905)
  • Der p-Wert ist 0.53. Hetereskedastizizät kann nicht nachgewiesen werden.
  • Wir fahren mit der ANOVA fort
In [4]:
stats.f_oneway(small_dose, medium_dose, large_dose)
Out[4]:
F_onewayResult(statistic=67.41573785674247, pvalue=9.532727011699946e-16)
  • Gabe von Vitamin C hat Einfluss auf das Zahnwachstum
  • Als nächstes würde man eine post-hoc Analyse machen

Beispiel: Barsche¶

In [5]:
df = pd.read_csv('barsche.csv')
df.head()
Out[5]:
Art Länge
0 gestreift 9.890006
1 gestreift 9.343944
2 gestreift 9.867069
3 gestreift 10.302781
4 gestreift 10.066964
In [6]:
sns.boxplot(data=df, x="Art", y="Länge");
No description has been provided for this image

Sieht heteroskedastisch aus

In [7]:
ds = df[df.Art=='gestreift'].Länge
dl = df[df.Art=='gefleckt'].Länge
db = df[df.Art=='blau'].Länge
dr = df[df.Art=='braun'].Länge
In [8]:
stats.levene(ds, dl, dr, db)
Out[8]:
LeveneResult(statistic=13.459492972830807, pvalue=1.3472893996510424e-07)

Probleme beim Test auf Heteroskedastizität¶

  • Die Nullhypothese beim Levene-Test ist

$H_0$: Die Daten sind homoskedastisch

  • Ein Hypothesentest "beweist" nie die Nullhypothese
    • bei starken Indizien dagegen lehnt er sie ab
    • bei starken Indizien dafür behält er sie bei
    • bei unklaren Indizien behält er sie auch bei
  • um zu erkennen, ob der Levene-Test Heteroskedastizität überhaupt erkennen kann, wäre eine Poweranalyse für den Levene-Test nötig, das ist aber unrealistisch
  • auch das andere Extrem ist möglich: Der Stichprobenumfang ist so groß, dass kleine Unterschiede schon signifikant werden

👁️eyeballing (scharfes Hinsehen)

Alexander-Govern-Test¶

Wenn die Daten heteroskedastisch, aber normalverteilt sind, dann rechnet man einen Alexander-Govern-Test

In [9]:
stats.alexandergovern(ds, dl, dr, db)
Out[9]:
AlexanderGovernResult(statistic=113.40810114676775, pvalue=2.02668339537414e-24)

Im homoskedastischen Fall ist der p-Wert des Alexander-Govern-Tests meist schlechter als der von f_oneway

In [10]:
stats.alexandergovern(small_dose, medium_dose, large_dose)
Out[10]:
AlexanderGovernResult(statistic=56.82538049315831, pvalue=4.57641509985116e-13)
In [11]:
stats.f_oneway(small_dose, medium_dose, large_dose)
Out[11]:
F_onewayResult(statistic=67.41573785674247, pvalue=9.532727011699946e-16)

Post-hoc Analyse¶

  • Der t-Test kann nur gerechnet werden, wenn die Varianzen der zu vergleichenden Datensätze übereinstimmen
  • Im heteroskedastischen Fall ist das nicht der Fall
  • Man rechnet dann einen Welch-Test
  • In scipy ist der Welch-Test wie folgt implementiert
In [12]:
stats.ttest_ind(db, dr, equal_var=False)
Out[12]:
TtestResult(statistic=9.647287139857793, pvalue=1.2289650206522807e-13, df=57.645418945809595)
  • Problem: Arbeitet nicht mit MultiComparison zusammen
In [13]:
# ganz kleines Programm

def welch(x, y):
    return stats.ttest_ind(x, y, equal_var=False)
In [14]:
from statsmodels.sandbox.stats.multicomp import MultiComparison
In [15]:
muc = MultiComparison(df.Länge, df.Art)
In [16]:
muc.allpairtest(welch, method='holm')[0]
Out[16]:
Test Multiple Comparison welch FWER=0.05 method=holm alphacSidak=0.01, alphacBonf=0.008
group1 group2 stat pval pval_corr reject
blau braun 9.6473 0.0 0.0 True
blau gefleckt 5.8735 0.0 0.0 True
blau gestreift 18.752 0.0 0.0 True
braun gefleckt 0.7068 0.484 0.484 False
braun gestreift 8.4956 0.0 0.0 True
gefleckt gestreift 3.3453 0.002 0.004 True

Normalverteilungsannahmen¶

  • Sowohl f_oneway als auch alexandergovern liefern nur für normalverteilte Daten richtige Ergebnisse
  • Normalverteilungsannahmen prüfen wir mit dem Q-Q-Plot
In [17]:
import statsmodels.api as sm
pp = sm.ProbPlot(db)
pp.qqplot();
No description has been provided for this image

ausreichende Übereinstimmung

Flügellängen von Libellen in mm (erfundene Daten)

In [18]:
df = pd.read_csv('libellen.csv')
df.head()
Out[18]:
Art Länge
0 graue 4.908840
1 graue 5.016692
2 graue 4.382700
3 graue 4.847548
4 graue 5.523503
In [19]:
df.Art.value_counts()
Out[19]:
Art
graue         60
grüne         60
ägyptische    60
Bilker        60
Name: count, dtype: int64
In [20]:
dg = df[df.Art=='graue'].Länge
du = df[df.Art=='grüne'].Länge
da = df[df.Art=='ägyptische'].Länge
dB = df[df.Art=='Bilker'].Länge
In [21]:
pp = sm.ProbPlot(dB)
pp.qqplot();
No description has been provided for this image

nicht normalverteilt

Kruskal-Wallis-Test¶

im Fall nicht normalverteilter Daten rechnet man den Kruskal-Wallis-Test

In [22]:
stats.kruskal(dg, du, da, dB)
Out[22]:
KruskalResult(statistic=16.028153526970982, pvalue=0.0011190121329907562)

Post-hoc Analyse¶

Das nicht-parametrische Analogon zum unverbundenen t-Test ist der Mann-Whitney-Test

In [23]:
muc = MultiComparison(df.Länge, df.Art)
In [24]:
muc.allpairtest(stats.mannwhitneyu, method='holm')[0]
Out[24]:
Test Multiple Comparison mannwhitneyu FWER=0.05 method=holm alphacSidak=0.01, alphacBonf=0.008
group1 group2 stat pval pval_corr reject
Bilker graue 2360.0 0.0033 0.0199 True
Bilker grüne 2167.0 0.0544 0.1088 False
Bilker ägyptische 2217.0 0.0288 0.0864 False
graue grüne 1260.0 0.0046 0.0228 True
graue ägyptische 1259.0 0.0046 0.0228 True
grüne ägyptische 1561.0 0.2106 0.2106 False
In [25]:
sns.boxplot(data=df, x='Art', y='Länge');
No description has been provided for this image