Ein wenig schneller zur simulierten Nullverteilung

Ein Nullhypothesentest ist schnell geschrieben. Will man den approximativen Weg gehen, so hilft R einem mit entsprechenden Tests mit einfachen Befehlen. Nimmt man MOSAIC dazu, so bekommt man u.a. für den Test auf Anteils- oder Mittelwerte sogar einen sehr einfachen, weil einheitlichen, Syntax.

Zwei Beispiele für approximative Hypothesentests mit MOSAIC

Laden wir unsere Testdaten, die tipping Daten wie folgt:

library(mosaic)
download.file("https://goo.gl/whKjnl", destfile = "tips.csv")
tips <- read.csv2("tips.csv")
set.seed(2009)

Dann erstellen wir zwei Forschungsfragen:

  1. Ist der mittlere Frauenanteil unter der Bezahler*innen zu den Zeitpunkten Lunch und Dinner gleich?
  2. Ist der mittlere Rechnungsbetrag zu den Zeitpunkten Lunch und Dinner gleich?

Im ersten Fall ist die Hypothese schnell geschrieben:

\[ H_0 : \pi_{\text{Lunch}} = \pi_{\text{Dinner}} \quad\text{vs.}\quad H_1 : \pi_{\text{Lunch}} \neq \pi_{\text{Dinner}} \] Der appoximative Test mit R und MOSAIC lautet nun:

prop.test(sex ~ time, success = "Female", data = tips)
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  tally(sex ~ time)
## X-squared = 9.3438, df = 1, p-value = 0.002237
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.36602563 -0.07247705
## sample estimates:
##    prop 1    prop 2 
## 0.2954545 0.5147059

Ähnlich sieht es für den zweiten Fall aus. Die Hypothese lautet hier:

\[ H_0 : \mu_{Lunch} = \mu_{Dinner} \quad\text{vs.}\quad H_1 : \mu_{Lunch} \neq \mu_{Dinner} \]

Der dazugehörige Test lautet dann:

t.test(total_bill ~ time, data = tips)
## 
##  Welch Two Sample t-test
## 
## data:  total_bill by time
## t = 3.123, df = 143.29, p-value = 0.002167
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.331877 5.925088
## sample estimates:
## mean in group Dinner  mean in group Lunch 
##             20.79716             17.16868

Simulation der Nullverteilung mit MOSAIC

Ein anderer Weg ist es die Stichprobe selber zu nutzen um daraus eine Verteilung der Nullhypothese (die Nullverteilung) ableiten zu können. Im ersten Fall schaut man sich die Anteilsunterschiede an, wenn man die (potentielle) Abhänigkeit von der Tageszeit (Lunch und Dinner) künstlich “abschaltet”:

set.seed(2009)
NullVtlgAntwert <- do(10000) * diffprop(sex ~ shuffle(time), 
    success = "Female", data = tips)
gf_histogram(~diffprop, nint = 25, data = NullVtlgAntwert)

Schaut man sich nun die Lage der Anteilsdifferenz der Stichprobe \(\hat\pi=`diffprop(sex ~ time, success="Female", data=tips)`\) in Bezug auf diese Nullverteilung geometrisch an, so kann man schon einen ersten Eindruck erlangen, ob die Nullhypothese abzulehnen ist oder nicht:

diffpropdach <- diffprop(sex ~ time, success = "Female", data = tips)
gf_histogram(~diffprop, nint = 25, data = NullVtlgAntwert) + 
    geom_vline(xintercept = diffpropdach, color = "blue")

Offenbar ist \(\hat\pi\) kein sehr häufiges Ereignis.

Der “p-Wert” ist auch leicht ermittelt:

pvalue_aw <- prop(~abs(diffprop) >= abs(diffpropdach), data = NullVtlgAntwert)
pvalue_aw
##   TRUE 
## 0.0027

Ähnlich sie die Situation im zweien Fall aus. Mittels weniger Befehler erzeugen wir die Verteilung.

set.seed(2009)
NullVtlgMittelwert <- do(10000) * diffmean(total_bill ~ shuffle(time), 
    data = tips)
gf_histogram(~diffmean, nint = 25, data = NullVtlgMittelwert)

Und können im Anschluss die Mittelwertsdifferent der Stichprobe geometrisch einordnen:

diffmeandach <- diffmean(total_bill ~ time, data = tips)
gf_histogram(~diffmean, nint = 25, data = NullVtlgMittelwert) + 
    geom_vline(xintercept = diffmeandach, color = "blue")

Auch den “p-Wert” können wir nun leicht bestimmen:

pvalue_mw <- prop(~abs(diffmean) >= abs(diffmeandach), data = NullVtlgMittelwert)
pvalue_mw
##   TRUE 
## 0.0039

Das Problem – Zeit

Das Problem bei der Simulation ist die Zeit, die R braucht um die Nullverteilungen zu generieren. Das liegt im wesendlichen an Mosaic. Mit den Routinen aus FastSimNullDistR lassen sich die Nullverteilungen deutlich schneller berechnen. Ein Vergleich:

library(mosaic)
library(mosaicCore)
source("https://raw.githubusercontent.com/NMarkgraf/FastSimNullDistR/master/R/fastSimNullDistRMean.R")
source("https://raw.githubusercontent.com/NMarkgraf/FastSimNullDistR/master/R/fastSimNullDistRProp.R")
source("https://raw.githubusercontent.com/NMarkgraf/FastSimNullDistR/master/R/fastSimNullDistR_work.R")
set.seed(2009)
system.time(NullDistMosaic_aw <- do(10000) * diffprop(sex ~ shuffle(time), 
    success = "Female", data = tips))
##        User      System verstrichen 
##      14.822       0.308      15.426
set.seed(2009)
system.time(NullDistFSNDR_aw <- fastSimNullDistRProp(sex ~ time, 
    success = "Female", data = tips))
##        User      System verstrichen 
##       1.415       0.034       1.476
set.seed(2009)
system.time(NullDistMosaic_mw <- do(10000) * diffmean(total_bill ~ 
    shuffle(time), data = tips))
##        User      System verstrichen 
##      15.943       0.377      16.683
set.seed(2009)
system.time(NullDistFSNDR_mw <- fastSimNullDistRMean(total_bill ~ 
    time, data = tips))
##        User      System verstrichen 
##       1.300       0.040       1.373

Das mit den beiden Routinen aus FastSimNullDistR die gleichen Ergebnisse zu erwarten sind, sie also ein “(quasi-)drop-in-replacements” der Mosaic Rouninen darstellen, kann man an den folgenden QQ-Plots erkennen:

qqplot(NullDistFSNDR_aw$diffprop, NullDistMosaic_aw$diffprop)

qqplot(NullDistFSNDR_mw$diffmean, NullDistMosaic_mw$diffmean)

Woher kommt die Geschwindigkeit?

Schaut man sich den Quellcode von Mosaic an, wird einem schnell klar, dass es zwar didaktisch sinnvoll ist die unabhängige Variable mit shuffle() zu bearbeiten, nicht aber programmiertechnisch. Und wenn, dann nicht in dem man die ganze Datenzeile für die Berechnung kopiert. Statt also \(10\,000\) mal die ganzen Daten im Speicher zu kopieren wäre es doch sinnvoll einen Index auf die unveränderten daten indirekt über einen Index zuzugreifen. Und genau da liegt der Zugang der Routinen. Nur dieser Zugriffsindex wird geshuffelt und das spart Speicherplatz und Rechenzeit.

comments powered by Disqus