Data Visualization in R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Data Visualization in R\n",
    "## by Diya Das\n",
    "\n",
    "### The goal\n",
    "Data visualization, like all visualization, is important because it tells a story. Take a moment to think about what stories you'd like to tell with your data. \n",
    "\n",
    "There's an important component of data visualization - deciding what kind of plot is appropriate to make a particular point - that we're not going to talk about today. In fact, today, we're going to make some visualizations that I consider to be pretty bad examples of data visualization, just to show you what is *technically* possible. However, as always, just because you *can* do something, it doesn't mean that you *should* do it. If you're interested in thinking about principles of good visualization, I encourage you to research that topic further. You could start with Edward Tufte.\n",
    "\n",
    "I'm going to be giving examples of how to make plots in R, using both base R and the occasional `ggplot` command, as well as a 3D plotting library called `rgl`. (Why? because *why should you not gif your data?* Actually, it's because it's sometimes easier to look at relationships among high-dimensional data in 3D, vs. 2D.)\n",
    "\n",
    "### The datasets\n",
    "We are going to be using the data from the R package [`nycflights13`](https://cran.r-project.org/web/packages/nycflights13/nycflights13.pdf). There are five datasets corresponding to flights departing NYC in 2013. We will **load directly into R from the library**, but the repository also includes TSV files we created for the purposes of the Python demo and can also be used to load the data into our R session."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "*** If you've never run Jupyter notebooks with R, please run `conda install -c r r-essentials`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repos=structure(c(CRAN=\"http://cran.cnr.berkeley.edu/\", \n",
    "BioCsoft=\"http://www.bioconductor.org/packages/release/bioc/\")))\n",
    "ipak <- function(pkg){\n",
    "     new.pkg <- pkg[!(pkg %in% installed.packages()[, \"Package\"])]\n",
    "     if (length(new.pkg))\n",
    "         install.packages(new.pkg, dependencies = TRUE)\n",
    "     sapply(pkg, require, character.only = TRUE)\n",
    " } #https://gist.github.com/stevenworthington/3178163\n",
    "\n",
    "pkgs <- c(\"nycflights13\", \"ggplot2\", \"rgl\", \"NMF\", \"dplyr\",\"tidyr\")\n",
    "ipak(pkgs)\n",
    "options(jupyter.plot_mimetypes = 'image/png',repr.plot.width=6, repr.plot.height=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#invisible(sapply(pkgs, library, character.only=TRUE )) \n",
    "    # if you knew you had installed the packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic scatterplots\n",
    "Let's run through an example using the `flights` dataset. This dataset includes...well what does it include? You could read the documentation, but let's take a look first. \n",
    "\n",
    "What are the dimensions of the flights data frame?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flights <- data.frame(flights) # dplyr has introduced a new data format that I am ignoring\n",
    "flights <- flights[complete.cases(flights),] # remove NAs to reduce some issues\n",
    "dim(flights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure out the categories of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "summary(flights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's just look at 1000 random flights in January."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "set.seed(5106)\n",
    "flights <- flights[flights$month==1,]\n",
    "flights <- flights[sample(1:nrow(flights),1000),]\n",
    "dim(flights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot(flights$dep_delay, flights$arr_delay)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing aesthetic parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if you prefer:\n",
    "- filled circles instead of hollow: adjust `pch`\n",
    "- different color points: adjust `col`\n",
    "- different sized points: adjust `cex`\n",
    "- different titles and axis labels: adjust `main` (title), `xlab`, `ylab`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot(flights$dep_delay, flights$arr_delay, cex=2, pch=19, \n",
    "     col=\"red\", main= \"Arrival vs. Departure Delay\", \n",
    "     xlab= \"Departure delay (minutes)\", ylab=\"Arrival delay (minutes)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coloring by one variable and adding a legend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Some tidying for plotting things\n",
    "flights$origin <- factor(flights$origin) # convert to factor to color by origin in plot \n",
    "colpal <- scales::alpha(c(\"red\",\"green\",\"blue\"),0.3) # how to make transparent colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4)\n",
    "plot(flights$dep_delay, flights$arr_delay, cex=0.5, pch=19,\n",
    "     col=colpal[flights$origin], main= \"Arrival vs. Departure Delay\")\n",
    "legend(\"bottomright\",levels(flights$origin),fill=colpal,\n",
    "       cex=0.75, y.intersp=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=3)\n",
    "hist(flights$dep_time, main=\"Flights by time of departure\", \n",
    "     xlab=\"Time\", ylab = \"# of flights in Jan 2013\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Change number of divisions in a histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hist(flights$dep_time, main=\"Flights by time of departure\", \n",
    "     xlab=\"Time\", ylab = \"# of flights in Jan 2013\",breaks=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot proportions, not frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hist(flights$dep_time, main=\"Flights by time of departure\", \n",
    "     xlab=\"Time\", ylab = \"Proportion of flights in Jan 2013\",breaks=24, freq=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot organization and margins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if you want to have two plots side by side, or change the margins?\n",
    "In R, this is controlled by `par`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4)\n",
    "par(mfrow=c(1,2)) # 1 row, 2 columns\n",
    "\n",
    "plot(flights$dep_delay, flights$arr_delay, cex=0.5, pch=19,\n",
    "     col=colpal[flights$origin], main= \"Arrival vs. Departure Delay\")\n",
    "legend(\"bottomright\",levels(flights$origin),fill=colpal,\n",
    "       cex=0.75, y.intersp=2, xpd=FALSE)\n",
    "\n",
    "hist(flights$dep_time, main=\"Flights by time of departure\", \n",
    "     xlab=\"Time\", ylab = \"# of flights in 2013\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "par.defaults <- par() # save default plot parameters before we mess around\n",
    "par.defaults"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dot plots with jitter\n",
    "## Turning off axes and labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=6)\n",
    "par(mfrow=c(2,1))\n",
    "hist(flights$dep_time, main=\"Flights by time of departure\", \n",
    "     xlab=\"Time\", ylab = \"# of flights in Jan 2013\", col=\"purple\")\n",
    "\n",
    "par(mar=c(5.1,1,1,1)) # change margins for second plot\n",
    "plot(flights$dep_time, jitter(rep(0, nrow(flights)), 0.3),\n",
    "     main=\"Flights by time of departure\", col=colpal[flights$origin],yaxt='n', ylab='',\n",
    "    xlab = \"Time\")\n",
    "\n",
    "\n",
    "# now reset to default margins and plot layout (1 x 1)\n",
    "par(par.defaults)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boxplot numerical variable by categorical variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4)\n",
    "boxplot(dep_time ~ origin, data=flights, col=colpal, ylab='Departure time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting multiple series\n",
    "### Base R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "delay_time <- flights %>% group_by(origin,time_hour) %>% \n",
    "    summarise(avg_dep_delay=mean(dep_delay), avg_arr_delay=mean(arr_delay))\n",
    "print(head(delay_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with(subset(delay_time, origin=='EWR'), plot(time_hour, avg_dep_delay, type=\"l\",col=colpal[1],\n",
    "                                            main=\"Average departure delay by airport\",\n",
    "                                            ylab = \"Departure delay\", xlab=\"Time\"))\n",
    "with(subset(delay_time, origin=='JFK'), lines(time_hour, avg_dep_delay,col=colpal[2]))\n",
    "with(subset(delay_time, origin=='LGA'), lines(time_hour, avg_dep_delay,col=colpal[3]))\n",
    "\n",
    "legend(\"top\",levels(flights$origin),fill=colpal,\n",
    "       cex=0.75, y.intersp=2, x.intersp=2, horiz=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `ggplot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ggplot(delay_time,aes(x=time_hour,y=avg_dep_delay,group=origin)) + geom_line(aes(colour = origin)) +\n",
    "ggtitle('Average departure delay vs. Time') +\n",
    "xlab('Time') + ylab('Arrival departure delay (minutes)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Bubble plots and diving into `ggplot`\n",
    "Bubble plots are useful for plotting three variables at once, where the third variable is some numerical quantity.\n",
    "\n",
    "Let's plot flights by origin and carrier. First, we need to construct a long-formatted data frame for input to `ggplot` (see Data Tidying and Manipulation notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "origin_carrier <- reshape2::melt(\n",
    "    prop.table(table(Origin=factor(flights$origin), \n",
    "                     Carrier=factor(flights$carrier)), 1), \n",
    "               value.name=\"Proportion\")\n",
    "\n",
    "head(origin_carrier)\n",
    "rowSums(prop.table(table(Origin=factor(flights$origin), \n",
    "                     Carrier=factor(flights$carrier)), 1))\n",
    "# Each carrier as a fraction of flights out of that airport"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `ggplot`: A conceptual introduction\n",
    "\n",
    "`ggplot` is based on the **grammar of graphics** - that is, a component-based approach to building graphics:\n",
    "\n",
    "- data\n",
    "- coordinate system\n",
    "- geoms: visual representation of data\n",
    "    - aesthetics: size, color, x, y locations\n",
    "\n",
    "The template is as follows, from the very helpful cheatsheet from RStudio:\n",
    "\n",
    "`ggplot(data = <DATA>) + \n",
    "  <GEOM_FUNCTION>(mapping=aes(<MAPPINGS>),stat= <STAT>,position=<POSITION>) +\n",
    "  <COORDINATE_FUNCTION> + <FACET_FUNCTION> + <SCALE_FUNCTION> + <THEME_FUNCTION>`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the base plot, describing what data are to be plotted. Note the size statement, with the `ifelse`: this is to alter the default behavior of `ggplot` which plots zeros as a tiny point, to distinguish them from `NA`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "or_car.plt <- ggplot(origin_carrier, \n",
    "                     aes(Carrier, Origin, size=ifelse(Proportion==0, NA, Proportion)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for plotting:\n",
    "\n",
    "- `geom_point` is a `GEOM_FUNCTION` that tells `ggplot` we want to plot points.\n",
    "- `scale_size_area` is a `SCALE_FUNCTION`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pdf(\"bubble.pdf\", width=6, height=4)\n",
    "plt2 <- or_car.plt + geom_point(color = \"blue\") + \n",
    "   scale_size_area(max_size=20, name='Proportion', breaks=c(0.1, 0.25, 0.5))\n",
    "plt2\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=3)\n",
    "theme1 <- theme(plot.background=element_blank(), panel.grid.minor=element_blank(), \n",
    "            panel.border=element_blank(), panel.background=element_blank(), \n",
    "            axis.line=element_blank(),axis.ticks=element_blank(), \n",
    "            legend.title=element_blank(), legend.background=element_blank(), \n",
    "            axis.text.x=element_blank(), axis.text.y=element_blank(),\n",
    "            legend.key= element_rect(fill=\"white\"))\n",
    "plt2 + theme1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "theme2 <- theme(plot.background=element_blank(), \n",
    "            panel.background=element_blank(), \n",
    "            legend.background=element_blank(), \n",
    "            legend.key= element_rect(fill=\"white\"))\n",
    "plt2 + theme2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stacked bar plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ggplot(origin_carrier, aes(x = Carrier, y = Proportion*100, fill = factor(Origin))) +\n",
    "  geom_bar(stat = \"identity\") + \n",
    "  scale_fill_manual(values=colpal, name=\"Origin\") + theme_bw() + \n",
    "  ggtitle('Proportion of flights from each airport by carrier') +\n",
    "  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),\n",
    "      panel.border=element_blank(),\n",
    "      axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +ylab(\"Percent\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More with `ggplot`: Going back over some of the graphics we've already made in base R\n",
    "\n",
    "### Basic scatterplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delay.plt <- ggplot(data=flights, aes(dep_delay, arr_delay,colour=origin)) \n",
    "delay.plt + geom_point(size=0.5) + \n",
    "ggtitle('Arrival delay vs. Departure delay') + theme(plot.title = element_text(hjust = 0.5)) +\n",
    "theme_bw() +\n",
    "xlab('Departure delay (minutes)') + ylab('Arrival delay (minutes)') +\n",
    "scale_colour_manual(values=colpal)+\n",
    "geom_smooth()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ggplot(data=flights, aes(dep_time)) + geom_histogram(binwidth=60) + \n",
    "  ggtitle(\"Flights by time of departure\")+ xlab(\"Time\") + ylab (\"# of flights in Jan 2013\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dot plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ggplot(data=flights, aes(dep_time, rep(0, nrow(flights))))+ geom_jitter(aes(colour = origin)) +\n",
    "  ggtitle('Flights by time of departure') + xlab('Time') +\n",
    "  scale_y_continuous(breaks=NULL) + ylab('')\n",
    "suppressWarnings(par(par.defaults))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Boxplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ggplot(data=flights, aes(origin,dep_time)) + geom_boxplot(fill=colpal) + ylab('Departure time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Faceting with `ggplot2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=3)\n",
    "delay.plt + geom_point() + geom_smooth() + facet_grid(. ~ origin) +\n",
    "ylab('Arrival delay') + xlab('Departure delay')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heatmaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "airtime <- left_join(flights, airports, by = c(\"dest\" = \"faa\")) %>% \n",
    "    select(origin, dest=name, air_time) %>% \n",
    "    group_by(origin, dest) %>% \n",
    "    summarize(avg_air_time = mean(air_time))\n",
    "avg_air <- data.frame(drop_na(spread(airtime, origin, avg_air_time)))\n",
    "rownames(avg_air) <- avg_air$dest\n",
    "head(avg_air)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=5)\n",
    "terminals <- c(\"3\",\"6\",\"4\") # Number of current terminals, as per Wikipedia\n",
    "aheatmap(avg_air[,2:ncol(avg_air)], annCol=data.frame(num_terminals=terminals), cexCol=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot3d(flights$distance, flights$air_time, flights$arr_delay)\n",
    "par3d(windowRect = c(20, 30, 800, 800))\n",
    "print(getwd())\n",
    "movie3d(spin3d(), 5, movie = \"movie\", frames = \"movie\", dir = getwd(), \n",
    "                    convert = TRUE, clean = TRUE, verbose=TRUE,\n",
    "                    top = TRUE, type = \"gif\") \n",
    "rgl.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pairs plots: 2D representations of high-dimensional data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pairs(flights[,c(\"distance\",\"air_time\", \"arr_delay\")], col=colpal[flights$origin])\n",
    "legend(\"bottomright\",levels(flights$origin),fill=colpal,cex=0.5, xpd=TRUE)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}