<style type="text/css">
p { text-align: left; }
</style>

<style>
    .reveal .slides {
        text-align: left;
    }
    .reveal .slides section>* {
        margin-left: 0;
        margin-right: 0;
    }
    
    .reveal .slide-number {
    font-size: 24pt;
    color: #ffffff;
</style>

### Data science fundamentals 03: 
Exploratory data analysis, Part 1 (Aug. 5th, 2019)

---

In this workshop we will do exploratory data analysis (EDA) with the HANES data set. Import the data set and change data types to reflect their correct annotation (housekeeping).


```{r eval=FALSE, message=FALSE, echo = TRUE}
  # Load the package RCurl
  library(RCurl)
  # Import the HANES data set from GitHub; break the string into two for readability
  # (Please note this readability aspect very carefully)
  URL_text_1 <- "https://raw.githubusercontent.com/kannan-kasthuri/kannan-kasthuri.github.io"
  URL_text_2 <- "/master/Datasets/HANES/NYC_HANES_DIAB.csv"
  # Paste it to constitute a single URL 
  URL <- paste(URL_text_1,URL_text_2, sep="")
  HANES <- read.csv(text=getURL(URL))
  # Rename the GENDER factor for identification
  HANES$GENDER <- factor(HANES$GENDER, labels=c("M","F"))
  # Rename the AGEGROUP factor for identification
  HANES$AGEGROUP <- factor(HANES$AGEGROUP, labels=c("20-39","40-59","60+"))
  # Rename the HSQ_1 factor for identification
  HANES$HSQ_1 <- factor(HANES$HSQ_1, labels=c("Excellent","Very Good","Good", "Fair", "Poor"))
  # Rename the DX_DBTS as a factor
  HANES$DX_DBTS <- factor(HANES$DX_DBTS, labels=c("DIAB","DIAB NO_DX","NO DIAB"))
  # Omit all NA from the data frame
  HANES <- na.omit(HANES)
  # Observe the structure
  str(HANES)
  # Load the tidyverse library
  library(tidyverse)
```

##

```{r eval=TRUE, message=FALSE, echo = FALSE}
  # Load the package RCurl
  library(RCurl)
  # Import the HANES data set from GitHub; break the string into two for readability
  # (Please note this readability aspect very carefully)
  URL_text_1 <- "https://raw.githubusercontent.com/kannan-kasthuri/kannan-kasthuri.github.io"
  URL_text_2 <- "/master/Datasets/HANES/NYC_HANES_DIAB.csv"
  # Paste it to constitute a single URL 
  URL <- paste(URL_text_1,URL_text_2, sep="")
  HANES <- read.csv(text=getURL(URL))
  # Rename the GENDER factor for identification
  HANES$GENDER <- factor(HANES$GENDER, labels=c("M","F"))
  # Rename the AGEGROUP factor for identification
  HANES$AGEGROUP <- factor(HANES$AGEGROUP, labels=c("20-39","40-59","60+"))
  # Rename the HSQ_1 factor for identification
  HANES$HSQ_1 <- factor(HANES$HSQ_1, labels=c("Excellent","Very Good","Good", "Fair", "Poor"))
  # Rename the DX_DBTS as a factor
  HANES$DX_DBTS <- factor(HANES$DX_DBTS, labels=c("DIAB","DIAB NO_DX","NO DIAB"))
  # Omit all NA from the data frame
  HANES <- na.omit(HANES)
  # Observe the structure
  str(HANES)
  # Load the tidyverse library
  library(tidyverse)
```

##

EDA is informally a set of procedures where we try to obtain directed insights about the data. They may comprise of the following steps:

1. Generate questions about the data.

2. Search for answers by visualising, transforming, and modelling the data.

3. Use the information to generate new questions.

##

#### Visualizing variations

It is natural for the same data to vary between measurements. 

There are several types of variations, including batch effects. Variations gives raise to distributions. 

There are ways to visualize distributions and they may depend on whether the distribution has categorical values or continuous. 

##

To examine the distribution of a categorical variable, we can use a `geom_bar()` function:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of people with diabetes diagnostic status 
  ggplot(data = HANES) +
    geom_bar(mapping = aes(x = DX_DBTS))
```
 
##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
   # Find the distribution of people with diabetes diagnostic status
  ggplot(data = HANES) +
    geom_bar(mapping = aes(x = DX_DBTS))
```

##

Manually, these values can be found with the `count()` function.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Manually find the counts
  HANES %>% count(DX_DBTS)
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
   # Manually find the counts
  HANES %>% count(DX_DBTS)
```

##

For continous variables (such as GLUCOSE, for example), we can use the histogram geom:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of the GLUCOSE levels in patients
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = GLUCOSE))
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
   # Find the distribution of the GLUCOSE levels in patients
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = GLUCOSE)) 
```

##

Multiple histogram plots can be overlayed by the geom `freqpoly()`.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of the GLUCOSE levels in different age groups
  ggplot(data = HANES, mapping = aes(x = GLUCOSE, colour = AGEGROUP)) +
    geom_freqpoly()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of the GLUCOSE levels in different age groups
  ggplot(data = HANES, mapping = aes(x = GLUCOSE, colour = AGEGROUP)) +
    geom_freqpoly() 
```

---

#### Typical values

In both bar charts and histograms, tall bars show the common values of a variable, and shorter bars show less-common values. Places that do not have bars reveal values that are not seen in this data. To turn this information into useful questions, we have to look for unexpected things:

1. Which values are the most common? Why? Are they typical?

2. Which values are rare? Why? Does that match your expectations?

3. Can we see any unusual patterns? What might explain them?

##

As an example, the histogram below suggests an interesting question:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of the HDL levels in different health status
  ggplot(data = HANES, mapping = aes(x = HDL, colour = HSQ_1)) +
    geom_freqpoly()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of the HDL levels in different health status
  ggplot(data = HANES, mapping = aes(x = HDL, colour = HSQ_1)) +
    geom_freqpoly() 
```

##

Why poor health has lesser mean HDL?

Consider the following histogram:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of log(UACR) for different age groups
  ggplot(data = HANES, mapping = aes(x = log(UACR), colour = AGEGROUP)) +
    geom_freqpoly()
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of log(UACR) for different age groups
  ggplot(data = HANES, mapping = aes(x = log(UACR), colour = AGEGROUP)) +
    geom_freqpoly() 
```

##

Bimodality suggests that subgroups exist in our data. To understand the subgroups, we need to ask:

1. How are the observations within each cluster similar to each other, in terms of other variables?

2. How are the observations in separate clusters different from each other, in terms of other variables?

3. Do the clusters have biological significance?

4. Why might the appearance of clusters be misleading?

---

#### Unusual values and outlier detection

Outliers can be extremely difficult to spot. 

Sometimes outliers are data entry errors; other times outliers suggest important new science. 

When we have a lot of data, outliers are sometimes difficult to see in a histogram. 

##

For example, take the distribution of the CREATININE variable from our HANES dataset. 

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of the CREATININE levels in patients
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = COTININE))
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of the CREATININE levels in patients
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = CREATININE))  
```

##

We can hardly spot any outlier. 

The only evidence of outliers is the unusually wide limits on the x-axis. 

However, zooming the x-axis by restricting the limits between 0 and 0.5 we see disconnected values suggesting the presense of outliers.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of the CREATININE levels in patients
  # restricted to the limits 0 to 0.5
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = COTININE)) + xlim(c(0,0.5)) 
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of the CREATININE levels in patients
  # restricted to the limits 0 to 0.5
  ggplot(data = HANES) +
    geom_histogram(mapping = aes(x = CREATININE)) + xlim(c(0,0.5)) 
```

##

Once we have visually inspected the outliers, we can isolate them out by data transformation tools discussed in the previous lecture.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # From the HANES data set
  unusual <- HANES %>% 
  # filter people with CREATININE levels between 0.42 and 0.37
  filter(CREATININE <= 0.42 & CREATININE > 0.375) %>% 
  # select only the people id (key) and gender
  select(CREATININE, KEY, GENDER) %>%
  # and order CREATININE as a first variable
  arrange(CREATININE)
  # And display it as a tibble
  as.tibble(unusual) 
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # From the HANES data set
  unusual <- HANES %>% 
  # filter people with CREATININE levels between 0.42 and 0.37
  filter(CREATININE <= 0.42 & CREATININE > 0.375) %>% 
  # select only the people id (key) and gender
  select(CREATININE, KEY, GENDER) %>%
  # and order CREATININE as a first variable
  arrange(CREATININE)
  # And display it as a tibble
  as.tibble(unusual)  
```

We see that all of these people are females.

---

### Exploratory data analysis, Part 2


#### Missing values

It may be tempting to ignore unsual values, but it is best to replace unusual values by `NA`. This can be done by `ifelse()` function:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE, na.rm=FALSE}
  # Input HANSES
  HANES_repl_unusual <- HANES %>% 
  # Replace using ifelse
  mutate(CREATININE = ifelse(CREATININE <= 0.42 & CREATININE > 0.375, NA, CREATININE))
```

##

`ifelse()` has three arguments. 

The first argument test should be a logical vector. 

The result will contain the value of the second argument, yes, when test is TRUE, and the value of the third argument, no, when it is false.

Sometimes `NA` could just mean no other option which we may want to compare with the rest of the options. 

For example, the HANES data set was collected at different times. Some of them who did not have their A1C levels measured may have GLUCOSE levels measured. 

We may want to compare, say, GLUCOSE levels of people with no A1C information with the levels for the people with A1C information. 

##

(Since we omited `NA` from HANES, we will re-import the data set):

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=TRUE}
# Load the package RCurl
  library(RCurl)
  # Import the HANES data set from GitHub; break the string into two for readability
  # (Please note this readability aspect very carefully)
  URL_text_1 <- "https://raw.githubusercontent.com/kannan-kasthuri/kannan-kasthuri.github.io"
  URL_text_2 <- "/master/Datasets/HANES/HANES.original.csv"
  # Paste it to constitute a single URL 
  URL <- paste(URL_text_1,URL_text_2, sep="")
  HANES_with_NA <- read.csv(text=getURL(URL))
```

##

And here is the code to do this -

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # From the HANES_with_NA dataset
  HANES_with_NA %>% 
  # Find the glucose levels of people with no A1C through the
  # logical operator/function is.na()
  mutate(no_A1C = is.na(A1C),
         glucose_level = GLUCOSE) %>% 
  # Now plot the result for the logical categories
  ggplot(mapping = aes(glucose_level)) + 
    geom_freqpoly(mapping = aes(colour = no_A1C), binwidth = 1/4)
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # From the HANES_with_NA dataset
  HANES_with_NA %>% 
  # Find the glucose levels of people with no A1C through the
  # logical operator/function is.na()
  mutate(no_A1C = is.na(A1C),
         glucose_level = GLUCOSE) %>% 
  # Now plot the result for the logical categories
  ggplot(mapping = aes(glucose_level)) + 
    geom_freqpoly(mapping = aes(colour = no_A1C), binwidth = 1/4)
```

We immediately infer that there are only a few such cases (to our delight)!

---

#### Covariation

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 visualise the relationship between two or more variables. 

First we will consider the variation between continuous and categorical variables.

##

##### Variation between continuous and categorical variable

One can plot the frequency distribution between continuous and categorical variable. For example, in the HANES data set, let us plot the frequency distribution of mercury in urine (MERCURYU variable) for different age groups. 

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot the frequency distribution of mercury for different age groups
  ggplot(data = HANES, mapping = aes(x = MERCURYU)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot the frequency distribution of mercury for different age groups
  ggplot(data = HANES, mapping = aes(x = MERCURYU)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

##

This plot that says people in age groups 20-39 is the largest group with average mercury in urine, followed by age group 40-59 and then by 60+, is misleading because the number of people in the group might differ. 

##

This is the case as we can see:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Find the distribution of people in each age group
  ggplot(HANES) + 
    geom_bar(mapping = aes(x = AGEGROUP))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Find the distribution of people in each age group
  ggplot(HANES) + 
    geom_bar(mapping = aes(x = AGEGROUP))
```

Number of people in 20-39 > # of people in 40-59 > # of people in 60+. 

##

To make the comparison easier we need to change what is displayed on the y-axis. 

Instead of displaying count, we’ll display **density**, which is the count standardised so that the area under each frequency polygon is one.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot density instead of count for MERCURYU variable
  ggplot(data = HANES, mapping = aes(x = MERCURYU, y = ..density..)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

##
  
```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot density instead of count for MERCURYU variable
  ggplot(data = HANES, mapping = aes(x = MERCURYU, y = ..density..)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

We can right away see that all groups have equal number of people with average mercury-in-urine levels. 

##

In contrast, if we plot the density of blood hemoglobin (A1C) levels for these age groups, we see:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot density instead of count for A1C variable
  ggplot(data = HANES, mapping = aes(x = A1C, y = ..density..)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot density instead of count for A1C variable
  ggplot(data = HANES, mapping = aes(x = A1C, y = ..density..)) + 
    geom_freqpoly(mapping = aes(colour = AGEGROUP), binwidth = 1)
```

the group 20-39 has the smallest variation in blood hemoglobin levels.

##

Another alternative to display the distribution of a continuous variable broken down by a categorical variable is the **boxplot**. 

A boxplot is a type of visual shorthand for a distribution of values that is popular among statisticians. 

##

Each boxplot consists of:

* A box that stretches from the 25th percentile of the distribution to the 75th percentile, a distance known as the interquartile range (IQR). 

In the middle of the box is a line that displays the median, i.e. 50th percentile, of the distribution. 

These three lines give you a sense of the spread of the distribution and whether or not the distribution is symmetric about the median or skewed to one side.

##

* Visual points that display observations that fall more than 1.5 times the IQR from either edge of the box. These outlying points are unusual so are plotted individually.

* A line (or whisker) that extends from each end of the box and goes to the farthest non-outlier point in the distribution.

---

![](EDA-boxplot.png)

---

##

We can look at the variation between A1C and AGEGROUP through a box plot.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Make a box plot for A1C and AGEGROUP
  ggplot(data = HANES, mapping = aes(x = AGEGROUP, y = A1C)) +
    geom_boxplot()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Make a box plot for A1C and AGEGROUP
  ggplot(data = HANES, mapping = aes(x = AGEGROUP, y = A1C)) +
    geom_boxplot()
```

Thus the box plot confirms our observation that the group 20-39 has the smallest variation in blood hemoglobin levels. 

##

We note that AGEGROUP is an ordinal (meaning which can be ordered). 

Some categorical variables such as health status (HSQ_1) cannot be ordered and understanding the trend might be difficult. 

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Make a box plot for HDL and HSQ_1
  ggplot(data = HANES, mapping = aes(x = HSQ_1, y = HDL)) +
    geom_boxplot()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Make a box plot for HDL and HSQ_1
  ggplot(data = HANES, mapping = aes(x = HSQ_1, y = HDL)) +
    geom_boxplot()
```

##

For such comparisons we can make use of `reorder` argument, say, based on median:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Make a box plot for HDL and HSQ_1 based on reordering by median values
  ggplot(data = HANES, mapping = aes(x = reorder(HSQ_1, HDL, FUN = median), y = HDL)) +
    geom_boxplot()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Make a box plot for HDL and HSQ_1 based on reordering by median values
  ggplot(data = HANES, mapping = aes(x = reorder(HSQ_1, HDL, FUN = median), y = HDL)) +
    geom_boxplot()
```

Such reordering helps us to see that average HDL increases as the health status becomes better.

##

If we have long variable names, `geom_boxplot()` will work better if you flip it 90°. We can do that with `coord_flip()`.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Make a box plot for HDL and HSQ_1 based on reordering by median values
  # and flip coordinates
  ggplot(data = HANES, mapping = aes(x = reorder(HSQ_1, HDL, FUN = median), y = HDL)) +
    geom_boxplot() + coord_flip()
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Make a box plot for HDL and HSQ_1 based on reordering by median values
  # and flip coordinates
  ggplot(data = HANES, mapping = aes(x = reorder(HSQ_1, HDL, FUN = median), y = HDL)) +
    geom_boxplot() + coord_flip()
```

##

##### Variation between two categorical variables

To visualise the covariation between categorical variables we can rely on the built-in `geom_count()` function:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot two categorical variables AGEGROUP and diabetes status DX_DBTS
  # using geom_count
  ggplot(data = HANES) +
    geom_count(mapping = aes(x = AGEGROUP, y = DX_DBTS))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot two categorical variables AGEGROUP and diabetes status DX_DBTS
  # using geom_count
  ggplot(data = HANES) +
    geom_count(mapping = aes(x = AGEGROUP, y = DX_DBTS))
```

##

The size of each circle in the plot displays how many observations occurred at each combination of values. 

Covariation will appear as a strong correlation between specific x values and specific y values.

We can make a more continuous map using the `geom_tile()` function:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Count the people in age group and health status
  HANES %>% 
  count(AGEGROUP, HSQ_1) %>%  
  ggplot(mapping = aes(x = AGEGROUP, y = HSQ_1)) +
  # and use the geom_tile to fill in the aesthetics
    geom_tile(mapping = aes(fill = n))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Count the people in age group and health status
  HANES %>% 
  count(AGEGROUP, HSQ_1) %>%  
  ggplot(mapping = aes(x = AGEGROUP, y = HSQ_1)) +
  # and use the geom_tile to fill in the aesthetics
    geom_tile(mapping = aes(fill = n))
```

##

##### Variation between two continuous variables

Covariation between two continuous variables can be viewed using `geom_point()` function: 

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL
  ggplot(data = HANES) +
    geom_point(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL
  ggplot(data = HANES) +
    geom_point(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))
```

We see the near linear relationship of LDL and cholesterol levels in people. 

##

Scatterplots become less useful as the size of the dataset grows, because points begin to overplot, and pile up into areas of uniform black (as above). 

One way to fix the problem is using the alpha aesthetic to add transparency.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # adding alpha aesthetics
  ggplot(data = HANES) +
    geom_point(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL), alpha = 1 / 5)
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # adding alpha aesthetics
  ggplot(data = HANES) +
    geom_point(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL), alpha = 1 / 5)
```

##

But using transparency can be challenging for very large datasets. Another solution is to use bin. 

`geom_bin2d()` and `geom_hex()` divide the coordinate plane into 2d bins and then use a fill color to display how many points fall into each bin. 

`geom_bin2d()` creates rectangular bins. 

`geom_hex()` creates hexagonal bins. 

##

**Note**: We need to install hexbin package to use `geom_hex()`.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # using geom_bin2d
  ggplot(data = HANES) +
    geom_bin2d(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))

  # install.packages("hexbin")
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # using geom_hex
  ggplot(data = HANES) +
    geom_hex(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # using geom_bin2d
  ggplot(data = HANES) +
    geom_bin2d(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # install.packages("hexbin")
  library(hexbin)
  # Plot the variation between LDLESTIMATE and CHOLESTEROLTOTAL 
  # using geom_hex
  ggplot(data = HANES) +
    geom_hex(mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL))
```

##

Another option is to bin one continuous variable so it acts like a categorical variable. 

Then we can use one of the techniques for visualising the combination of a categorical and a continuous variable that we learned. 

##

For example, we could bin LDLESTIMATE and then for each group, display a boxplot with variable width representing the number of points in the category. 

This can be accomplished by the argument `cut_width()` available as an aesthetic along with the `varwidth` argument available in the geom `geom_boxplot`. 

`cut_width(x, width)` divides x into bins of width `width` while `varwidth = TRUE` defines variable width for the box plot depending on the count of the points in the bin.

##

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Bin LDLESTIMATE using cut_width argument and do a box plot like we 
  # did between categorical and continuous variable allowing the size 
  # of the box plot represent the number of points - varwidth argument
  ggplot(data = HANES, mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL)) + 
    geom_boxplot(mapping = aes(group = cut_width(LDLESTIMATE, 20)), varwidth = TRUE)
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Bin LDLESTIMATE using cut_width argument and do a box plot like we 
  # did between categorical and continuous variable allowing the size 
  # of the box plot represent the number of points - varwidth argument
  ggplot(data = HANES, mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL)) + 
    geom_boxplot(mapping = aes(group = cut_width(LDLESTIMATE, 20)), varwidth = TRUE)
```

##

Another approach would be to display the same number of points for each bin which can be defined by the argument `cut_number`.

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Use the same number of points for LDLESTIMATE using cut_number argument 
  # and do a box plot like we did between categorical and continuous variable 
  ggplot(data = HANES, mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL)) + 
    geom_boxplot(mapping = aes(group = cut_number(LDLESTIMATE, 20)))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Use the same number of points for LDLESTIMATE using cut_number argument 
  # and do a box plot like we did between categorical and continuous variable 
  ggplot(data = HANES, mapping = aes(x = LDLESTIMATE, y = CHOLESTEROLTOTAL)) + 
    geom_boxplot(mapping = aes(group = cut_number(LDLESTIMATE, 20)))
```

---

#### Patterns and models

Relationships in data can be determined when we observe patterns. Systematic relationship between two variables will result in a pattern. The following are key questions to ask when we observe patterns:

1. Could the pattern be due to random chance?

2. If not how can we describe the relationship implied by the pattern?

3. How strong is the relationship implied by the pattern?

4. What other variables might affect the relationship?

5. Does the relationship change if we look at individual subgroups of the data?

##

For example, we already saw:

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Make a ggplot out of log(A1C) and log(UACR) variables
  ggplot(data = HANES) + 
    geom_point(mapping = aes(x = log(A1C), y = log(UACR)))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Make a ggplot out of log(A1C) and log(UACR) variables
  ggplot(data = HANES) + 
    geom_point(mapping = aes(x = log(A1C), y = log(UACR)))
```

##
Patterns provide one of the most useful tools for data scientists because they reveal covariation. 

If we think of variation as a phenomenon that creates uncertainty, covariation is a phenomenon that reduces it. 

If two variables covary, we 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 we can use the value of one variable to control the value of the second.

Models are a tool for extracting patterns out of data which we will cover in more detail later.

##

For example, in HANES, it’s hard to understand the relationship between health status (HSQ_1) and total cholesterol, because health status and LDL cholesterol, and LDL cholesterol and total cholesterol are tightly related. 

It’s possible to use a model to remove the very strong relationship between total cholesterol and LDL cholesterol so we can explore the subtleties that remain. The following code fits a model that predicts total cholesterol from LDL cholesterol and then computes the residuals (the difference between the predicted value and the actual value). 

The residuals give us a view of the total cholesterol of the patients, once the effect of LDL cholesterol has been removed.

##

```{r eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
  # Load the modelr library
  library(modelr)
  # Compute the model using the `lm` - linear regression function
  mod <- lm(CHOLESTEROLTOTAL ~ LDLESTIMATE , data = HANES)
  # add the residuals to the data and name it as resid
  HANES_mod <- HANES %>% 
    add_residuals(mod) %>% 
      mutate(resid = resid)
  # Plot the total cholesterol and residuals
  ggplot(data = HANES_mod) + 
    geom_point(mapping = aes(x = CHOLESTEROLTOTAL, y = resid))
  # Boxplot the health status and residuals that gives the 
  # effect of total cholesterol removing the effects of LDL
  # and reorder by the median
  ggplot(data = HANES_mod) + 
    geom_boxplot(mapping = aes(x = reorder(HSQ_1, resid, FUN = median), y = resid))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}
  # Load the modelr library
  library(modelr)
  # Compute the model using the `lm` - linear regression function
  mod <- lm(CHOLESTEROLTOTAL ~ LDLESTIMATE , data = HANES)
  # add the residuals to the data and name it as resid
  HANES_mod <- HANES %>% 
    add_residuals(mod) %>% 
      mutate(resid = resid)
  # Plot the total cholesterol and residuals
  ggplot(data = HANES_mod) + 
    geom_point(mapping = aes(x = CHOLESTEROLTOTAL, y = resid))
```

##

```{r eval=TRUE, message=FALSE, warning=FALSE, echo=FALSE}  
  # Boxplot the health status and residuals that gives the 
  # effect of total cholesterol removing the effects of LDL
  # and reorder by the median
  ggplot(data = HANES_mod) + 
    geom_boxplot(mapping = aes(x = reorder(HSQ_1, resid, FUN = median), y = resid))
```

##

This exercise shows people with "Good" and "Very Good" health have low total cholesterol compared to people with poor health.

---

#### Selected materials and references

[R for Data Science - Exploratory Data Analysis Part](http://r4ds.had.co.nz/exploratory-data-analysis.html#questions)

---
