r_lectures/w04-R_lecture.Rmd

   1 ---
   2 title: "Week 4 R lecture"
   3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
   4 author: "Aaron Shaw"
   5 date: "April 18, 2019"
   6 output: html_document
   7 ---
   8
   9 ```{r setup, include=FALSE}
  10 knitr::opts_chunk$set(echo = TRUE)
  11 ```
  12
  13 This week, we'll focus on one more way to manage date-time objects and some ways to generate distributions.
  14
  15 ## as.Date
  16
  17 First, something I meant to include in last week's materials. The `as.Date()` function provides an alternative to `as.POSIX()` that is far more memorable and readable, but far less precise. Note that it truncates the time of day and the timezone from the ouput
  18
  19 ```{r}
  20 m <- "2019-02-21 04:35:00"
  21 class(m)
  22
  23 a.good.time <- as.Date(m, format="%Y-%m-%d %H:%M:%S", tz="CDT")
  24 class(a.good.time)
  25 a.good.time
  26 ```
  27
  28 ## Binomial and factorial functions
  29
  30 In Chapter 3 (and in last week's 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.
  31
  32 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!$
  33 ```{r}
  34 choose(10,2)
  35
  36 factorial(10)
  37 ```
  38
  39 ## Distribution functions
  40
  41 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).
  42
  43 The key to this is that R has a set of 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. 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.
  44
  45 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 pages and the internet can be helpful here).
  46
  47 Onwards to the example code, which looks at a uniform distribution between 0 and 3:
  48
  49 ```{r}
  50 dunif(x=1, min=0, max=3) # What proportion of the area is the to the left of 1?
  51
  52 punif(q=1, min=0, max=3) # Same as the prior example in this case.
  53
  54 qunif(p=0.5, min=0, max=3) # 50th percentile
  55
  56 runif(n=4, min=0, max=3) # Random values in [0,3]
  57 ```
  58 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?).
  59
  60 ### A quick simulation (using a for-loop!)
  61
  62 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.
  63
  64 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. 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:
  65
  66 ```{r}
  67 d <- rlnorm(10000, meanlog=0, sdlog=1)
  68
  69 head(d)
  70 mean(d)
  71 sd(d)
  72 hist(d)
  73 ```
  74
  75 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.
  76
  77 ```{r}
  78 sample.mean <- function(pop, n){
  79   s <- sample(pop, n)
  80   return(mean(s))
  81 }
  82
  83 ## Run it once to see how it goes:
  84 sample.mean(d, 100)
  85 ```
  86 Next step: let's run that 500 times. Here's where the for-loop comes in handy. A couple of things about the syntax of for-loops in R:  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 they work. The loop iterates through the integers between 1-10 and prints the square of each value:
  87 ```{r}
  88 for(x in c(1:10)){
  89   print(x^2)
  90 }
  91 ```
  92
  93 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 in a moment.
  94 ```{r}
  95 s.means <- 0
  96 ```
  97 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 call the current index value `i` during each loop. Each time through the loop, the value of `i` advances through  the index (in this case, it goes up by 1). The result is that each time through it 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.
  98
  99 ```{r}
 100 for(i in c(1:500)){
 101   s.means[i] <- sample.mean(d, 100)
 102   next
 103 }
 104 ```
 105 The `s.means` variable now contains a distribution of sample means! What are the characteristics of the distribution? You already know how to do that.
 106
 107 ```{r}
 108 summary(s.means)
 109 ```
 110 Let's plot it too:
 111 ```{r}
 112 library(ggplot2)
 113 qplot(s.means, geom="density")
 114 ```
 115
 116 That looks pretty "normal."
 117
 118 Experiment with this example by changing the size of the sample and/or the number of samples we draw.
 119
 120 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`?
 121
 122 ## Quantile quantile plots
 123
 124 Last, but not least, you might have admired the quantile-quantile plots presented in some of the examples in *OpenIntro*. 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 distribution. You too can create these plots!
 125
 126 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.
 127
 128 ```{r}
 129 s.means <- data.frame(s.means)
 130 ggplot(s.means, aes(sample=s.means)) + geom_qq() + geom_qq_line(color="red")
 131
 132 ```
 133
 134
 135 And/or (finally) we could even standardize the values of `s.means` as z-scores using the `scale()` function:
 136
 137 ```{r}
 138 s.z <- data.frame(scale(s.means)); names(s.z) <- "z"
 139 ggplot(s.z, aes(sample=z)) + geom_qq() + geom_qq_line(color="red")
 140 ```
 141