Lab08: Time

Today, we’re going to focus on visualizing data over time. This is an area where we can’t use linear regression – if our data is time-dependent, the observations are not independent from one another.

While regression can be adapted for forecasting in time-series (which is called ARIMA: Auto-Regressive Integrated Moving Average), we’ll focus on graphical methods for understanding time series, rather than predicting new observations.

Part 1: Time Series

A time series graph is usually a connected scatterplot. The first part of the lab will walk you through how to make some of the time series plots in Chapter 9 of The Truthful Art. We’ll return to the covid data from the Delphi group for this lab.

covid = read_csv("data/covid_data.csv")
covid

## # A tibble: 65,688 x 6
##    signal       geo_value time_value value stderr sample_size
##    <chr>        <chr>     <date>     <dbl>  <dbl>       <dbl>
##  1 smoothed_cli ak        2020-10-01 1.59  0.429         760 
##  2 smoothed_cli al        2020-10-01 0.983 0.141        3742.
##  3 smoothed_cli ar        2020-10-01 1.06  0.174        2587.
##  4 smoothed_cli az        2020-10-01 0.597 0.0905       5682.
##  5 smoothed_cli ca        2020-10-01 0.450 0.0399      21930.
##  6 smoothed_cli co        2020-10-01 0.561 0.0917       5137.
##  7 smoothed_cli ct        2020-10-01 0.444 0.0978       3866 
##  8 smoothed_cli dc        2020-10-01 0.268 0.224         533 
##  9 smoothed_cli de        2020-10-01 0.203 0.134        1129.
## 10 smoothed_cli fl        2020-10-01 0.516 0.0472      17767.
## # … with 65,678 more rows

This time, we want to examine this data over time. Since there’s a huge number of signals (see options with unique(covid$signal)), and 50 states + DC, we’ll restrict this analysis to Pennsylvania, and focus on the signal smoothed_tested_14d. This signal contains a smoothed estimate of the percentage of people who have been tested for COVID in the last 14 days.

pa_tested = covid %>%
  filter(signal == "smoothed_tested_14d" & geo_value == "pa")
pa_tested

## # A tibble: 92 x 6
##    signal              geo_value time_value value stderr sample_size
##    <chr>               <chr>     <date>     <dbl>  <dbl>       <dbl>
##  1 smoothed_tested_14d pa        2020-10-01  6.17  0.224      11505.
##  2 smoothed_tested_14d pa        2020-10-02  6.17  0.223      11622.
##  3 smoothed_tested_14d pa        2020-10-03  6.12  0.223      11558.
##  4 smoothed_tested_14d pa        2020-10-04  6.09  0.222      11564.
##  5 smoothed_tested_14d pa        2020-10-05  6.11  0.223      11588.
##  6 smoothed_tested_14d pa        2020-10-06  6.04  0.219      11817.
##  7 smoothed_tested_14d pa        2020-10-07  6.14  0.217      12279.
##  8 smoothed_tested_14d pa        2020-10-08  6.09  0.215      12423.
##  9 smoothed_tested_14d pa        2020-10-09  6.11  0.215      12443.
## 10 smoothed_tested_14d pa        2020-10-10  6.16  0.216      12379.
## # … with 82 more rows

Let’s make a basic time series line graph:

ggplot(pa_tested, aes(x = time_value, y = value, ymin = value - 2*stderr, ymax = value + 2*stderr)) +
  geom_line() +
  geom_ribbon(alpha = .5) +
  lab08_theme +
  labs(
    x = "",
    y = "",
    title = "% of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Smoothed Estimate",
    caption = "Source: Delphi Symptom Survey"
  )

I’ve included a geom_ribbon to emphasize that this is an estimate and so there is some uncertainty associated with it. In order to calculate the ymin and ymax values for the geom_ribbon, I’ve used the formula \(\hat{y} \pm 2 \times SE(\hat{y})\) for a 95% confidence interval.

Now, note that R recognizes that the time_value variable is a <date>. Date variable type comes in handy when making axes. In the graph below, I’ve changed the way that the date is formatted on the x-axis to %Y %b %d, which corresponds to “Year”, “abbreviated month”, “day”. Pull up the help page for scale_x_date, or the table here to see other options.

ggplot(pa_tested, aes(x = time_value, y = value, ymin = value - 2*stderr, ymax = value + 2*stderr)) +
  geom_line() +
  geom_ribbon(alpha = .5) + 
  scale_x_date(date_labels = "%Y %b %d") +
  lab08_theme +
  labs(
    x = "",
    y = "",
    title = "% of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Smoothed Estimate",
    caption = "Source: Delphi Symptom Survey"
  )

Now, let’s adjust the y-axis of the graph to show the deviation from the average rather than the raw percentage. This does not add much to this graph, and is mostly useful for looking at trends over many years where we expect a lot of seasonality. It’s hard to find seasonality in a 3 month period, but let’s do it anyways! Note that the graph looks exactly the same, except now the y-axis is centered at zero.

pa_tested %>%
  mutate(
    value_dev = value - mean(value)
  ) %>%
  ggplot(., aes(x = time_value, y = value_dev, ymin = value_dev - 2*stderr, ymax = value_dev + 2*stderr)) +
  geom_line() +
  geom_ribbon(alpha = .5) + 
  labs(
    x = "",
    y = "",
    title = "% of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Deviation of Smoothed Estimate from average over Oct 1 - Dec 31",
    caption = "Source: Delphi Symptom Survey"
  ) + 
  lab08_theme

Even though we shouldn’t fit a linear regression line, we can still add the loess smoother to the graph. All this does is compute a moving average, so it still works pretty well here:

pa_tested %>%
  mutate(
    value_dev = value - mean(value)
  ) %>%
  ggplot(., aes(x = time_value, y = value_dev, ymin = value_dev - 2*stderr, ymax = value_dev + 2*stderr)) +
  geom_line() +
  geom_ribbon(alpha = .5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(
    x = "",
    y = "",
    title = "% of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Deviation from average; includes LOESS trend",
    caption = "Source: Delphi Symptom Survey"
  ) +
  lab08_theme 

To annotate this graph, we can use the as.Date() function to pass x values to it:

pa_tested %>%
  mutate(
    value_dev = value - mean(value)
  ) %>%
  ggplot(., aes(x = time_value, y = value_dev, ymin = value_dev - 2*stderr, ymax = value_dev + 2*stderr)) +
  geom_line() +
  geom_ribbon(alpha = .5) +
  annotate("rect", xmin = as.Date("2020-11-26") - 14, xmax = as.Date("2020-11-26"), ymin = -4, ymax = 4, fill = "darkorange2", alpha = .3) +
  annotate("text", x = as.Date("2020-11-26") - 7, y = -3, size = 2.5, label = "2 Week Period\n before Thanksgiving", col = "darkorange4") +
  geom_smooth(method = "loess", se = FALSE) +
  labs(
    title = "% Of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Deviation from average; includes LOESS trend",
    x = "",
    y = "",
    caption = "Source: Delphi Symptom Survey"
  ) +
  lab08_theme

Finally, we might be interested in calculating the loess trend beforehand, and plotting the deviations from the moving average. In the graph above, this is equivalent to plotting the difference between the black line and the blue line.

We first have to calculate the loess model, and then we can augment it the same way we did with lm and glm. Note that we have to use as.numeric(time_value) in order to get the model to work. I’m also cheating here a little bit and adding a new column to the augmented dataset that contains the original, non-numeric time_value variable so that we can actually plot the results on the same graph.

loess_model = loess(value ~ as.numeric(time_value), data = pa_tested)
loess_aug = augment(loess_model) %>%
  mutate(time_value = pa_tested$time_value)

To plot the difference between the fitted and observed values, we can just make a ggplot using the results from augment:

ggplot(loess_aug) +
  geom_line(aes(x = time_value, y = .fitted - value)) +
  labs(
    title = "% Of Pennsylvanians who were tested for COVID over previous 14 days",
    subtitle = "Deviation from LOESS trend",
    x = "",
    y = "",
    caption = "Source: Delphi Symptom Survey"
  ) +
  lab08_theme

(One thing to be careful of with this data is that we already have a smoothed estimate – we’re basically “double smoothing” when we use LOESS here.)

Decomposing Time Series

Chapter 9 of The Truthful Art spent some time discussing seasonality and decomposing time series into the seasonal, trend, and random components. This is super easy in R! But it doesn’t work for this data (since we only have 3 months of data). Below is some code in “base R” that does not work, but shows the sequence of commands that would make such a plot for you.

pa_tested_ts = ts(pa_tested$value)
pa_tested_decomp = decompose(pa_tested_ts)
plot(pa_testeed_decomp)

You’ll notice this isn’t ggplot2! I don’t work with time series very often, but from what I can tell the ggseas package is the best bet for doing a similar thing in ggplot. Here’s a blog post with a quick example if this is something you’re interested in for any of your projects.

Part 2: Alluvial Diagrams

In order to make this graph, you’ll need to install the ggalluvial package from CRAN. You can do so either by running install.packages('ggalluvial') in the console, or by clicking “Tools –> Install Packages…” in the toolbar and then typing ggalluvial in the search bar. The data we’re using here is built into the package.

data(vaccinations)

First, we’ll do some cleaning so that our data displays nicely:

vaccinations = vaccinations %>%
  mutate(
    response = factor(response, rev(levels(response))),
    survey = case_when(survey == "ms153_NSA" ~ "First Survey",
                      survey == "ms432_NSA" ~ "Second Survey", 
                      TRUE ~ "Third Survey")) 
head(vaccinations)

##         survey freq subject response start_date   end_date
## 1 First Survey   48       1   Always 2010-09-22 2010-10-25
## 2 First Survey    9       2   Always 2010-09-22 2010-10-25
## 3 First Survey   66       3   Always 2010-09-22 2010-10-25
## 4 First Survey    1       4   Always 2010-09-22 2010-10-25
## 5 First Survey   11       5   Always 2010-09-22 2010-10-25
## 6 First Survey    1       6   Always 2010-09-22 2010-10-25

Examine the vaccinations data before the cleaning step, and see if you can figure out what these steps did. The case_when function is especially useful for dealing with factors/text (it’s essentially a nested if-else statement).

The syntax for ggalluvial is different than a lot of other geoms we’ve seen, and frankly it’s still confusing to me. I do think these plots are beautiful and useful, so I wanted to include them. If you want to learn more about how to make them, I recommend the package vignette.

Below is the code for a very basic alluvial diagram. Try commenting out (one line at a time) the geom_flow(), geom_stratum(alpha = .5), and geom_text(stat = "stratum", size = 3) to get a feel for what each of them does.

ggplot(vaccinations,
       aes(x = survey, 
           stratum = response, 
           alluvium = subject,
           y = freq,
           fill = response, 
           label = response)) +
  geom_flow() +
  geom_stratum(alpha = .5) +
  geom_text(stat = "stratum", size = 3) +
  theme(legend.position = "none") +
  ggtitle("vaccination survey responses at three points in time")

This graph is often used to display relationships between categorical data, but since it “feels” like the data is flowing from left to right, I find it especially useful for time-related data. It would of course be too busy for a usual time series, but it’s a good option if you have a few “snapshots” of a variable at different periods of time and want to visualize how the distribution changed.

To clean this up a bit, try some of the following to get as close to the below graph as you can. (I used the Optima font and scale_color_viridis(end = .75 ) color scale, but you are welcome to choose others!)

• geom_flow(width = 1/6)
• scale_x_discrete(expand = c(.1, .1))
• geom_stratum(col = "black", fill = "black", width = 1/6, alpha = .8)
• geom_text(stat = "stratum", size = 3, col = "white", angle = 45)

Part 3: Open-Ended

Using a different signal from the covid data (or more!), and/or restricting to any states of your choosing, make any of the time-related plots above, or something related to them. It may be helpful to return to the description of Delphi’s signal variables when choosing a signal and creating your titles, axes labels, and legend.

You should make sure to incorporate the following custom elements in your graph:

• X-axis labels using scale_x_time if making a time series
• Color palette
• Theme
• Informative title/subtitle/axes labels/legends, etc. (If you find you do not need one or more of these elements to fully explain what is going on in your graph, you can leave them out if you prefer!)