January 11, 2019

Pythagorean win expectation (a simple application of the tidyverse)

At the All Star break, a little over mid-way through the 2018 Major League Baseball (MLB) season, the Seattle Mariners were sitting in second place in the American League West, with a record of 58 wins and 39 losses, a winning percentage of .598. This success had been in spite of a negative run differential; they had scored 2 fewer runs than they had allowed over the 97 games played to that point. They had been losing some games as blowouts, and had been winning a lot of close games.
The Mariners’ success had been noted throughout the season;
Of course, regression toward the mean is a thing, so we might have anticipated that by the season’s end the Mariners’ win/loss ratio would more closely reflect their run differential. But did it?

I’ve written before about run scoring and prevention (index here). This time, I will look at the simplest of the approaches to calculating “win expectation” that have burbled up in the sabermetric community over the years; the other approaches may be worthy of consideration for a future post. This exercise will also give us a way, in subsequent posts, to look at the ways that the statistical programming language R works with regression models.

Pythagorean win ratio

Bill James, the godfather of sabermetrics, developed the Pythagorean win expectation model (wikipedia page). The basic idea is that there is a relationship between the runs a team scores (\(RS\)) and allows (\(RA\)), and the proportion of the games that they can be expected to win (\(WE\)). The equation is expressed thus:
\[ WE = RS^2 / (RS^2 + RA^2)\]

an R function

Let’s write a little R function for this equation…in that way, we can save some typing later.
winexp_fun <- function(RS, RA) {
  RS^2 / (RS^2 + RA^2)
}

The data

First, we’ll load the packages we need. Note that tidyverse contains multiple packages, including the graphing package ggplot2 and the data wrangling package dplyr.
For this analysis, we’ll use the Major League Baseball data package. To get the data, we’ll rely on the CRAN version of the Lahman package, which will (at this writing; an update is pending) take us through the 2016 season.
library(tidyverse)

library(Lahman)
The code chunk below accesses the Teams table from the Lahman database, and wrangles it a bit, starting with filtering the series to only those years from 1961 (the start of the expansion era) to the most recent season in the data package.
The code then calculates and adds to the data table (though the dplyr::mutate function) two new variables: the team’s winning percentage, and using the winexp_fun function we wrote above, the win expectation.
data(Teams)

Teams_sel <- Teams %>%
  filter(yearID >= 1961) %>%
  rename(RS = R) %>%
  mutate(winpct = W / G, 
         winexp = winexp_fun(RS, RA))

plot

Now we’ll use ggplot2 to look at the relationship between the Pythagorean estimate of win expectation and the actual value of winning percentage. We can do this in a couple of ways: one is to overlay the density plots of the two variables, and the other is an X-Y scatterplot.
First the density plot.
plot_winexp_density <- ggplot(data = Teams_sel) +
  geom_density(aes(x = winexp, colour = "winexp"), show.legend = FALSE) +
  stat_density(aes(x = winpct, colour = "winpct"),
               geom = "line", position = "identity", size = 0) +
  guides(colour = guide_legend(override.aes=list(size=1)))


# tips from https://stackoverflow.com/questions/29563375/plotting-line-legend-for-two-density-curves-with-ggplot2

plot_winexp_density

In the above plot, we can see that there’s not a perfect match between the two lines. First, there are gaps between the two lines at either tail of the curve. But more prominently, the blue line (representing the actual winning percentage) isn’t a smooth curve at the top–there’s a hollow around .500, and increased proportions on either side. Something to investigate another day!
Next, the scatter plot. Because we are going to return to the foundations of this plot (i.e. the calculated win expectancy winexp as the X axis and the end-of-season final winning percentage winpct plotted on the Y axis), we’ll create a blank frame in an object called we_scatterplot. Once this object is created, we can build a variety of plots by simply overlaying different data representations. (It’s not lazy, it’s efficient.)
Note that there’s a few things going on here:
  • the use of the geom_blank function; usually, we would call geom_point for a scatter plot, but in this case we don’t want to see the data points.
  • the coord_fixed means that the X and Y scales have the units represented by equal length on both (one tenth of a point is the same length on both axes).
  • the scale_x_continuous function and its equivalent for Y set the grid marks and length of the two axes.
In this approach to plotting, the initial chunk of code creates the underlying framework using geom_blank(). The data is in the object but not rendered yet. This will work effectively for our project, since we are going to be plotting different representations of the single data set. The resulting we_scatterplot object contains the winexp and winpct data points, which we will summon by using different geom_ functions.
plot_we_scatterplot <- ggplot(data = Teams_sel, aes(x = winexp, y = winpct)) +
  geom_blank() +
  coord_fixed() +
  scale_x_continuous(breaks = seq(0.2, 0.8, by = 0.1),
                     limits = c(0.2, 0.8)) +
  scale_y_continuous(breaks = seq(0.2, 0.8, by = 0.1),
                     limits = c(0.2, 0.8)) 

plot_we_scatterplot

Now, we’ll render that object but add the geom_point so we can see the winexp and winpct values on an X-Y scatterplot.
plot_we_scatterplot_point <- plot_we_scatterplot +
  geom_point()
In the above plot it’s easy to see the strong relationship between the win expectation (the Pythagorean estimate, winexp on the X axis) and the winning percentage (winpct, on the Y).
(Yes, this looks a lot like the “Winning Percentage vs Run Differential” plot that appears in Jeff Sullivan’s July 3rd article. That’s because the winexp variable above is a permutation of run differential. Same values, different equation.)
To the above plot, let’s now add a red line showing where the win expectation, based on the Pythagorean equation, equals the winning percentage recorded by the team (i.e. where the value on the X axis equals the value on the Y.) The individual data points will be dialed back using a shade of grey (another option would be to use the alpha aesthetic to make the points somewhat transparent).
plot_we_scatterplot_line <-
plot_we_scatterplot +
  geom_point(colour = "grey75") +
  geom_segment(aes(x = 0.25, xend = 0.75, y = 0.25, yend = 0.75), colour = "red", size = 1.5)

plot_we_scatterplot_line

The individual data points above the red line are where teams have outperformed their win expectancy, and those below the line have failed to win as many games as the Pythagorean model would predict.
While in general the trend is clear, it’s not a perfect relationship. Over a 162 game season, there is still plenty of variation, with some teams above the line (that is, winning more games than the Pythagorean model would predict) and other teams losing more games than the model prediction appearing below the line.

Seattle Mariners, 2018

So how did the 2018 season end for the Seattle Mariners? Did they regress to the mean, or end up one of the clutchiest teams on record?
They ended with a record of 89 wins and 73 loses, a 0.549 record.
Mariner_winpct <- (89 / sum(89 + 73))

Mariner_winpct
## [1] 0.5493827
But on the run differential front, they allowed 34 more runs than they scored (677 scored vs. 711 allowed.) Let’s plug those numbers into the winexp_fun():
RS= 677
RA = 711

Mariner_winexp <- winexp_fun(RS, RA)

Mariner_winexp
## [1] 0.475519
Mariner_winexp * 162
## [1] 77.03408
The Mariners’ predicted winning percentage for the season, based on the Pythagorean model, is 0.475519, well below their final result. In terms of the number of games the Pythagorean model would predict they’d win in a 162 game season would be 77 … far fewer than the 89 wins they actually registered.
Finally, let’s add the point (0.475519 , 0.5493827) to our X-Y scatterplot:
plot_we_scatterplot_SM18 <-
plot_we_scatterplot_line +
  geom_point(x = Mariner_winexp, y = Mariner_winpct, size = 3, colour = "#005C5C")


plot_we_scatterplot_SM18

That Northwest Green (hex code #005C5C) dot well above the line? That’s the 2018 Seattle Mariners. They started the season over-performing relative to their run differential, and finished that way…virtually no regression to the mean.
In my next post, I’ll use linear regression–that workhorse of statistics, machine learning, artificial intelligence, econometrics, etc.–to look more deeply at the relationship between run differential and winning percentage. As part of this, I’ll use the broom package to dig into the regression model, and quantify the 2018 Seattle Mariners relative to other teams.

Further reading

FanGraphs “BaseRuns” page
Jay Heumann, “An improvement to the baseball statistic ‘Pythagorean Wins’”, Journal of Sports Analytics 2 (2016) 49–59
-30-

November 23, 2018

EARL conference recap: Seattle 2018

I had the pleasure of attending the EARL (Enterprise Applications of the R Language) Conference held in Seattle on 2018-11-07, and the honour of being one of the speakers. The EARL conferences occupy a unique niche in the R conference universe, bringing together the I-use-it-at-work contingent of the R community. The Seattle event was, from my perspective (I use R at work, and lead a team of data scientists that uses R) a fantastic conference. Full marks to the folks from Mango Solutions for organizing it!

The conference started with a keynote, “Text Mining with Tidy Data Principles”, from the always-brilliant Julia Silge. She’s an undisputed leader in the field of text analysis with R (the book she co-authored with David Robinson, Text Mining with R: A Tidy Approach, is already a cornerstone resource), and although I’d heard her deliver some of the same material at the Joint Statistical Meetings in July, this talk
  1. was longer and
  2. introduced some of her thinking about problems she’s tackling at her job at Stack Overflow.
It was fascinating to see where the utility of R as a text analysis tool is going, and Julia’s engaging manner, energy, and enthusiasm was a great start to the day.

Next up was a panel of leaders in the R community, on “Examining the future of R in industry”. The panelists were:
  • the aforementioned Julia Silge,
  • David Smith from Microsoft (he has the title “Cloud Developer Advocate–AI & Data Science”, but he’s also famous in the R community for his editorship and contributions to the Revolutions blog), and
  • Joe Cheng (the creator of Shiny and the CTO and Shiny team lead at RStudio).
With a trio of this calibre it was no surprise they had a wide-ranging and thoughtful discussion of the questions from the floor, covering everything from the pros and cons of different open source licenses to implementing R into production environments. The panel seemed, in my opinion, to land on a consensus that the future of R is bright, and that we will continue to see it remain specialized as a data science tool, and that we will continue to see integration with other tools.

The rest of the day was dedicated to the presentations, which covered a wide range of topics from modeling the relationship between roadway speed (from Joonbum Lee at Battelle Memorial Institute) and quantitative risk assessment at Starbucks (David Severski) to using deep learning on satellite images (Damian Rodziewicz of Appsilon Data Science). All of the speakers were engaging, had a great perspective on their topics, and only one (full disclosure: me) nattered on and didn’t leave any time for questions from the floor.

Intending no slight to the other speakers, three presentations really struck a nerve with me.
Eina Ooka from The Energy Authority spoke about her experience moving to R (and all of the benefits, from reproducibility to accuracy) in what she termed an “Excel-pervasive” environment. The space I work in is much the same; Excel is a workhorse for a lot of numeric analysis, and it is a go-to tool for many people in the clients we serve. Some of those clients expect delivery of their data tables in an Excel file. Eina’s success tackling the transition, in spite of the hurdles she faced, was inspiring.

Stephanie Kirmer from Uptake delivered what was, to me, perhaps the most immediately relevant talk: “The case for R packages as team collaboration tools”. I particularly liked the matrix showing the “Progression of Team Collaboration Infrastructure”, with version control, code sharing, and code storage and dissemination at four levels of sophistication. I was struck by how far my colleagues and I have to go to move up the ladder, but immediately recognized at least one project where a package would be an ideal way for us to start to collaborate more effectively.

And finally Aimee Gott from Mango Solutions, whose closing talk “Building a data science teams with R” was the perfect summary of everything that had preceded it. Again, it was a typology that stuck with me–in this case, types of R users, from the Super Users to the Cut & Paste Tweakers.
In short, the conference was a great way to hear from and meet R users who are finding applications for it in a business (or in my case, government) setting. Thanks again to Mango Solutions.

The 2018 EARL road show continued on to Houston (2018-11-09) and Boston (2018-11-13), each with different slates of speakers.

My only hope is that next year’s EARL road show makes a stop in Canada!

Note: looking for the slides and full narrative of my talk?
Bonus note: this post can be found in the B.C. Government GitHub repo dedicated to public presentations on the topic of R.


-30-

August 31, 2018

Smoke from a distant fire

Forest fires and air quality

August 31, 2018


It was recently announced that during 2018, British Columbia has seen the most extensive forest fire season on record. As I write this (2018-08-31) there are currently 442 wildfires burning in British Columbia. These fires have a significant impact on people’s lives–many areas are under evacuation order and evacuation alert, and there are reports that homes have been destroyed by the blazes.

The fires also create a significant amount of smoke, which has been pushed great distances by the shifting winds. This includes the large population centres of Vancouver and Victoria in British Columbia, as well as the Seattle metropolitan region and elsewhere in Washington. (Clifford Mass, Professor of Atmospheric Sciences at the Universtiy of Washington in Seattle, has written extensively about the smoke events in the region; see for example Western Washington Smoke: Darkest Before the Dawn from 2018-08-22.)

The Province of British Columbia has many air quality monitoring stations around the province, and makes the data available. The measure most used for monitoring the effects on health is PM25 or PM2.5, for fine particles with a diameter of 2.5 microns (millionths of a metre). The B.C. government has a Current Particulate Matter map that colour codes the one hour average measures for all the testing stations around the province.

The data file and a simple plot


The DataBC Catalogue provides access to air quality data. There’s “verified” to the end of 2017, and “unverified” for the past 30 days. Since we want to see what happened this month, it’s the latter we want. (The page with the links to the raw files is here.)

The files are arranged by particulate or gas type; there’s a table for ozone and another for sulpher dioxide, and others for the particulate matter. Note that the data are made available under the Province of B.C.’s Open Data license, and are in nice tidy form. And the date format is ISO 8601, which makes me happy.

To make sure we’ve got a reproducible version, I’ve saved the file I downloaded early this morning to my google drive. The link to the folder is here.

For the first plot, let’s look at the PM2.5 level for my hometown of Victoria, B.C. The code below loads the R packages we'll use, reads the data, and generates the plot.


# tidyverse packages
library(tidyverse)
library(glue)

PM25_data <- readr::read_csv("PM25_2018-08-31.csv")
filter(STATION_NAME == "Victoria Topaz") %>% ggplot() + geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) + labs(x = "date", title = glue("Air quality: Victoria Topaz"), subtitle = "one hour average, µg/m3 of PM2.5", caption = "data: B.C. Ministry of Environment and Climate Change Strategy")



There are 61 air quality monitoring stations around British Columbia. It would be interesting to see how the air quality was in other parts of the region–and since over half (54% in 2017) of the province’s population lives in the Vancouver Census Metropolitan Area (CMA), let’s plot the air quality there. There are multiple stations in the Vancouver CMA, so I chose the one at Burnaby South…it’s fairly central in the region.We can run this line of code to see a listing of all 61 stations (but we won’t do that now…)

# list all the air quality stations for which there is PM2.5 data
unique(PM25_data$STATION_NAME)

And since we’re going to be doing this often, let’s wrap the code that filters for the location we want and runs the plot in a function. Note that we’ll create a new variable station_name so all we need to do to change the plot is assign the name of the station we want, and off we go. Not only does this simplify our lives now, but is all-but-essential for a Shiny application.

# the air quality plot
PM25_plot <- function(datafile, station_name){
  datafile %>%
  filter(STATION_NAME == station_name) %>%
  ggplot() +
  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +
  labs(x = "date",
       title = glue("Air quality: ", station_name),
       subtitle = "one hour average, µg/m3 of PM2.5", 
       caption = "data: B.C. Ministry of Environment and Climate Change Strategy")
}

Now that we've got the function, the code to create the plot for the Burnaby South station is significantly simplified: assign the station name, and call the function.

# our Burnaby plot
station_name <- "Burnaby South"

PM25_plot(PM25_data, station_name)




And what about the towns that are the closest to the fires? While there are fires burning across the province, the fires that are burning the forests of the Nechako Plateau have understandably received a lot of attention. You may have seen the news stories and images from Prince George like this and this, or the images of the smoke plume from the NASA Worldview site.

Prince George is east of many major fires, downwind of the prevailing westerly winds. So what has the air quality in Prince George been like?

station_name <- "Prince George Plaza 400"
PM25_plot(PM25_data, station_name)




Or still closer to the fires, the town of Burns Lake.

station_name <- "Burns Lake Fire Centre"
PM25_plot(PM25_data, station_name)



The town of Smithers is west of the fires that are burning on the Nechako Plateau and producing all the smoke experienced in Burns Lake and Prince George. The residents of Smithers have had a very different experience, only seeing smoke in the sky when the winds shifted to become easterly.

station_name <- "Smithers St Josephs" 
PM25_plot(PM25_data, station_name)




multiple stations in one plot

You may have noticed that the Y axis on the plots can be quite different–for example, Victoria reaches 300, Smithers gets to 400, and Prince George is double that at 800, and Burns Lake is more than double again. There are two ways we can compare multiple stations: a single plot, or faceted plots.

a line plot with four stations

station_name <- c("Burns Lake Fire Centre", "Prince George Plaza 400", 
                  "Smithers St Josephs", "Victoria Topaz")

PM25_data %>%
  filter(STATION_NAME %in% station_name) %>%
  ggplot() +
  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE, colour = STATION_NAME)) +
  labs(x = "date",
       title = glue("Air quality: Burns Lake, Prince George, Smithers, Victoria"),
       subtitle = "one hour average, µg/m3 of PM2.5", 
       caption = "data: Ministry of Environment and Climate Change Strategy")




With four complex lines as we have here, it can be hard to discern which line is which. 

Use facets to plot the four stations separately


Facets give us another way to view the comparisons. In the first version, with the facets stacked vertically, it emphasizes comparisons on the X axis–that is, over time. In this way, we can see that the four locations have had smoke events that have occurred at different times.

station_name <- c("Burns Lake Fire Centre", "Prince George Plaza 400", 
                  "Smithers St Josephs", "Victoria Topaz")


PM25_data %>%
  filter(STATION_NAME %in% station_name) %>%
  ggplot() +
  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +
  facet_grid(STATION_NAME ~ .) +
  labs(title = glue("Air quality: Burns Lake, Prince George, Smithers, Victoria"),
       subtitle = "one hour average, µg/m3 of PM2.5", 
       caption = "data: B.C. Ministry of Environment and Climate Change Strategy") +
  theme(axis.text.x=element_text(size=rel(0.75), angle=90),
        axis.title = element_blank())



In the second version, the facets are placed horizontally, making comparisons on the Y axis clear. The smoke events in the four locations have been of very different magnitudes.
PM25_data %>%
  filter(STATION_NAME %in% station_name) %>%
  ggplot() +
  geom_line(aes(x = DATE_PST, y = REPORTED_VALUE)) +
  facet_grid(. ~ STATION_NAME) +
  labs(title = glue("Air quality: Burns Lake, Prince George, Smithers, Victoria"),
       subtitle = "one hour average, µg/m3 of PM2.5", 
       caption = "data: B.C. Ministry of Environment and Climate Change Strategy") +
  theme(axis.text.x=element_text(size=rel(0.75), angle=90),
        axis.title = element_blank())




These two plots show not only that Burns Lake and Prince George have had the most extreme smoke events, but that they have had sustained periods of poor air quality through the whole month. While the most extreme event in Smithers exceeds that of Victoria, there hasn’t been a prolonged period of smoke in the air like the other three locations.

-30-






April 3, 2017

Storytelling with Data: consumer alert

In March 2016 I gave a favourable review to a newly published book, Storytelling with Data by Cole Knaflic. Today (2017-04-03) she posted a new blog entry, "the book you're holding might be a fake!", in response to the discovery that poor-quality pirate editions of her books are available for purchase.

Knaflic's response has been exemplary: notifying people who bought the book what to look for in the knock-offs, and how to exchange the book for a proper copy. If you have a copy, go to the linked site above and check your copy against the description. (Fortunately, my copy is legit.)

This problem has alerted me to something I neglected to mention in my original review: the high quality of the physical book. The paper has a soft sheen which allows the print and graphics to stand out, the colours are sharp and consistent, and printing is clear.

A year later, my high opinion of this book has not shifted.

-30-

March 26, 2017

Updated Shiny app

A short post to alert the world that my modest Shiny application, showing Major League Baseball run scoring trends since 1901, has been updated to include the 2016 season. The application can be found here:
https://monkmanmh.shinyapps.io/MLBrunscoring_shiny/.

In addition to the underlying data, the update removed some of the processing that was happening inside the application, and put it into the pre-processing stage. This processing needs to happen only the once, and is not related to the reactivity of the application. This will improve the speed of the application; in addition to reducing the processing, it also shrinks the size of the data table loaded into the application.

The third set of changes were a consequence of the updates to the Shiny and ggplot2 packages in the two years that have passed since I built the app. In Shiny, there was a deprecation for "format" in the sliderInput widget. And in ggplot2, it was a change in the quotes around the "method" specification in stat_smooth(). A little thing that took a few minutes to debug! Next up will be some formatting changes, and a different approach to one of the visualizations.

 -30-

November 15, 2016

Subtitles and captions with ggplot2 v.2.2.0

Back in March 2016, I wrote about an extension to the R package ggplot2 that allowed subtitles to be added to charts. The process took a bit of fiddling and futzing, but now, with the release of ggplot2 version 2.2.0, it’s easy.

Let’s retrace the steps, and create a chart with a subtitle and a caption, the other nifty feature that has been added.

First, let’s read the packages we’ll be using, ggplot2 and the data carpentry package dplyr:

# package load 
library(ggplot2)
library(dplyr)

Read and summarize the data

For this example, we’ll use the baseball data package Lahman (bundling the Lahman database for R users), and the data table ‘Teams’ in it.

Once it’s loaded, the data are filtered and summarized using dplyr.
  • filter from 1901 [the establishment of the American League] to the most recent year,
  • filter out the Federal League
  • summarise the total number of runs scored, runs allowed, and games played
  • using `mutate`, calculate the league runs (leagueRPG) and runs allowed (leagueRAPG) per game
library(Lahman)
data(Teams)

MLB_RPG <- Teams %>%
  filter(yearID > 1900, lgID != "FL") %>%
  group_by(yearID) %>%
  summarise(R=sum(R), RA=sum(RA), G=sum(G)) %>%
  mutate(leagueRPG=R/G, leagueRAPG=RA/G)

A basic plot

You may have heard that run scoring in Major League Baseball has been down in recent years…but what better way to see if that’s true than by plotting the data?

For the first version of the plot, we’ll make a basic X-Y plot, where the X axis has the years and the Y axis has the average number of runs scored. With ggplot2, it’s easy to add a trend line (the geom_smooth option).

The scale_x_continuous options set the limits and breaks of the axes.

MLBRPGplot <- ggplot(MLB_RPG, aes(x=yearID, y=leagueRPG)) +
  geom_point() +
  geom_smooth(span = 0.25) +
  scale_x_continuous(breaks = seq(1900, 2015, by = 20)) +
  scale_y_continuous(limits = c(3, 6), breaks = seq(3, 6, by = 1))

MLBRPGplot




So now we have a nice looking dot plot showing the average number of runs scored per game for the years 1901-2015. (The data for the 2016 season, recently concluded, has not yet been added to the Lahman database.)

With the basic plot object now created, we can make the changes in the format.  In the past, the way we would set the title, along with X and Y axis labels, would be something like this.

MLBRPGplot +
  ggtitle("MLB run scoring, 1901-2014") +
  theme(plot.title = element_text(hjust=0, size=16)) +
  xlab("year") +
  ylab("team runs per game")


Adding a subtitle and a caption: the function

A popular feature of charts–particularly in magazines–is a subtitle that has a summary of what the chart shows and/or what the author wants to emphasize.

In this case, we could legitimately say something like any of the following:
  • The peak of run scoring in the 2000 season has been followed by a steady drop
  • Teams scored 20% fewer runs in 2015 than in 2000
  • Team run scoring has fallen to just over 4 runs per game from the 2000 peak of 5 runs
  • Run scoring has been falling for 15 years, reversing a 30 year upward trend
I like this last one, drawing attention not only to the recent decline but also the longer trend that started with the low-scoring environment of 1968.

How can we add a subtitle to our chart that does that, as well as a caption that acknowledges the source of the data? The new labs function, available in ggplot2 version 2.2.0, lets us do that.

Note that labs contains the title, subtitle, caption, as well as the X and Y axis labels.

MLBRPGplot +
  labs(title = "MLB run scoring, 1901-2015",
       subtitle = "Run scoring has been falling for 15 years, reversing a 30 year upward trend",
       caption = "Source: the Lahman baseball database", 
       x = "year", y = "team runs per game") 




Easy.

Thanks to everyone involved with ggplot2 who made this possible.

The code for this post (as an R markdown file) can be found in my Bayesball github repo.


-30-