---
title: |
  ![](images/Logo_Use.png){width=10in}
  
  "R for Data Science Exercises: Exploratory Data Analysis"
format: html
editor: source
author: Lorcán Mason
bibliography: reference.bib
---

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

```

# R for Data Science 2nd Edition Exercises [@wickham2023r]

## Exploratory Data Analysis

### Run the code in your script for the answers! I'm just exploring as I go.

### Packages to load

```{r}
library(tidyverse)
library(nycflights13)
library(ggbeeswarm)
library(lvplot)
library(ggstance)
```

### Introduction

Exploratory Data Analysis (EDA) is an repetitive cycle where you:

1.  Generate questions about your data.

2.  Search for answers by visualizing, transforming, and modelling your data.

3.  Use what you learn to refine your questions and/or generate new questions.

The goal during EDA is to develop an understanding of your data using questions as tools to guide your exploration. Two types of questions will always be useful for making discoveries within your data:

1.  What type of variation occurs within my variables?

2.  What type of covariation occurs between my variables?

By considering these questions, the objectives of EDA can be funnelled down to the following:

-   Explore the *variation* within the variables of your observations.
-   Deal with *outliers* and *missing values* in your data.
-   Explore the *covariation* between the variables of your observations.
-   Recognise how *models* can be used to explore *patterns* in your data.

## Vocabulary

-   **Variable:** a quantity, quality, or property that you can measure.
-   **Value:** the state of a variable when you measure it. This can change.
-   **Observation:** a set of measurements made under similar conditions. One value per variable.
-   **Tabular Data:** observations of variables.
-   **Tidy Data:** 1 observation per row, 1 variable per column, 1 value per cell. Definition of "tidy" for a dataset can depend on the questions you're trying to answer.

## Variation

-   **variation:** the tendency of values of a variable to change between measurements.
-   **categorical variable:** can only take certain values. Visualize variation with bar chart.
-   **continuous variables:** can take on infinite set of ordered values. Visualize variation with histogram.

**Variation** is the tendency of the values of a variable to change from measurement to measurement. You can see variation easily in real life; if you measure any continuous variable twice, you will get two different results. Error is inherent in any data collection process. Variables can vary if you measure across different subjects or at different times. Every variable has its own pattern of variation, which gives interesting information about how much it varies between measurements in the same observation as well as across observations.

### Typical values

Typical values are the most common values of a variable. Think of a bell-curve - the higher the curve the more common a value is; the lower the curve the less-common the values are. The areas not under the curve show values that were not seen in your data. To turn this information into useful questions, look for anything unexpected:

-   Which values are the most common? Why?

-   Which values are rare? Why? Does that match your expectations?

-   Can you see any unusual patterns? What might explain them?

A good starting point for EDA are **Distributions.** **Distributions** can answer some of these questions within the data, while some will require domain expertise. Many of these questions will prompt you to explore a relationship *between* variables, to see if the values of one variable can explain the behaviour of another variable.

### Unusual values

Outliers are observations that are unusual; data points that don't seem to fit the general pattern of the data which are difficult to observe within distributions. Outliers can sometimes be data entry errors, sometimes simply values at the extremes, and other times they suggest important new discoveries.

If you encounter unusual values in your dataset, and simply want to move on to the rest of your analysis, you have two options:

1.  Drop the entire row with the strange values

2.  Replace the unusual values with missing values (`NA`).

    The former is not recommended! A single invalid value doesn't imply that all other values for that observation are also invalid. Additionally, if you have low quality data, applying this approach to every variable might leave you without any data!

You may also want to make distinctions between outliers to discover a pattern, whereby the outlier or missing data actually gives important information.

### Questions

1.  Explore the distribution of each of the `x`, `y`, and `z` variables in `diamonds`. What do you learn? Think about a diamond and how you might decide which dimension is the length, width, and depth.

2.  Explore the distribution of `price`. Do you discover anything unusual or surprising? (Hint: Carefully think about the `binwidth` and make sure you try a wide range of values.)

3.  How many diamonds are 0.99 carat? How many are 1 carat? What do you think is the cause of the difference?

4.  Compare and contrast `coord_cartesian()` vs. `xlim()` or `ylim()` when zooming in on a histogram. What happens if you leave `binwidth` unset? What happens if you try and zoom so only half a bar shows?

5.  What happens to missing values in a histogram? What happens to missing values in a bar chart? Why is there a difference in how missing values are handled in histograms and bar charts?

6.  What does `na.rm = TRUE` do in `mean()` and `sum()`?

7.  Using the `nycflights13::flights`[^1] recreate the frequency plot of `scheduled_dep_time` coloured by whether the flight was cancelled or not.

[^1]: Remember that when we need to be explicit about where a function (or dataset) comes from, we'll use the special form `package::function()` or `package::dataset`.

### Answers

Solution 1:

Calculate summary statistics for these variables and plot their distributions.

```{r}
summary(select(diamonds, x, y, z))

ggplot(diamonds) +
  geom_histogram(mapping = aes(x = x), binwidth = 0.01) # x

ggplot(diamonds) +
  geom_histogram(mapping = aes(x = y), binwidth = 0.01) # y

ggplot(diamonds) +
  geom_histogram(mapping = aes(x = z), binwidth = 0.01) # z
```

There's a few noticeable features of the distributions:

1.  `x` and `y` are larger than `z`,
2.  There are outliers,
3.  They are all right skewed

There are two types of outliers in this data. Some diamonds have values of zero and some have abnormally large values of `x`, `y`, or `z`.

```{r}
summary(select(diamonds, x, y, z))
```

These appear to be either data entry errors, or an undocumented convention in the dataset for indicating missing values. Alternatively it could be that values of zero are the result of rounding values minuscule down, but since there are no diamonds with values of 0.01, that does not seem to be the case.

```{r}
filter(diamonds, x == 0 | y == 0 | z == 0)
```

There are also some diamonds with values of `y` and `z` that are abnormally large. There are diamonds with `y == 58.9` and `y == 31.8`, and one with `z == 31.8`. These are probably data entry errors as they don’t seem in line with the values of the other variables.

```{r}
diamonds %>%
  arrange(desc(y)) %>%
  head() # y
diamonds %>%
  arrange(desc(z)) %>%
  head() # z
```

According to the `diamonds` documentation, `x` is length, `y` is width, and `z` is depth in mm. If documentation were unavailable, you would compare the values of the variables to match them to the length, width, and depth. A quick google search for the definitions of length, width, and depth with respect to diamond cuts will suffice. [Depth](https://en.wikipedia.org/wiki/Diamond_cut) can be expressed as a percentage of the length/width of the diamond, which means it should be less than both the length and the width.

```{r}
summarise(diamonds, mean(x > y), mean(x > z), mean(y > z))
```

It appears that depth (`z`) is always smaller than length (`x`) or width (`y`). Length is greater than width in less than half the observations.

Solution 2:

-   The price data has many spikes, but it is difficult to tell what each spike corresponds to. The following plots don't show much difference in the distributions in the last one or two digits.
-   There are no diamonds with a price of \$1,500 (between \$1,455 and \$1,545, inclusive).
-   There's a swell in the distribution around \$750.

```{r}
ggplot(filter(diamonds, price < 2500), aes(x = price)) +
  geom_histogram(binwidth = 10, center = 0)

ggplot(filter(diamonds), aes(x = price)) +
  geom_histogram(binwidth = 100, center = 0)
```

The last digits of prices are often not uniformly distributed. They are often round, ending in 0 or 5. Another common pattern is in values ending in 99.

```{r}
diamonds %>%
  mutate(ending = price %% 10) %>%
  ggplot(aes(x = ending)) +
  geom_histogram(binwidth = 1, center = 0)

diamonds %>%
  mutate(ending = price %% 100) %>%
  ggplot(aes(x = ending)) +
  geom_histogram(binwidth = 1)

diamonds %>%
  mutate(ending = price %% 1000) %>%
  filter(ending >= 500, ending <= 800) %>%
  ggplot(aes(x = ending)) +
  geom_histogram(binwidth = 1)
```

Solution 3:

There are more than 70 times as many 1 carat diamonds as 0.99 carat diamond.

```{r}
diamonds %>%
  filter(carat >= 0.99, carat <= 1) %>%
  count(carat)
```

Presumably there is a premium for a 1 carat diamond vs. a 0.99 carat diamond beyond the expected increase in price due to a 0.01 carat increase.[^2]

[^2]: Prices for diamonds rise at key weights, such as 0.8, 0.9, and 1 carat. See <https://www.hpdiamonds.com/en-us/extra/136/education-undercarat.htm>.

To check this, look at the number of diamonds in each carat range to see if there is an unusually low number of 0.99 carat diamonds, and an abnormally large number of 1 carat diamonds.

```{r}
diamonds %>%
  filter(carat >= 0.9, carat <= 1.1) %>%
  count(carat) %>%
  print(n = Inf)
```

Solution 4:

-   The `coord_cartesian()` function zooms in on the area specified by the limits, after having calculated and drawn the geoms.
-   The `xlim()` and `ylim()` functions influence actions before the calculation of the stats related to the histogram. Thus, any values outside the x- and y-limits are dropped before calculating bin widths and counts. This can influence how the histogram looks.

```{r}
ggplot(diamonds) +
  geom_histogram(mapping = aes(x = price)) +
  coord_cartesian(xlim = c(100, 5000), ylim = c(0, 3000)) # coord_cartesian

ggplot(diamonds) +
  geom_histogram(mapping = aes(x = price)) +
  xlim(100, 5000) + # xlim
  ylim(0, 3000) # ylim
```

Solution 5: Missing values are removed when the number of observations in each bin are calculated. See the warning message: `Removed 9 rows containing non-finite values (stat_bin)`

```{r}
diamonds2 <- diamonds %>%
  mutate(y = ifelse(y < 3 | y > 20, NA, y))

ggplot(diamonds2, aes(x = y)) +
  geom_histogram()
```

In the `geom_bar()` function, `NA` is treated as another category. The `x` aesthetic in `geom_bar()` requires a discrete (categorical) variable, and missing values act like another category.

```{r}
diamonds %>%
  mutate(cut = if_else(runif(n()) < 0.1, NA_character_, as.character(cut))) %>%
  ggplot() +
  geom_bar(mapping = aes(x = cut))
```

In a histogram, the `x` aesthetic variable needs to be numeric, and `stat_bin()` groups the observations by ranges into bins. Since the numeric value of the `NA` observations is **unknown**, they cannot be placed in a particular bin, and are dropped.

Solution 6:

This option removes `NA` values from the vector prior to calculating the mean and sum.

```{r}
mean(c(0, 1, 2, NA), na.rm = TRUE)
sum(c(0, 1, 2, NA), na.rm = TRUE)
```

Solution 7:

-   Frequency Plot

```{r}
nycflights13::flights |> 
  mutate(
    cancelled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dep_time = sched_hour + (sched_min / 60)
  ) |> 
  ggplot(aes(x = sched_dep_time)) + 
  geom_freqpoly(aes(color = cancelled), binwidth = 1/4)
```

-   Box Plot

```{r}
nycflights13::flights %>%
  mutate(
    cancelled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dep_time = sched_hour + sched_min / 60
  ) %>%
  ggplot() +
  geom_boxplot(mapping = aes(y = sched_dep_time, x = cancelled))
```

## Covariation

If variation describes the behavior *within* a variable, covariation describes the behavior *between* variables. **Covariation** is the tendency for the values of two or more variables to vary together in a related way. The best way to spot covariation is to visualize the relationship between two or more variables.

### Questions

1. Based on EDA, what variable in the diamonds dataset appears to be most important for predicting the price of a diamond? How is that variable correlated with cut? Why does the combination of those two relationships lead to lower quality diamonds being more expensive?

2.  Instead of exchanging the x and y variables, add `coord_flip()` as a new layer to the vertical boxplot to create a horizontal one. How does this compare to exchanging the variables?

3.  One problem with boxplots is that they were developed in an era of much smaller datasets and tend to display a prohibitively large number of "outlying values". One approach to remedy this problem is the letter value plot. Install the lvplot package, and try using `geom_lv()` to display the distribution of price vs. cut. What do you learn? How do you interpret the plots?

4.  Create a visualization of diamond prices vs. a categorical variable from the `diamonds` dataset using `geom_violin()`, then a faceted `geom_histogram()`, then a colored `geom_freqpoly()`, and then a colored `geom_density()`. Compare and contrast the four plots. What are the pros and cons of each method of visualizing the distribution of a numerical variable based on the levels of a categorical variable?

5.  If you have a small dataset, it's sometimes useful to use `geom_jitter()` to avoid overplotting to more easily see the relationship between a continuous and categorical variable. The `ggbeeswarm` package provides a number of methods similar to `geom_jitter()`. List them and briefly describe what each one does.

### Answers

Solution 1:

For investigating the general relationships of each variable with the price of the diamonds, first consider the variables: `carat`, `clarity`, `color`, and `cut`. Ignore the dimensions of the diamond since `carat` measures size, and thus incorporates most of the information contained in these variables.

Since both `price` and `carat` are continuous variables, use a scatter plot to visualise their relationship.

```{r}
ggplot(diamonds, aes(x = carat, y = price)) +
  geom_point()
```

However, since there is a large number of points in the data, use a boxplot by binning `carat`, as suggested in the chapter:

```{r}
ggplot(data = diamonds, mapping = aes(x = carat, y = price)) +
  geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)), orientation = "x")
```

-   Note that the choice of the binning width is important, as if it were too large it would obscure any relationship, and if it were too small, the values in the bins could be too variable to reveal underlying patterns.

The variables `color` and `clarity` are ordered categorical variables. The chapter suggests visualizing a categorical and continuous variable using frequency polygons or boxplots. There is a weak negative relationship between `color` and `price`. The scale of diamond color goes from D (best) to J (worst). Before plotting, reverse the order of the `color` levels so they will be in increasing order of quality on the x-axis.

```{r}
diamonds %>%
  mutate(color = fct_rev(color)) %>%
  ggplot(aes(x = color, y = price)) +
  geom_boxplot()
```

There is also weak negative relationship between `clarity` and `price`. The scale of clarity goes from `I1` (worst) to `IF` (best).

```{r}
ggplot(data = diamonds) +
  geom_boxplot(mapping = aes(x = clarity, y = price))
```

For both `clarity` and `color`, there is a much larger amount of variation within each category than between categories. `Carat` is clearly the single best predictor of diamond prices. Since this is an example of a continuous (`carat`) and categorical (`cut`) variable, it can be visualized with a box plot.

```{r}
ggplot(diamonds, aes(x = cut, y = carat)) +
  geom_boxplot()
```

There is a lot of variability in the distribution of `carat` sizes within each cut category. There is a slight negative relationship between `carat` and `cut.` Noticeably, the largest `carat` diamonds have a cut of "Fair" (the lowest).

This negative relationship can be due to the way in which diamonds are selected for sale. A larger diamond can be profitably sold with a lower quality `cut`, while a smaller diamond requires a better `cut.`

Solution 2: 

First, recreate this horizontal box plot of the distribution `hwy` by `class`, using `geom_boxplot()` and `coord_flip()`:   
```{r}
ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y = hwy)) +
  coord_flip()
```

In this case the output looks the same, but `x` and `y` aesthetics are flipped.
```{r}
library("ggstance")

ggplot(data = mpg) +
  geom_boxploth(mapping = aes(y = reorder(class, hwy, FUN = median), x = hwy))
```

In `ggplot2` (since [version 3.3.0](https://ggplot2.tidyverse.org/news/index.html#new-features)) all geoms can choose the direction. The direction is inferred from the aesthetic mapping. In this case, switching `x` and `y` produces a horizontal boxplot.
```{r}
ggplot(data = mpg) +
  geom_boxplot(mapping = aes(y = reorder(class, hwy, FUN = median), x = hwy))
```

The `orientation` argument is used to explicitly specify the axis orientation of the plot.
```{r}
ggplot(data = mpg) +
  geom_boxplot(mapping = aes(y = reorder(class, hwy, FUN = median), x = hwy), orientation = "y")
```

Solution 3: 
Like box-plots, the boxes of the letter-value plot correspond to quantiles. However, they incorporate
far more quantiles than box-plots. They are useful for larger datasets because,

1.  larger datasets can give precise estimates of quantiles beyond the quartiles
2.  larger datasets should have more outliers (in absolute numbers).

```{r}
ggplot(diamonds, aes(x = cut, y = price)) +
  geom_lv()
```

Solution 4: 

The `geom_freqpoly()` is better for look-up: meaning that given a price, it is easy to tell which `cut` has the highest density. However, the overlapping lines makes it difficult to distinguish how the overall distributions relate to each other. The `geom_violin()` and faceted `geom_histogram()` have similar strengths and weaknesses. It is easy to visually distinguish differences in the overall shape of the distributions (skewness, central values, variance, etc). However, since we can't easily compare the vertical values of the distribution, it is difficult to look up which category has the highest density for a given price. All of these methods depend on tuning parameters to determine the level of smoothness of the distribution.

```{r}
ggplot(data = diamonds, mapping = aes(x = price, y = ..density..)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500) # geom_freqpoly

ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_histogram() + # geom_histogram
  facet_wrap(~cut, ncol = 1, scales = "free_y")

ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
  geom_violin() + # geom_violin
  coord_flip()
```


Solution 5:

There are two methods:

-   `geom_quasirandom()` produces plots that are a mix of jitter and violin plots. There are several different methods that determine exactly how the random location of the points is generated.
-   `geom_beeswarm()` produces a plot similar to a violin plot, but by offsetting the points.

Use the `mpg`  box plot example since these methods display individual points, they are better suited for smaller datasets.

```{r}
ggplot(data = mpg) +
  geom_quasirandom(mapping = aes(
    x = reorder(class, hwy, FUN = median),
    y = hwy
  ))

ggplot(data = mpg) +
  geom_quasirandom(
    mapping = aes(
      x = reorder(class, hwy, FUN = median),
      y = hwy
    ),
    method = "tukey"
  )

ggplot(data = mpg) +
  geom_quasirandom(
    mapping = aes(
      x = reorder(class, hwy, FUN = median),
      y = hwy
    ),
    method = "tukeyDense"
  )

ggplot(data = mpg) +
  geom_quasirandom(
    mapping = aes(
      x = reorder(class, hwy, FUN = median),
      y = hwy
    ),
    method = "frowney"
  )

ggplot(data = mpg) +
  geom_quasirandom(
    mapping = aes(
      x = reorder(class, hwy, FUN = median),
      y = hwy
    ),
    method = "smiley"
  )

ggplot(data = mpg) +
  geom_beeswarm(mapping = aes(
    x = reorder(class, hwy, FUN = median),
    y = hwy
  ))
```


## Two Categorical Variables

### Questions

1.  How could you rescale the count dataset above to more clearly show the distribution of cut within color, or color within cut?

2.  What different data insights do you get with a segmented bar chart if color is mapped to the `x` aesthetic and `cut` is mapped to the `fill` aesthetic?
    Calculate the counts that fall into each of the segments.

3.  Use `geom_tile()` together with dplyr to explore how average flight departure delays vary by destination and month of year.
    What makes the plot difficult to read?
    How could you improve it?

### Answers

Solution 1:

To clearly show the distribution of `cut` within `color`, calculate a new variable `prop` which is the proportion of each cut within a `color`. This is done using a grouped mutate.
```{r}

diamonds %>%
  count(color, cut) %>%
  group_by(color) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(mapping = aes(fill = prop))
```

Similarly, to scale by the distribution of `color` within `cut`,
```{r}
diamonds %>%
  count(color, cut) %>%
  group_by(cut) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(mapping = aes(fill = prop))
```

Add `limit = c(0, 1)` to put the color scale between (0, 1). These are the logical boundaries of proportions. This makes it possible to compare each cell to its actual value, and would improve comparisons across multiple plots. However, it ends up limiting the colors and makes it harder to compare within the dataset. However, using the default limits of the minimum and maximum values makes it easier to compare within the dataset the emphasizing relative differences, but harder to compare across datasets.

Solution 2:

The scale of diamond color goes from D (best) to J (worst). Colourless diamonds are considered better than the ones with yellowish-brownish tint. This dataset contains diamonds of 7 colours- “D” to “J”. D, E, F are colourless & G-J have a very faint colour. The code below calculates the `count` of each `cut` within each `color`.

```{r}
diamonds %>%
  count(color, cut) %>%
  group_by(color)
```

To show this visually, execute the below code. This will show the proportion of each `cut` in each `color`. `position="dodge"` creates the segmented bars to display the `count` of each `cut` in each `color` category.
```{r}
diamonds %>%
  ggplot(mapping = aes(fill=cut, x=color)) +
  geom_bar(position="dodge")
```

Solution 3:

```{r}
flights %>%
  group_by(month, dest) %>%
  summarise(dep_delay = mean(dep_delay, na.rm = TRUE)) %>%
  ggplot(aes(x = factor(month), y = dest, fill = dep_delay)) +
  geom_tile() +
  labs(x = "Month", y = "Destination", fill = "Departure Delay")
```

There are several things that could be done to improve it,

-   sort destinations by a meaningful quantity (distance, number of flights, average delay)
-   remove missing values

How to treat missing values is difficult. In this case, missing values correspond to airports which don't have regular flights (at least one flight each month) from NYC. These are likely smaller airports (with higher variance in their average due to fewer observations). When we group all pairs of (`month`, `dest`) again by destination, we should have a total count of 12 (one for each month) per group (`dest`). This makes it easy to filter.

```{r}
flights %>%
  group_by(month, dest) %>%                                 # This gives (month, dest) pairs
  summarise(dep_delay = mean(dep_delay, na.rm = TRUE)) %>%
  group_by(dest) %>%                                        # group all (month, dest) pairs by destination ..
  filter(n() == 12) %>%                                     # and only select those that have one entry per month 
  ungroup() %>%
  mutate(dest = reorder(dest, dep_delay)) %>%
  ggplot(aes(x = factor(month), y = dest, fill = dep_delay)) +
  geom_tile() +
  labs(x = "Month", y = "Destination", fill = "Departure Delay")
```


## Two Numerical Variables

### Questions

1.  Instead of summarizing the conditional distribution with a boxplot, you could use a frequency polygon.
    What do you need to consider when using `cut_width()` vs. `cut_number()`?
    How does that impact a visualization of the 2d distribution of `carat` and `price`?

2.  Visualize the distribution of `carat`, partitioned by `price`.

3.  How does the price distribution of very large diamonds compare to small diamonds?
    Is it as you expect, or does it surprise you?

4.  Combine two of the techniques you've learned to visualize the combined distribution of cut, carat, and price.

5.  Two dimensional plots reveal outliers that are not visible in one dimensional plots.
    For example, some points in the following plot have an unusual combination of `x` and `y` values, which makes the points outliers even though their `x` and `y` values appear normal when examined separately.
    Why is a scatterplot a better display than a binned plot for this case?

    ```{r}
    #| eval: false
    diamonds |> 
      filter(x >= 4) |> 
      ggplot(aes(x = x, y = y)) +
      geom_point() +
      coord_cartesian(xlim = c(4, 11), ylim = c(4, 11))
    ```

6.  Instead of creating boxes of equal width with `cut_width()`, we could create boxes that contain roughly equal number of points with `cut_number()`.
    What are the advantages and disadvantages of this approach?

    ```{r}
    #| eval: false
    ggplot(diamonds, aes(x = carat, y = price)) + 
      geom_boxplot(aes(group = cut_number(carat, 20)))
    ```

### Answers

Solution 1:

Both `cut_width()` and `cut_number()` split a variable into groups. When using `cut_width()`, we need to choose the width, and the number of bins will be calculated automatically. When using `cut_number()`, we need to specify the number of bins, and the widths will be calculated automatically.

In either case, choose the bin widths and number large enough to aggregate observations and remove noise, but not so large as to remove all the signal/data.

If categorical colors are used, no more than eight colors should be used in order to keep them distinct. Using `cut_number`, split carats into quantiles (five groups).

```{r message=FALSE}
ggplot(
  data = diamonds,
  mapping = aes(color = cut_number(carat, 5), x = price)
) +
  geom_freqpoly() +
  labs(x = "Price", y = "Count", color = "Carat")
```

Alternatively, you can use `cut_width` to specify widths at which to cut. Choose 1-carat widths. Since there are very few diamonds larger than 2-carats, this is not as informative. However, using a width of 0.5 carats creates too many groups, and splitting at non-whole numbers is unappealing.

```{r message=FALSE}
ggplot(
  data = diamonds,
  mapping = aes(color = cut_width(carat, 1, boundary = 0), x = price)
) +
  geom_freqpoly() +
  labs(x = "Price", y = "Count", color = "Carat")
```

Solution 2:

Plotted with a box plot with 10 bins with an equal number of observations, and the width determined by the number of observations.
```{r}
ggplot(diamonds, aes(x = cut_number(price, 10), y = carat)) +
  geom_boxplot() +
  coord_flip() +
  xlab("Price")
```

Plotted with a box plot with 10 equal-width bins of \$2,000. The argument `boundary = 0` ensures that first bin is \$0--\$2,000.
```{r}
ggplot(diamonds, aes(x = cut_width(price, 2000, boundary = 0), y = carat)) +
  geom_boxplot(varwidth = TRUE) +
  coord_flip() +
  xlab("Price")
```

Solution 3:

The distribution of very large diamonds is more variable. This may be due to the way in which diamonds are selected for retail sales. For example, someone selling a diamond only finds it profitable to sell it if some combination size, cut, clarity, and color are above a certain threshold. The smallest diamonds are only profitable to sell if they are exceptional in all the other factors (cut, clarity, and color), so the small diamonds sold have similar characteristics. However, larger diamonds may be profitable regardless of the values of the other factors. Thus the data will show large diamonds with a wider variety of cut, clarity, and color and thus more variability in prices.

Solution 4:

There are many options to try, so solutions may vary. Here are a few options that work.

```{r}
ggplot(diamonds, aes(x = carat, y = price)) +
  geom_hex() +
  facet_wrap(~cut, ncol = 1) # 1 Facets by cut
```

```{r}
ggplot(diamonds, aes(x = cut_number(carat, 5), y = price, colour = cut)) +
  geom_boxplot() # 2 Price of cuts using boxplots
```

```{r}
ggplot(diamonds, aes(colour = cut_number(carat, 5), y = price, x = cut)) +
  geom_boxplot() # 3 Price of carats using boxplots
```

Solution 5:

```{r}
ggplot(data = diamonds) +
  geom_point(mapping = aes(x = x, y = y)) +
  coord_cartesian(xlim = c(4, 11), ylim = c(4, 11))
```

In this case, there is a strong relationship between `x` and `y`. The outliers in this case are not extreme in either `x` or `y`. A binned plot would not reveal these outliers, and may lead us to conclude that the largest value of `x` was an outlier even though it appears to fit the bivariate pattern well.

Solution 6:

`cut_width(x, width)`, as used above in Solution 1, divides `x` (the dataset) into bins of width `width`.
By default, boxplots look roughly the same (apart from number of outliers) regardless of how many observations there are, so it's difficult to tell that each boxplot summaries a different number of points. One way to show that is to make the width of the boxplot proportional to the number of points with `varwidth = TRUE`.

## Patterns and Models

If a systematic relationship exists between two variables it will appear as a pattern in the data. If you spot a pattern, ask yourself:

-   Could this pattern be due to coincidence (i.e. random chance)?

-   How can you describe the relationship implied by the pattern?

-   How strong is the relationship implied by the pattern?

-   What other variables might affect the relationship?

-   Does the relationship change if you look at individual subgroups of the data?

Patterns in your data provide clues about relationships, i.e., they reveal covariation. If you think of variation as a phenomenon that creates uncertainty, covariation is a phenomenon that reduces it. If two variables covary, you can use the values of one variable to make better predictions about the values of the second. If the covariation is due to a causal relationship (a special case), then you can use the value of one variable to control the value of the second. Models are a tool for extracting patterns out of data. For example, consider the diamonds data. It's hard to understand the relationship between cut and price, because cut and carat, and carat and price are tightly related. It's possible to use a model to remove the very strong relationship between price and carat so we can explore the subtleties that remain.