]> code.communitydata.science - stats_class_2020.git/blob - r_tutorials/w04-R_tutorial.rmd
typo
[stats_class_2020.git] / r_tutorials / w04-R_tutorial.rmd
1 ---
2 title: "Week 4 R tutorial"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw"
5 date: "October 6, 2020"
6 output:
7   html_document:
8     toc: yes
9     number_sections: true
10     toc_depth: 3
11     toc_float:
12       collapsed: false
13       smooth_scroll: true
14     theme: readable
15   pdf_document:
16     toc: yes
17     toc_depth: '3'
18 ---
19
20 ```{r setup, include=FALSE}
21 knitr::opts_chunk$set(echo = TRUE)
22 ```
23
24 Welcome back! This week we'll get into some more advanced fundamentals to help you import and tidy/manage data, define and run your own functions, generate distributions and samples, and manage dates. Some of these topics are more useful for future problem sets (e.g., the stuff about distributions and dates), but I've included them here to start introducing the ideas.
25
26
27 # Importing (yet more) data
28
29 So far, we have imported datasets that come installed with R packages with the `data()` function as well as the `load()` function to read R data files (.RData, .rda, etc.). For better and worse, you should be able to import data in other formats too.
30
31 Before I get any further, I want to note that my approach here is to use R commands directly that you can type in RMarkdown scripts or run at the console yourself. RStudio also provides a number of handy data import tools through the graphical interface and drop-down menus. [This how-to article](https://support.rstudio.com/hc/en-us/articles/218611977-Importing-Data-with-RStudio) introduces some of these resources.
32
33 Tabular (rows and columns) data files formatted as plain text with "comma-separated values" (".csv's") are quite common, so we'll look at those. R comes with a handy `read.csv()` command that does exactly what you'd expect. Here's an example using a csv file I created from one of R's built-in datasets called `mtcars`, which has old data about fuel consumption cars. Run `help(mtcars)` to learn more about where it comes from and to read the variable descriptions. Since it's built-in, you can import it using `data(mtcars)`, but I also posted it to the course data repository so we can use the `url()` command to point `read.csv()` to download it:
34
35 ```{r}
36
37 data.url <- url("https://communitydata.science/~ads/teaching/2020/stats/data/week_04/mtcars.csv")
38
39 my.mtcars <- read.csv(data.url)
40
41 head(my.mtcars)
42 ```
43
44 As always, take a look at the help documentation for `read.csv()` to learn more about some of the arguments that you can use. Because data comes in so many (weird) formats, there are many possible arguments!
45
46 You might notice that the documentation for `read.csv()` is actually part of the documentation for another command called `read.delim()`. Turns out `read.delim()` is just a more general-purpose way to read in tabular data and that `read.csv()` is short-hand for `read.delim()` with some default values that make sense for csv files. Here is a command that produces identical output to the previous one
47
48 ```{r}
49
50 more.cars <- read.delim(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/mtcars.csv"), sep=",")
51
52 head(more.cars)
53
54 table(more.cars == my.mtcars)
55 ```
56
57 When find yourself trying to load a tabular data file that consists of plain text, but has some idiosyncratic difference from a csv (e.g., it is tab-separated instead of comma-separated), you may want to use `read.delim()`.
58
59 Note that the Tidyverse has some related functions, namely `read_csv()` that can also import csv's very efficiently and with helpful defaults to try and guess what kinds of variables you're working with. The guesses are usually pretty good! Here's an example using `read_csv()`
60
61 ```{r}
62 library(tidyverse)
63
64 yet.more.cars <- read_csv(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/mtcars.csv"))
65
66 yet.more.cars
67
68 table(yet.more.cars == my.mtcars)
69 ```
70
71 How will you know what to use? The best practice is always to get to know your data first! Seriously, try opening it the file (or at least opening up part of it) using a text editor and/or spreadsheet software. Looking at the "raw" plain text can help you figure out what arguments you need to use to make the data load up exactly the way you want it.
72
73 For example, you might notice that my import of the mtcars.csv file introduces an important difference from the original `mtcars` dataset. In the original `mtcars`, the car model names are `row.names` attributes of the dataframe instead of a variable. 
74
75 ```{r}
76 head(mtcars)
77 row.names(mtcars)
78 ```
79
80 Since it's the first column of the raw data, you can fix this with an additional argument to `read.csv` (I'm sure there's also a way to do this with `read_csv()` but I didn't look it up in time to include in this tutorial):
81
82 ```{r}
83
84 my.mtcars <- read.csv(url("https://communitydata.science/~ads/teaching/2020/stats/data/week_04/mtcars.csv"), 
85                       row.names=1)
86 head(my.mtcars)
87
88 table(my.mtcars == mtcars)
89 ```
90
91 This illustrates a common issue that relates back to variable types (classes). Most of the commands in R that import data try to "guess" what class is appropriate for each column of your dataset. Surprise, surprise, these guesses are sometimes not so great and often quite different from what you might guess. As a result, it's a great idea to inspect the classes of every column of a dataset after you import it (review last week's R lecture materials for more on this). 
92
93 ## Importing proprietary data formats
94
95 R has libraries that can read (and write) many proprietary data file formats, including files from Stata, SAS, MS Excel, and SPSS (among others). 
96
97 That same [helpful Rstudio how-to data import article](https://support.rstudio.com/hc/en-us/articles/218611977-Importing-Data-with-RStudio) includes several examples of these. You can find other suggestions and examples online.
98
99
100 # Tidy and manage data  
101
102 Last tutorial introduced the "*apply" family of functions (e.g., `sapply`, `lapply`, and `tapply`). This time around, I want to introduce an example of how to use them repeatedly over a set of objects. 
103
104 ## Conditional means with nested `*apply` functions
105
106 I used `mtcars` for this last time around, so let's start off with it again. Before I used `tapply()` to calculate conditional means of miles per gallon for each number of cylinders. What if I wanted to do calculate conditional means for a bunch of the of the other variables in the `mtcars` dataset? You *could* just write the same line of `tapply` code over and over again until you had one line for each variable. However, there's a much more efficient way:
107
108 ```{r}
109
110 variables <- c("mpg", "disp", "hp", "wt")
111
112 sapply(mtcars[variables], function(v){
113   tapply(v, mtcars$cyl, mean) 
114   } 
115   )
116
117 ```
118
119 Do you follow the example? Try reading it "inside-out" to understand out what it does. The innermost part is the call to `tapply()` and that should look pretty familiar from the previous tutorial example. The difference now is that it contains this weird call to some variable `v`. Working backwards, the previous line of code defines `v` as an argument to a function. The function is, in turn, being passed as an argument to `sapply()`, which it is applying to the columns in `mtcars` indexed by my vector of variable names.
120
121 Put another way: I have used `sapply` to call `tapply` to calculate conditional means (by cylinder number) for each item in the `variables` vector. My example only included four variables, but you could do the same thing for an arbitrarily large set of variables and it would work just fine. Indeed, it would be a very fast, efficient way of solving a seemingly complicated data processing task with very little code.
122
123 ## Conditional means in Tidyverse code 
124 I should reiterate here that the `*apply` functions are part of Base R. Other functions to do similar things exist and you may find these other functions more intelligible or useful for you. In particular, this is the sort of task that the Tidyverse handles in an arguably more intuitive fashion than Base R because it allows you to organize functions as a sequence of actions with function names that are verbs and this generally leads to more readable code. Here's an example snippet that replicates the same output. :
125
126 ```{r}
127 library(tidyverse)
128
129 mtcars %>% 
130   group_by(cyl) %>%
131   summarize( across(
132     .cols=c(mpg, disp, hp, wt),
133     .fns=list(
134       avg=mean) ## writing this in this way allows the summarize function to create a useful variable name
135     ))
136
137 ```
138
139 A little bit more verbose, but still quite concise![^1] Let's unpack it. The first thing you'll notice is the `%>%` symbol on the first line. In Tidyverse-speak `%>%` is referred to as a "pipe." It tells R to take the thing to the left and use it to perform the action that follows on the right (or below in this example). The other thing to notice is that when I'm working inside of these "piped" commands, the Tidyverse doesn't require me to use quotations to indicate variable/column names.
140
141 [^1]: Indeed, this command could be yet more concise, but I've intentionally left it a little "wordy" to facilitate readability and allow the example to develop further.
142
143 Once you know those Tidyverse fundamentals, the rest is nearly readable: you take the `mtcars` data set, group it by the `cyl` variable, and then summarize across a series of columns (`.cols`) with some function(s) (`.fns`). 
144
145 Since the Tidyverse code takes a grammatical structure that nearly maps to English, adding another action (like calculating an additional summary statistic such as the conditional variance) can be done in a fairly straightforward way:
146
147 ```{r}
148
149 mtcars %>% 
150   group_by(cyl) %>%
151   summarize( across(
152     .cols=c(mpg, disp, hp, wt),
153     .fns=list(avg = mean,
154               sigma.sq = var)
155     ))
156
157 ```
158
159
160 Don't worry if you didn't follow that yet. The [Tidyverse documentation](https://tidyverse.org/) is extensive. Indeed, my example looks a lot like some of the examples from [the documentation for the `summarize()` function](https://dplyr.tidyverse.org/reference/summarise.html) and [the `across()` function](https://dplyr.tidyverse.org/reference/across.html) (Hadley Wickham, the author of most of the Tidyverse uses `mtcars` a **lot** in his example code). For now, I mostly wanted to introduce a Tidyverse data cleanup/management example because it's an area where the Tidyverse really shines.
161
162 ## Merging data frames
163
164 Combining various chunks of data together poses another family of common data management and processing challenges.
165
166 One way to do this when you have multiple dataframes is to use the `merge()` function. A merge might involve dataframes with mutually exclusive rows/columns or it might involve dataframes that already have some things in common. I'll present an example here using two example dataframes that I create by hand and that share a column:
167
168 ```{r}
169 spiders <- c("black widow", "brown recluse", "daddy longlegs")
170 poisonous <- c(TRUE, TRUE, FALSE)
171 region <- c("south", "central", "all")
172 genus <- c("Latrodectus", "Locosceles", "many")
173 number <- c(2, 1, 250)
174
175
176 # Notice that the datasets here share one column in common:
177 s1 <- data.frame(spiders, poisonous, region)
178 s2 <- data.frame(spiders, genus, number)
179
180 s.combined <- merge(s1, s2, by="spiders")
181 s.combined
182
183 ```
184 The `by="spiders"` part of the call to `merge()` tells R to use the shared column to produce a combined dataframe. 
185
186 Take a look at the documentation for `merge()` and search for additional examples on StackOverflow and beyond if/when you want to pursue more complicated tasks along these lines. The arguments you can pass via `by` and `all`  can be especially useful.
187
188 ## Cleaning and organizing data
189
190 You already saw several ways to cleanup data in my initial presentation of dataframes last week. Let's return to that and build out some additional examples here. A lot of what follows is about taking the results of various data management steps and making sure they are arranged so as to make subsequent steps more intelligible. While it may seem tedious to spend time writing or revising your code to do this kind of thing, the payoff is code and data that are more human-readable. When your code and data are more human-readable, it's easier to develop intuitions and find mistakes. 
191
192 Take a look at the class of the object returned by a previous example:
193
194 ```{r}
195 class(
196   sapply(mtcars[variables], function(v){
197     tapply(v, mtcars$cyl, mean)
198 }) )
199
200 ```
201
202 Hmm, what if you wanted a data frame? You already know how to do that! I'll call `data.frame()` and assign the output to a new object:
203
204 ```{r}
205 cyl.conditional.means <- data.frame(
206   lapply(mtcars[variables], function(v){
207     tapply(v, mtcars$cyl, mean)
208 }) )
209
210 head(cyl.conditional.means)
211 ```
212
213 We can even do some more cleanup to give the cylinder count its own variable and rename the other columns to reflect that they consist of averages:
214
215 ```{r}
216 names(cyl.conditional.means) <- c("mean.mpg", "mean.disp", "mean.hp", "mean.wt")
217
218 # Notice that I'm converting the values back to numeric
219 cyl.conditional.means$cyl.number <- as.numeric(row.names(cyl.conditional.means)) 
220
221 # Maybe look at the row names before you run this next line so that you can see what it does.
222 row.names(cyl.conditional.means) <- NULL 
223
224 head(cyl.conditional.means)
225
226 ```
227
228 This data should be clearer and easier to work with now.
229
230 What if, instead of having that dataframe sorted by the number of cylinders, you wanted to sort it by the value of `mean.mpg`? The `sort.list` function can help:
231
232 ```{r}
233 cyl.conditional.means[sort.list(cyl.conditional.means$mean.mpg),]
234 ```
235
236 Again, you could run `row.names(cyl.conditional.means)` <- NULL` to reset the row numbers.
237
238 By now you might have noticed that a some of my code in this section replicates parts of the output that the Tidyverse example above created by default. As I said earlier, the Tidyverse really shines when it comes to tidying data (shocking, I know). There are also functions ("verbs") in the Tidyverse libraries to perform all of the additional steps (such as sorting data by the values of a column). If you are interested in doing so, check out more of the [Tidyverse documentation](https://tidyverse.org/) and the [Wickham and Grolemund *R for Data Science* book](https://r4ds.had.co.nz/) listed among the course resources. 
239
240 ## Working with dates (in Base R)
241
242 Date and time objects are another more advanced data class in R. Managing date and time data can be a headache. This is because dates and times tend to be formatted inconsistently and are usually given to you as character variables, so you should know how to transform them into a format that R can "understand" as dates. There are many packages for working with dates and times, but for now I'll introduce you to the Base R way of doing so. This uses a data type formally called "POSIX" (no need to worry about why it's called that).
243
244 To build up an example, I'll create some date-time values, add a little noise, and convert them into a character vector:
245
246 ```{r}
247 add.an.hour <- seq(0, 3600*24, by=3600)
248 some.hours <- as.character(Sys.time() + add.an.hour) ## Look up Sys.time() to see what it does.
249 ```
250 Now, let's take a look at the results of that:
251 ```{r}
252 class(some.hours)
253 head(some.hours)
254 ```
255 These are beautifully formatted timestamps, but R will not understand them as such. This is often how you might receive data in, for example, a dataset you import from Qualtrics, scrape from the web, or elsehwere. You can convert the `some.hours` vector into an object class that R will recognize as a time object using the `as.POSIXct()` function. Notice that it even adds a timezone back in!
256 ```{r}
257 as.POSIXct(some.hours)
258 ```
259
260 If things aren't formatted in quite the way R expects, you can also tell it how to parse a character string as a POSIXct object:
261
262 ```{r}
263 sometime <- "2019-02-21 04:35:00"
264 class(sometime)
265
266 a.good.time <- as.POSIXct(sometime, format="%Y-%m-%d %H:%M:%S", tz="CDT")
267 class(a.good.time)
268 ```
269 Once you have a time object, you can even do date arithmetic with `difftime()` (but watch out as this can get complicated):
270
271 ```{r}
272 difftime(Sys.time(), a.good.time, units="weeks")
273 ```
274 This calculated the number of weeks elapsed between the current time and an example date/time I created above.
275
276 ### `as.Date` vs. `as.POSIX()` 
277
278 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. I'll use the same 
279
280 ```{r}
281 a.good.time <- as.Date(sometime, format="%Y-%m-%d %H:%M:%S", tz="CDT")
282 class(a.good.time)
283 a.good.time
284 ```
285
286 ### An easier way (most of the time) in the Tidyverse
287
288 The Tidyverse (via the `lubridate` package) usually handles dates and times quite well. When you need to import and work with date and time objects, you may find it best to try Tidyverse data import tools (e.g., `read_csv()`) as a starting point for this reason. 
289
290 I *highly recommend* reading [this chapter](https://r4ds.had.co.nz/dates-and-times.html) of Wickham and Grolemund on dates and times as they introduce a number of challenges and nuances that can befuddle even the most brilliant statistical programmers.
291
292 # Creating (repeated) sequences, distributions, and samples
293
294 The *OpenIntro* reading this week introduces some core concepts and tools for probability. R has a number of functions that you can use to illustrate these concepts, replicate examples, and generate simulations of your own. This section of the tutorial introduces a few useful starting points. You may also find these snippets and functions valuable for managing real data.
295
296 First, let's create two sequences of numbers using the `seq()` function:
297
298 ```{r}
299 odds <- seq(from=1, to=100, by=2)
300 evens <- seq(from=2, to=100, by=2)
301
302 head(odds)
303 head(evens)
304 ```
305 Or maybe (for some reason) I want to have a vector that repeats the odd integers between 1 and 100 five times? Try `rep()`:
306
307 ```{r}
308 more.odds <- rep(seq(from=1, to=100, by=2), 5)
309 ```
310 You can use `sample()` to draw samples from an object:
311
312 ```{r}
313 sample(x=odds, size=3)
314 sample(x=evens, size=3)
315 ```
316
317 You can also sample "with replacement." Here I take 100 random draws from the binomial distribution, which is analogous to 100 independent trials of, say, a coin flip:
318 ```{r}
319 draws <- sample(x=c(0,1), size=100, replace=TRUE)
320
321 table(draws)
322 ```
323
324 What if you wanted to take a random set of 10 observations from a dataframe? You can use `sample` on an index of the rows:
325
326 ```{r}
327 odds.n.evens <- data.frame(odds, evens)
328
329 odds.n.evens[ sample(row.names(odds.n.evens), 10), ]
330
331 ```
332
333 ## Managing randomness
334
335 Try running one of the `sample` commands above again. You will (probably!) get a different result because `sample` makes a random draw each time it runs.
336
337 Even in statistics, you sometimes want things to be a little *less random*. For example, let's say you want to reproduce a specific random draw. To do this, you use the `set.seed()` command, which takes any integer as its argument. Setting the seed value in this way ensures that the next random draw will be the same each time:
338
339 ```{r}
340 set.seed(42) 
341 x <- sample(odds, 10)
342
343 y <- sample(odds, 10) # different
344
345 set.seed(42)
346 z <- sample(odds, 10) # 
347
348 head(x)
349 head(y)
350 head(z)
351
352 table(x==y) ##  One in ten!
353 table(x==z)
354 ```
355 Note that the first random draw after I set the seed to the same value is the same both times (see the results of `table(x==z)`), but not the same when I create `y`.
356
357 # Creating functions
358
359 You should learn how to define and run your own functions. This is a big part of what makes statistical **programming** a thing in R and what has made it possible for people to develop such a wide range of packages on top of the R Base. 
360
361 Defining a function in R requires that you use the `function()` command to let R know what your function can do as well as what arguments you will pass to it. Here is an example of a function that I'll call `times.two()`. It takes a number and multiplies it by two:
362
363 ```{r}
364 times.two <- function(x){
365   twicex <- x * 2
366   return(twicex)
367 }
368
369 times.two(21)
370 ``` 
371
372 A few things to notice here. First, you can see that I assigned my function to an object/variable (in this case `times.two`) just like we've done with other kinds of objects. (**Important**: You need to execute the code that defines the function in order to be able to use it!)
373
374 I use curly braces (`{ }`) to contain the content of the function. Everything inside the curly braces is the "action" that the function executes when I call it later (in this case, the function takes an arbitrary value, multiplies it by two, and returns the result). 
375
376 The `return()` syntax at the end of the function tells R exactly what output I want my function to "return" when I execute it. As usual in R, you can leave out the explicit call to something like `return()` and R will try to guess what you want based on the last thing you asked it to do inside the function. That guess will usually be bad, so you're almost always better off calling `return()`.
377
378 Any function needs something inside the parentheses, known as an "argument." In this case the argument, `x` is included in the definition and content of my function. The important thing to know is that this `x` will really only have meaning *inside* the context of the function. If another `x` is defined elsewhere in my R environment at the time I run my function, R won't necessarily "see" it within the function.
379
380 Finally, in the last line of that code block, you can see that I have "called" my function with the argument "21." Since I execute the part of my code that defines the function, R now knows what `two.times()` means, so when I call it and provide the argument `21`, R goes ahead and tries to execute the code I wrote substituting the argument for `x`. 
381
382 ## Example: Creating my own `mean()` function:
383 Here is a slightly more complicated function that takes a numeric vector, removes any `NA` (missing) values, and calculates the arithmetic mean:
384
385 ```{r}
386 my.mean <- function(z){
387   z <- z[!is.na(z)] 
388   sigma <- sum(z)
389   n <- length(z)
390   out.value <- sigma/n
391   return(out.value)
392 }
393 ```
394 You can try running this on the built-in `rivers` dataset (run `help(rivers)` to see what it's about). Here I'll add a missing value and then I'll compare the output of `my.mean()` to the output of the `mean()` function included with R:
395
396 ```{r}
397 data(rivers)
398 rivers[7] <- NA
399
400 my.mean(rivers)
401
402 mean(rivers, na.rm=TRUE)
403 ```
404
405 A key reminder: in order to access a user-defined function you must execute the code that defines the function first. If you don't do this, R won't know what you're talking about when you ask it to run the new function.
406
407 ## Two tips for writing and debugging your functions
408
409 When you're creating a function of you're own it's common to have things not-quite-work at some point along the way (especially while you're just learning the syntax of the R language!). In other words, your code will have bugs. What to do?!
410
411 Debugging is a big topic and a reality of programming in any language (just ask someone you know with previous programming experience!). Without going into much depth, I offer two very simple tips for avoiding, identifying, and fixing bugs in your R code:
412
413 (@) **Work small.** 
414
415 You may have big plans for your function, but it's often helpful to build towards those big plans in small chunks. For example, notice that the `my.mean()` function I defined above only does one thing on each line. This is intentional. It is possible to write an elaborate single line of code to do all of the things in that function. However, doing so increases the likelihood that there will be an error somewhere and that the error will be much, much harder to see. Shorter, simpler lines of code (especially while you're learning) are often better.
416
417 (@) **Test often.**
418
419 If you have a function that involves multiple steps, any one of them may go wrong. A great way to evaluate/prevent this is to test each step as you go to make sure it's producing the output you expect. One way to do this is by building up your function adding one piece at a time and inspecting the output. For example, you might create some test data and run each line of your code individually on the test data to make sure it does what you expect. Here's an example that test the first line of `my.mean()`:
420
421 ```{r}
422 y <- c(2,12,31,26,NA,25)
423
424 y
425
426 y[!is.na(y)]
427 ```
428 This line-by-line approach is good, but eventually you need to also make sure things are working *inside* the function. A basic way to do this is using the `print()` command to print the results of intermediate steps in your function. 
429
430 Here is an example of that using the `y` vector I just defined and a new version of the `my.mean` function with a debugging (`print()`) line included:
431
432 ```{r}
433 debug.my.mean <- function(z){
434   z <- z[!is.na(z)] 
435   sigma <- sum(z)
436   print(sigma) ## here's my debugging line.
437   n <- length(z)
438   out.value <- sigma/n
439   return(out.value)
440
441
442 debug.my.mean(y)
443
444 sum(y, na.rm = TRUE) ## Checks my code against a built-in function
445
446 ```
447
448 Notice that the call to `debug.my.mean()` produces two pieces of output: the thing I asked to be printed in the middle of the function as well as the returned value of the function itself. 
449

Community Data Science Collective || Want to submit a patch?