Mathematik für Biologiestudierende II¶

Sommersemester 2024

16.04.2024

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

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

Normalverteilungsannahmen¶

konservative Tests¶

  • Der $t$-Test verwendet eine Verteilungsannahme: Daten müssen normalverteilt sein
  • Tests, die auch bei Verletzung der Verteilungsannahmen noch gute Ergebnisse liefern, heißen konservativ
  • Der $t$-Test ist konservativ

Q-Q-Plot¶

  • Mit dem Quantil-Quantil-Plot kann man auf graphischem Wege beurteilen, ob Messwerte Realisierungen einer normalverteilten Zufallsvariablen sind
  • Man trägt dazu auf der $x$-Achse die Quantile der Standardnormalverteilung und auf der $y$-Achse die Quantile der Beobachtungsdaten auf
  • Wenn diese Punkte annähernd auf einer Geraden liegen, sind die Daten näherungsweise normalverteilt, ansonsten nicht
  • es gibt auch Testverfahren, um auf Normalverteilungsannahmen zu testen
  • nicht ganz klar, wie sinnvoll das ist

Beispieldaten aus Lektion 12

In [2]:
u = "https://www.math.uni-duesseldorf.de/~braun/bio2324/data/schadstoffe.csv"
df = pd.read_csv(u, index_col=0)
df
Out[2]:
Messstelle Konzentration
0 5 0.000867
1 3 0.000490
2 1 0.000589
3 1 0.000950
4 4 0.001152
... ... ...
75 5 0.000918
76 3 0.000528
77 3 0.000961
78 4 0.001272
79 3 0.001012

80 rows × 2 columns

In [3]:
import statsmodels.api as sm
In [4]:
pp = sm.ProbPlot(df.Konzentration)
pp.qqplot();
No description has been provided for this image

Wunderbar normalverteilt

Die Daten aus dem synthetischen Medikamentenexperiment aus Lektion 13

In [5]:
df = pd.read_csv('treatment.csv', index_col=0)
In [6]:
pp = sm.ProbPlot(df.t0)
pp.qqplot();
No description has been provided for this image

Die Daten sind nicht normalverteilt, weil ich oben bei 100 abgeschnitten hatte

Beispiel Galapagos Inseln¶

Ein Datensatz zum Buch "Linear Models with Python" von Faraway

In [7]:
df = pd.read_csv("galapagos.csv", index_col=0)
df
Out[7]:
Species Area Elevation Nearest Scruz Adjacent
Baltra 58 25.09 346 0.6 0.6 1.84
Bartolome 31 1.24 109 0.6 26.3 572.33
Caldwell 3 0.21 114 2.8 58.7 0.78
Champion 25 0.10 46 1.9 47.4 0.18
Coamano 2 0.05 77 1.9 1.9 903.82
Daphne.Major 18 0.34 119 8.0 8.0 1.84
Daphne.Minor 24 0.08 93 6.0 12.0 0.34
Darwin 10 2.33 168 34.1 290.2 2.85
Eden 8 0.03 71 0.4 0.4 17.95
Enderby 2 0.18 112 2.6 50.2 0.10
Espanola 97 58.27 198 1.1 88.3 0.57
Fernandina 93 634.49 1494 4.3 95.3 4669.32
Gardner1 58 0.57 49 1.1 93.1 58.27
Gardner2 5 0.78 227 4.6 62.2 0.21
Genovesa 40 17.35 76 47.4 92.2 129.49
Isabela 347 4669.32 1707 0.7 28.1 634.49
Marchena 51 129.49 343 29.1 85.9 59.56
Onslow 2 0.01 25 3.3 45.9 0.10
Pinta 104 59.56 777 29.1 119.6 129.49
Pinzon 108 17.95 458 10.7 10.7 0.03
Las.Plazas 12 0.23 94 0.5 0.6 25.09
Rabida 70 4.89 367 4.4 24.4 572.33
SanCristobal 280 551.62 716 45.2 66.6 0.57
SanSalvador 237 572.33 906 0.2 19.8 4.89
SantaCruz 444 903.82 864 0.6 0.0 0.52
SantaFe 62 24.08 259 16.5 16.5 0.52
SantaMaria 285 170.92 640 2.6 49.2 0.10
Seymour 44 1.84 147 0.6 9.6 25.09
Tortuga 16 1.24 186 6.8 50.9 17.95
Wolf 21 2.85 253 34.1 254.7 2.33
In [8]:
pp = sm.ProbPlot(df.Area)
pp.qqplot();
No description has been provided for this image

Nicht-parametrische Tests¶

Beispiel für Situationen, in denen man nicht-parametrische Tests macht:

  • Wenn die Verteilungsannahmen nicht erfüllt sind
  • Wenn die Stichprobenumfänge zu klein sind
  • Wenn die Zufallsvariable diskret ist

Vergleich zweier Erwartungswerte bzw. zweier Mediane¶

Vergeich parametrisch nicht-parametrisch
mit Referenzwert t-Test für verbundene Stichproben Wilcoxon-Test
vorher-nachher t-Test für verbundene Stichproben Wilcoxon-Test
verschiedene Populationen t-Test für unverbundene Stichproben Mann-Whitney-U-Test

Wilcoxon-Signed-Rank-Test¶

Den Wilcoxon Test verwendet man zum Vergleich der Mediane verbundener Datensätze, wenn die Normalverteilungsannahme nicht gesichert ist.

Er ist ein Rangtest: Das bedeutet, dass nur eine Rolle spielt, ob Werte größer sind als andere, aber nicht, um wie viel.

In [9]:
from scipy import stats
In [10]:
df = pd.read_csv(u, index_col=0)
df['Referenz'] = 0.08 / 100
  • Wir vergleichen die Konzentration mit dem Referenzwert 0.08
  • Die Nummer der Messstelle benötigen wir erst in einer späteren Auswertung
In [11]:
res = stats.wilcoxon(df.Konzentration, df.Referenz, alternative="greater")
res
Out[11]:
WilcoxonResult(statistic=2169.0, pvalue=0.004229703509534525)

Zum Vergleich:

In [12]:
stats.ttest_rel(df.Konzentration, df.Referenz, alternative="greater")
Out[12]:
TtestResult(statistic=2.768040010585661, pvalue=0.0035114445640696246, df=79)
  • Der p-Wert ist ein kleines bisschen schlechter
  • Dieser Unterschied ist unerheblich
  • Im allgemeinen ist nicht klar, ob die unberechtigte Nutzung eines parametrischen Tests den p-Wert verbessert oder verschlechtert

Wir bestimmen jetzt auch noch die Effektstärke gemäß Cohen's r

Dazu ist es erforderlich, den Test noch einmal mit method="approx" zu rechnen, um die z-Statistik zu bekommen

In [13]:
res = stats.wilcoxon(df.Konzentration, df.Referenz, alternative="greater", method="approx")
res
Out[13]:
WilcoxonResult(statistic=2169.0, pvalue=0.004229703509534525)
In [14]:
res.zstatistic
Out[14]:
2.6331616685404655
In [15]:
n = df.Konzentration.count()
n
Out[15]:
80
In [16]:
r = abs(res.zstatistic / np.sqrt(n))
r
Out[16]:
0.2943964243301625

Interpretation der Effektstärke¶

r-Wert Interpretation
0.1 geringer Effekt
0.3 mittlerer Effekt
0.5 starker Effekt

Wir haben also einen mittleren Effekt beobachtet

Mann-Whitney-U-Test¶

Den Mann-Whitney Test verwendet man zum Vergleich der Mediane unverbundener Datensätze, wenn die Normalverteilungsannahme nicht gesichert ist

In [17]:
df = pd.read_csv('treatment.csv', index_col=0)
dfv = df[df.Treatment=='Verum']
dfp = df[df.Treatment=='Placebo']
In [18]:
res = stats.mannwhitneyu(dfv.Difference, dfp.Difference, alternative='greater')
res
Out[18]:
MannwhitneyuResult(statistic=416814.5, pvalue=0.008213617991035723)

zum Vergleich

In [19]:
stats.ttest_ind(dfv.Difference, dfp.Difference, alternative='greater')
Out[19]:
TtestResult(statistic=2.314493969317715, pvalue=0.010377396661800722, df=1767.0)
In [20]:
pp = sm.ProbPlot(dfv.Difference)
pp.qqplot();
No description has been provided for this image

Auch hier müssen wir für die Effektstärke die Statistik noch einmal rechnen, und zwar mit vertauschten Datensätzen

In [21]:
res2 = stats.mannwhitneyu(dfp.Difference, dfv.Difference, alternative='less')
res2
Out[21]:
MannwhitneyuResult(statistic=365519.5, pvalue=0.008213617991035723)

Der $p$-Wert ist derselbe, aber die Statistik ist eine andere

Die Formel für die Effektstärke nach Cohen's r ist $$ r = 1 - \frac{2U}{n_1 \cdot n_2} $$

wobei $U$ die kleinere der beiden Statistiken ist und $n_1$ und $n_2$ die Stichprobenumfänge in den beiden Populationen sind

In [22]:
U = res2.statistic
n1 = dfp.Difference.count()
n2 = dfv.Difference.count()
In [23]:
r = 1 - 2*U/(n1*n2)
r
Out[23]:
0.06556662499648491

Ein sehr geringer Effekt