initial commit of this file. r lecture materials.
[stats_class_2019.git] / r_lectures / w03-R_lecture.Rmd
1 ---
2 title: "Week 3 R lecture"
3 subtitle: "Statistics and statistical programming  \nNorthwestern University  \nMTS 525"
4 author: "Aaron Shaw"
5 date: "April 11, 2019"
6 output: html_document
7 ---
8
9 ```{r setup, include=FALSE}
10 knitr::opts_chunk$set(echo = TRUE)
11 ```
12
13 Welcome back! This week we'll get into some more advanced topics to help you import and manage data, define and run your own functions, generate and sample distributions, and manage dates. Some of these topics are more useful for next week's problem set, but I've included them here to start introducing the ideas.
14
15
16 ## Importing data
17
18 So far, you have used datasets that come installed with R as well as the `load()` function to read `.RData` files. For better and for worse, you also need to be able to import data in other formats and from the web!
19
20 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 today. R comes with a handy `read.csv()` command that does exactly what you'd expect. Here's an example using a sample csv file I created from that same `mtcars` dataset we looked at last week. I uploaded it to the internet so I also use the `url()` command to get it:
21
22 ```{r}
23
24 data.url <- url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/mtcars.csv")
25
26 my.mtcars <- read.csv(data.url)
27
28 head(my.mtcars)
29 ```
30
31 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!
32
33 Also, 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 does produces identical output to the previous one
34
35 ```{r}
36
37 more.cars <- read.delim(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/mtcars.csv"), sep=",")
38
39 head(more.cars)
40
41 table(more.cars == my.mtcars)
42 ```
43
44 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 should use `read.delim()`.
45
46 How will you know what to use? 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.
47
48 For example, you might have noticed that my import of the mtcars.csv file had one difference from the original `mtcars` dataset. In the original `mtcars`, the car model names are "row.name" attributes of the dataframe instead of a whole separate variable. Since it's the first column of the raw data, you can fix this with an additional argument to `read.csv`:
49
50 ```{r}
51
52 my.mtcars <- read.csv(url("https://communitydata.cc/~ads/teaching/2019/stats/data/week_03/mtcars.csv"), 
53                       row.names=1)
54 head(my.mtcars)
55
56 table(my.mtcars == mtcars)
57 ```
58
59 Another common issue 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. 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). 
60
61 Finally (for now), you might also find it useful to know that R has a library called `foreign` that can read and write many proprietary data file formats, including files from Stata, SAS, and SPSS (among others). You can read the [documentation for the package](https://cran.r-project.org/web/packages/foreign/foreign.pdf) to learn more.
62
63 ## Defining functions
64
65 You also need to learn how to define and run your own functions. Doing this is a big part of what makes it possible to do statistical **programming** with R and what has made it possible for people to develop such a wide range of packages on top of the R "base." 
66
67 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 `two.times()`. It takes a number and multiplies it by two:
68
69 ```{r}
70 two.times <- function(x){
71   twicex <- x * 2
72   return(twicex)
73 }
74
75 two.times(21)
76 ``` 
77
78 A few things to notice here. First, you can see that I assigned my function to an object/variable (in this case `two.times`) just like we've done with other kinds of objects. Second is the use of the 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). Third is the use of the `return()` syntax at the end of the function to tell R exactly what output I want my function to "return" when I execute it. Fourth, there is the argument, `x` that I have 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 environment at the time I run my function, R won't "see"" it within the function. Finally, in the last line of that code block, you can see that I have "called" my function with the argument "21." Since R now knows what `two.times()` means, it goes ahead and executes the code I wrote substituting the argument for `x`. 
79
80 Here is a slightly more complicated function that takes a numeric vector, removes any `NA` (missing) values, and calculates the arithmetic mean:
81
82 ```{r}
83
84 my.mean <- function(z){
85   z <- z[!is.na(z)] 
86   sigma <- sum(z)
87   n <- length(z)
88   out.value <- sigma/n
89   return(out.value)
90 }
91
92
93 ```
94 You can try running this on the `rivers` dataset that we have used in the past. Here I'll add a missing value first just to test that aspect of the function and I'll compare my output against the `mean()` function included with R:
95
96 ```{r}
97 my.rivers <- rivers
98 my.rivers[7] <- NA
99
100 my.mean(my.rivers)
101
102 mean(my.rivers, na.rm=TRUE)
103 ```
104
105 ### Two tips for writing and debugging your functions
106
107 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?!
108
109 Debugging is a big topic and a reality of programming in any language (just ask the people in our class 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:
110
111 (@) **Work small.** 
112
113 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.
114
115 (@) **Test often.**
116
117 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()`:
118
119 ```{r}
120 y <- c(2,12,31,26,NA,25)
121
122 y
123
124 y[!is.na(y)]
125 ```
126 This line-by-line approach is good, but eventually you need to also make sure things are working *inside* the function. A great way to do this is using the `print()` command to print the results of intermediate steps in your function. 
127
128 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 line included:
129
130 ```{r}
131 debug.my.mean <- function(z){
132   z <- z[!is.na(z)] 
133   sigma <- sum(z)
134   print(sigma) ## here's the debugging line.
135   n <- length(z)
136   out.value <- sigma/n
137   return(out.value)
138
139
140 debug.my.mean(y)
141
142
143 sum(y, na.rm = TRUE) ## Checks my code against a built-in function
144
145 ```
146
147 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. 
148
149 ## Still more data processing and management
150
151 Last week's R lecture 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. 
152
153 ### Using "*apply" to call functions repeatedly
154
155 We 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 better way:
156
157 ```{r}
158
159 variables <- c("mpg", "disp", "hp", "wt")
160
161 sapply(mtcars[variables], function(v){
162   tapply(v, mtcars$cyl, mean) 
163   } 
164   )
165
166 ```
167
168 Do you follow the example? Try reading it "inside-out" to understand out what it does. The innermost part is the call to `tapply()`. That should look pretty familiar from the Week 2 R lecture example, except that it now 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.
169
170 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.
171
172 I should point out that using the `*apply` functions are part of the base R packages. Many other functions to do similar things exist and you may find these other functions more intelligible or useful for you. For example, if you like the `ggplot2` approach to the world might want to take a look at the other [Tidyverse](https://www.tidyverse.org/) packages, many of which are widely adopted at this point and *very* well-documented/supported by the same people who bring you RStudio. Tidyverse functions are also used in example code throughout the Kieran Healy *Data Visualization* book. Others of you who plan to work with extremely large datasets and/or have a background working with databases might find the `DataTable` package more powerful and/or intuitive.
173
174 Whatever the case, the base R `*apply` functions are widely used as well. Learning to work with them to call other functions will make you a much more effective R programmer.
175
176 ### Merging data frames
177
178 Combining various chunks of data together poses another family of common data management and processing challenges.
179
180 One common 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 a pretty straightforward example here using two example dataframes that I create by hand:
181
182 ```{r}
183 spiders <- c("black widow", "brown recluse", "daddy longlegs")
184 poisonous <- c(TRUE, TRUE, FALSE)
185 region <- c("south", "central", "all")
186 genus <- c("Latrodectus", "Locosceles", "many")
187 number <- c(2, 1, 250)
188
189
190 # Notice that they share one column in common:
191 s1 <- data.frame(spiders, poisonous, region)
192 s2 <- data.frame(spiders, genus, number)
193
194 s.combined <- merge(s1, s2, by="spiders")
195 s.combined
196
197 ```
198
199 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.
200
201 ### Cleaning up and organizing data
202
203 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. 
204
205 Take a look at the class of the object returned by a previous example:
206
207 ```{r}
208 class(
209   sapply(mtcars[variables], function(v){
210     tapply(v, mtcars$cyl, mean)
211 }) )
212
213 ```
214
215 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:
216
217 ```{r}
218 cyl.conditional.means <- data.frame(
219   lapply(mtcars[variables], function(v){
220     tapply(v, mtcars$cyl, mean)
221 }) )
222
223 head(cyl.conditional.means)
224 ```
225
226 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:
227
228 ```{r}
229 names(cyl.conditional.means) <- c("mean.mpg", "mean.disp", "mean.hp", "mean.wt")
230
231 # Notice that I'm converting the values back to numeric
232 cyl.conditional.means$cyl.number <- as.numeric(row.names(cyl.conditional.means)) 
233
234 # Maybe look at the row names before you run this next line so that you can see what it does.
235 row.names(cyl.conditional.means) <- NULL 
236
237 head(cyl.conditional.means)
238
239 ```
240
241 This data should be clearer and easier to work with now.
242
243 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:
244
245 ```{r}
246 cyl.conditional.means[sort.list(cyl.conditional.means$mean.mpg),]
247 ```
248
249 Again, you could run `row.names(cyl.conditional.means)` <- NULL` to reset the row numbers.
250
251
252 ## Creating distributions and sampling
253
254 The *OpenIntro* reading this week includes a lot of material about distributions, so I thought it would be good to also give you some R code that you can use to create and sample from distributions. I don't think you'll need this for a problem set just yet, but you'll have access to it for later.
255
256 First, let's create two sequences of numbers using the `seq` function:
257
258 ```{r}
259 odds <- seq(from=1, to=100, by=2)
260 evens <- seq(from=2, to=100, by=2)
261
262 head(odds)
263 head(evens)
264 ```
265 Or maybe (for some reason) I want to have a vector that repeats the odd integers between 1 and 100 five times? Try `rep()`:
266
267 ```{r}
268 more.odds <- rep(seq(from=1, to=100, by=2), 5)
269
270 ```
271 You can use `sample()` to draw samples from an object:
272
273 ```{r}
274 sample(x=odds, size=3)
275 sample(x=evens, size=3)
276 ```
277
278 You can also sample "with replacement." Here I take 100 random draws from the binomial distribution:
279 ```{r}
280 draws <- sample(x=c(0,1), size=100, replace=TRUE)
281
282 table(draws)
283 ```
284
285 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:
286
287 ```{r}
288 odds.n.evens <- data.frame(odds, evens)
289
290 odds.n.evens[ sample(row.names(odds.n.evens), 10), ]
291
292 ```
293
294
295 ### Managing randomness
296
297 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.
298
299 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:
300
301 ```{r}
302 set.seed(42) 
303 x <- sample(odds, 10)
304
305 y <- sample(odds, 10) # different
306
307 set.seed(42)
308 z <- sample(odds, 10) # 
309
310 head(x)
311 head(y)
312 head(z)
313
314 table(x==y) ##  One in ten!
315 table(x==z)
316
317
318 ```
319
320 ## Dates 
321
322 Date and time objects are another 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 will need 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).
323
324 To build up an example, I'll create some date-time values, add a little noise, and convert them into a character vector:
325
326 ```{r}
327 add.an.hour <- seq(0, 3600*24, by=3600)
328 some.hours <- as.character(Sys.time() + add.an.hour) ## Look up Sys.time() to see what it does.
329
330 class(some.hours); head(some.hours)
331 ```
332
333 These are beautifully formatted timestamps, but R will not see them as such. You can convert them back 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!
334
335 ```{r}
336 as.POSIXct(some.hours)
337 ```
338
339 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:
340
341 ```{r}
342 m <- "2019-02-21 04:35:00"
343 class(m)
344
345 a.good.time <- as.POSIXct(m, format="%Y-%m-%d %H:%M:%S", tz="CDT")
346 class(a.good.time)
347 ```
348
349 Once you have a time object, you can even do date arithmetic with `difftime()` (but watch out as this can get complicated):
350
351 ```{r}
352 difftime(Sys.time(), a.good.time, units="weeks")
353 ```
354
355

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