r_tutorials/w06-R_tutorial.rmd

   1 ---
   2 title: "Week 6 R Lecture"
   3 author: "Jeremy Foote and Aaron Shaw"
   4 date: "October 20, 2020"
   5 output:
   6   html_document:
   7     toc: yes
   8     toc_depth: 3
   9     toc_float:
  10       collapsed: true
  11       smooth_scroll: true
  12     theme: readable
  13   pdf_document:
  14     toc: yes
  15     toc_depth: '3'
  16   header-includes:
  17   - \newcommand{\lt}{<}
  18   - \newcommand{\gt}{>}
  19 ---
  20
  21 ```{r setup, include=FALSE}
  22 knitr::opts_chunk$set(echo = TRUE)
  23 ```
  24
  25 ## Analyzing categorical data in R
  26
  27 The goal of this script is to help you think about analyzing categorical data, including proportions, tables, chi-squared tests, and simulation.
  28
  29 ### Estimating proportions
  30
  31 If a survey of 50 randomly sampled Chicagoans found that 45% of them thought that Giordano's made the best deep dish pizza, what would be the 95% confidence interval for the true proportion of Chicagoans who prefer Giordano's?
  32
  33 Can we reject the hypothesis that 50% of Chicagoans prefer Giordano's?
  34
  35
  36 ```{r}
  37 est = .45
  38 sample_size = 50
  39 SE = sqrt(est*(1-est)/sample_size)
  40
  41 conf_int = c(est - 1.96 * SE, est + 1.96 * SE)
  42 conf_int
  43 ```
  44
  45 What if we had the same result but had sampled 500 people?
  46
  47
  48 ```{r}
  49 est = .45
  50 sample_size = 500
  51 SE = sqrt(est*(1-est)/sample_size)
  52
  53 conf_int = c(est - 1.96 * SE, est + 1.96 * SE)
  54 conf_int
  55 ```
  56
  57 ### Tabular Data
  58
  59 The Iris dataset is composed of measurements of flower dimensions. It comes packaged with R and is often used in examples. Here we make a table of how often each species in the dataset has a sepal width greater than 3.
  60
  61 ```{r}
  62 data(iris)
  63 table(iris$Species, iris$Sepal.Width > 3)
  64
  65 ```
  66
  67
  68 As described in the most recent *OpenIntro* reading, the chi-squared ($\chi^2$) test is a test of whether the frequencies we see in a table differ from what we would expect if there was no difference between the groups. Base-R provides a `chisq.test()` function to conduct these tests. Here's an example of the test evaluating whether the frequency of Iris sepal widths greater than 3 differs by species (the null hypothesis is that there is no difference in the frequencies and the groups are drawn from the same population). It's always good to examine the contingency table before conducting the test!
  69
  70 ```{r}
  71 table(iris$Species, iris$Sepal.Width > 3)
  72
  73 chisq.test(table(iris$Species, iris$Sepal.Width > 3))
  74 ```
  75
  76 While we know little about the dataset, sample, or measurement, the incredibly low p-value means that the observed frequencies likely did not come from the same distribution and that the frequency of sepal widths more/less than 3 differs by species.
  77
  78 Note that you can create tables to conduct $\chi^2$ tests by hand. The `chisq.test()` function will need the object to be of "class" `table` or `matrix` in order to run the test.
  79
  80 ## BONUS: Using simulation to test hypotheses and calculate "exact" p-values
  81
  82 Section 6.1.4 of *OpenIntro* mentions a simulation-based approach to inference when assumptions of $\chi^2$ tests aren't met. You possess all the tools to implement something like this in R, so let's do it. Basically, the idea is that we use simulated data to construct a null distribution of outcomes and then compare our observed test statistic against the null to estimate an "exact" p-value corresponding to the proportion of samples drawn from the null would produce statistics at least as large as our observed data.
  83
  84 The earlier (3rd) edition of Chapter 6 of the *OpenIntro* textbook uses an example of a medical practitioner who has 3 complications out of 62 procedures, while the typical rate is 10%. We can use that information to create an example here based on the proportions. The null hypothesis is that this practitioner's true rate is also 10%, so we're trying to figure out how rare it would be to have 3 or fewer complications, if the true rate is 10%.
  85
  86 ```{r}
  87 # We write a function that we are going to replicate
  88 simulation <- function(rate = .1, n = 62){
  89   # Draw n random numbers from a uniform distribution from 0 to 1
  90   draws = runif(n)
  91   # If rate = .4, on average, .4 of the draws will be less than .4
  92   # So, we consider those draws where the value is less than `rate` as complications
  93   complication_count = sum(draws < rate)
  94   # Then, we return the total count
  95   return(complication_count)
  96 }
  97
  98 # The replicate function runs a function many times
  99 simulated_complications <- replicate(5000, simulation())
 100
 101 ```
 102
 103 We can look at our simulated complications
 104
 105 ```{r}
 106 hist(simulated_complications)
 107 ```
 108
 109 And determine how many of them are as extreme or more extreme than the value we saw. This is the "exact" p-value.
 110 ```{r}
 111
 112 sum(simulated_complications <= 3)/length(simulated_complications)
 113 ```
 114