TidyTuesday 2022, Week 32

Let’s keep this TidyTuesday rolling by loading packages and reading the data:

library(tidyverse)
library(lubridate)
library(lme4)
library(arm)
library(ggeffects)

wheels <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2022/2022-08-09/wheels.csv')

head(wheels)

## # A tibble: 6 x 22
##      X1 name  height diameter opened     closed     country location
##   <dbl> <chr>  <dbl>    <dbl> <date>     <date>     <chr>   <chr>   
## 1     1 360 ~   200       NA  2012-07-03 2013-01-01 USA     Pensaco~
## 2     2 Amur~   303      200. 2004-01-01 NA         Japan   Kagoshi~
## 3     3 Asia~   200      200  2012-12-15 NA         Tailand Asiatiq~
## 4     4 Auro~   295      272  NA         NA         Japan   Nagashi~
## 5     5 Bagh~   180       NA  2011-01-01 NA         Iraq    Al-Zawr~
## 6     6 Beij~   693.     643. NA         NA         China   Chaoyan~
## # ... with 14 more variables: number_of_cabins <dbl>,
## #   passengers_per_cabin <dbl>, seating_capacity <dbl>, hourly_capacity <dbl>,
## #   ride_duration_minutes <dbl>, climate_controlled <chr>,
## #   construction_cost <chr>, status <chr>, design_manufacturer <chr>,
## #   type <chr>, vip_area <chr>, ticket_cost_to_ride <chr>,
## #   official_website <chr>, turns <dbl>

Things are relatively clean this week! We are going to focus on replicating the second image provided on GitHub. To do so, let’s steal some secrets from Stack Overflow

circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){
    r = diameter / 2
    tt <- seq(0,2*pi,length.out = npoints)
    xx <- center[1] + r * cos(tt)
    yy <- center[2] + r * sin(tt)
    return(data.frame(x = xx, y = yy))
}

dat <- circleFun(c(1,-1),2.3,npoints = 100)
#geom_path will do open circles, geom_polygon will do filled circles
ggplot(dat,aes(x,y)) + 
  geom_path() + 
  coord_fixed()

We got a circle! Now let’s think about how to replicate the image for one wheel. We need one for which we have all of (1) name, (2) height, (3) diameter, and (4) number_of_cabins.

wheels %>% 
  dplyr::select(name, height, diameter, number_of_cabins) %>% 
  na.omit(.) -> full_wheels 

We may as well calculate the centers right away:

full_wheels %>% 
  mutate(center = height - diameter/2) -> full_wheels

Let’s start with the first one.

arbitrary_wheel <- full_wheels[1,]

Let’s start by plotting the basic circle and cleaning things up:

  with(arbitrary_wheel,
       circleFun(center = c(0,center), 
                 diameter = diameter,
                 npoints = number_of_cabins)) %>% 
  ggplot(aes(x,y)) +
    geom_path() + 
    coord_fixed() +
    theme_bw() +
    xlab("") +
    ylab("") +
    ggtitle(arbitrary_wheel$name) +
    scale_y_continuous(expand = c(0, 0), 
                       limits = c(0, arbitrary_wheel$diameter/2 + 
                                     arbitrary_wheel$center + 10))-> basic

basic

Sweet sweet victory! Now let’s add the legs. A feature which is a bit annoying about this approach is that we have to re-introduce our underlying dataset.

basic +
  geom_segment(data=arbitrary_wheel,aes(x = 0 - center/3,
                                        y = 0,
                                        xend = 0,
                                        yend = center)) +
  geom_segment(data=arbitrary_wheel,aes(x = 0 + center/3,
                                        y = 0,
                                        xend = 0,
                                        yend = center)) -> basic2

basic2

This one is a bit strange due to its height, but the basic idea is there.

Now let’s add on some carts. We need however many carts are indicated, broken into three groups, and arrayed across the wheel. Luckily, we have those points already in our underlying circle dataframe!

basic2 +
  geom_point(pch=24, 
             position = position_nudge(y=-5),
             size = 2,
             aes(fill = as.factor(rep(1:3,length = length(x))))) +
  theme(legend.position = "none")

Cute. Now that we have the basic structure down, let’s go ahead and get things ready to make these in bulk. The strategy will be to generate a new dataset which includes the above plotting fields for all of the wheels in our set. Here is a quick approach, subsetted down to the few names in the

wheelzzz <- full_wheels$name

lapply(1:length(wheelzzz), function(i){
  with(full_wheels[i,],
       circleFun(center = c(0,center), 
                 diameter = diameter,
                 npoints = number_of_cabins)
       ) -> out
  out$name <- wheelzzz[i]
  out$diameter <- rep(full_wheels[i,"diameter"],nrow(out))
  out$center <- rep(full_wheels[i,"center"],nrow(out))
  out
}) -> wheelzzz

wheelzzz <- do.call("rbind",wheelzzz)

Now let’s plot them all as a facet wrap, subsetting down to just those six wheels included in the graphic:

which_wheels <- c("Beijing Great Wheel",
                  "Golden Gate Flyer",
                  "High Roller",
                  "London Eye",
                  "Singapore Flyer",
                  "Wiener Riesenrad")

wheelzzz %>% 
  filter(name %in% which_wheels) -> less_wheelzzz

less_wheelzzz %>% 
  dplyr::select(name, diameter, center) %>% 
  group_by(name) %>% 
  summarise(name = unique(name),
            diameter = as.numeric(unique(diameter)),
            center = as.numeric(unique(center)),
            .groups = "keep") -> legzzz

less_wheelzzz %>% 
  ggplot(aes(x,y)) +
    facet_wrap(~name, ncol = 3) +
    geom_path() + 
    coord_fixed() +
    theme_bw() +
    xlab("") +
    ylab("") +
    geom_point(pch=24, 
               position = position_nudge(y=-20),
               size = 1.5,
               aes(fill = as.factor(rep(1:3,length = length(x))))) +
    theme(legend.position = "none") +
    geom_segment(data=legzzz,aes(x = 0 - center/3,
                                          y = 0,
                                          xend = 0,
                                          yend = center)) +
    geom_segment(data=legzzz,aes(x = 0 + center/3,
                                          y = 0,
                                          xend = 0,
                                          yend = center)) +
    scale_y_continuous(expand = c(0, 0), 
                       limits = c(-9, legzzz$diameter/2 + 
                                      legzzz$center + 50))

Hey, good enough for government work! Now that the replication is complete, is there anything more interesting we can look at with this data?

One opportunity to practice some Tidy skillz might be to look at the cost of building these wheely-boiis, fitting a model, and making a nice graphical illustration of the model. To start, we need to clean the cost vector, which is super easy with Tidy:

head(wheels$construction_cost)

## [1] "Unknown"          "Unknown"          "Unknown"          "Unknown"         
## [5] "$6 million USD"   "$290 million USD"

Gross. Luckily there is a common structure for all of these observations, except for the Golden Gate Flyer, which reports a range of estimates for the construction cost. Regardless, everything is listed in millions of USD so we might just …

wheels$cost <- parse_number(wheels$construction_cost)

## Warning: 40 parsing failures.
## row col expected  actual
##   1  -- a number Unknown
##   2  -- a number Unknown
##   3  -- a number Unknown
##   4  -- a number Unknown
##   8  -- a number Unknown
## ... ... ........ .......
## See problems(...) for more details.

The parsing errors are fine, they just mean “could not find a number” here. Now let’s run a model after a little more cleaning. Note that there are not many observations, and there is a good amount of missingness which guides covariate selection, and the resulting model will be almost by definition quite crummy.

wheels$country <- as.factor(wheels$country)

wheels %>% 
  dplyr::select(cost, height, country) %>% 
  na.omit(.) -> wheel_dat

m <- lmer(cost ~ poly(height,2) + (1 |country), 
            data = wheel_dat)

pr <- ggpredict(m, "height [all]")
plot(pr) +
  xlab("Wheel Height") +
  ylab("Wheel Cost")

The above is the predicted cost as a function of height. As we would expect, it’s increasing at an increasing rate. Now let’s think about plotting the BLUPs.

ranefs <- data.frame(rownames(ranef(m)$country),
                     ranef(m)$country,
                     se.ranef(m)$country)
colnames(ranefs) <- c("country", "intercept", "se")
ranefs$title <- "Random Intercept Estimates"

ranefs$country <- factor(as.character(ranefs$country),
                         levels = ranefs$country[order(ranefs$intercept)])

ggplot(ranefs, aes(x = country, y = intercept)) +
  geom_point() +
  geom_errorbar(aes(ymin = intercept - 2*se,
                    ymax = intercept + 2*se),width=0) +
  coord_flip() +
  theme_bw() +
  facet_grid(~title) +
  geom_hline(yintercept = 0,
             lty = 2) +
  ylab("Expenditure Deviation from Average") +
  xlab("Country") +
  theme(strip.text.x = element_text(size = 12))

So it looks like, taking height into account, the USA spends about 48 million USD more than average on their wheels whereas China spends about 76 million USD less. As the saying goes Chabuduo! Close enough …

p.s. Farewell to the ICPSR 2022 Summer Program!