--- title: Exploratory Data Analysis / Plotting date: 3 September 2020 (Lecture 2) author: 36-467/667, Fall 2020 output: slidy_presentation bibliography: locusts.bib --- ```{r, include=FALSE} # General set-up options library(knitr) opts_chunk$set(size="small",background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, tidy=TRUE, warning=FALSE, message=FALSE, echo=FALSE) ``` ## House-keeping - Signed up for Canvas and Piazza? - Read the policies and Turabian? Done HW 0? - Downloaded HW 1 and the data set? ## Kinds of spatio-temporal data {.smaller} 1. Data are some kind of measurement for regions of space or intervals of time ("region" data, "areal" data) + Often called "time series" for temporal data + Often equally-sized time intervals or areas (every day or week, every square mile) + But not always (every county) 2. Data are the specific coordinates where some event occurs ("point" data) + E.g., times and locations of earthquakes + E.g., times at which a nerve cell "spiked" + E.g., locations of trees + E.g., times and locations of (reported) crimes + Point data can come with "marks", qualitative or quantitative: * Magnitude of earthquake * Species of tree * Type of crime * (Spikes are pretty much identical so doesn't make sense to mark them) ## Plotting over time: Time series - Use a scatterplot: time on horizontal axis, value on vertical (always) + Often a good idea to draw lines connecting values * Be careful when there are NAs! - Worry about the vertical scale + Are you exggerating small changes? + Are you downplaying meaningful changes? + "Adjust the vertical scale to the range of the data" and "Make sure the vertical scale includes zero" are both good rules of thumb but can conflict + $\therefore$ You need to think - Often nice: rug plot on vertical axis to show distribution + Sometimes nice to add to horizontal to show where data is present ## Back to Kyoto {.smaller} ```{r} kyoto <- read.csv("http://www.stat.cmu.edu/~cshalizi/dst/20/data/kyoto-2020.csv") plot(Flowering.DOY ~ Year.AD, data=kyoto, type="l", ylab="Day in year of full flowering", xlab="Year (CE)", main="Cherry blossoms in Kyoto") ``` - Vertical axis doesn't go to 0 + cherry trees have never flowered in Kyoto on January 1st, and never will) - Used `type="l"` so it _only_ draws connecting lines, with gaps for NAs + No marking for points with NAs on both sides! ## Back to Kyoto {.smaller} ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="o", cex=0.2, ylab="Day in year of full flowering", xlab="Year (CE)", main="Cherry blossoms in Kyoto") rug(x=kyoto$Flowering.DOY, side=2, col="pink") ``` - Now we see all the data values, even the isolated ones - Many repeated values, so rug gives more a sense of the _range_ than of the _distribution_ + Add a little visual noise through `jitter()` ## Back to Kyoto ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="o", cex=0.2, ylab="Day in year of full flowering", xlab="Year (CE)", main="Cherry blossoms in Kyoto") rug(x=jitter(kyoto$Flowering.DOY), side=2, col="pink") ``` ## Plotting over time: point data / events - Use a timeline + Spread out the events + If marked, use shape, color or size to indicate the mark + Continuous color and size better for quantitative marks than shape because there's a more intuitive mapping ## Plotting over time: events Earthquakes, in a rectangle that's (roughly) Japan, with a magnitude of at least 5.5 on the Richter scale, since January 1 2000 ```{r} quakes <- read.csv("https://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&minlatitude=31.5&maxlatitude=44.5&minlongitude=130.5&maxlongitude=145&starttime=2000-01-01&minmagnitude=5.5") # R needs some help realizing the first column is dates quakes$time <- as.Date(quakes$time) plot(x=quakes$time, y=rep(0, nrow(quakes)), type="o", ylab="", yaxt="n", xlab="Date", main="Earthquakes around Japan, 2000--present", sub="Magnitude 5.0 or bigger") ``` (This has a lot of wasted vertical space!) ## Plotting over time: events - Alternative: plot the _cumulative_ number of events from the beginning of the data up to time $t$ + Sometimes called the **counting process** ```{r} plot(x=quakes$time, y=sapply(quakes$time, function(t) { sum(quakes$time <= t) }), type="l", xlab="time", ylab="Cumulative number of earthquakes", main="Earthquakes around Japan, 2000--present") ``` - There's a jump at each event + There's a big jump if many events happen at once ## Plotting over time: tricks: index numbers - Have measurements $X(t_1), X(t_2), \ldots X(t_n)$ - Pick _one_ time, $t^*$, and then plot $X(t_i)/X(t^*)$ + Often $t^*=t_1$ but not necessarily - These ratios are called **index numbers**[^indexnumbers] + Often multiplied by 100 - Lets you avoid thinking about the actual units and look at just magnitude relative to some baseline + Often useful when the actual magnitude isn't illuminating _or_ we want to make comparisons across different time series [^indexnumbers]: To be precise, this is the simplest kind of index number. When measuring the economy, we often want to come up with a sort of summary of multiple pieces of information which will give us an over-all sense of change of time, and these are also called "index numbers". For instance, if we want to measure inflation, we have to look at the prices of, in principle, all the different sorts of goods and services that people buy, and how the prices of (comparable) goods change over time. But different goods will have prices changing at different rates, so we need to some how combine that into _one_ number. The "consumer price index" (CPI), which is what people usually mean when they talk about inflation, is a weighted average of the index numbers for a "basket" of goods and services, with the weights reflecting how consumers spent their money at the start of the time period. ("Chained" indices adjust the weights every so often.) Other choices of weights would give different indices --- there are actually _two_ CPIs, one for urban consumers and one rural consumers; a "producer price index" geared to what businesses rather than households spend money on; and so on. Sometimes, when econometricians talk about the "index number problem", they mean the problem of how to pick the weights; sometimes the mean the fact that no set of weights is ideal for all purposes. ## Plotting over time: tricks: indexing ![](https://live.staticflickr.com/3131/2295650288_a8f7a9a51b_c.jpg) ["Yukon Lynx Family" by David Cartier, Sr.](https://www.flickr.com/photos/dcartiersr/2295650288/) ## Plotting over time: tricks: indexing ```{r, echo=TRUE} library(datasets) data(lynx) ``` - A famous data set showing the _number_ of lynxes trapped each year in the Mackenzie River region of Canada - Ecologists don't actually care about this number, but about how many lynxes there were - Assume total number of lynxes $\propto$ number of trapped lynxes + But we don't know the proportion! + (Yes, people have looked in to whether the assumption makes sense here) - Use an index number, basing it off say the first year in the data: ## Plotting over time: tricks: indexing ```{r, echo=TRUE} lynx.index.numbers <- 100*lynx/lynx[1] plot(lynx.index.numbers, main="Relative abundance of lynxes", ylab="Index (1821==100)") ``` ## Plotting over time: tricks: differencing - Rate of change is \[ \Delta X(t_i) = \frac{X(t_i) - X(t_{i-1})}{t_i - t_{i-1}} \] - The numerator is the **first difference** of the series + What'd the **second difference** be? - If observations are equally spaced in time, can just look at $X(t) - X(t-1)$ - This is like looking at the derivative, $\frac{dX}{dt}$ - In R, we take differences with `diff()` (and it gives us somethign 1 shorter than the original vector) ## Differencing - Start with the population of the US over time: ```{r} library(pdfetch) # Total US population in 1000s pop.fred <- pdfetch_FRED("POP") # This comes as a complicated data type, so break it down to a simple data # frame, converting dates to a year with a decimal fraction library(xts) # Functions needed for pdfetch's preferred format library(lubridate) # Provides useful date-conversion functions pop <- data.frame(year=decimal_date(index(pop.fred)), y=1000*as.numeric(pop.fred)) plot(pop, type="l",ylab="US population", ylim=c(0, max(pop))) ``` ## Differencing: {.smaller} ```{r, echo=TRUE} nrow(pop) delta.pop <- diff(pop$y) length(delta.pop) pop$rate.of.change <- 12*c(delta.pop, NA) ``` (Why `12`? Why pad with the `NA`?) ## Differencing: {.smaller} ```{r, echo=TRUE} plot(rate.of.change ~ year, data=pop, type="l", xlab="year", main="US population rate of change", ylab="people/yr") ``` What do you think happened in 2010? ## Plotting over time: tricks: accumulating - Go the other way: add up rates of change - Like taking the integral - R: `cumsum()` (for `cum`ulative `sum`mation) ```{r} plot(cumsum(rate.of.change/12) ~ year, data=pop, type="l", ylab="Cumulative change in population") ``` ## Plotting over time: tricks: relative time - Instead of making one plot with one line over time, make multiple lines where we re-start when something happens, and plot time _relative to_ those events + Times at which those events happen are sometimes called "epochs" ## Relative time - Make the epoch-defining events earthquakes with magnitude at least 7 - Look $\pm 30$ days on either side ```{r} # Find all the big quakes, both in the data set and in time big.quake.rows <- which(quakes$mag > 7) big.quake.dates <- with(quakes, time[big.quake.rows]) # Make an empty plot plot(0, xlim=c(-30, 30), ylim=range(big.quake.dates), xlab="Days relative to first big quake", ylab="", yaxt="n", type="n", main="Earthquakes before and after big earthquakes", sub="Big = magnitude 7 or more") # Put nicely-formated values on the vertical axis axis(side=2, at=big.quake.dates, labels=as.character(big.quake.dates), las=2, cex.axis=0.5) # Now look for 30 days on either side of a 7.0 quake, and plot the time # _relative to_ this big quake for (big.one in big.quake.dates) { abline(h=big.one, col="grey") early <- big.one - 30 late <- big.one + 30 relevant.rows <- with(quakes, which(time >= early & time < late)) points(x=quakes$time[relevant.rows]-big.one, y=rep(big.one, length(relevant.rows))) } ``` ## Plotting over time: tricks: cycle time - Like relative time, but the "event" happens on a regular cycle, say annually ## Cycle time ```{r} plot(rate.of.change ~ year, data=pop, type="l", xlab="year", main="US population rate of change", ylab="people/yr") ``` ## Cycle time ```{r} # Make an empty plot plot(0, xlim=c(1, 12), ylim=range(na.omit(pop$rate.of.change)), type="n", ylab="Rate of change (people/yr)", xlab="Month", xaxt="n") axis(side=1, at=1:12, label=month.name, las=2, cex.axis=0.5) for (yr in 1:floor(nrow(pop)/12)) { rows <- (1:12) + (yr-1)*12 lines(x=1:12, y=pop$rate.of.change[rows], lwd=0.2) } ``` What's the difference between winter and summer? ## Dynamics: scatterplots - Plot $X(t+1)$ vs $X(t)$ ## Dynamics ```{r} plot(x=head(lynx.index.numbers, -1), y=tail(lynx.index.numbers, -1), xlab="Relative abundance of lynxes this year", ylab="Relative abundance of lynxes next year", sub="Indexed so 1821==100") ``` - Of course you might also try plotting $X(t+1)$ vs $((X(t), X(t-1))$... ```{r} library(scatterplot3d) n <- length(lynx.index.numbers) scatterplot3d(x=lynx.index.numbers[1:(n-2)], y=lynx.index.numbers[2:(n-1)], z=lynx.index.numbers[3:n], xlab="X(t-2)", ylab="X(t-1)", zlab="X(t)", main="Dynamics of lynx relative abundance", sub="scaled so 1821 abundance = 100", type="l") ``` ## Relationships between two variables {.smaller} - You often see plots with two variables with different vertical scales ```{r, fig.retina=NULL, out.width=600, echo=FALSE} knitr::include_graphics("cage-vs-drowning.jpg") ``` [(Spurious Correlations)](https://www.tylervigen.com/spurious-correlations) - This can be OK but it's often not a great idea, because you can make the relationship look very different by changing the axes + Imagine having the red vertical axis go down to 0! - Usually a better idea to use a scatterplot of one variable against another + More obvious if you're pulling a fast one with the axes + Use lines to connect successive measurements * Or use numbers as plotting symbols, etc. ## Relationships between two variables {.smaller} - There are also issues when each variable is correlated _with itself_ + Say $X(t)$ is $+$vely correlated with $X(t+1)$ and so are $Y(t)$ and $Y(t+1)$ + Then if $X(t)$ is a above average/below, we expect $X(t+1)$ to also be above/below average + The same for $Y(t)$ and $Y(t+1)$ + $\therefore$ knowing if $X(t)$ is above average helps us guess if $Y(t)$ is above average - We will come back to dealing with _auto_-correlation later in the course ## Plotting over space: maps - Coordinates of the plot match the coordinates in space - For point data: end of story + Unless there are marks; use shape, color or size - For region data: you need to indicate the value of $X(r)$ + Put a point at the coordinate: use shape, color or size + Show the whole region: use color * (or lines for fill-in; more old-fashioned) * The jargon name for a map that shades each region to represent a value is a **choropleth** (Greek _choros_ "region" + _plethos_ "many") (I tend to mis-spell it as ch_l_oropleth...) ## Plotting over space: maps: complications - If you have data on a regular grid, especially a _rectangular_ grid, R has many nice tools for making your plots + E.g., `image()` - If you have data for irregularly-shaped areas (counties, states, countries, ...), then you need: 1. A computer-friendly way of naming each area + For the US, the are standardized "FIPS" codes giving numerical identifiers for states and counties 2. A definition for the shape and location of each area, usually a **polygon** + The `maps` package may offer what you need 3. A good way to fill in each polygon with the right color or texture from the data set - There are very elaborate geographic information system (GIS) software packages for doing stuff like this - For straightforward examples I recommend looking at @Healy-on-visualization, ch. 7 [(online)](https://socviz.co/maps.html) ## Mapping earthquakes ```{r} library(maps) map("world", xlim=range(quakes$longitude), ylim=range(quakes$latitude)) points(latitude ~ longitude, data=quakes, pch=16) ``` Exercise: Modify color or size of the points to indicate magnitude ## Plotting over space: maps: projections {.smaller} - The Earth is round - The page (or screen) is flat - _No_ coordinate system for a sphere goes on to a flat page without distortion and discontinuity[^geometry] - Not a big deal on a small[^curvature] scale - On a large scale, you need to pick how to "project" the round Earth on to the flat page - Different projection schemes all have different virtues/drawbacks + Using longitude and latitude as $x$ and $y$ is easy and familiar + Other projections preserve area + Yet other projections preserve orientation * the "Mercator" projection that makes Greenland look bigger than Africa (unless you flip the north and south poles) + _Nothing_ preserves distance - If you want to calculate distance between points using latitude and longitude, you've got to do some trigonometry + Many R packages have functions to do the trig for you [^geometry]: This is actually a kind of profound mathematical fact. Two sets of points $A$ and $B$ have the same dimension if, and only if, we can find a continuous function $f$ from $A$ to $B$ with a continuous inverse $f^{-1}$. The sphere is three dimensional and the page is two dimensional, so we can't. If you want to follow up on this thought, and be convinced of the theorem, I strongly recommend @McCleary-topology-first and @Schutz-geometrical (which isn't just for physicists). [^curvature]: "Small" here means "small compared to the curvature of the Earth". Pretend that the Earth is perfectly spherical. (It isn't; working out exactly how it departs from being a sphere ["geodesy"] was very important to the development of statistics [@Stigler-history].) Pick two points on the surface of the Earth, $A$ and $B$, and now imagine drawing lines from the center of the Earth $C$ to those points. The angle $\angle ACB$ has a certain magnitude, call it $\theta$ in radians. (Remember $2\pi$ radians = $360$ degrees.) The distance between $A$ and $B$ if we measure _along_ the surface of the Earth ("as the crow flies") is $r\theta$, where $r = \overline{AC} = \overline{BC}$ is the radius of the Earth. But the distance between $A$ and $B$ on a straight line, $\overline{AB}$, which would cut through the surface of the Earth, is $2 r \sin{(\theta/2)}$. (Find the mid-point between $A$ and $B$ on the line between them, say $D$. Draw the triangles $\triangle ACD$ and $\triangle BCD$. These are both right triangles with hypotenuse $r$, and the angles $\angle ACD$ and $\angle BCD$ are both $\theta/2$, so the sides opposite those angles ($\overline{AD}$ and $\overline{BD}$) have length $r \sin{(\theta/2)}$. For small angles $x$, $\sin{x} \approx x$, so for small angles, the distance along the surface of the Earth, $r\theta$, is approximately the same as the straight-line distance, $\approx 2r \theta/2 = r\theta$, and we don't really care that the Earth is round. When it starts to matter depends on how much inaccuracy we're willing to tolerate. If we're OK with a 1% error distances, we can go out to a $\theta$ where $2\sin{(\theta/2)}/\theta = 0.99$, which numerically is about `r signif(uniroot(function(theta) { 2*sin(theta/2)/theta - 0.99}, lower=1e-6, upper=1)$root,4)` radians. Since the radius of the Earth is 6378 kilometers, this corresponds to a surface distance of `r signif(6378 * uniroot(function(theta) { 2*sin(theta/2)/theta - 0.99}, lower=1e-6, upper=1)$root,4)` kilometers. But, as the proverb says, getting 99% of the way across the ocean still means drowning; if we want an accuracy of $0.01$%, then we're looking at an angle of only `r signif(uniroot(function(theta) { 2*sin(theta/2)/theta - 0.9999}, lower=1e-6, upper=1)$root,4)` radians, and so only `r signif(6378 * uniroot(function(theta) { 2*sin(theta/2)/theta - 0.9999}, lower=1e-6, upper=1)$root,4)` kilometers along the surface. ## Plotting over space: relationships between variables - You'll often see side-by-side maps - Usually not a great way to establish a relationship + Partly because of the multiple-scales issue (as with time series) + Partly because of auto-correlation issues (as with time series) + Partly because you're often just showing they're both tracking population (or whatever) - Recommendation: make a scatter plot + And consider regressing-out things like population, or at least normalizing by population + Again, we'll come back to how to deal with auto-correlation ## Spatial relationships ```{r, fig.retina=NULL, out.height=450, echo=FALSE} knitr::include_graphics("https://imgs.xkcd.com/comics/heatmap_2x.png") ``` [XKCD no. 1138, "Heatmaps"](https://xkcd.com/1138/) ## Spatial Relationships ```{r, fig.retina=NULL, out.height=450, echo=FALSE} knitr::include_graphics("https://kieranhealy.org/files/misc/us-pop-density-2014.png") knitr::include_graphics("https://kieranhealy.org/files/misc/us-pct-black-2013.png") ``` [Kieran Healy, "America's Ur-Choropleths"](https://kieranhealy.org/blog/archives/2015/06/12/americas-ur-choropleths/) ## Spatio-temporal data - Make a 3D plot, one axis being time + Often useful for point data + _Or_ when there's some sort of "background" you can make invisible + Not so good for coloring every area over time, usually - Make a movie + Lots of packages for animations - Make lots of side-by-side maps + Arrange them in order of time! - Make lots of side-by-side time series plots + Might try putting each plot in the right spatial location but that can get crowded, especially if areas have different sizes ## Other forms of EDA - Everything we knew how to do with non-spatial non-temporal data is still valid: boxplots, histograms, tables, scatterplots showing relations between variables - Often useful to sub-divide into different regions of time or space, and do the EDA for each one ## Conventional EDA, Sub-Divided ```{r} par(mfrow=c(1,2)) with(kyoto, plot(hist(Flowering.DOY[Year.AD <= 1750], plot=FALSE), main="Day of full bloom, up to 1750")) with(kyoto, plot(hist(Flowering.DOY[Year.AD > 1750], plot=FALSE), main="Day of full bloom, after 1750")) ``` (_You_ can persuade R to show both of these in one plot) ## Final thought - Nearby values are usually similar and differences are noisy and accidental - We can get rid of noise by averaging - Plots are often improved by _smoothing_ - We'll spend next week looking at that ## Take-aways - Know your data type (points vs. regions) - Make plots - Plots over time are pretty simple: + Horizontal coordinate shows time, vertical coordinate shows value + Some tricks can highlight relative change and/or repeating patterns + Use a scatter-plot rather than two vertical axes - Plots over space are usually maps: + 2D plot coordinates show spatial coordinates + Values are shown by a 3rd axis or by color + You may need to remember that the Earth is round + You may need to remember that most of the Earth has few people on it - Don't forget your old EDAs, but do try sub-dividing the data to look for variation in time or space ## After-notes - On the way earthquakes cluster in time and space but are still hard to predict, see @Hough-on-earthquakes + So far as I know, the best way to predict earthquakes is still the rule proposed by @Luen-Stark, which is basically "they either happen at random, or near a big quake, shortly after a big quake" + On the 2011 Tohoku earthquake, and the subsequent disaster at the Fukushima nuclear power plant, the [Wikipedia article](https://en.wikipedia.org/wiki/2011_T%C5%8Dhoku_earthquake_and_tsunami) is surprisingly decent + It's also worth tracking down Peter Galison's film [_Containment_](http://www.containmentmovie.com/) - The lynxes: we'll be revisiting them multiple times over the course, but if you're impatient, look at @Sigmund-games-of-life, ch. 3 - On the "Great Migration" of blacks from the South to cities in the North (and West), @Wilkerson-other-suns is a great popular history + Some of what I was saying about the geographic roots of the "black belt" across the South, and its continuing influence, comes via @Deep-Roots ## References