{"id":894,"date":"2015-06-24T13:26:28","date_gmt":"2015-06-24T17:26:28","guid":{"rendered":"https:\/\/www.causeweb.org\/sbi\/?p=894"},"modified":"2015-06-24T18:33:06","modified_gmt":"2015-06-24T22:33:06","slug":"how-i-teach-sbi-using-r-2","status":"publish","type":"post","link":"https:\/\/www.causeweb.org\/sbi\/?p=894","title":{"rendered":"How I teach SBI using R"},"content":{"rendered":"

\"RPruim\"<\/a><\/em>R Pruim,\u00a0Calvin College<\/strong><\/p>\n

I\u2019m not writing to convince you that you should<\/em> use R<\/em> to teach simulation-based inference (SBI). My goal is to convince you that you can<\/em> use R<\/em> for SBI, even with students (and instructors) who have never used R<\/em> before. Along the way I\u2019ll mention some guiding principles and illustrate some tools that my colleagues Danny Kaplan and Nick Horton and I have assembled in the mosaic<\/code> R package<\/a> to make SBI (and EDA and traditional inference procedures) much easier.[pullquote]The biggest key to using R<\/em> well is to provide a lot of creative opportunity with as little R<\/em> as possible.[\/pullquote]<\/p>\n

<\/p>\n

\n

If you are unfamiliar with R<\/em>, this post will be too short to give you a good introduction, but I hope it will make you curious enough to find out more. See the references for places to learn more about R<\/em> and teaching with R<\/em>.<\/p>\n

\n

Less Volume, More Creativity<\/h3>\n

[pullquote]I find that using the mosaic<\/code> 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<\/em> objects, etc.). [\/pullquote]The biggest key to using R<\/em> well is to provide a lot of creative opportunity with as little R<\/em> as possible. There can be a lot of ways to skin a cat in R<\/em>, 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<\/em> before you show them, then you will know you are on the right track.
\n<\/p>\n<\/div>\n

\n

An example: The Lady Tasting Tea<\/h3>\n

I typically introduce inference on the first day doing a hypothesis test for a proportion \u201cinformally\u201d. I don\u2019t use terms like null hypothesis, or p-value, but all the big ideas are introduced.<\/p>\n

Often I use some variation on Fisher\u2019s Lady Tasting Tea<\/a> 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.<\/p>\n

Let\u2019s 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<\/code> provides a coin flipper:<\/p>\n

require(mosaic)    # load the mosaic package\r\nrflip()            # think: random flip<\/code><\/pre>\n
## \r\n## Flipping 1 coin [ Prob(Heads) = 0.5 ] ...\r\n## \r\n## H\r\n## \r\n## Number of Heads: 1 [Proportion Heads: 1]<\/code><\/pre>\n
Since our design calls for 10 cups of tea, we need to flip 10 coins. Fortunately, we don\u2019t have to use rflip()<\/code> ten times and record the results manually. We just ask for 10 flips.<\/p>\n
rflip(n = 10)<\/code><\/pre>\n
## \r\n## Flipping 10 coins [ Prob(Heads) = 0.5 ] ...\r\n## \r\n## H H H T T T H H H T\r\n## \r\n## Number of Heads: 6 [Proportion Heads: 0.6]<\/code><\/pre>\n
Now we need to do this a lot of times. The do()<\/code> 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.<\/p>\n
do(1000) * rflip(n = 10)   # do 1000 times <\/code><\/pre>\n
##    n heads tails prop\r\n## 1 10     6     4  0.6\r\n## 2 10     5     5  0.5\r\n## 3 10     4     6  0.4<\/code><\/pre>\n
## <997 lines omitted><\/code><\/pre>\n
Of course, we should look at this some other way than by having it all scroll past on our screen. Let\u2019s save the results and look at numerical and graphical summaries.<\/p>\n
GuessingLadies <-\r\n  do(1000) * rflip(n = 10)  # simulate 1000 guessing ladies\r\ntally(~heads, data = GuessingLadies)<\/code><\/pre>\n
## \r\n##   0   1   2   3   4   5   6   7   8   9  10 \r\n##   2   5  50 111 202 241 216 105  57  10   1<\/code><\/pre>\n
histogram(~heads, data = GuessingLadies, width = 1)<\/code><\/pre>\n
Based on our simulation, it appears that someone can get 9 or 10 just by guessing only about 1.1% of the time.<\/p>\n
If we want to have a more precise estimate of this proportion (the p-value), we can increase the number of times we do()<\/code> 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<\/a>.<\/p>\n<\/div>\n
\n
Grand Tour of Randomization Distributions<\/h3>\n
If we want to test for a different proportion (e.g., 0.25), we simply tell rflip()<\/code> the probability of obtaining heads (p<\/code>) and the size of our sample (n = 100<\/code> in the example below):<\/p>\n
Proportion.null <- do(1000) * rflip(n = 100, p = 0.25) <\/code><\/pre>\n
##     n heads tails prop\r\n## 1 100    31    69 0.31\r\n## 2 100    25    75 0.25<\/code><\/pre>\n
<998 rows omitted><\/code><\/pre>\n
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:<\/p>\n
Galton2 <- Galton %>% group_by(family) %>% sample_n(1)  # one from each family\r\nNull1 <- do(1000) *  diffprop(  homeless ~ shuffle(sex),       data = HELPrct)\r\nNull2 <- do(1000) *  diffmean(  height ~ shuffle(sex),       data = Galton2)\r\nNull3 <- do(1000) *  lm(  age ~ shuffle(substance), data = HELPrct)  # ANOVA\r\nNull4 <- do(1000) *  lm(  father ~ shuffle(mother),    data = Galton2)  # regression<\/code><\/pre>\n
Looking at the first few rows of Null4<\/code>, we see that do()<\/code> extracts and records several useful bits of information about a linear model:<\/p>\n
##   Intercept  mother sigma r.squared     F\r\n## 1      79.8 -0.1641  2.60   0.02173 4.332\r\n## 2      74.1 -0.0736  2.62   0.00437 0.856\r\n## 3      63.9  0.0849  2.62   0.00582 1.141<\/code><\/pre>\n
We can use the slope (labeled mother<\/code>), r2<\/sup>, or F<\/em> as our test statistic.<\/p>\n
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<\/em>0<\/sub>: \u03bc = 98.6 into H<\/em>0<\/sub>: \u03bc – 98.6 = 0 and sampling with replacement from the sample:<\/p>\n
require(Lock5withR)\r\nNull <- do(1000) * mean( ~ (BodyTemp - 98.6), data = resample(BodyTemp50))<\/code><\/pre>\n
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.<\/p>\n<\/div>\n
\n
Bootstrap, too<\/h3>\n
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\u2019s data, resampling separately from the males and females:<\/p>\n
HeightBySex.boot <- \r\n  do(1000) * diffmean( height ~ sex, data = resample(Galton, groups = sex) )\r\nhistogram(~diffmean, data = HeightBySex.boot)<\/code><\/pre>\n
Simple confidence intervals can be computed using percentiles<\/p>\n
cdata(0.95, ~diffmean, data = HeightBySex.boot)  # central 95%<\/code><\/pre>\n
##       low        hi central.p \r\n##      4.78      5.46      0.95<\/code><\/pre>\n
or the bootstrap standard error, which is just the standard deviation of the bootstrap distribution.<\/p>\n
sd(~diffmean, data=HeightBySex.boot)<\/code><\/pre>\n
## [1] 0.166<\/code><\/pre>\n
Eventually, I may show my students how to automate generating the confidence intervals from the boostrap distribution.<\/p>\n
confint(HeightBySex.boot, method=c(\"percentile\", \"se\"))<\/code><\/pre>\n
##       name lower upper level   method estimate margin.of.error  df\r\n## 1 diffmean  4.78  5.46  0.95 quantile       NA              NA  NA\r\n## 2 diffmean  4.80  5.45  0.95   stderr     5.12           0.326 897<\/code><\/pre>\n<\/div>\n
\n
Striking a Balance<\/h3>\n
I find that using the mosaic<\/code> 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<\/em> objects, etc.). There are other R<\/em> packages (resample<\/code><\/a> and boot<\/code><\/a>, 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.<\/p>\n
The R<\/em>-related cognitive load for students to do SBI using mosaic<\/code> amounts to<\/p>\n