By State

penis %>% 
  gather(key = state, value = volume, volume_flaccid, volume_erect) %>%
  mutate(state = factor(state, 
                        levels = c("volume_flaccid", "volume_erect"), 
                        labels =  c("Flaccid", "Erect"))) %>%
  {
    ggplot(., aes(x = volume, fill = state, color = state)) +
    geom_histogram(binwidth = 5, alpha = .7, position = "dodge") +
    scale_x_continuous(breaks = seq(0, 500, 25), 
                       minor_breaks = seq(0, 500, 5),
                       labels = math_format(.x ~cm^3),
                       sec.axis = sec_axis(trans = ~./2.54^3, 
                                           labels = math_format(.x ~inch^3))) +
    scale_fill_brewer(palette = "Paired") +
    scale_color_brewer(palette = "Paired") +
    labs(title = "World Penis Data", subtitle = "Volume by State",
         x = "Volume", y = "Frequency", 
         fill = "State", color = "State", caption = plot_caption) +
    theme(legend.position = "top")
  }

By State and Region

penis %>% 
  gather(key = state, value = volume, volume_flaccid, volume_erect) %>%
  mutate(state = factor(state, 
                        levels = c("volume_flaccid", "volume_erect"), 
                        labels =  c("Flaccid", "Erect"))) %>%
  group_by(Region) %>%
  mutate(order_volume = mean(volume)) %>%
  {
    ggplot(., aes(x = reorder(Region, order_volume), y = volume, 
                  fill = state, color = state)) +
    geom_boxplot(alpha = .7) +
    coord_flip() +
    scale_y_continuous(breaks = seq(0, 500, 25), 
                       minor_breaks = seq(0, 500, 5),
                       labels = math_format(.x ~cm^3),
                       sec.axis = sec_axis(trans = ~./2.54^3, 
                                           labels = math_format(.x ~inch^3))) +
    scale_fill_brewer(palette = "Paired") +
    scale_color_brewer(palette = "Paired") +
    labs(title = "World Penis Data", subtitle = "Volume by State",
         x = "", y = "Volume", 
         fill = "State", color = "State", caption = plot_caption) +
    theme(axis.text.y = element_text(size = rel(1.2)),
          legend.position = "top")
  }

# CIs
penis_long %>%
  group_by(Region) %>%
  mutate(order_volume = mean(volume)) %>%
  {
    ggplot(., aes(x = reorder(Region, order_volume), y = volume,
                  fill = state, color = state)) +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar") +
    stat_summary(fun.y = mean, geom = "point", size = 2) +
    coord_flip() +
    scale_y_continuous(breaks = seq(0, 500, 25),
                       minor_breaks = seq(0, 500, 5),
                       labels = math_format(.x ~cm^3),
                       sec.axis = sec_axis(trans = ~./2.54^3, 
                                           labels = math_format(.x ~inch^3))) +
    scale_color_brewer(palette = "Paired") +
    labs(title = "World Penis Data", subtitle = "Volume by State",
         x = "", y = "Volume",
         fill = "State", color = "State", caption = plot_caption) +
    theme(axis.text.y = element_text(size = rel(1.2)),
          legend.position = "top")
  }

By Length

Since volume is calculated by circumference and length, the relationship between volume and length is pretty self-explanatory. We’re still going to plot it, of course.

ggplot(data = penis_long, aes(x = length, y = volume, color = state)) +
  geom_smooth(method = lm, se = F, color = "gray") +
  geom_smooth(method = lm, se = F) +
  geom_point(size = 3, color = "black") +
  geom_point(size = 2) +
  scale_x_continuous(labels = label_cm,
                     sec.axis = sec_axis(trans = ~./2.54,
                                         labels = label_in)) +
  scale_y_continuous(labels = math_format(.x ~cm^3),
                     sec.axis = sec_axis(trans = ~./2.54^3,
                                         labels = math_format(.x ~inch^3))) +
  scale_color_brewer(palette = "Paired") +
  labs(title = "World Penis Data",
       subtitle = "Volume by Length",
       x = "Length", y = "Volume", color = "State",
       caption = plot_caption) +
  theme(legend.position = "top")

Choropleth-Penismap

A choropleth map displaying the worldwide distribution of the volume of erect penisses. Enjoy.

highchart() %>% 
  hc_add_series_map(worldgeojson, map, value = "volume_erect", joinBy = "iso3",
                    name = "Penis Volume (erect)") %>% 
  hc_title(text = "Erect Penis Volume by Country") %>% 
  hc_tooltip(valueDecimals = 2, valuePrefix = "<b>", valueSuffix = " cm<sup>3</sup></b>", 
             useHTML = TRUE) %>% 
  hc_colorAxis(stops = idk, min = 50, max = 250) %>% 
  hc_legend(enabled = TRUE) %>% 
  hc_mapNavigation(enabled = TRUE)