Exploratory Data Analysis (Completed Notebook)

R for Data Science

Set-up

Load the tidyverse:

library(tidyverse)

Acknowledgements

This notebook is based on Chapter 7 in R for Data Science (Grolemund and Wickham 2020). It does a nice treatment of the diamonds data set. So we will look at some other examples.

1 Model building

Models are about compressing information: take many data points and compress them into a summary. Throw out information to gain information!

The estimator (i.e. model) for the unconditional population mean is a form of compression:

\[ E[X] = \frac{1}{n} \sum_{i=1}^n x_i \]

Take \(n\) observations and compress them into one summarization of central tendency.

The same principle applies if we run a fancier estimator (e.g., regression, estimator for the conditional population mean). Models produce “insight” because they compress information.

1.1 Building up to building models

Before we can build models we first need to figure out what should go in them.

Say we have some outcome of interest \(y\) (e.g., number of website visitors). We think the number of unique visitors is a function of some variables:

\[ y = f(\cdot) \]

and if we new a) what those variables were and b) how they influence the outcome, we would have a good idea how to increase readership.

This is exploratory data analysis (EDA).

At a shallow level EDA boils down to calculating sumamry statistics and making plots.

But the bigger picture behind EDA is figuring out what goes inside \(f(\cdot)\) so you can build (and estimate the parameters of) a model.

1.2 Types of models

Remember there are roughly two types of models:

• models that predict (“what” will people buy?")
• models that infer (“why” will people buy it?")

Prediction is the domain of traditional machine learning. You don’t necessarily care about why – you just want to generate the most accurate prediction possible.

Inference is about hypothesis testing: the “how” or “why” some variables affect others.

In this notebook we will focus on inference. Specifically, we will use EDA to generate testable hypotheses.

2 Height

Let’s start with a trivial example: variation in height across males and females.

Re-load the nhanes survey data:

nhanes = read_csv("https://raw.githubusercontent.com/lrdegeest/r_for_data_science/main/data/nhanes.csv")

We have basic medical data of over 10,000 randomly surveyed individuals (the survey was carried out by the CDC). One of those variables is height.

How is height distributed in our sample? Let’s plot a histogram

ggplot(data = nhanes, aes(x = height)) + # fill this in!
  geom_histogram() 

Height is clearly normally distributed around it’s mean:

nhanes %>% 
  summarise(mean(height))

So if we were to randomly draw somebody from the sample using slice_sample (another of the slice family of operators):

# try playing this chunk many times!!
nhanes %>% # take the data
  slice_sample(n=1) %>% # draw one random person
  select(height) # view their height

more often than not we will draw somebody whose height is close to the mean height.

2.1 Height by sex

But what if we were to randomly draw two people: one male, one female? Play this chunk a few times:

# try playing this chunk many times!!
nhanes %>% # take the data
  group_by(sex) %>% # group observations by sex
  slice_sample(n=1) %>% # draw one random person (from each sex)
  select(height) # view their height

More often than not it seems like we are drawing a male taller than the female. This would suggest that men are systematically taller than women.

We could plot the distributions of height by sex (in one panel, using fill):

ggplot(data = nhanes, aes(x = height, fill = sex)) + 
  geom_histogram()

OK, looks like there is a clear difference in the distributions. But this plot covers up the right half of the female distribution. Can we clean it up?

Let’s use a 75% transparency (alpha = 0.75) and position the histograms “side-by-side” (position = "identity"):

ggplot(data = nhanes, aes(x = height, fill = sex)) + 
  geom_histogram(alpha = 0.75, position = "identity")

2.2 Generating a hypothesis

So yes, it looks like the distributions are different. We could calculate the averages by sex:

nhanes %>% 
  group_by(sex) %>% 
  summarise(mean(height))

and now think about a hypothesis: on average, men are taller than women – not just in this sample, but in the population from which this sample was drawn.

More specifically we can formulate a null and alternative hypothesis:

• Null: no relationship between sex and height (hence “null”)
• Alternative: some relationship between sex and height

We want to see if our data provides evidence to support the alternative hypothesis. This is inference in a nutshell.

2.3 Testing the hypothesis

We can test our hypothese with a model where height is a function of sex plus some random error (\(\epsilon\)):

\[ \text{height} = f(\text{sex}) + \epsilon \]

Now we have to choose the functional form of \(f()\). The simplest is a line:

\[ \begin{aligned} \text{height} &= f(\text{sex}) + \epsilon \\ &= \beta_0 + \beta_1(\text{sex}) + \epsilon \end{aligned} \] and use lm() (for “linear model”) to estimate the parameters \(\beta_0\) (the intercept) and more importantly \(\beta_1\) (the relationship between sex and height) and then summary() to view the results (coefficient estimates and inference):

lm(formula = height ~ sex, data = nhanes) %>% # estimate the model
  summary() # view the coefficients AND hypothesis tests

Call:
lm(formula = height ~ sex, data = nhanes)

Residuals:
     Min       1Q   Median       3Q      Max 
-29.9451  -4.4451  -0.0403   4.6559  27.3587 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 161.23933    0.09374 1719.99   <2e-16 ***
sexMale      13.50274    0.13604   99.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.912 on 10349 degrees of freedom
Multiple R-squared:  0.4877,    Adjusted R-squared:  0.4876 
F-statistic:  9851 on 1 and 10349 DF,  p-value: < 2.2e-16

Ample evidence to support the alternative hypothesis. On average, men are taller than women.

2.4 lm()

In another notebook we’ll talk more about linear models. For now, just note that lm() takes two arguments:

• formula =: the model formula. y ~ x means \(y = f(x)\).
• data =: the data that house the variables used in formula

3 Avocados

Let’s consider a more interesting example: demand for organic avocados.

Here is some data on avocados Kaggle thanks to the Hass Avocado Board:

avocados = read_csv("https://raw.githubusercontent.com/lrdegeest/r_for_data_science/main/data/avocados.csv")

3.1 Checkpoint

Plot a histogram of avocado prices (AveragePrice).

ggplot(data = avocados, aes(x = AveragePrice)) + 
  geom_histogram()

3.2 Checkpoint

Now plot a histogram of avocado prices by organice/not organic (type). Make sure to handle overlapping distributions.

ggplot(data = avocados, aes(x = AveragePrice, fill=type)) + 
  geom_histogram(alpha = 0.75, position = "identity")

3.3 Checkpoint

Use lm() and summary() to estimate a model of avocado prices as a linear function of avocado type

lm(formula = AveragePrice ~ type, data = avocados) %>% 
  summary()

Call:
lm(formula = AveragePrice ~ type, data = avocados)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.21400 -0.20400 -0.02804  0.18600  1.59600 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.158040   0.003321   348.7   <2e-16 ***
typeorganic 0.495959   0.004697   105.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3173 on 18247 degrees of freedom
Multiple R-squared:  0.3793,    Adjusted R-squared:  0.3792 
F-statistic: 1.115e+04 on 1 and 18247 DF,  p-value: < 2.2e-16

3.4 Demand shifts

OK, but does this price shift lead to a demand shift?

If so, we should the demand curve shift down for organic avocados.

So let’s look at sales of avocados and see if there is a relationship between sales and avocado type.

3.4.1 Renaming columns

Unfortunately some of the column names have spaces in them:

colnames(avocados)

 [1] "X1"           "Date"         "AveragePrice" "Total Volume" "4046"        
 [6] "4225"         "4770"         "Total Bags"   "Small Bags"   "Large Bags"  
[11] "XLarge Bags"  "type"         "year"         "region"      

including the column on sales, “Total Volume”. R will throw an error if we try to plot its distribution because it will think “Total” and “Volume” are separate objects:

ggplot(data = avocados, aes(x = Total Volume)) + 

Error: unexpected symbol in "ggplot(data = avocados, aes(x = Total Volume"

We can use rename() to rename that column:

avocados = avocados %>% # make sure you reassign the data!
  rename(TotalVolume = Total Volume)

Error: unexpected symbol in:
"avocados = avocados %>% # make sure you reassign the data!
  rename(TotalVolume = Total Volume"

Argh. Why won’t this work? Now you see why spaces are so annoying!

(You might say it’s the “pits”…https://tinyurl.com/yyawvnwn)

To “escape” the space we have to wrap the old name in quotes:

avocados = avocados %>% # make sure you reassign the data!
  rename(TotalVolume = "Total Volume")

3.5 Demand shifts (continued)

OK, now we’re in business:

ggplot(data = avocados, aes(x = TotalVolume)) + 
  geom_histogram()

Sales are highly skewed. That means the average is not representative. So if we build a model about the conditional average (e.g., average sales conditional on avocado type), it won’t be useful.

Fortunately we know that we can transform data onto a new scale without changing it’s underlying information (we saw this when we did z-scores).

A common method to transform a non-normal variable into a normal variable is the log transformation. (Just make sure you don’t have negative or zero values! To see why, try typing log(0) and log(-1) in your console.)

3.5.1 Checkpoint

Use mutate() and log() to create a new variable called “log_sales” that log transforms TotalVolume.

avocados = avocados %>% 
  mutate(log_sales = log(TotalVolume))

3.5.2 Checkpoint

Plot the distribution of log_sales:

ggplot(data = avocados, aes(x = log_sales)) + 
  geom_histogram()

(More) normally distributed!

3.5.3 Checkpoint

Make a scatter plot of log sales as a function of price. Use a transparency (try different values of alpha!) Hint: alpha can go inside geom_point(), just like we did with geom_histogram().

ggplot(data = avocados, aes(x = AveragePrice, y = log_sales)) + 
  geom_point(alpha = 0.20)

3.5.4 Checkpoint

Now re-do the plot and this time add geom_smooth(method = "lm"). This will add a regression line to the plot.

ggplot(data = avocados, aes(x = AveragePrice, y = log_sales)) + 
  geom_point(alpha = 0.20) + 
  geom_smooth(method = "lm")

3.5.5 Checkpoint

Estimate the parameters of that line with lm() and summary():

lm(formula = log_sales ~ AveragePrice, data = avocados) %>% 
  summary()

Call:
lm(formula = log_sales ~ AveragePrice, data = avocados)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2673 -1.2907  0.0552  1.1642  6.5046 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  15.98827    0.05008  319.23   <2e-16 ***
AveragePrice -3.32293    0.03425  -97.03   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.863 on 18247 degrees of freedom
Multiple R-squared:  0.3404,    Adjusted R-squared:  0.3403 
F-statistic:  9415 on 1 and 18247 DF,  p-value: < 2.2e-16

3.5.6 Checkpoint

Now, re-plot the scatterplot of log sales and avocado type with a regression line, but this time use facet_wrap() to make separate panels for conventional and organic avocados:

ggplot(data = avocados, aes(x = AveragePrice, y = log_sales)) + 
  geom_point(alpha = 0.20) +
  geom_smooth(method = "lm") + 
  facet_wrap(~type) 

3.6 Demand shifts: bringing it all together

Notice two key things in that last plot:

1. the intercept of the regression line for organic avocados shifted down relative to the intercept for conventional avocados
2. the slope of the regression lines look different; the relationship between price and sales varies across avocado types

Now we have a good idea of a model we can build that summarizes the demand for avocados. The model should:

• allow the intercept to vary by avocado type, and;
• allow the slope to vary by avocado type

Here is such a model that includes a dummy variable (to capture the changing intercept) and an interaction effect (to capture the changing slope):

\[ \text{log(sales)} = \beta_0 + \beta_1(\text{price}) + \beta_2(\text{organic}) + \beta_3(\text{price} \times \text{organic}) + \epsilon \]

Let’s estimate the parameters with lm():

lm(formula = log_sales ~ AveragePrice + type + AveragePrice*type, data = avocados) %>% 
  summary()

Call:
lm(formula = log_sales ~ AveragePrice + type + AveragePrice * 
    type, data = avocados)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0958 -0.9735 -0.1847  0.7429  4.8303 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              14.49652    0.06420  225.81  < 2e-16 ***
AveragePrice             -1.16907    0.05406  -21.62  < 2e-16 ***
typeorganic              -3.87157    0.09226  -41.96  < 2e-16 ***
AveragePrice:typeorganic  0.48246    0.06673    7.23 5.03e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.358 on 18245 degrees of freedom
Multiple R-squared:  0.6493,    Adjusted R-squared:  0.6492 
F-statistic: 1.126e+04 on 3 and 18245 DF,  p-value: < 2.2e-16

What does this model say?

• The intercept for conventional avocados is 14.49652
• The shift in the intercept from conventional to organic is -1.16907
• The effect of price on sales (the slope) for coventional is -1.16907
• The shift in the slope from conventional to organic is 0.48246

So it seems that people who buy organic avocados are less sensitive to changes in price!

Of course, take this analysis with a grain of salt.

There are probably several confounders to address, like variation in prices over time and across regions.

Also, who is buying these avocados? This might explain why we see a right shift in volume for conventional avos. Is it folks like us at the store? Restaurants buying in bulk? We don’t know.

4 Blood pressure

Let’s go back to nhanes and consider a more interesting example: variation in blood pressure over time and by sex.

When your heart beats, it creates pressure to circulated oxygenated blood throughout your body. That pressure has two parts, systolic and diastolic.

That is why when they take your blood pressure they say “120 over 80” (units of millimeters of mercury, or mmHg), meaning 120 systolic over 80 diastolic. 120/80 is considered a normal blood pressure. Anything above 130 systolic is considered high blood pressure.

The harder your blood travels the more it deteriorates your cardiovascular system. This excess force is hypertension or high blood pressure. Over time hypertension can lead to a number of health issues.

In nhanes the variable bpsystol captures systolic blood pressure.

4.1 Checkpoint

Calculate avereage systolic blood pressure by sex:

nhanes %>% 
  group_by(sex) %>% 
  summarise(mean(bpsystol))

4.2 Checkpoint

Plot the distribution of systolic blood pressure by sex:

ggplot(data = nhanes, aes(x = bpsystol, fill = sex)) + 
  geom_histogram(alpha = 0.75, position = "identity")

4.3 Blood pressure over time

It doesn’t look like there is much a of difference in the unconditional average blood pressure between men and women.

But what if we need to look at the conditional average? For instance, what if blood pressure varies over time (age)?

4.3.1 Checkpoint

First plot the bpystol as a function of age. Include a regression line.

ggplot(data = nhanes, aes(x = age, y = bpsystol)) + 
  geom_point() + 
  geom_smooth(method = "lm")

4.3.2 Checkpoint

That plot is a bit messy. Let’s simplify it. Now plot average blood pressure by age. Include a regression line. To do this, first use summarise and create a variable called “mean_bpsystol” by age, then pipe (%>%) the output to ggplot. (You might find it helpful to first build the summary table). Don’t include a regression line.

nhanes %>% 
  group_by(age) %>% 
  summarise(mean_bpsystol = mean(bpsystol)) %>% 
  ggplot(., aes(age, mean_bpsystol)) + 
  geom_point() + 
  geom_line() 

4.3.3 Checkpoint

Does this relationship vary by sex?

Re-do that plot, but this time by age and sex. Instead of fill= use color= to color the plot by sex:

nhanes %>% 
  group_by(age, sex) %>% 
  summarise(mean_bpsystol = mean(bpsystol)) %>% 
  ggplot(., aes(age, mean_bpsystol, color = sex)) + 
  geom_point() + 
  geom_line() 

4.4 Modeling the relationship

When we plot average systolic blood pressure by age for men and women, we see an upward trend in both sexes, but not necessarily the same trend. The lines do not look parallel.

Another way to think about is that the slope of the lines for age are different for men and women. In the graph it looks like the effect of age is more pronounced – the slope is steeper – among women.

Based on our exploratory analysis we might consider the following model:

\[ \text{blood pressure} = \beta_0 + \beta_1 \text{age} + \beta_2 \text{sex} + \beta_3 (\text{sex} \times \text{age}) + \epsilon \]

4.4.1 Checkpoint

Modify the model we estimated for avocado sales to estimate the model above. Recall that inside lm() we can interact two variables x1 and x2 by typeing x1*x2.

lm(formula = bpsystol ~ sex + age + sex*age, data = nhanes) %>% 
  summary()

Call:
lm(formula = bpsystol ~ sex + age + sex * age, data = nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-65.990 -13.214  -2.101  10.606 156.826 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 90.11098    0.80336  112.17   <2e-16 ***
sexMale     20.45813    1.16526   17.56   <2e-16 ***
age          0.81635    0.01583   51.56   <2e-16 ***
sexMale:age -0.34573    0.02304  -15.01   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.14 on 10347 degrees of freedom
Multiple R-squared:  0.2551,    Adjusted R-squared:  0.2548 
F-statistic:  1181 on 3 and 10347 DF,  p-value: < 2.2e-16

Interestingly, the effect of age on blood pressure among men, sexMale:age, is negative.

Meaning that as men age, their blood pressure goes up (because the average effect of age is positive), but slower than for women (because the interaction effect is negative).

The general takeaway is that a patient should be compared against the right benchmark. For instance, we can’t just compare a male patient’s blood pressure against the overall average. They need to be compared to the average male at the same age.

5 Big picture

Don’t worry if you haven’t followed every step, especially the stuff with lm(). We’ll do more of that later.

For now, just notice how we were able to use summary statistics and plots to generate hypotheses, build models, and then estimate those models to compress the data and create insights.