r_tutorials/w05-R_tutorial.rmd

   1 ---
   2 title: "Week 5 R tutorial"
   3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
   4 author: "Aaron Shaw"
   5 date: "October 13, 3030"
   6 output:
   7   pdf_document:
   8     toc: yes
   9     toc_depth: '3'
  10   html_document:
  11     toc: yes
  12     number_sections: true
  13     toc_depth: 3
  14     toc_float:
  15       collapsed: false
  16       smooth_scroll: true
  17     theme: readable
  18 ---
  19
  20 ```{r setup, include=FALSE}
  21 knitr::opts_chunk$set(echo = TRUE)
  22 ```
  23
  24 This week, we'll introduce some ways to do some the discrete math, generate distributions, and use a for-loop to run a small simulation.
  25
  26 # Binomial and factorial functions
  27
  28 In Chapter 3 (and in the previous problem set), you needed to calculate some binomial choice arithmetic and/or factorials. They weren't absolutely necessary for the problem set, but here are the corresponding functions in R.
  29
  30 Let's say we want to calculate how many possible pairs you can draw from a population of ten individuals, a.k.a., $10 \choose 2$ or, instead you wanted to calculate $10!$
  31 ```{r}
  32 choose(10,2)
  33
  34 factorial(10)
  35 ```
  36 Note that factorial arithmetic can get quite large quite fast, so consider your processor/memory constraints and options before you try to calculate something truly large like $365!$.
  37
  38 # Distribution functions
  39
  40 R has a number of built-in functions to help you work with distributions in various ways that also started to come up in *OpenIntro* Chapter 3. I will introduce a couple of points about them here, but I also highly recommend you look at the relevant section of the Verzani *Using R Introductory Statistics* book (pp 222-229) for more on this (and, honestly, for more on most of the topics we're covering in R).
  41
  42 The key to using R to analyze distributions is that R has a set of built-in distributions (e.g. uniform, normal, binomial, log-normal, etc.) and a set of functions (`d`, `p`, `q`, and `r`) that can be run for each distribution. In the example that follows, I'll use a uniform distribuition (a distribution between any two values (*min*, *max*) where the values may occur with uniform probability) for my example below. Verzani has others for when you need them.
  43
  44 The `d` function gets you information about the *density function* of the distribution. The `p` function works with the *cumulative probabilities*. The `q` function gets you *quantiles* from the distribution. The `r` function allows you to generate *random samples* from the distribution. As you can see, the letters corresponding to each function *almost* make sense...<*sigh*>. They also each take specific arguments that can vary a bit depending on which kind of distribution you are using them with (as always, the help documentation and the internet can be helpful here).
  45
  46 Onwards to the example code, which uses the different functions to calculate information about a uniform distribution between 0 and 3 (take a moment to think about what that would look like in terms of raw data and/or a plot):
  47
  48 ```{r}
  49 dunif(x=1, min=0, max=3) # What proportion of the area is the to the left of 1?
  50
  51 punif(q=1, min=0, max=3) # Same as the prior example in this case.
  52
  53 qunif(p=0.5, min=0, max=3) # 50th percentile
  54
  55 runif(n=4, min=0, max=3) # Random values in [0,3]
  56 ```
  57 Look at the Verzani text for additional examples, including several that can solve binomial probability calculations (e.g., if you flip a fair coin 100 times, what are the odds of observing heads 60 or more times?).
  58
  59 # For-loops and a quick simulation
  60
  61 Beyond proving invaluable for rapid calculations of solutions to problem set questions, the distribution functions are very, very useful for running simulations. We won't really spend a lot of time on simulations in class, but I'll give you an example here that can generalize to more complicated problems. I also use a programming technique we haven't talked about yet called a for-loop to help repeat the sampling process multiple times. For-loops feature prominently in some programming tasks/languages, but I encourage you to minimize your use of them in R for reasons that are sort of beyond the scope of the course. That said, it's still super important to learn how they work!
  62
  63 For my simulation let's say that I want to repeatedly draw random samples from a distribution and examine the distribution of the resulting sample means (this example is going to feature prominently in Chapter 5 of *OpenIntro*). I'll start by generating a vector of 10,000 random data points drawn from a log-normal distribution where the mean and standard deviation of the log-transformed values are 0 and 1 respectively:
  64
  65 ```{r}
  66 d <- rlnorm(10000, meanlog=0, sdlog=1)
  67
  68 head(d)
  69 mean(d)
  70 sd(d)
  71 hist(d)
  72 ```
  73
  74 That sure does look like a logarithmic distribution!
  75
  76 Okay, now, I want to draw 500 samples of 100 observations from this population and take the mean of each sample. Time to write a function! Notice that I require two inputs into my function: the population data and the sample size.
  77
  78 ```{r}
  79 sample.mean <- function(pop, n){
  80   s <- sample(pop, n)
  81   return(mean(s))
  82 }
  83
  84 ## Run it once to see how it goes:
  85 sample.mean(d, 100)
  86 ```
  87 Next step: let's run that 500 times. Here's where the for-loop comes in handy. In somewhat abstract terms, a for-loop allows me to define an index and then repeat an operation for some range of values that can assume a function of that index. The simplest example is to take a range of integers as my index (e.g., 1-10) and then I run my loop ten times (once for each integer). You can also do more sophisticated things such as using the index to substitute and/or subset the objects in your loop along the way.
  88
  89 We can break that down further through an example. The basic syntax of a for-loop is that you call some operation to occur over some index. Here's a simple example that illustrates how that works. The loop iterates through the integers between 1-10 and prints the square of each value:
  90 ```{r}
  91 for(x in c(1:10)){
  92   print(x^2)
  93 }
  94 ```
  95 Notice that the initial statement defines a variable (`x` in this case) within the loop (that's a way you can think about the `for(x in...))` part of the first line of code. The part that comes after `in` defines the index (in this case, the integers between 1 and 10). The part between the curly braces defines what R will do for each iteration of the loop, in which it will substitute the corresponding value from the index for variable.
  96
  97 Now let's go back to my sample means example. Since I want to store the output of my sample means loop, I will first create an object `s.means` that is a numeric vector with one value (0) that will be replaced as I run through the loop.
  98 ```{r}
  99 s.means <- 0
 100 ```
 101 Onwards to the loop itself. In the block of code below, you'll notice that I once again declare an index over which to iterate. That's what happens inside that first set of parentheses where I have `i in c(1:30)`. That's telling R to go through the loop for each value from 1:30 and to assign the current index value to `i` during that iteration of the loop. Then, each time through the loop, the value of `i` advances through  the index (in this case, it just goes up by 1). The result is that each iteration will take the output of my `sample.mean` function and append it as the $i^{th}$ value of `s.means`. The `next` call at the end is optional, but can be important sometimes to help you keep track of what's going on.
 102
 103 ```{r}
 104 for(i in c(1:500)){
 105   s.means[i] <- sample.mean(d, 100)
 106   next
 107 }
 108 ```
 109 In case you're coming to R from other programming languages that index from 0, you should note that this example very much takes advantage of the fact that R indexes from 1 (i.e., the first value of some vector `v` can be returned by `v[1]`).
 110
 111 Once we've run that loop the `s.means` variable now contains a distribution of sample means. What are the characteristics of the distribution? In other words, how would you summarize the distribution? Well, you already know how to do that.
 112
 113 ```{r}
 114 summary(s.means)
 115 ```
 116 Let's plot it too:
 117 ```{r}
 118 library(ggplot2)
 119 qplot(s.means, geom="density")
 120 ```
 121
 122 That looks pretty "normal."
 123
 124 Experiment with this example by changing the size of the sample and/or the number of samples we draw.
 125
 126 Now, think back to the original vector `d.` Can you explain what fundamental statistical principle is illustrated in this example? Why do the values in `s.means` fluctuate so much? What is the relationship of `s.means` to `d`?
 127
 128 # Quantile quantile plots
 129
 130 Last, but not least, you might have admired the quantile-quantile plots (a.k.a. "Q-Q plots") presented in some of the examples in the most recent *OpenIntro* chapter. The usual idea with Q-Q- plots is that you want to see how the observed (empirical) quantiles of some data compare against the theoretical quantiles of a normal (or other) distribution. You too can create these plots!
 131
 132 Here's an example that visualizes the result of our simulation (labeled "sample") against a normal distribution with the same mean and standard deviation (labeled "theoretical"). Notice that to accommodate `ggplot2` I have to turn `s.means` into a data frame first.
 133
 134 ```{r}
 135 s.means <- data.frame(s.means)
 136 ggplot(s.means, aes(sample=s.means)) + geom_qq() + geom_qq_line(color="red")
 137
 138 ```
 139
 140 And/or (finally) we could even standardize the values of `s.means` as z-scores using the `scale()` function:
 141
 142 ```{r}
 143 s.z <- data.frame(scale(s.means)); names(s.z) <- "z"
 144 ggplot(s.z, aes(sample=z)) + geom_qq() + geom_qq_line(color="red")
 145 ```
 146