Lab06: Nonlinear Models

Today is more of a hodgepodge of tips, tricks, and tools rather than a coherent “lesson” or data analysis.

I’ll illustrate techniques using simulated data, and you’ll use gapminder data for your submission. (You’ll need to install the package and run library(gapminder) in order to access the data). The gapminder data was compiled by Hans Rosling, another household name in data viz, and there’s a super cool online tool for exploring the data.

Spend a little bit of time using the online tool to explore the variables and visualizations, thinking about what geoms the variables could be mapped to in order to recreate the visualizations.

First, we’re going to define a ggplot theme for the entire document. This is helpful when you want multiple graphs to share the same theme, but don’t want to change the theme every single time. In the chunk of code below, I’m defining a variable called lab06_theme that I’ll use as my theme throughout this lab.

lab06_theme = theme_minimal() + 
  theme(text = element_text(family = "serif", size = 12))

Part 1: Nonlinear regression

First, we’ll look at some data which is clearly non-linear:

data = read_csv("data/swarthmore.csv", col_names = FALSE)
ggplot(data, aes(x = X1, y = X2)) + 
  geom_point() +
  lab06_theme

We can add an lm() line to this graph:

ggplot(data, aes(x = X1, y = X2)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  lab06_theme

As expected, it adds absolutely no value! Let’s try adding a polynomial - the poly(x,2) creates an orthogonal polynomial of degree 2 (you should rarely fit these models in practice, so it’s OK to skip understanding the details! Think of it as a quadratic equation best fit to the data)

ggplot(data, aes(x = X1, y = X2)) + 
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE) +
  geom_smooth(method = "lm", se = FALSE, col = "lightgray") +
  lab06_theme

Now, let’s try a better statistical method. Add geom_smooth(method = "loess", se = FALSE) to the graph above - your output should look something like the graph below (you can delete the other lines if you like). LOESS stands for Locally Estimated Scatterplot Smoothing. The basic idea is that it takes a moving window and draws the best fit line for that window, and then smooths all the windows together.

We can change the size of the moving window by changing the span argument. Smaller spans lead to bumpier curves, and larger spans lead to smoother fits. The plot below uses span = 0.2, but you should try a few different values of span:

ggplot(data, aes(x = X1, y = X2)) + 
  geom_point() +
  geom_smooth(method = "loess", span = .2, se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, col = "lightgray") +
  geom_smooth(method = "lm", se = FALSE, alpha =.2, col = "lightgray") +
  lab06_theme

Another option in geom_smooth is gam, which stands for Generalized Additive Model. GAMs are super flexible but still interpretable (as compared to more “black box” methods often used in machine learning). There’s a great free interactive course on GAMs in R which I highly recommend if you want to learn more about them. Add geom_smooth(method = "gam", se = FALSE) to see what a GAM looks like on this graph.

If we’re adding a bunch of different fits to the data, it would be helpful if our legend showed what they were! Here’s a sneaky way to define a color aesthetic as we go, where each of the lines is assigned a color that’s equal to the model name, and then we rename the aesthetic “Model” at the end in order to make the legend.

ggplot(data, aes(x = X1, y = X2)) + 
  geom_point() +
  geom_smooth(method = "gam", se = FALSE, aes(col = "GAM")) +
  geom_smooth(method = "loess", span = .2, se = FALSE, aes(col = "Loess")) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, aes(col = "Quadratic")) +
  geom_smooth(method = "lm", se = FALSE, alpha =.2, aes(col = "OLS")) +
  lab06_theme + 
  scale_color_viridis_d("Model", end = .75, option = "C") 

Part 2: Density Plots

Like we saw in class today, density plots are closely related to histograms, and you’ve all seen a histogram! There’s a great Interactive Essay you can check out if you have time to see how complicated putting together histograms really is under the hood.

One of the main problems with histograms is the bumpiness of the bars. This is often because we’re plotting a sample, while a density estimate provides an estimate for the population. Here’s an example of how different samples can look when they’re drawn from the exact same distribution. The set.seed function makes sure that if you run this same code on your computer, you’ll get the exact same set of random samples. The rnorm function tells R to draw 1000 random samples of a normally distributed variable, and the first few lines of code saves that data in a tibble with 10 columns (corresponding to 10 samples) and 1000 rows (corresponding to observations).

set.seed(011221)
samples = data.frame(matrix(rnorm(1000), ncol = 10))
colnames(samples) = paste0("sample", 1:10)
samples = tibble(samples)
samples %>%
  pivot_longer(1:10) %>%
  ggplot(., aes(x = value, fill = name)) + 
  geom_histogram() + 
  facet_wrap(vars(name)) + 
  scale_fill_viridis_d(end = .75, option = "C") + 
  lab06_theme + 
  theme(legend.position = "none")

The code below adds a density estimate to each of the histograms. Note that I had to change the y aesthetic in geom_histogram to map to ..density.. rather than ..count... This is for the probability reasons we talked about in class today, so that the area under the density curve is equal to one. Here, I’ve added bw = .1 to the geom_density, which results in a pretty bumpy curve. Try playing around with different values of bw to get a smoother curve that looks roughly the same across all of our samples.

samples %>%
  pivot_longer(1:10) %>%
  ggplot(., aes(x = value, fill = name)) + 
  geom_histogram(aes(y = ..density..)) + 
  geom_density(alpha = .3, bw = .1, aes(col = name)) +
  facet_wrap(vars(name)) + 
  scale_fill_viridis_d(end = .75, option = "C") + 
  scale_color_viridis_d(end = .75, option = "C") + 
  lab06_theme + 
  theme(legend.position = "none")

Just like a density plot is related to a histogram, a violin plot is related to a box plot. Below is the same set of samples, but plotted using boxplots instead of histograms.

samples %>%
  pivot_longer(1:10) %>%
  ggplot(., aes(y = value, x = name, fill = name)) + 
  geom_boxplot(alpha = .7) + 
  scale_fill_viridis_d(end = .75, option = "C") + 
  scale_color_viridis_d(end = .75, option = "C") + 
  lab06_theme + 
  theme(legend.position = "none")

Now, add geom_jitter(alpha = X, aes(col = name)) to the code above to replicate the graph below. Play around with different alpha values to get something that looks nice. This is a great example of not hiding the raw data behind a model.

Finally, let’s replace geom_boxplot with geom_violin. A violin plot is just a density estimate rotated to be on the other axis, and then mirrored to create a “violin” shape. This graph is a little “busy” and I’d want to tweak it before publishing, but it does a nice job of showing the raw data (geom_jitter) alongside the statistical model (geom_violin).

Finally, if we have many distributions that we want to compare to each other, sometimes it’s easiest if we have a single axes to deal with. Enter ggridges! (Yes, this is yet another package)

library(ggridges)
samples %>%
  pivot_longer(1:10) %>%
  ggplot(., aes(y = name, x = value, fill = name)) + 
  geom_density_ridges(alpha = .8) +
  scale_fill_viridis_d(end = .75, option = "C") + 
  scale_color_viridis_d(end = .75, option = "C") + 
  lab06_theme + 
  theme(legend.position = "none")

Part 3: Open-ended

Reminder: here you’ll use the gapminder data:

library(gapminder)
head(gapminder)

## # A tibble: 6 x 6
##   country     continent  year lifeExp      pop gdpPercap
##   <fct>       <fct>     <int>   <dbl>    <int>     <dbl>
## 1 Afghanistan Asia       1952    28.8  8425333      779.
## 2 Afghanistan Asia       1957    30.3  9240934      821.
## 3 Afghanistan Asia       1962    32.0 10267083      853.
## 4 Afghanistan Asia       1967    34.0 11537966      836.
## 5 Afghanistan Asia       1972    36.1 13079460      740.
## 6 Afghanistan Asia       1977    38.4 14880372      786.

Nonlinear smoothers

Find two variables from the gapminder, make a scatterplot, add a geom_smooth layer with a linear relationship and a geom_smooth layer with a nonlinear relationship (your choice!). Which fits better? Add colors and a legend.

Density

Choose a variable (or combination of variables) from the gapminder dataset and make either a kernel density plot, violin plot, or density ridges plot. Write 1-2 sentences summarizing your findings.