University of Iowa

College of Liberal Arts and Sciences

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

Time Series Plots

Time Series Data

A lot of data is collected over time:

• economic performace measures;

• weather, climate data;

• medical information.

The objectives can be

• understanding the past – how we got where we are;

• understanding patterns in variation;

• forecasting what will happen next.

Even if time is not of primary interest, data are often collected with a time stamp.

Examining the behavior over time can reveal interesting features that may hint at confounding variables:

• temperature varying with time of day;

• humidity varying over the year.

Time series data can be collected at regularly spaced intervals:

• yearly, e.g. the gapminder data

• monthly, e.g. the river data

• daily, e.g. high temperatures.

Data might be missing

• irregularly, e.g. faulty instrument, holidays;

• regularly, e.g. weekends.

Time series can also be recorded at irregular times:

• blood pressure at doctor visits;

Variables recorded over time can be numerical or categorical.

Time Series Objects

There are a number of specialized object classes for dealing with time series.

A number of packages provide plot or autoplot methods and other utilities for these objects.

Some of the most useful:

Class	Package	Features
ts, mts	base	regular series; methods for plot and lines.
	forecast	methods for autoplot and autolayer.
zoo	zoo	irregular and regular series; autoplot support.
	broom	tidy for converting to tidy data frames.

autoplot and autolayer

Many data structures provide methods for autoplot and autolayer.

From the autoplot help page:

autoplot uses ggplot2 to draw a particular plot for an object of a particular class in a single command. This defines the S3 generic that other classes and packages can extend.

From the autolayer help page:

autolayer uses ggplot2 to draw a particular layer for an object of a particular class in a single command. This defines the S3 generic that other classes and packages can extend.

Example: River Data

Creating a ts object for the river data:

river <- scan("https://www.stat.uiowa.edu/~luke/data/river.dat")
riverTS <- ts(river)
riverTS
## Time Series:
## Start = 1 
## End = 128 
## Frequency = 1 
##  [1]  8.95361  9.49409 10.19430 10.95660 11.07770 10.98170 10.48970 10.27540
##  [9] 10.16620  9.23484  8.56537  8.51064
##  [ reached getOption("max.print") -- omitted 116 entries ]

Using tidy from package broom to create a data frame from the ts object:

riverDF <- broom::tidy(riverTS)
head(riverDF, 3)
## # A tibble: 3 × 2
##   index value
##   <dbl> <dbl>
## 1     1  8.95
## 2     2  9.49
## 3     3 10.2

Change to more convenient names:

riverDF <- rename(riverDF, month = index, flow = value)
head(riverDF, 3)
## # A tibble: 3 × 2
##   month  flow
##   <dbl> <dbl>
## 1     1  8.95
## 2     2  9.49
## 3     3 10.2

Basic Plots for a Single Numeric Time Series

Numeric time series are usually plotted as a line chart.

This is what the autoplot method for ts objects provided by package forecast does.

For the river flow data:

library(forecast)
p <- autoplot(riverTS)
p

Using symbols in addition to lines is sometimes useful:

p + geom_point()

Aspect Ratio

Aspect ratio has a significant effect on how well we perceive slope differences.

• Banking to 45 degrees: choose an aspect ratio so the magnitude of the important slopes is close to 45 degrees.

• This usually requires the horizontal axis to be longer than the vertical axis.

Controlling the aspect ratio:

• For base graphics use the asp argument to plot.

• lattice uses the aspect argument.

• For ggplot use the ratio argument to coord_fixed.

p + coord_fixed(ratio = 4)

Lattice supports several aspect ratio strategies:

• "fill" to avoid unused white space;

• "xy" to use a 45 degree banking rule;

• "iso" for isometric scales, e.g. distances.

Some work may be needed to avoid excess white space around figures.

Using multiple panels can help.

In an interactive setting, zoom and pan support can help.

Multiple panels:

mutate(riverDF,
       period = cut_width(month,
                          center = 21.5,
                          width = 43)) |>
    ggplot(aes(month, flow)) +
    geom_line() +
    facet_wrap(~ period,
               scales = "free_x",
               ncol = 1)

Adding zoom and pan support using plotly:

library(plotly)
pp <- p + geom_line() + coord_fixed(ratio = 4)
ggplotly(pp)

There is not always a single best aspect ratio.

The co2 data set in the datasets package contains monthly concentrations of CO2 for the years 1959 to 1997 recorded at the Mauna Loa observatory.

The co2 data is stored as an object of class ts:

str(co2)
##  Time-Series [1:468] from 1959 to 1998: 315 316 316 318 318 ...

plot and xyplot have methods for ts objects that simplify time series plotting.

Two aspect ratios:

xyplot(co2)

xyplot(co2, aspect = "xy")

Is there a slowdown in the increase?

A graph with more recent data is available.

The full data is available as well.

autoplot does not yet provide an easy way to bank to 45 degrees; you need to experiment or calculate an appropriate aspect ratio yourself.

The default autoplot:

library(forecast)
autoplot(co2)

Aspect ratio adjusted with coord_fixed:

r45 <- median(diff(time(co2))) / median(abs(diff(co2)))
autoplot(co2) + coord_fixed(ratio = r45)

Updated Mauna Loa CO2 Data

Downloading and reading in the data:

library(lubridate)
url <- "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt"
if (! file.exists("co2new.dat") ||
    now() > file.mtime("co2new.dat") + weeks(4))
    download.file(url, "co2new.dat")
co2new <- read.table("co2new.dat")

Clean up the variable names:

head(co2new)
##     V1 V2       V3     V4     V5 V6    V7    V8
## 1 1958  3 1958.203 315.70 314.43 -1 -9.99 -0.99
## 2 1958  4 1958.288 317.45 315.16 -1 -9.99 -0.99
## 3 1958  5 1958.370 317.51 314.71 -1 -9.99 -0.99
## 4 1958  6 1958.455 317.24 315.14 -1 -9.99 -0.99
## 5 1958  7 1958.537 315.86 315.18 -1 -9.99 -0.99
## 6 1958  8 1958.622 314.93 316.18 -1 -9.99 -0.99
names(co2new) <- c("year", "month", "decimal_date", "average",
                   "deseason", "ndays", "stdd", "umm")

An initial look:

ggplot(co2new, aes(decimal_date, average)) + geom_line()

Missing values are encoded as -9.99 or -99.99.

If there are any such values they need to be replaced with R’s missing value representation:

filter(co2new, average < -9)
## [1] year         month        decimal_date average      deseason    
## [6] ndays        stdd         umm         
## <0 rows> (or 0-length row.names)
co2new <- mutate(co2new, average = ifelse(average < -9, NA, average))

This produces a reasonable plot:

ggplot(co2new, aes(decimal_date, average)) +
    geom_line()

The plot would show gaps for the missing values.

Missing values could also be interpolated with na.interp from package forecast.

Converting the new data to a ts object makes adding the original data a little easier:

co2newTS <- ts(co2new$average, start = c(1958, 3), frequency = 12)
autoplot(co2newTS) + autolayer(co2)

Seasonal Decomposition of Time Series by Loess (STL)

For time series with a strong seasonal component it can be useful to look at a Seasonal Decomposition of Time Series by Loess, or (STL).

The methodology was suggested by Clevaland and coworkers.

The stl function in the base package computes such a decomposition.

It requires a series without missing values.

co2newTS_no_NA <- ts(co2new$average, start = c(1958, 3), frequency = 12)
co2stl <- stl(co2newTS_no_NA, s.window = 21)
head(co2stl$time.series)
##            seasonal    trend  remainder
## Mar 1958  1.1151418 314.9331 -0.3482841
## Apr 1958  2.2812065 315.0057  0.1630786
## May 1958  2.7700298 315.0783 -0.3383173
## Jun 1958  2.2252349 315.1509 -0.1360949
## Jul 1958  0.8621909 315.2178 -0.2199862
## Aug 1958 -1.0841455 315.2847  0.7294149

Package forecast provides an autoplot method for stl objects that shows the original data and the three components in separate panels:

autoplot(co2stl)

Some Alternative Visualizations

Bar charts are sometimes used when comparison to a base line is important; for example for illustrating climate change.

Calendar plots, or calendar heatmaps, are sometimes useful for displaying daily data where day of the week/weekend effects might be important. One example is available here.

A blog post shows how to build a calendar plot with ggplot.

GDP Growth Data

Bar charts are sometimes used for economic time series, in particular growth rate data.

Some GDP growth rate data is available on this web page:

Read the data from the web page:

library(rvest)
gdpurl <- "http://www.multpl.com/us-real-gdp-growth-rate/table/by-quarter"
gdppage <- read_html(gdpurl)
gdptbl <- html_table(gdppage)[[1]]
head(gdptbl)
## # A tibble: 6 × 2
##   Date         Value
##   <chr>        <chr>
## 1 Dec 31, 2023 3.09%
## 2 Sep 30, 2023 2.93%
## 3 Jun 30, 2023 2.38%
## 4 Mar 31, 2023 1.72%
## 5 Dec 31, 2022 0.65%
## 6 Sep 30, 2022 1.71%

Clean the data:

library(dplyr)
library(lubridate)
names(gdptbl) <- c("Date", "Value")
gdptbl <- mutate(gdptbl,
                 Value = as.numeric(sub("%", "", Value)),
                 Date = mdy(Date))
head(gdptbl)
## # A tibble: 6 × 2
##   Date       Value
##   <date>     <dbl>
## 1 2023-12-31  3.09
## 2 2023-09-30  2.93
## 3 2023-06-30  2.38
## 4 2023-03-31  1.72
## 5 2022-12-31  0.65
## 6 2022-09-30  1.71

Data prior to 1947 is annual, so drop earlier data:

gdptbl <- filter(gdptbl,
                 Date >= make_date(1948, 1, 1))

Rearrange to put Date in increasing order:

gdptbl <- arrange(gdptbl, Date)

Create and plot a time series:

library(ggplot2)
library(forecast)
gdpts <- ts(gdptbl$Value,
            start = 1948, frequency = 4)
autoplot(gdpts) +
    geom_hline(yintercept = 0, linetype = 2)

A bar chart with color used to distinguish positive and negative growth can be created from the original table gdptbl:

library(ggplot2)
mutate(gdptbl,
       growth = ifelse(Value >= 0,
                       "pos",
                       "neg")) |>
    ggplot(aes(x = Date,
               y = Value,
               fill = growth)) +
    geom_col() +
    theme(legend.position = "top") +
    scale_fill_manual(
        values = c(pos = "blue", neg = "red"))

The bar chart can also be created from the time series object:

data.frame(rate = as.numeric(gdpts),
           time = as.numeric(time(gdpts))) |>
    mutate(growth = ifelse(rate >= 0,
                           "pos",
                           "neg")) |>
    ggplot(aes(time, rate, fill = growth)) +
    geom_col() +
    theme(legend.position = "top") +
    scale_fill_manual(
        values = c(pos = "blue", neg = "red"))

Airline Flights From NYC

The nycflights13 data set contains data on all flight leaving New York City in 2013.

Similar data is available for other cities and years, e.g. from the anyflights package.

We can look at the number of departures for each day:

library(nycflights13)
library(ggplot2)
nf <- count(flights, year, month, day)

It will be useful to have a combined date field in our data frame.

The lubridate function make_date function can be used for this:

library(lubridate)
nf <- mutate(nf, date = make_date(year, month, day))

A line plot of the number of flights:

ggplot(nf, aes(date, n)) +
    geom_line()

Weekday variation is clearly visible.

Breaking out the days of the week into separate facets may help:

nnf <-
    mutate(nf, wd = wday(date, label = TRUE))

nnf_thm <-
    theme(axis.text.x =
              element_text(angle = -45,
                           hjust = 0),
          plot.margin = margin(r = 30))

ggplot(nnf, aes(date, n)) +
    geom_line() +
    facet_wrap(~ wd) +
    nnf_thm

July 4 in 2013:

as.character(wday(ymd("2013-07-04"),
                  label = TRUE))
## [1] "Thu"

Identifying when some major holidays occur may also be useful:

holidays <-
    data.frame(Holiday = c("4th of July",
                           "Thanksgiving"),
               date = ymd(c("2013-07-04",
                            "2013-11-28")))
hlayer <-
    geom_vline(aes(xintercept = date,
                   color = Holiday),
               linetype = 2,
               linewidth = 1,
               data = holidays)
ggplot(nf, aes(date, n)) +
    hlayer +
    geom_line()

Holidays in the faceted plot:

ggplot(nnf, aes(date, n)) +
    hlayer +
    geom_line() +
    facet_wrap(~ wd) +
    nnf_thm

This marks the holiday weeks.

Holidays in the faceted plot:

This shows only the holiday days.

hlayer2 <-
    geom_vline(aes(xintercept = date,
                   color = Holiday),
               linewidth = 1,
               data = mutate(holidays,
                             wd = wday(date,
                                       TRUE)))
ggplot(nnf, aes(date, n)) +
    hlayer2 +
    geom_line() +
    facet_wrap(~ wd) +
    nnf_thm

Holidays in the faceted plot:

Show holidays as solid lines, days in holiday weeks as dashed lines.

ggplot(nnf, aes(date, n)) +
    hlayer +
    hlayer2 +
    geom_line() +
    facet_wrap(~ wd) +
    nnf_thm

A simple calendar plot:

This uses geom_tile and faceting.

monthweek <- function(d, w)
    ceiling((d - w) / 7) + 1

nf <- mutate(nf,
             wd = wday(date, label = TRUE),
             wd = fct_rev(wd),
             mw = monthweek(day, wday(date)))

p <- ggplot(nf, aes(x = as.character(mw),
                    y = wd,
                    fill = n)) +
    geom_tile(color = "white") +
    facet_wrap(~ month(date, label = TRUE))
p

Some theme and color adjustments:

p <- p +
    scale_fill_gradient(low = "white") +
    ylab("") + xlab("Week of Month") +
    theme(panel.grid.major = element_blank(),
          panel.border = element_blank())
p

This layout with days on the vertical axis works well with a single row, which could be extended to a multi-year plot:

p + facet_wrap(~ month(date, TRUE), nrow = 1)

A more standard layout:

nf2 <-
    mutate(nf,
           wd = wday(date, label = TRUE),
           mw = factor(monthweek(day,
                                 wday(date))),
           mw = fct_rev(mw))

ggplot(nf2, aes(x = wd, y = mw, fill = n)) +
    geom_tile(color = "white") +
    facet_wrap(~ month(date, TRUE)) +
    scale_fill_gradient(low = "white") +
    ylab("") + xlab("Week of Month") +
    theme(panel.grid.major = element_blank(),
          axis.text.x =
              element_text(angle = 45,
                           hjust = 1),
          panel.border = element_blank())

Reading

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

Exercises

1. A useful feature of the lattice package is the ability to easily request an aspect ratio that is banked to 45 degrees by specifying aspect = "xy" in the arguments to lattice plotting functions. For example, This code produces a reasonable aspect ratio for a plot of an unemployment time series:

   library(ggplot2)
   library(forecast)
   library(lattice)
   unemp_ts <- ts(economics$unemploy, start = c(1967, 7), frequency = 12)
   xyplot(unemp_ts, aspect = "xy")

   

   Which of the following plots most closely matches the aspect ratio in the lattice plot?

   1. autoplot(unemp_ts, aspect = 200)
   2. autoplot(unemp_ts) + coord_fixed(ratio = 0.0007)
   3. autoplot(unemp_ts) + coord_fixed(asp = 0.0007)
   4. autoplot(unemp_ts) + coord_fixed(70)

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