University of Iowa

College of Liberal Arts and Sciences

• STAT:4580
• Syllabus
• Outline
• Homework
• Notes
• Resources
-->

Plots for Multiple Time Series

Multiple Time Series

Some questions for multiple time series:

• How do the series behave individually?

• How are the series related?

• Does one series lead or lag the others?

Some options for plotting multiple series:

• separate panels in a trellis display;

• multiple series in a single plot - this will require standardizing if the scales vary substantially;

• a multivatiate plot with time represented by connecting line segments or animation.

Another option for two time series:

• use a separate y axis for each.

This is usually a bad idea.

Some data structures for multiple time series:

• mts objects created by ts().

  • All series have to be numeric.
  • autoplot methods may be sufficient.

• Regular data frames.

  • The series will have to be plotted in separate layers,
  • or the data frames will need converting into long form.

Crimean War Casualties

The HistData package contains the data used by Florence Nightingale to call attention to terrible sanitary conditions for the Britich army in Crimea.

data(Nightingale, package = "HistData")
head(select(Nightingale, -Date))
##   Month Year  Army Disease Wounds Other Disease.rate Wounds.rate Other.rate
## 1   Apr 1854  8571       1      0     5          1.4         0.0        7.0
## 2   May 1854 23333      12      0     9          6.2         0.0        4.6
## 3   Jun 1854 28333      11      0     6          4.7         0.0        2.5
## 4   Jul 1854 28722     359      0    23        150.0         0.0        9.6
## 5   Aug 1854 30246     828      1    30        328.5         0.4       11.9
## 6   Sep 1854 30290     788     81    70        312.2        32.1       27.7

A plot of the Disease, Wounds, and Other:

nDWO <- select(Nightingale,
               Disease, Wounds, Other) |>
    ts(start = c(1854, 4), frequency = 12)

autoplot(nDWO)

Using the rate variables produces a very similar pattern:

nDWO.rate <- select(Nightingale,
                    ends_with("rate")) |>
    ts(start = c(1854, 4), frequency = 12)

autoplot(nDWO.rate)

The army size does not vary enough over the time where the casualties are high to affect the pattern much:

nArmy <- select(Nightingale, Army) |>
    ts(start = c(1854, 4), frequency = 12)

autoplot(nArmy)

Instead of showing three series in one plot you can also use faceting.

autoplot(nDWO, facets = TRUE)

To create the plot of Disease, Wounds, and Other using basic ggplot it is best to first convert the data into long form.

You can create a long version of the data, using either pivot_longer or broom::tidy.

Using broom:tidy on the mts object:

broom::tidy(nDWO) |>
    rename(Date = index,
           Deaths = value,
           Cause = series) |>
    ggplot(aes(x = Date,
               y = Deaths,
               color = Cause)) +
    geom_line()

Using pivot_longer on the original data frame:

select(Nightingale,
       Date, Disease, Wounds, Other) |>
    pivot_longer(-Date,
                 names_to = "Cause",
                 values_to = "Deaths") |>
    ggplot(aes(x = Date,
               y = Deaths,
               color = Cause)) +
    geom_line()

Using faceting instead of a single plot:

select(Nightingale,
       Date, Disease, Wounds, Other) |>
    pivot_longer(-Date,
                 names_to = "Cause",
                 values_to = "Deaths") |>
    ggplot(aes(x = Date, y = Deaths)) +
    geom_line() +
    facet_wrap(~ Cause, ncol = 1)

To focus on the shape of the curves it is often useful to allow the vertical axis ranges to adjust to the data:

select(Nightingale,
       Date, Disease, Wounds, Other) |>
    pivot_longer(-Date,
                 names_to = "Cause",
                 values_to = "Deaths") |>
    ggplot(aes(x = Date, y = Deaths)) +
    geom_line() +
    facet_wrap(~ Cause, ncol = 1,
               scales = "free_y")

This is particularly useful when the ranges are very different, for example if the Army series is included:

select(Nightingale,
       Army, Date, Disease, Wounds, Other) |>
    pivot_longer(-Date,
                 names_to = "series",
                 values_to = "Deaths") |>
    ggplot(aes(Date, Deaths)) +
    geom_line() +
    facet_wrap(~ series, ncol = 1,
               scales = "free_y")

Some Notes

Using autoplot is often a quick way to get useful plots.

Using basic ggplot graphics is a bit more work but provides more flexibility.

Many Time Series

Annual flow for dams in Italy:

Clearly showing all series in a single plot is not helpful!

data(hydroSIMN, package = "nsRFA")
annualflows <- mutate(annualflows,
                      cod = factor(cod))

ggplot(annualflows,
       aes(x = anno, y = dato,
           group = cod,
           color = factor(cod))) +
    geom_line()

With small multiples some patterns can be seen:

ggplot(annualflows,
       aes(x = anno, y = dato,
           group = cod)) +
    geom_line() +
    facet_wrap(~ cod) +
    theme(axis.text.x =
              element_text(angle = 45,
                           hjust = 1))

Small multiples can also be arranged in ways that convey spatial information:

library(geofacet)
ggplot(state_unemp, aes(year, rate)) +
    geom_line() +
    facet_geo(~ state,
              grid = "us_state_grid2",
              label = "name") +
    scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
    labs(title = "Seasonally Adjusted US Unemployment Rate 2000-2016",
         caption = "Data Source: bls.gov",
         x = "Year",
         y = "Unemployment Rate (%)") +
    theme(strip.text.x = element_text(size = 6),
          axis.text.y = element_text(size = 6),
          axis.text.x = element_blank())

The Fall and Rise of Income Inequality

A 2015 post on the NPR Planet Money site uses two graphs to show transitions in income inequality.

The page includes a link to the data source.

This has been unavailable for some time.

I tried to obtain the same data but may not have done it quite right.

The closest source I could find is https://wid.world/data/.

The data is available locally:

income <- read.table("https://www.stat.uiowa.edu/~luke/data/nprincome.dat")

Plots over time:

select(income, -year) |>
    ts(start = 1913) |>
    autoplot()

The ranges are too different for this to be effective.

One alternative: standardize the series to one, or 100%, at some point.

Here the initial value in 1913 is used:

select(income, -year) |>
    transmute(L90rel = L90 / L90[1],
              U01rel = U01 / U01[1]) |>
    ts(start = 1913) |>
    autoplot()

Faceted plots, with free \(y\) axes, are another option:

select(income, -year) |>
    ts(start = 1913) |>
    autoplot(facets = TRUE)

With two series another option is a connected scatter plot.

ggplot(income, aes(U01, L90)) +
    geom_path(linewidth = 2) +
    labs(x = "Average income for the top 1%",
         y = "Average income for the bottom 90%")

This simple form does not show the direction or scale of time.

Color or annotation can be used to show direction and scale of time:

library(ggrepel)
ggplot(income, aes(U01, L90, color = year)) +
    geom_path(size = 2) +
    geom_text_repel(
        aes(label =
                ifelse(year %% 5 == 0 |
                       year == max(year),
                       year,
                       "")),
        color = "grey50") +
    labs(x = "Average income for the top 1%",
         y = "Average income for the bottom 90%")

Animated is another option:

p <- ggplot(income,
            aes(x = U01, y = L90,
                color = year)) +
    geom_path(size = 2,
              color = "lightgrey",
              data = mutate(income,
                            year = NULL)) +
    geom_path(size = 2) +
    labs(title = "Year: {frame_along}",
         x = "Average income for the top 1%",
         y = "Average income for the bottom 90%")

library(gganimate)
animate(p + transition_reveal(year),
        fps = 5, end_pause = 10)

Stack Loss Data

Operational data obtained from 21 days of operation of a plant for the oxidation of ammonia (NH3) to nitric acid (HNO3).

The dependent variable, stack.loss is 10 times the percentage of the ingoing ammonia to the plant that escapes from the absorption column unabsorbed.

library(lattice)
splom(stackloss)

The data are collected over time:

library(dplyr)
sl <- mutate(stackloss,
             intex = seq_len(nrow(stackloss)))
cpal <- colorRampPalette(c("#132B43", "#56B1F7"))
cols <- cpal(nrow(stackloss) - 1)
pfun <- function(x, y, ...) {
    n <- length(x)
    lsegments(x[-n], y[-n], x[-1], y[-1],
              col = cols, lwd = 3)
}
splom(stackloss, panel = pfun)

Randomly reordering the indices can help calibrate whether there is a time effect:

splom(sample_n(stackloss, nrow(stackloss)),
      panel = pfun)

This is sometimes called graphical inference.

Hadley Wickham, Dianne Cook, Heike Hofmann, and Andreas Buja, “Graphical Inference for Infovis,” IEEE Transactions on Visualization and Computer Graphics, Vol. 16, No. 6, November/December 2010.

Picking one plot and embedding the true order among 19 random orderings:

True order: plot 20.

The code:

N <- 20
k <- sample(N, 1)

ord <- seq_len(nrow(stackloss))
ss <- function(i) {
    if (i == k)
        sl <- stackloss
    else
        sl <- sample_n(stackloss, nrow(stackloss))
    mutate(sl, order = ord, sample = i)
}

d <- lapply(1 : N, ss)
d <- do.call(rbind, d)

library(ggplot2)
ggplot(d, aes(x = Acid.Conc., y = stack.loss, color = order)) +
    geom_path(size = 0.8) + facet_wrap(~ sample) +
    guides(color = "none")

Lynx-Hare Data

Trapping data on Canadian lynx and snowshoe hare:

Cycle patterns are similar.

lynxhare <- read.table("https://www.stat.uiowa.edu/~luke/data/LynxHare.txt")
names(lynxhare) <- c("year", "hare", "lynx")
lhts <- ts(select(lynxhare, -year), start = 1845)
autoplot(lhts, facets = TRUE)

Lynx cycles often lag behind:

autoplot(lhts)

Things to Look For

Changes in data definition.

• Administrative definitions, such as employment/unemployment can change over time.

• Administrative entities, such as countries and states, can change over time.

• Instruments can change (satelite imagery replacing surface measurements).

Outliers need not be globally outlying: departures from local patterns may also be important.

Level shifts and changes in volatility.

Shocks affecting level or volatility.

Australian Airline Data

Weekly passenger loads, in thousands, on Ansett Airlines between Australia’s two largest cities (available as melsyd in package fpp2):

Some features:

• Strike in 1989.

• Increase in passenger load starting in second half of 1991.

• Dip in economy, rise in business in 1992 corresponding to an experiment with fewer economy/more business seats.

data(melsyd, package = "fpp2")
autoplot(melsyd, facets = TRUE)

Weekly SP 500 Closing Returns

Weekly closing returns of the SP 500 from 2003 to September, 2012:

Some features:

• Increased volatility at the start of the financial crisis.

• Volatility seems to be decreasing, but not very fast.

library(xts)
library(astsa)
autoplot(sp500w)

Forecasting

Studying the past may help predict the future.

Forecasting a time series usually involves choosing a model and running the model forward.

Accuracy of forecasts decreases rapidly the farther ahead the forecast is made.

An example from the forecast package: WWWusage is a time series of the numbers of users connected to the Internet.

library(forecast)
fit <- ets(WWWusage)
fc <- forecast(WWWusage, model = fit)
autoplot(fc)

Labeling Trend Lines

For trend lines using a standard guide can be a little awkward.

Average life expectancy against year for the continents in the gapminder data:

library(gapminder)
d <- group_by(gapminder,
              continent, year) |>
    summarize(avgLifeExp = mean(lifeExp)) |>
    ungroup()

p <- ggplot(d,
            aes(x = year,
                y = avgLifeExp,
                color = continent)) +
    geom_line(size = 1)
p

With the legend on the right, a small improvement is to order the levels to match the order of the right hand curve values:

## Warning: <ggplot> %+% x was deprecated in ggplot2 4.0.0.
## ℹ Please use <ggplot> + x instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

d_last <- filter(d, year == max(year)) |>
    mutate(continent =
               reorder(continent, -avgLifeExp))

p %+%
    mutate(d,
           continent =
               factor(continent,
                      levels(d_last$continent)))

Direct labeling is often a better choice.

One option is to use geom_text at the ends of the lines.

This may need some additional adjustments to margins and clipping behavior:

p <- p + guides(color = "none")

d_last <- filter(d, year == max(year))

p + geom_text(aes(label = continent),
              data = d_last,
              size = 5,
              hjust = -0.1) +
    coord_cartesian(clip = "off",
                    expand = FALSE) +
    theme(plot.margin =
              margin(30, 60, 10, 10))

Another option is to use a second axis:

p + scale_y_continuous(
        sec.axis = dup_axis(
            breaks = d_last$avgLifeExp,
            labels = d_last$continent,
            name = NULL)) +
    scale_x_continuous(expand = c(0, 0))

It is also possible to place the labels along the lines.

The geomtextpath package provides support for this.

library(geomtextpath)
ggplot(d, aes(x = year,
              y = avgLifeExp,
              color = continent)) +
    geom_textline(aes(label = continent),
                  linewidth = 1,
                  vjust = -0.5,
                  hjust = 0.95) +
    guides(color = "none")

Bump Charts

Top 10 countries in terms of GDP per capita in 2007:

top_gdp <- filter(gapminder,
                  continent == "Europe",
                  year == max(year)) |>
    slice_max(gdpPercap, n = 10) |>
    pull(country)
d <- filter(gapminder,
            country %in% top_gdp) |>
    group_by(year) |>
    mutate(
        rank = min_rank(desc(gdpPercap))) |>
    ungroup()
d_last <- filter(d, year == max(year))

ggplot(d, aes(year,
              gdpPercap,
              color = country)) +
    geom_line()

The dominant feature is the overall trend.

One option: show deviations from the overall trend.

Another option: show ranks within years.

This is called a bump chart.

A simple bump chart:

ggplot(d, aes(year, rank, color = country)) +
    geom_line(size = 1) +
    scale_y_reverse(breaks = 1:10) +
    geom_text(aes(label = country),
              data = d_last,
              size = 5,
              hjust = -0.1) +
    coord_cartesian(clip = "off",
                    expand = FALSE) +
    theme(plot.margin =
              margin(30, 90, 10, 10)) +
    guides(color = "none")

The ggbump package provides geom_bump for drawing smooth bump chart lines:

library(ggbump)
ggplot(d, aes(year, rank, color = country)) +
    geom_bump(size = 1) +
    scale_y_reverse(breaks = 1:10) +
    geom_text(aes(label = country),
              data = d_last,
              size = 5,
              hjust = -0.1) +
    coord_cartesian(clip = "off",
                    expand = FALSE) +
    theme(plot.margin =
              margin(30, 90, 10, 10)) +
    guides(color = "none")

Other Approaches

A talk from 2019 describes a wide range of options for visualizing temporal data.

Reading

Chapters Visualizing time series and other functions of an independent variable and Visualizing trends in Fundamentals of Data Visualization.

Exercises

1. Creating a time series or multiple time series object and using autoplot can be an efficient way to obtain line plots of time series data. The following code uses the weather table in the nycflights13 package to create a data frame containing minimum and maximum temperatures for each month in 2013:

   library(tidyverse)
   library(nycflights13)
   library(lubridate)
   library(forecast)
   
   temps <- group_by(weather, year, month) |>
       summarize(max_temp = max(temp, na.rm = TRUE),
                 min_temp = min(temp, na.rm = TRUE)) |>
       ungroup()

   Which of the following expressions produces a line plot showing the minimum and maximum temperature series (and no other series)?

   1. autoplot(temps)
   2. autoplot(temps$max_temp, temps$min_temp)
   3. ts(temps) |> autoplot()
   4. select(temps, max_temp, min_temp) |> ts() |> autoplot()

This work is licensed under the CC BY-NC 3.0 Creative Commons License.