Ein wenig schneller zur simulierten Nullverteilung

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 approximative 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 between group Dinner and group Lunch 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ängigkeit 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=0.2192513\) 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 ebenfalls leicht zu ermitteln:

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

## prop_TRUE 
##    0.0018

Mit einem Anteilswert (p-Wert) von 0.0018 zweigen wir wie selten das Ereignis unter der \(H_0\) ist.

Ähnlich sieht die Situation im zweien Fall aus. Mit Hilfe weniger Befehle erzeugen wir die Nullverteilung:

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 Mittelwertsdifferenz 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 wieder leicht bestimmen:

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

## prop_TRUE 
##    0.0047

Das Problem – Zeit

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

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

df.diffprop <- data_frame(diffprop = c(NullDistFSNDR_aw$diffprop,
    NullDistMosaic_aw$diffprop), type = c(rep("FSNDR", 10000),
    rep("mosaic", 10000)))

## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

gf_qq(~diffprop, color = ~type, data = df.diffprop)

df.diffmean <- data_frame(diffmean = c(NullDistFSNDR_mw$diffmean,
    NullDistMosaic_mw$diffmean), type = c(rep("FSNDR", 10000),
    rep("mosaic", 10000)))
gf_qq(~diffmean, color = ~type, data = df.diffmean)

# qqplot(NullDistFSNDR_aw\(diffprop, NullDistMosaic_aw\)diffprop) gf_qq(FSNDR ~ Mosaic, data=df) # qqplot(NullDistFSNDR_mw\(diffmean, NullDistMosaic_mw\)diffmean) gf_qq(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 sinnvoller mit Hilfe eines Index auf die unveränderten Daten zuzugreifen. Und genau das machen die zwei Routinen. Es wird also nur dieser Zugriffsindex wird geshuffelt und das spart Speicherplatz und deutlich auch Rechenzeit.