Less Volume, More Creativity

[pullquote]I find that using the mosaic package keeps the focus on what is important (thinking about where randomization enters and why, and what to do with the distribution once we have it) while hiding distracting details (looping structures, extracting the relevant information from R objects, etc.). [/pullquote]The biggest key to using R well is to provide a lot of creative opportunity with as little R as possible. There can be a lot of ways to skin a cat in R, what you need is a systematic approach that allows you to economically achieve your goals. If technology is the most difficult thing in your course, then your technology is too hard or your questions are too simple (probably both). On the other hand, if your students are able to guess how to do new things in R before you show them, then you will know you are on the right track.

An example: The Lady Tasting Tea

I typically introduce inference on the first day doing a hypothesis test for a proportion “informally”. I don’t use terms like null hypothesis, or p-value, but all the big ideas are introduced.

Often I use some variation on Fisher’s Lady Tasting Tea example because I find that students remember it well and I can refer back to it as a motivating example for the rest of the semester. I find the example works equally well in an Intro Stats course as in my upper level courses, and I like that it connects them in a small way to the history of the discipline.

Let’s see how we can test whether the lady is just guessing given that she correctly identifies 9 of 10 randomly prepared cups of tea. Guessing is like flipping a coin, so mosaic provides a coin flipper:

require(mosaic)    # load the mosaic package
rflip()            # think: random flip

## 
## Flipping 1 coin [ Prob(Heads) = 0.5 ] ...
## 
## H
## 
## Number of Heads: 1 [Proportion Heads: 1]

Since our design calls for 10 cups of tea, we need to flip 10 coins. Fortunately, we don’t have to use rflip() ten times and record the results manually. We just ask for 10 flips.

rflip(n = 10)

## 
## Flipping 10 coins [ Prob(Heads) = 0.5 ] ...
## 
## H H H T T T H H H T
## 
## Number of Heads: 6 [Proportion Heads: 0.6]

Now we need to do this a lot of times. The do() function provides an easy syntax and is clever about how it stores the results, in this case storing the number of flips, number of heads, number of tails, and proportion that are heads.

do(1000) * rflip(n = 10)   # do 1000 times 

##    n heads tails prop
## 1 10     6     4  0.6
## 2 10     5     5  0.5
## 3 10     4     6  0.4

## <997 lines omitted>

Of course, we should look at this some other way than by having it all scroll past on our screen. Let’s save the results and look at numerical and graphical summaries.

GuessingLadies <-
  do(1000) * rflip(n = 10)  # simulate 1000 guessing ladies
tally(~heads, data = GuessingLadies)

## 
##   0   1   2   3   4   5   6   7   8   9  10 
##   2   5  50 111 202 241 216 105  57  10   1

histogram(~heads, data = GuessingLadies, width = 1)

Based on our simulation, it appears that someone can get 9 or 10 just by guessing only about 1.1% of the time.

If we want to have a more precise estimate of this proportion (the p-value), we can increase the number of times we do() things. We could also 999 (or 9999) instead of 1000 and include the observed statistic in the null distribution as recommended in the article Tim Hesterberg previews in this SBI blog post.

Grand Tour of Randomization Distributions

If we want to test for a different proportion (e.g., 0.25), we simply tell rflip() the probability of obtaining heads (p) and the size of our sample (n = 100 in the example below):

Proportion.null <- do(1000) * rflip(n = 100, p = 0.25) 

##     n heads tails prop
## 1 100    31    69 0.31
## 2 100    25    75 0.25

<998 rows omitted>

Null distributions for permutation tests can be simulated by shuffling the appropriate labels. Here are examples for tests involving the difference between two proportions, the difference between two means, or a linear model:

Galton2 <- Galton %>% group_by(family) %>% sample_n(1)  # one from each family
Null1 <- do(1000) *  diffprop(  homeless ~ shuffle(sex),       data = HELPrct)
Null2 <- do(1000) *  diffmean(  height ~ shuffle(sex),       data = Galton2)
Null3 <- do(1000) *  lm(  age ~ shuffle(substance), data = HELPrct)  # ANOVA
Null4 <- do(1000) *  lm(  father ~ shuffle(mother),    data = Galton2)  # regression

Looking at the first few rows of Null4, we see that do() extracts and records several useful bits of information about a linear model:

##   Intercept  mother sigma r.squared     F
## 1      79.8 -0.1641  2.60   0.02173 4.332
## 2      74.1 -0.0736  2.62   0.00437 0.856
## 3      63.9  0.0849  2.62   0.00582 1.141

We can use the slope (labeled mother), r², or F as our test statistic.

One null distribution that is somewhat more challenging is the null distribution for a test of a single mean. The challenge is to determine how to simulate data with a mean equal to that of the null hypothesis. This can be done non-parametrically by resampling from a shifted version of the empirical distribution. I find it simplest to do this by changing H₀: μ = 98.6 into H₀: μ – 98.6 = 0 and sampling with replacement from the sample:

require(Lock5withR)
Null <- do(1000) * mean( ~ (BodyTemp - 98.6), data = resample(BodyTemp50))

Both the shifting and the introduction of sampling with replacement in the context of hypothesis testing may be more distracting than they are worth since a confidence interval for a mean is generally more useful anyway, but it rounds out our tour.

Bootstrap, too

Generating bootstrap distributions to form confidence intervals is similar. Typically, instead of shuffling a variable we resample the entire data set (possibly within groups) to produce a bootstrap distribution. Here is a bootstrap distribution for the difference in mean heights of men and women using Galton’s data, resampling separately from the males and females:

HeightBySex.boot <- 
  do(1000) * diffmean( height ~ sex, data = resample(Galton, groups = sex) )
histogram(~diffmean, data = HeightBySex.boot)

Simple confidence intervals can be computed using percentiles

cdata(0.95, ~diffmean, data = HeightBySex.boot)  # central 95%

##       low        hi central.p 
##      4.78      5.46      0.95

or the bootstrap standard error, which is just the standard deviation of the bootstrap distribution.

sd(~diffmean, data=HeightBySex.boot)

## [1] 0.166

Eventually, I may show my students how to automate generating the confidence intervals from the boostrap distribution.

confint(HeightBySex.boot, method=c("percentile", "se"))

##       name lower upper level   method estimate margin.of.error  df
## 1 diffmean  4.78  5.46  0.95 quantile       NA              NA  NA
## 2 diffmean  4.80  5.45  0.95   stderr     5.12           0.326 897

Striking a Balance

I find that using the mosaic package keeps the focus on what is important (thinking about where randomization enters and why, and what to do with the distribution once we have it) while hiding distracting details (looping structures, extracting the relevant information from R objects, etc.). There are other R packages (resample and boot, for example) that are faster and use more sophisticated methods, but hide some of the things I want students to see and think about when they are learning how simulation-based inference works.

The R-related cognitive load for students to do SBI using mosaic amounts to

• a common template for graphical and numerical summaries, which is needed for non-SBI things as well

• rflip(), shuffle() and resample() for generating random data
• do() for repeating things

These can be combined to handle a wide variety of situations. That’s a lot of SBI creativity for little volume.

Additional Resources

The following resources are provided by Project MOSAIC (Danny Kaplan, Nick Horton, and Randy Pruim) for those interested in learning more about using R to teach statistics.

• eCOTS 2014 Workshop: “Effective Teaching using R, RStudio, and the MOSAIC Package”.
• Less Volume, More Creativity – an introduction to the mosaic package emphasizing the use of a common template for numerical summaries, graphical summaries, and modeling.
• Randomization based Inference [view] [download] – examples of simulation-based inference with mosaic based on the 2013 USCOTS Randomization Bake Off
• Start Teaching with R
• A Student’s Guide to R
• Lock5withR – a companion to Unlocking the Power of Data by Five Locks (with Lana Park)
• ISIwithR – a companion to Introduction to Statistical Inference by Nathan Tintle et al. (with Lana Park)