]> code.communitydata.science - stats_class_2020.git/blob - r_tutorials/w05a-R_tutorial.rmd
705978490b063dce94e18901311851e5ed09a0d4
[stats_class_2020.git] / r_tutorials / w05a-R_tutorial.rmd
1 ---
2 title: "Week 5 R tutorial (supplement)"
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 library(formatR)
22 knitr::opts_chunk$set(echo = TRUE, tidy='styler', message = FALSE)
23 ```
24
25 # Getting started (more better plots)
26 This is a supplement to the Week 5 R tutorial focused on elaborating some examples of time series plots and more polished plots using [`ggplot2`](https://ggplot2.tidyverse.org). I'll work some data on state-level COVID-19 in the United States published by *The New York Times* (*NYT*). You can access the data an details about the sources, measurement, and different datasets available via [the *NYT* github repository](https://github.com/nytimes/covid-19-data). 
27
28 To start, I'll load up the `tidyverse` library and also attach the `lubridate` package to help handle dates and times. Then I'll import the "raw csv" from the web, and take a look at the dataset:
29
30 ```{r}
31 library(tidyverse)
32 library(lubridate)
33
34 data_url <- url("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv")
35
36 d <- read_csv(data_url)
37
38 d
39 ```
40
41 For the sake of my examples, I'm planning to work with the `date`, `state`, `cases`, and `deaths` variables. Notice that by using the `read_csv()` function to import the data, R already recognizes the `date` column as dates. It looks like I need to convert the state variable to a factor, however. After I do that I can get a quick sense of how much data I have for each state with a univariate table that just counts the number of observations (rows) for each value of `state`.
42 ```{r}
43 d$state <- factor(d$state)
44 table(d$state)
45 ```
46
47
48
49 # Plotting a univariate time series
50
51 I recommend using [`geom_path()`](https://ggplot2.tidyverse.org/reference/geom_path.html) to create univariate time series plots. Specifically, I'll call `geom_line()`, which is a specialized version of `geom_path()` that connects observations in order according to the values of variable that is mapped to the x-axis. By convention, a univariate time series maps dates to the x-axis, so this will just plot a line connecting the dots over time.
52
53 For my first example, I want to build up a plot of weekly case counts in Illinois. I can start off by just plotting the cumulative cases for all of the states and work my way towards the specific plot I want from there:
54
55 ```{r}
56 ggplot(data=d, aes(date, cases)) + geom_line()
57 ```
58
59 Notice that ggplot handles the `date` variable quite well by default! It recognizes the units of time and generates axis labels in terms of months. Also notice that ggplot handles the axis labels for the `cases` variable...less well. I don't know about you, but my brain doesn't parse scientific notation quickly/easily.
60
61 # Tidying timeseries data for better plots
62
63 Okay, let's get to work cleaning all this up. At this point, my next steps are to (1) restrict the data to the Illinois cases; (2) reorganize the *cumulative* daily case counts into weekly counts; and (3) plot it again with better axis labels and a nice title. 
64
65 I can restrict the data to Illinois in a few ways. Since I'm using ggplot, I'll work with Tidyverse "pipes" (`%>%`) and "verbs" (in this case, `filter`):
66 ```{r}
67 d %>%
68   filter(state == "Illinois") %>%
69   ggplot(aes(date, cases)) + geom_line()
70 ```
71
72 That's already much less cluttered. Inserting a call to the Tidyverse  `mutate`, `group_by`, and `summarize` verbs can help me generate the weekly counts I'm looking for. Here's the code to produce a new object. I'll walk through it below:
73 ```{r}
74 il_weekly_cases <- d %>%  
75   filter(state == "Illinois") %>%  
76   mutate(diff_cases = c(cases[1], diff(cases, lag = 1)),
77          weekdate = cut(date, "week")) %>%
78   group_by(weekdate) %>%
79   summarize(new_cases = sum(diff_cases, na.rm = T),)
80
81 il_weekly_cases
82 ```
83 There's quite a lot happening there. I'll go through it verb-by-verb.
84
85 First, I use `mutate` to create a `diff_cases` variable that disaggregates the cumulative values of `cases` (read the documentation for `diff` to learn more about this one). Differenced values alone wouldn't produce the same number of items (try running `length(1:10)` and compare that with `length(diff(1:10, 1))` to see what I mean), so I stores the first value of my `cases` variable and then append the differenced values after that. Within the same call to mutate I also create a new variable `weekdate` that collapses the dates into weeks (see the documentation for `cut.Date`) and stores the resulting strings as factors (e.g., a factor where the levels correspond to a series of Mondays: "2020-01-20", "2020-01-27"...). Hopefully, so far so good?
86
87 Next, I use `group_by` to aggregate everything by my `weekdate` factor values. 
88
89 Finally I use `summarize` to reshape my data and collapse everything into weekly counts of new cases (notice that I use `sum` inside the `summarize` call to add up the case counts within the grouping variable). 
90 Okay, let's see about plotting this now:
91
92 Hmm. looks like I have a problem with my dates. Let's troubleshoot this:
93
94 ```{r}
95 class(il_weekly_cases$weekdate)
96 ```
97 Whoops. It looks like I need to convert that `weekdate` variable into an object of class "date" so that it will work with ggplot. There are a number of ways I could do this, but I'll just make a new variable by first converting `weekdate` to a character vector and then converting that into a date using `as.Date` (and remember that it is sometimes easier to read these "nested" commands from the inside-out).
98 ```{r}
99 il_weekly_cases$date = as.Date(as.character((il_weekly_cases$weekdate)))
100 il_weekly_cases
101 ```
102 That ought to work now:
103 ```{r}
104 plot1 <- il_weekly_cases %>%
105   ggplot(aes(date, new_cases)) + geom_line()
106
107 plot1
108 ```
109
110 Much better! Notice that the final week of the data appears to fall off a cliff. That's just an artifact of the way that the *NYT* has published the data for part of the most recent week. Once it updates, the case count probably won't drop like that (yikes). Anyhow, onwards to cleaning things up and adding a title.
111
112 # Working on ggplot axis labels, titles, and scales
113 As I mentioned briefly in class `ggplot2` treats labels, titles, and scales as "layers" within it's "grammar of graphics" (and yes, I'm rolling my eyes as I type those scare-quotes).  For the purposes of our example here I'm going to use `scale_date` to work with the x-axis, `scale_continuous` to work with the y-axis, and `labs` to clean up the title and axis labels.
114
115 For starters, let's see whether there might be any way I want to improve the axis labels. The ggplot defaults for my `date` variable are pretty good already, but maybe I want to incorporate a label/break for each month as well as a more granular grid in the background that shows the weeks? Here's what all of that looks like:
116
117 ```{r}
118 plot2 <- plot1 + scale_x_date(date_labels = "%b", date_breaks= "1 month", date_minor_breaks = "1 week")
119 plot2
120 ```
121
122 The ggplot documentation for [`scale_date`](https://ggplot2.tidyverse.org/reference/scale_date.html) can give you some other examples and ideas. Also, notice how I appended the `scale_date` layer to my existing plot and stored it as a new object? This can make it easier to work iteratively without losing any of my earlier layers along the way. 
123
124 Now I can fix up the y-axis labels a bit using a call to the `labels` argument after I load the `scales` package.
125 ```{r}
126 library(scales)
127 plot3 <- plot2 + scale_y_continuous(label=comma)
128 plot3
129 ```
130
131 Nearly done. All that's left is a title and better axis names. I'll do that with yet another layer.
132 ```{r}
133 plot4 <- plot3 + labs(x="Week (in 2020)", y="New cases", title="COVID-19 cases in Illinois")
134 plot4
135 ```
136
137 Last, but not least, I mentioned in our class session that ggplot also has "themes" that can be useful for styling plots. One I have used for publications is the "light" theme. Here's how to apply that:
138 ```{r}
139 plot4 + theme_light()
140 ```
141
142 That's looking much better than when we started! If you wanted to export it as a standalone file (e.g., .png, .pdf, or whatever), I recommend looking at the documentation for the `ggsave()` function, which is available via ggplot2. Base R also has a `save()` function that you can work with, although it can be a bit more complicated to get comfortable with.
143
144 # Long versus wide data (and why long data is often helpful)
145
146 So what if you wanted to plot a multivariate time series (e.g., the same plot for more than one state and/or for more than one measure)? As always, you have a number of options, but the most effective way to achieve this with ggplot involves learning to work with "long" format data.
147
148 Thus far, we have worked mostly with "wide" format data where (nearly) every row corresponds to a single unit/observation and every column corresponds to a variable (for which we usually have no more than one value attributed to any unit/observation). Wide format data is great for many things, but it turns out that learning to work with long format data can be super helpful for a number of purposes. Producing richer, multidimensional ggplot visualizations is one of them.
149
150 Consider the format of my tidied dataframe that I used for plotting:
151 ```{r}
152 il_weekly_cases
153 ```
154 This dataframe is in a "wide" format. Each row is a week and each column is a variable unique to that week.
155
156 Our original dataframe was a bit "longer":
157 ```{r}
158 d
159 ```
160 We see multiple observations per state (I think I would say the units or rows correspond to "state-dates" or something like that). It's not completely "long" however, because we also have multiple columns corresponding to the two variables of interest: `cases` and `deaths`. The point I want to make is that there are a number of ways we can make this data "longer." For the purposes of producing a multi-state plot like the one above, the most important of these is going to involve dropping the step where I filtered by `state=="Illinois"` and replacing by a `group_by` step before I create my `weekdate` variable. I'm also going to go ahead and drop the `date` and `fips` variables because they're just getting in my way at this point. I'll start there
161 ```{r}
162 weekly <- d %>%
163   group_by(state) %>%
164   mutate(
165     weekdate = cut(date, "week"),
166   ) %>% select(state, cases, deaths, weekdate)
167 weekly
168 ```
169 I'm getting somewhere with this, I promise. One of the principles of "tidy" data is to make it so that every variable has a column, every observation has a row, and every value has a cell. Right now, I've got multiple observations for each state-week spread across multiple rows. Remember that my `cases` and `deaths` variables are actually cumulative counts, so I really only need to store the maximum value for each state-week in order to calculate the new cases per state-week. Let's see what to do about that:
170
171 ```{r}
172 tidy_weekly <- weekly %>% 
173   group_by(state, weekdate) %>%
174   summarize(
175     cum_cases = max(cases, na.rm=T),
176     cum_deaths = max(deaths, na.rm=T)
177     )
178 ```
179
180 ```{r}
181 tidy_weekly$weekdate <- as.Date(as.character(tidy_weekly$weekdate))
182
183 tidy_weekly <- tidy_weekly %>% 
184   group_by(state) %>%
185   arrange(-desc(weekdate)) %>%
186   mutate(
187     new_cases = c(cum_cases[1], diff(cum_cases, lag=1)),
188     new_deaths = c(cum_deaths[1], diff(cum_deaths, lag=1)),
189   ) 
190   
191 tidy_weekly
192 ```
193 This is headed in the right direction. For some purposes, though, it's still not quite "long" enough For starters, I can drop the cumulative cases and deaths columns. The other thing I can do is "pivot" the data to organize the `new_cases` and `new_deaths` measures a little differently. To manage this, I'll use the `pivot_longer()` function (part of the `tidyr` package from the tidyverse). I will also go ahead and coerce my `weekdate` into a Date object again:
194 ```{r}
195 long_weekly <- tidy_weekly %>%
196   select(state, weekdate, new_cases, new_deaths) %>%
197   pivot_longer(
198     cols = starts_with("new"),
199     names_to = "variable",
200     values_to = "value"
201   )
202
203 long_weekly
204 ```
205 Can you see what that did? I now have two rows of data for every state-week. One that contains a value for `new_cases` and one that contains a value for `new_deaths`. Both of those variables have been "pivoted" into a single `variable` column. 
206
207 Before we move forward I'm going to clean up the values of `variable`.
208 ```{r}
209 long_weekly <- long_weekly %>%
210   mutate(
211     variable = recode(variable, new_cases = "new cases", new_deaths = "new deaths")
212   )
213
214 ```
215 Okay, prepared with my `tidy_weekly` and my `long_weekly` tibbles, I'm now ready to generate some more interesting multidimensional plots. Let's start with the same sort of time series of new cases we made for Illinois before so we can see how to replicate that with this new data structure:
216 ```{r}
217 long_weekly %>% filter(
218   state == "Illinois" & variable == "new cases"
219 )  %>% ggplot(aes(weekdate, value)) + geom_line()
220
221 ```
222
223 Now we can easily plot Illinois cases against deaths from the same tibble:
224 ```{r}
225 long_weekly %>% filter(state == "Illinois") %>% 
226   ggplot(aes(weekdate, value, color=variable)) + geom_line()
227 ```
228
229 That plot isn't so great because the death counts are dwarfed by the case counts. Thank goodness!
230
231 Now let's compare Illinois case counts against some its neighbors in the upper midwest:
232 ```{r}
233 upper_midwest = c("Illinois", "Michigan", "Wisconsin", "Iowa", "Minnesota")
234
235 long_weekly %>% 
236   filter(state %in% upper_midwest & variable == "new cases") %>% 
237   ggplot(aes(weekdate, value, color=state)) + geom_line()
238 ```
239
240 Now that's getting a bit more interesting.
241
242 What about finding some way to also incorporate the death counts? Well, ggplot has another layer option called "facets" that can help produce multiple plots and present them alongside each other (or in a grid). Here's an example that creates a faceted "grid" (really just a side-by-side comparison) of case counts and deaths for the same five states.
243 ```{r}
244 midwest_plot <- long_weekly %>% 
245   filter(state %in% upper_midwest) %>% 
246   ggplot(aes(weekdate, value, color=state)) + geom_line() + facet_grid(rows=vars(variable), scales="free_y")
247
248 midwest_plot
249 ```
250
251 Now we can clean up some of the other elements we worked on with the original plot (axes, title, etc.). I'll bake that into a single chunk below.
252
253 ```{r}
254 midwest_plot + scale_x_date(date_labels = "%b", date_breaks= "1 month", date_minor_breaks = "1 week") + scale_y_continuous(label=comma) + labs(x="Week (in 2020)", y="", title="COVID-19 cases in the Upper Midwest") + theme_light()
255 ```
256

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