suppressPackageStartupMessages(library("modelr"))
package ă¤¼ă¸±modelră¤¼ă¸² was built under R version 3.6.3
suppressPackageStartupMessages(library("tidyverse"))
package ă¤¼ă¸±tidyverseă¤¼ă¸² was built under R version 3.6.3
suppressPackageStartupMessages(library("gapminder"))
package ă¤¼ă¸±gapminderă¤¼ă¸² was built under R version 3.6.3

1. A linear trend seems to be slightly too simple for the overall trend. Can you do better with a quadratic polynomial? How can you interpret the coefficients of the quadratic? Hint you might want to transform year so that it has mean zero.)

The following code replicates the analysis in the chapter but replaces the function country_model() with a regression that includes the year squared.

lifeExp ~ poly(year, 2)
lifeExp ~ poly(year, 2)
country_model <- function(df) {
  lm(lifeExp ~ poly(year - median(year), 2), data = df)
}

by_country <- gapminder %>%
  group_by(country, continent) %>%
  nest()

by_country <- by_country %>%
  mutate(model = map(data, country_model))
by_country <- by_country %>%
  mutate(
    resids = map2(data, model, add_residuals)
  )
by_country

unnest(by_country, resids) %>%
  ggplot(aes(year, resid)) +
  geom_line(aes(group = country), alpha = 1 / 3) +
  geom_smooth(se = FALSE)


by_country %>%
  mutate(glance = map(model, broom::glance)) %>%
  unnest(glance, .drop = TRUE) %>%
  ggplot(aes(continent, r.squared)) +
  geom_jitter(width = 0.5)
The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
All list-columns are now preserved.
This warning is displayed once per session.
Call `lifecycle::last_warnings()` to see where this warning was generated.

2. Explore other methods for visualizing the distribution of \(R^2\) per continent. You might want to try the ggbeeswarm package, which provides similar methods for avoiding overlaps as jitter, but uses deterministic methods.

suppressPackageStartupMessages(library("ggbeeswarm"))
by_country %>%
  mutate(glance = map(model, broom::glance)) %>%
  unnest(glance, .drop = TRUE) %>%
  ggplot(aes(continent, r.squared)) +
  geom_beeswarm()

3. To create the last plot (showing the data for the countries with the worst model fits), we needed two steps: we created a data frame with one row per country and then semi-joined it to the original dataset. It’s possible to avoid this join if we use unnest() instead of unnest(.drop = TRUE). How?

gapminder %>%
  group_by(country, continent) %>%
  nest() %>%
  mutate(model = map(data, ~ lm(lifeExp ~ year, .))) %>%
  mutate(glance = map(model, broom::glance)) %>%
  unnest(glance) %>%
  unnest(data) %>%
  filter(r.squared < 0.25) %>%
  ggplot(aes(year, lifeExp)) +
  geom_line(aes(color = country))

LS0tDQp0aXRsZTogImdhcG1pbmRlciINCm91dHB1dDogDQogIGh0bWxfbm90ZWJvb2s6DQogICAgdG9jOiB0cnVlDQogICAgdG9jX2Zsb2F0OiB0cnVlDQotLS0NCg0KYGBge3J9DQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeSgibW9kZWxyIikpDQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeSgidGlkeXZlcnNlIikpDQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeSgiZ2FwbWluZGVyIikpDQpgYGANCg0KIyMjIDEuIEEgbGluZWFyIHRyZW5kIHNlZW1zIHRvIGJlIHNsaWdodGx5IHRvbyBzaW1wbGUgZm9yIHRoZSBvdmVyYWxsIHRyZW5kLiBDYW4geW91IGRvIGJldHRlciB3aXRoIGEgcXVhZHJhdGljIHBvbHlub21pYWw/IEhvdyBjYW4geW91IGludGVycHJldCB0aGUgY29lZmZpY2llbnRzIG9mIHRoZSBxdWFkcmF0aWM/IEhpbnQgeW91IG1pZ2h0IHdhbnQgdG8gdHJhbnNmb3JtIHllYXIgc28gdGhhdCBpdCBoYXMgbWVhbiB6ZXJvLikNCg0KVGhlIGZvbGxvd2luZyBjb2RlIHJlcGxpY2F0ZXMgdGhlIGFuYWx5c2lzIGluIHRoZSBjaGFwdGVyIGJ1dCByZXBsYWNlcyB0aGUgZnVuY3Rpb24gYGNvdW50cnlfbW9kZWwoKWAgd2l0aCBhIHJlZ3Jlc3Npb24gdGhhdCBpbmNsdWRlcyB0aGUgeWVhciBzcXVhcmVkLg0KDQpgYGB7cn0NCmxpZmVFeHAgfiBwb2x5KHllYXIsIDIpDQpjb3VudHJ5X21vZGVsIDwtIGZ1bmN0aW9uKGRmKSB7DQogIGxtKGxpZmVFeHAgfiBwb2x5KHllYXIgLSBtZWRpYW4oeWVhciksIDIpLCBkYXRhID0gZGYpDQp9DQoNCmJ5X2NvdW50cnkgPC0gZ2FwbWluZGVyICU+JQ0KICBncm91cF9ieShjb3VudHJ5LCBjb250aW5lbnQpICU+JQ0KICBuZXN0KCkNCg0KYnlfY291bnRyeSA8LSBieV9jb3VudHJ5ICU+JQ0KICBtdXRhdGUobW9kZWwgPSBtYXAoZGF0YSwgY291bnRyeV9tb2RlbCkpDQpieV9jb3VudHJ5IDwtIGJ5X2NvdW50cnkgJT4lDQogIG11dGF0ZSgNCiAgICByZXNpZHMgPSBtYXAyKGRhdGEsIG1vZGVsLCBhZGRfcmVzaWR1YWxzKQ0KICApDQpieV9jb3VudHJ5DQoNCnVubmVzdChieV9jb3VudHJ5LCByZXNpZHMpICU+JQ0KICBnZ3Bsb3QoYWVzKHllYXIsIHJlc2lkKSkgKw0KICBnZW9tX2xpbmUoYWVzKGdyb3VwID0gY291bnRyeSksIGFscGhhID0gMSAvIDMpICsNCiAgZ2VvbV9zbW9vdGgoc2UgPSBGQUxTRSkNCg0KYnlfY291bnRyeSAlPiUNCiAgbXV0YXRlKGdsYW5jZSA9IG1hcChtb2RlbCwgYnJvb206OmdsYW5jZSkpICU+JQ0KICB1bm5lc3QoZ2xhbmNlLCAuZHJvcCA9IFRSVUUpICU+JQ0KICBnZ3Bsb3QoYWVzKGNvbnRpbmVudCwgci5zcXVhcmVkKSkgKw0KICBnZW9tX2ppdHRlcih3aWR0aCA9IDAuNSkNCmBgYA0KDQojIyMgMi4gRXhwbG9yZSBvdGhlciBtZXRob2RzIGZvciB2aXN1YWxpemluZyB0aGUgZGlzdHJpYnV0aW9uIG9mICRSXjIkICBwZXIgY29udGluZW50LiBZb3UgbWlnaHQgd2FudCB0byB0cnkgdGhlIGBnZ2JlZXN3YXJtYCBwYWNrYWdlLCB3aGljaCBwcm92aWRlcyBzaW1pbGFyIG1ldGhvZHMgZm9yIGF2b2lkaW5nIG92ZXJsYXBzIGFzIGppdHRlciwgYnV0IHVzZXMgZGV0ZXJtaW5pc3RpYyBtZXRob2RzLg0KDQpgYGB7cn0NCnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KCJnZ2JlZXN3YXJtIikpDQpieV9jb3VudHJ5ICU+JQ0KICBtdXRhdGUoZ2xhbmNlID0gbWFwKG1vZGVsLCBicm9vbTo6Z2xhbmNlKSkgJT4lDQogIHVubmVzdChnbGFuY2UsIC5kcm9wID0gVFJVRSkgJT4lDQogIGdncGxvdChhZXMoY29udGluZW50LCByLnNxdWFyZWQpKSArDQogIGdlb21fYmVlc3dhcm0oKQ0KYGBgDQoNCiMjIyAzLiBUbyBjcmVhdGUgdGhlIGxhc3QgcGxvdCAoc2hvd2luZyB0aGUgZGF0YSBmb3IgdGhlIGNvdW50cmllcyB3aXRoIHRoZSB3b3JzdCBtb2RlbCBmaXRzKSwgd2UgbmVlZGVkIHR3byBzdGVwczogd2UgY3JlYXRlZCBhIGRhdGEgZnJhbWUgd2l0aCBvbmUgcm93IHBlciBjb3VudHJ5IGFuZCB0aGVuIHNlbWktam9pbmVkIGl0IHRvIHRoZSBvcmlnaW5hbCBkYXRhc2V0LiBJdOKAmXMgcG9zc2libGUgdG8gYXZvaWQgdGhpcyBqb2luIGlmIHdlIHVzZSBgdW5uZXN0KClgIGluc3RlYWQgb2YgYHVubmVzdCguZHJvcCA9IFRSVUUpYC4gSG93Pw0KDQpgYGB7cn0NCmdhcG1pbmRlciAlPiUNCiAgZ3JvdXBfYnkoY291bnRyeSwgY29udGluZW50KSAlPiUNCiAgbmVzdCgpICU+JQ0KICBtdXRhdGUobW9kZWwgPSBtYXAoZGF0YSwgfiBsbShsaWZlRXhwIH4geWVhciwgLikpKSAlPiUNCiAgbXV0YXRlKGdsYW5jZSA9IG1hcChtb2RlbCwgYnJvb206OmdsYW5jZSkpICU+JQ0KICB1bm5lc3QoZ2xhbmNlKSAlPiUNCiAgdW5uZXN0KGRhdGEpICU+JQ0KICBmaWx0ZXIoci5zcXVhcmVkIDwgMC4yNSkgJT4lDQogIGdncGxvdChhZXMoeWVhciwgbGlmZUV4cCkpICsNCiAgZ2VvbV9saW5lKGFlcyhjb2xvciA9IGNvdW50cnkpKQ0KYGBgDQo=