---
title: "Summarizing Data with Group-By Analyses"
author: "Rian DeFaccio , for BH Analytics LLC"
output: 
  html_document:
    toc: true
    toc_depth: 3
---

```{r setup, include=FALSE}
knitr::opts_knit$set(warning = FALSE, message = FALSE)
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```

```{r}
# load libraries and adjust settings
library(tidyverse)
```

## Introduction

------------------------------------------------------------------------

In this lesson we will take a deep dive into using group-by analyses and do some more different kinds of multivariate exploratory analysis. We will cover the following topics:

-   Group-By Analysis with dplyr
-   Multiple Group-Bys
-   Pivot Tables
-   Cross Tabulations
-   Exploring Categorical vs. Continuous Variables
-   Exploring Two Categorical Variables

## Gapminder Data

------------------------------------------------------------------------

This notebook uses the gapminder data. This is an interesting data set containing life expectancy, GDP, and population for all countries on the planet in 5 year increments from 1952 to 2007. You can learn more about the data at the gapminder website:

🌐 <https://www.gapminder.org/data/>

```{r}
df_gapminder <- read.csv('../data/gapminder.csv')

# create a subset of data for just the year 2007
df_2007 <- df_gapminder %>%
  filter(year == 2007 & continent != 'Oceania')
head(df_2007)
```

Now let's examine the structure of the data frame:

```{r}
str(df_2007)
```

## Adding Group-By to a Data Pipe

------------------------------------------------------------------------

One way to identify if two variables are related to each other is to use a group-by analysis, also called 'split-apply-combine'. The basic idea is as follows:

1.  **Split** the data set apart into groups (based on a categorical variable).
2.  **Apply** a summary function to those groups.
3.  **Combine** the results into a final table for display.

To see how this works, first consider how you would get the overall mean life expectancy in 2007 using dplyr.

```{r}
## overall mean
df_2007 %>%
  summarise(mean_life_exp = mean(lifeExp))
```

To get this same summary statistic by groups, we use the `group_by()` function as part of the pipeline.

```{r}
## mean conditional on continent
df_2007 %>%
  group_by(continent) %>%
  summarise(mean_life_exp = mean(lifeExp))
```

If you want to get multiple summary statistics, add an argument to `summarise()` for each statistic.

```{r}
## calculate multiple statistics on each continent
df_2007 %>%
  group_by(continent) %>%
  summarise(mean_life_exp   = mean(lifeExp),
            sd_life_exp     = sd(lifeExp),
            median_life_exp = median(lifeExp))
```

## Multiple Group-By

We can have multiple variables define the groups by simply adding them into the `group_by()` function. Here we return to the original gapminder data, because we want to also group by year.

```{r}
## group by multiple columns
df_gapminder %>%
  group_by(continent, year) %>%
  summarise(mean_life_exp = mean(lifeExp))
```

When you do a group-by analysis with multiple grouping variables, only the last group is stripped away from the result, meaning you still have a grouped data frame. This can be problematic to analyses that are downstream of this result.

### ✍️ Your Turn! Identify the Grouping Variable

Can you identify the remaining grouping variable in the above result?

**Answer here:**

There are two ways to remove left-over grouping variables after a group-by analysis, see the next two examples for both methods.

```{r}
## remove left-over grouping variable (method 1)
df_gapminder %>%
  group_by(continent, year) %>%
  summarise(mean_life_exp = mean(lifeExp)) %>%
  ungroup()

## remove left-over grouping variable (method 2)
df_gapminder %>%
  group_by(continent, year) %>%
  summarise(mean_life_exp = mean(lifeExp),
            .groups = "drop")
```

Now we have a regular data frame with no grouping. However, the output is in a long format. We may want to pivot these results to a wider format for easy viewing. We will use the `pivot_wider()` function from the tidyr package.

### 🛠️Coding Practice: Pivot Data Wider

Finish the following code so that you get a data frame where each row is a continent, and each column is a year.

```{r}
# FINISH THIS CODE
df_gapminder %>%
  group_by(continent, year) %>%
  summarise(mean_life_exp = mean(lifeExp),
            .groups = "drop")
# more code here ...
```

**END OF EXERCISE**

------------------------------------------------------------------------

## A Pivot Table is Like a Group-By

------------------------------------------------------------------------

In Excel, SQL, SAS, and many other data analysis software there is a tool called a pivot table which summarizes several variables at one time.

> ❓**QUICK POLL:** Who has used pivot tables previously in some other software?

We can see in the following picture how the usual arguments to a pivot table procedure map very conveniently to a double group-by procedure using tidyverse tools.

> `dplyr::group_by(var1, var2)` ➜ `dplyr::summarise(...)` ➜ `tidyr::pivot_wider(...)`

Let's consider our prior procedure. Here's the arguments we might expect to give to a pivot table function in some other software:

![](images/pivot_table-01.png){width="1000"}

And for reference, here's the same procedure implemented with tidyverse tools:

```{r}
df_gapminder %>%
  group_by(continent, year) %>%
  summarise(mean_life_exp = mean(lifeExp),
            .groups = "drop") %>%
  pivot_wider(names_from = year, values_from = mean_life_exp)
```

## A Cross-Tabulation is Like a Pivot Table

------------------------------------------------------------------------

Finally, let's consider the case when we are doing a pivot table, but the aggregating function is just a count of rows. There is actually a name for this type of procedure, it's called a **cross tabulation**. There are a couple of ways to build these in R.

First, the quick and easy base R way, just use the `table()` function!

```{r}
# first, making a variable of high vs. low population
df_2007$life_exp_binary <- ifelse(df_2007$lifeExp >= 70, "High", "Low")

# count of high vs. low life expectancy by continent in 2007
table(df_2007$continent, df_2007$life_exp_binary)
```

We want to show you a trick now which is the base R `with()` function. Check it out in the next example:

```{r}
# same table but using with()
with(df_2007, table(continent, life_exp_binary))
```

Now, we can also do this with a familiar double group-by procedure, this is a kind of pivot table.

```{r}
df_2007 %>%
  group_by(continent, life_exp_binary) %>%
  summarise(num_countries = n(),
            .groups = "drop") %>%
  pivot_wider(names_from = life_exp_binary, values_from = num_countries)
```

### 🎯Quick Check

What was the aggregating function in the last code cell?

**Answer here:**

Performing the pivot to the wider form introduced some `NA`s. We can fill those in with the `values_fill=` argument to `pivot_wider()`.

```{r}
df_2007 %>%
  group_by(continent, life_exp_binary) %>%
  summarise(num_countries = n(),
            .groups = "drop") %>%
  pivot_wider(names_from = life_exp_binary, 
              values_from = num_countries,
              values_fill = 0)
```

And finally, the creators of the tidyverse knew this would be a common task, so they created the `count()` function to simplify things a bit:

```{r}
df_2007 %>%
  group_by(continent, life_exp_binary) %>%
  count() %>%
  pivot_wider(names_from = life_exp_binary, 
              values_from = n,
              values_fill = 0)
```

In fact, we can even just supply the grouping variables as arguments to `count()`!

```{r}
df_2007 %>%
  count(life_exp_binary, continent) %>%
  pivot_wider(names_from = life_exp_binary, 
              values_from = n,
              values_fill = 0)
```

## 📈 EDA: Categorical by Continuous Data

------------------------------------------------------------------------

Let's get back into some statistics now. Or should we call it data science 😜?

We'll start by making some appropriate plots for when you have one categorical variable and one continuous variable.

### Box Plots

Let's check out life expectancy across the continents in 2007.

The question we are asking is this: *Is the distribution of life expectancy different conditional on continent?*

```{r}
## boxplot
ggplot(df_2007, aes(x = continent, y = lifeExp)) + 
  geom_boxplot() +
  ggtitle("Life Expectancy in 2007")
```

### Colored Density Plot

Another type of plot that would be appropriate here would be a density plot, with the levels of the categorical variable mapped to different colors.

```{r}
ggplot(df_2007, aes(x = lifeExp, fill = continent)) + 
  geom_density(alpha = 0.5) +
  ggtitle("Life Expectancy in 2007")
```

### ✍️ Your Turn! Pick a Plot

Which plot do you prefer for viewing the relationship between life expectancy and continent in 2007? Why?

**Answer here:**

### BONUS: Ridge Plots

One last type of plot that can be very satisfying is a ridge plot. This is great when you want to view the distribution of a continuous variable as it varies over an ordinal categorical variable.

```{r}
library(ggridges)
library(viridis)

# make year into a factor variable
df_gapminder$year_f <- factor(df_gapminder$year)

# build the ridge plot
ggplot(df_gapminder, aes(x = lifeExp, y = year_f, fill = ..x..)) +
  geom_density_ridges_gradient(scale = 3, rel_min_height = 0.001) +
  scale_fill_viridis(name = "Temp. [F]", option = "C") +
  labs(title = 'Life Expectancy on Planet Earth', y = NULL, 
       x = 'Life Expectancy (Years)', y = NULL) +
  theme_bw() +
  theme(
    legend.position="none",
    panel.spacing = unit(0.1, "lines"),
    strip.text.x = element_text(size = 8),
    plot.title = element_text(hjust = 0.5)
  )
```

### 🛠️Coding Practice: Box Plot Version of Ridge Plot

Let's compare this to a more traditional plot. Can you make a box-plot that shows how the distribution of life expectancy has changed over the years?

🗝️**HINT:** You may need to specify the `group =` argument inside your call to `aes()`.

```{r}
# ANSWER HERE

```

**END OF EXERCISE**

------------------------------------------------------------------------

### Group-By Analyses

Then, the final way to examine a categorical by continuous relationship is to do what we just learned about, a group by analysis! We can calculate statistics stratified by year in this case.

### 🛠️Coding Practice: Life Expectancy Grouped by Year

Create a table that summarizes the minimum, median, mean, standard deviation, and maximum values of life expectancy across each year.

-   round your results to 2 decimal places.

-   Pipe your result through the function `knitr::kable()` and see if you like the output.

-   Add the argument `caption = "Life Expectancy"` to the `kable()` function if you'd like.

```{r}
# ANSWER HERE

```

**END OF EXERCISE**

------------------------------------------------------------------------

## 📈 EDA: Categorical by Categorical Data

------------------------------------------------------------------------

The last set of variable types we can consider in our exploratory analysis is when you have a pair of categorical variables. We are considering this from the exploratory analysis angle now, but we just considered the mechanics of such a procedure earlier.

First let's review how to summarize a single categorical variable.

```{r}
# base R tools
table(df_2007$continent)
table(df_2007$continent) %>% prop.table()
```

```{r}
## now with dplyr
df_2007 %>%
  count(continent)
```

```{r}
# the count and proportion
df_2007 %>%
  count(continent) %>%
  mutate(prop = n / sum(n))
```

Now for when we are counting observations across two categorical variables...

```{r}
# the count and proportion
df_2007 %>%
  count(continent, life_exp_binary) %>%
  pivot_wider(names_from = life_exp_binary, 
              values_from = n,
              values_fill = 0)
```

Getting row or column percents...

```{r}
# getting row percents
df_2007 %>%
  count(continent, life_exp_binary) %>%
  pivot_wider(names_from = life_exp_binary, 
              values_from = n,
              values_fill = 0) %>%
  mutate(perc_high = High / (High + Low),
         perc_low  = Low / (High + Low)) %>%
  select(continent, contains('high'), contains('low'))
```

Getting column proportions is not so easy.

In practice, base R tools are actually very convenient for these types of tasks:

```{r}
# row proportions
table(df_2007$continent, df_2007$life_exp_binary) %>%
  prop.table(margin=1)
```

```{r}
# column proportions
table(df_2007$continent, df_2007$life_exp_binary) %>%
  prop.table(margin=2)
```

### 🛠️Coding Practice: Adding `with()`

Can you adjust the code from the previous cell to use the `with()` function?

```{r}
# ANSWER HERE
```

**END OF EXERCISE**

------------------------------------------------------------------------

## Conclusion

------------------------------------------------------------------------

This concludes our deep dive into group-by analyses as well as performing multivariate exploratory analyses on continuous and categorical variables.

## Answers to Questions

------------------------------------------------------------------------

### 🛠️Coding Practice: Pivot Data Wider

```{r}
df_gapminder %>%
  group_by(continent, year) %>%   
  summarise(mean_life_exp = mean(lifeExp),             
            .groups = "drop") %>%   
  pivot_wider(names_from = year, values_from = mean_life_exp) 
```

### 🛠️Coding Practice: Box Plot Version of Ridge Plot

```{r}
ggplot(df_gapminder, aes(x = lifeExp, y = year, group = year)) +   
  geom_boxplot() 
```

### 🛠️Coding Practice: Life Expectancy Grouped by Year

```{r}
df_gapminder %>%   
  group_by(year) %>%   
  summarise(Minimum = min(lifeExp) %>% round(2),             
            Median  = mean(lifeExp) %>% round(2),             
            Mean    = median(lifeExp) %>% round(2),             
            Maximum = max(lifeExp) %>% round(2)) %>%   
  rename(Year = year) %>%   
  knitr::kable(caption = "Life Expectancy") 
```

### 🛠️Coding Practice: Adding `with()`

```{r}
with(df_2007, table(continent, life_exp_binary)) %>%   
  prop.table(margin=2) 
```
