Sampling distributions

A sampling distribution answers a question about what else could have happened. If the empty model were true — if condition made no difference to tips at all — how big a b1 would we see just from the luck of which tables happened to land in which group?

do() answers that by running the same estimate many times over data that has been shuffled, so any relationship between the two variables is broken by construction.

set.seed(42)
sdob1 <- do(1000) * b1(shuffle(Tip) ~ Condition, data = TipExperiment)
head(sdob1)
#>           b1
#> 1 -5.0454545
#> 2 -4.0454545
#> 3 -7.9545455
#> 4 -3.5000000
#> 5 -1.7727273
#> 6  0.2272727

Each row is one b1 from one shuffle. Plotted, they are the range of estimates the empty model produces.

gf_histogram(~b1, data = sdob1, binwidth = 1)

Marking off regions

The distribution part functions turn a region of that distribution into a fill aesthetic. middle() marks the values we would expect to see often.

gf_histogram(~b1, data = sdob1, fill = ~ middle(b1, .95), binwidth = 1)

tails() marks the same cutoffs with the opposite coloring: the 5% most extreme values, the ones that would be surprising under the empty model.

gf_histogram(~b1, data = sdob1, fill = ~ tails(b1, .95), binwidth = 1)

outer() marks the same region but takes the tail proportion directly, which reads more naturally when you are thinking in terms of how much is out rather than how much is in.

gf_histogram(~b1, data = sdob1, fill = ~ outer(b1, .05), binwidth = 1)

For a directional hypothesis, all of the 5% goes in one tail.

gf_histogram(~b1, data = sdob1, fill = ~ upper(b1, .05), binwidth = 1)

gf_histogram(~b1, data = sdob1, fill = ~ lower(b1, .05), binwidth = 1)

Marking the cutoffs themselves

Filling a region colors it. show_cutoffs() marks where it ends, and it takes the region directly, so a plain histogram can have cutoffs without being colored.

gf_histogram(~b1, data = sdob1, binwidth = 1) %>%
  show_cutoffs(middle(b1, .95))

Because the region is named rather than read off the fill, the marked region and the colored one can differ: color the middle 95%, mark where the middle 99% ends.

gf_histogram(~b1, data = sdob1, fill = ~ middle(b1, .95), binwidth = 1) %>%
  show_cutoffs(middle(b1, .99))

Two calls put two sets on one plot, which is how you see what changing alpha costs you.

gf_histogram(~b1, data = sdob1, binwidth = 1) %>%
  show_cutoffs(middle(b1, .95)) %>%
  show_cutoffs(middle(b1, .99), color = "firebrick")

The wider pair is harder to fall outside, so a result that clears it is rarer.

Bootstrapping instead of shuffling

Swap shuffle() for resample() and the same pattern produces a bootstrap confidence interval rather than a null distribution. Nothing is shuffled, so the relationship stays intact; what varies is which observations are drawn.

set.seed(42)
sdob1_boot <- do(1000) * b1(Tip ~ Condition, data = resample(TipExperiment))
gf_histogram(~b1, data = sdob1_boot, fill = ~ middle(b1, .95), bins = 100)

Framing the distribution with its DGP

gf_squareplot() draws the same distribution as countable squares, and show_dgp() frames it with the data generating process: the population model on the top axis, the sample estimate on the bottom, and a marker at the null hypothesis. Both overlays take a plot and give one back, so the same frame goes just as readily on a gf_histogram() of the same shuffles.

With only ten shuffles the mean of the distribution can land well away from the null.

set.seed(42)
small <- do(10) * b1(shuffle(Tip) ~ Condition, data = TipExperiment)
gf_squareplot(~b1, data = small, binwidth = 2) %>%
  gf_lims(x = c(-30, 30)) %>%
  gf_refine(expand_limits(y = 10)) %>%
  show_mean() %>%
  show_dgp()

With a hundred it settles close to the null, which is what the empty model predicts. expand_limits(y = ) sets a floor under the count axis, so a distribution with a low peak still gets a full-height panel; a hundred shuffles stack higher than ten and the axis grows to fit them.

set.seed(42)
larger <- do(100) * b1(shuffle(Tip) ~ Condition, data = TipExperiment)
gf_squareplot(~b1, data = larger, binwidth = 2) %>%
  gf_lims(x = c(-30, 30)) %>%
  gf_refine(expand_limits(y = 10)) %>%
  show_mean() %>%
  show_dgp()