Exploratory Data Analysis

Introduction

This chapter will show you how to use visualisation and transformation to explore your data in a systematic way, a task that statisticians call exploratory data analysis, or EDA for short. EDA is an iterative cycle. You:

1. Generate questions about your data.
2. Search for answers by visualising, transforming, and modelling your data.
3. Use what you learn to refine your questions and/or generate new questions.

EDA is not a formal process with a strict set of rules. More than anything, EDA is a state of mind. During the initial phases of EDA you should feel free to investigate every idea that occurs to you. Some of these ideas will pan out, and some will be dead ends. As your exploration continues, you will home in on a few particularly productive areas that you'll eventually write up and communicate to others.

EDA is an important part of any data analysis, even if the questions are handed to you on a platter, because you always need to investigate the quality of your data. Data cleaning is just one application of EDA: you ask questions about whether your data meets your expectations or not. To do data cleaning, you'll need to deploy all the tools of EDA: visualisation, transformation, and modelling.

Prerequisites

In this chapter we'll combine what you've learned about dplyr and ggplot2 to interactively ask questions, answer them with data, and then ask new questions.

You will use the datasetindia04that you created as part of Assignment 2, together with a template for Assignment 3,Assignment3_Start.ipynb. Both files will be provided to you on theAssignment3.zipzipped file uploaded on Canvas.

1. Download and extract the contents of theAssignment3.zipfile in your preferred working directory (folder). This is an archived folder, so you will need to extract it.

Once succesfully extracted, you should be able to see in your working directory both theAssignment3_Start.ipynbnotebook and theindia04.Rdatadata file. Make sure you keep theindia04data file in the same folder as the Jupyter Notebook.

1. Using Jupyter Notebook, open the providedAssignment3_Start.ipynbtemplate notebook, and rename itFirstName_LastName_A3. Next, type and run the following lines of code at the beginning of your assignment:

 library(tidyverse) library(nycflights13) data(diamonds) load("india04.Rdata") 

to load the dataset into your R session. Make sure to include these lines at the top of the notebook you are going to submit for the assignment, preferrably where you are prompted to do so in the template, under theLoad packages and data needed for this assignmentcell.

Submission

As before, you should submit a zip-folder with your assignment. The zip folder should be namedFirstName_LastName_A3.zip. It should contain

• The Jupyter Notebook with your code and formatted markdown answers to the assignment questions, "FirstName_LastName_A3.ipynb"
• Theindia04.Rdatafile.
• Any html files that get created using stargazer in the exercises that follow.

You are responsible for naming your files appropriately, and that the code should run without errors.Please also make sure to include your name within your Jupyter Notebook at the top of your script.

Questions

"There are no routine statistical questions, only questionable statistical routines." Sir David Cox

"Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise." John Tukey

Your goal during EDA is to develop an understanding of your data. The easiest way to do this is to use questions as tools to guide your investigation. When you ask a question, the question focuses your attention on a specific part of your dataset and helps you decide which graphs, models, or transformations to make.

EDA is fundamentally a creative process. And like most creative processes, the key to askingqualityquestions is to generate a largequantityof questions. It is difficult to ask revealing questions at the start of your analysis because you do not know what insights are contained in your dataset. On the other hand, each new question that you ask will expose you to a new aspect of your data and increase your chance of making a discovery. You can quickly drill down into the most interesting parts of your dataand develop a set of thought-provoking questionsif you follow up each question with a new question based on what you find.

There is no rule about which questions you should ask to guide your research. However, two types of questions will always be useful for making discoveries within your data. You can loosely word these questions as:

1. What type of variation occurs within my variables?
2. What type of covariation occurs between my variables?

The rest of this chapter will look at these two questions. I'll explain what variation and covariation are, and I'll show you several ways to answer each question. To make the discussion easier, let's define some terms:

• Avariableis a quantity, quality, or property that you can measure.
• Avalueis the state of a variable when you measure it. The value of a variable may change from measurement to measurement.
• Anobservationis a set of measurements made under similar conditions (you usually make all of the measurements in an observation at the same time and on the same object). An observation will contain several values, each associated with a different variable. I'll sometimes refer to an observation as a data point.
• Tabular datais a set of values, each associated with a variable and an observation. Tabular data istidyif each value is placed in its own "cell", each variable in its own column, and each observation in its own row.

So far, all of the data that you've seen has been tidy. In real-life, most data isn't tidy, so we'll come back to these ideas again in [tidy data].

Variation

Variationis 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. This is true even if you measure quantities that are constant, like the speed of light. Each of your measurements will include a small amount of error that varies from measurement to measurement. Categorical variables can also vary if you measure across different subjects (e.g.the eye colors of different people), or different times (e.g.the energy levels of an electron at different moments). Every variable has its own pattern of variation, which can reveal interesting information. The best way to understand that pattern is to visualise the distribution of the variable's values.

Visualising distributions

How you visualise the distribution of a variable will depend on whether the variable is categorical or continuous. A variable iscategoricalif it can only take one of a small set of values. In R, categorical variables are usually saved as factors or character vectors. To examine the distribution of a categorical variable, use a bar chart:

ggplot(data = india04) + geom_bar(mapping = aes(x = edattain)) 

The height of the bars displays how many observations occurred with each x value. You can compute these values manually withdplyr::count():

india04 %>% count(edattain) ## # A tibble: 5 x 2 ## edattain n ##   ## 1  ## 2 primary 52381 ## 3 secondary 35147 ## 4 university 11723 ## 5 Unknown 85 

A variable iscontinuousif it can take any of an infinite set of ordered values. Numbers and date-times are two examples of continuous variables. To examine the distribution of a continuous variable, use a histogram:

ggplot(data = india04) + geom_histogram(mapping = aes(x = incwage), binwidth = 500) 

You can compute this by hand by combiningdplyr::count()andggplot2::cut_width():

india04 %>% count(cut_width(incwage, 500)) ## # A tibble: 35 x 2 ## `cut_width(incwage, 500)` n ##   ## 1 [-250,250] 10596 ## 2 (250,750] 14681 ## 3 (750,1.25e+03] 3206 ## 4 (1.25e+03,1.75e+03] 2556 ## 5 (1.75e+03,2.25e+03] 1905 ## 6 (2.25e+03,2.75e+03] 978 ## 7 (2.75e+03,3.25e+03] 762 ## 8 (3.25e+03,3.75e+03] 443 ## 9 (3.75e+03,4.25e+03] 259 ## 10 (4.25e+03,4.75e+03] 127 ## # ... with 25 more rows 

A histogram divides the x-axis into equally spaced bins and then uses the height of a bar to display the number of observations that fall in each bin. In the graph above, the tallest bar shows that 35,100 individuals have anincwagevalue between 250 and 750, which are the left and right edges of the bar.

You can set the width of the intervals in a histogram with thebinwidthargument, which is measured in the units of thexvariable. You should always explore a variety of binwidths when working with histograms, as different binwidths can reveal different patterns. For example, here is how the graph above looks when we zoom into just the individuals with a wage of less than 1000 rupees and choose a smaller binwidth.

smaller <- india04 %>% filter(incwage < 1000) ggplot(data = smaller, mapping = aes(x = incwage)) + geom_histogram(binwidth = 100) 

If you wish to overlay multiple histograms in the same plot, I recommend usinggeom_freqpoly()instead ofgeom_histogram().geom_freqpoly()performs the same calculation asgeom_histogram(), but instead of displaying the counts with bars, uses lines instead. It's much easier to understand overlapping lines than bars.

ggplot(data = smaller, mapping = aes(x = incwage, colour = edattain)) + geom_freqpoly(binwidth = 100) 

There are a few challenges with this type of plot, which we will come back to invisualising a categorical and a continuous variable.

Now that you can visualise variation, what should you look for in your plots? And what type of follow-up questions should you ask? I've put together a list below of the most useful types of information that you will find in your graphs, along with some follow-up questions for each type of information. The key to asking good follow-up questions will be to rely on your curiosity (What do you want to learn more about?) as well as your skepticism (How could this be misleading?).

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 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?

Clusters of similar values suggest that subgroups might exist in your data. To understand the subgroups, ask:

• How are the observations within each cluster similar to each other?
• How are the observations in separate clusters different from each other?
• How can you explain or describe the clusters?
• Why might the appearance of clusters be misleading?

The histogram below shows the length (in minutes) of 272 eruptions of the Old Faithful Geyser in Yellowstone National Park. Eruption times appear to be clustered into two groups: there are short eruptions (of around 2 minutes) and long eruptions (4-5 minutes), but little in between.

ggplot(data = faithful, mapping = aes(x = eruptions)) + geom_histogram(binwidth = 0.25) 

Many of the questions above will prompt you to explore a relationshipbetweenvariables, for example, to see if the values of one variable can explain the behavior of another variable. We'll get to that shortly.

Unusual values

Outliers are observations that are unusual; data points that don't seem to fit the pattern. Sometimes outliers are data entry errors; other times outliers suggest important new science. When you have a lot of data, outliers are sometimes difficult to see in a histogram. For example, take the distribution of theincwagevariable we saw before. The only evidence of outliers is the unusually wide limits on the x-axis.

ggplot(india04) + geom_histogram(mapping = aes(x = incwage), binwidth = 500) 

There are so many observations in the common bins that the rare bins are so short that you can't see them (although maybe if you stare intently at 0 you'll spot something). To make it easy to see the unusual values, we need to zoom to small values of the axes withx_lim():

ggplot(india04) + geom_histogram(mapping = aes(x = incwage), binwidth = 1000) + xlim(10000, 30000) 

This allows us to see that there are several unusual values, individuals that earn more than 10000 rupees. It is natural to wonder if these observations correspond to true outliers or if these values are the result of data entry errors.

It's good practice to repeat your analysis with and without the outliers. If they have minimal effect on the results, and you can't figure out why they're there, it's reasonable to replace them with missing values, and move on. However, if they have a substantial effect on your results, you shouldn't drop them without justification. You'll need to figure out what caused them (e.g.a data entry error) and disclose that you removed them in your write-up.

Section 1 Exercises

1. (3 pts)Filter out theUnknowncategory from theedattainvariable in theindia04dataframe. Save the resulting dataframe asindia04_new.

Use theindia04_newdataframe to answer the remainder of the questions in this section.

1. (7 pts)Show one example of how you could explore the distribution of each of theage,edattain, andincwagevariables inindia04_new. What do you learn?
2. (5 pts)Compare and contrastcoord_cartesian()vsxlim()when zooming in on a histogram.
3. (5 pts)What doesna.rm = TRUEdo inmax()andmean()functions? Test each command on the variableincwageinindia04_new.
4. (5 pts)What is the highest income level observed in the data?
5. (5 pts)How many individuals earn exactly 200 rupees? How many earn 201? What are your thoughts on this difference?

Missing values

If you've encountered 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:

india04_1 <- india04 %>% filter(incwage < 10000) 

1. I don't recommend this option because just because one measurement is invalid, doesn't mean all the measurements are. Additionally, if you have low quality data, by time that you've applied this approach to every variable you might find that you don't have any data left!
2. Instead, I recommend replacing the unusual values with missing values. The easiest way to do this is to usemutate()to replace the variable with a modified copy. You can use theifelse()function to replace unusual values withNA:

india04_2 <- india04 %>% mutate(incwage_adj = ifelse(incwage > 10000, NA, incwage)) 

1. ifelse()has three arguments. The first argumenttestshould be a logical vector. The result will contain the value of the second argument,yes, whentestisTRUE, and the value of the third argument,no, when it is false. Alternatively to ifelse, usedplyr::case_when().case_when()is particularly useful inside mutate when you want to create a new variable that relies on a complex combination of existing variables.

Like R, ggplot2 subscribes to the philosophy that missing values should never silently go missing. It's not obvious where you should plot missing values, so ggplot2 doesn't include them in the plot, but returns a warning that they've been removed.

To suppress that warning, you can setna.rm = TRUE:

ggplot(india04) + geom_histogram(mapping = aes(x = incwage), binwidth = 1000, na.rm = TRUE) 

Other times you want to understand what makes observations with missing values different to observations with recorded values. For example, innycflights13::flights, missing values in thedep_timevariable indicate that the flight was cancelled. So you might want to compare the scheduled departure times for cancelled and non-cancelled times. You can do this by making a new variable withis.na().

cancelled <-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 ) cancelled %>% ggplot(mapping = aes(sched_dep_time)) + geom_freqpoly(mapping = aes(colour = cancelled), binwidth = 1/4) 

However this plot isn't great because there are many more non-cancelled flights than cancelled flights. In the next section we'll explore some techniques for improving this comparison.

Section 2 Exercises

Start from theindia04_newdataframe you created in Section 1 Exercise 1:

1. (5 pts)Use theifelse()function to create a new variable calledincwage_adjin a new dataframe calledindia04_new2, where all the income values that arenotbetween the 1st and the 99th percentile are substituted with missing values (NA).
2. (5 pts)Usingindia04_new2, create a new variable calledlog_incwage_adj, that is the natural logarithm ofincwage_adj. Save the resulting dataframe as a new objectindia04_new3.
3. (5 pts)Plot a histogram oflog_incwage_adjfromindia04_new3constructed above. How does it look different from the same histogram constructed using the natural log transformation of theincwagevariable from theindia04_newdataframe?

Covariation

If variation describes the behaviorwithina variable, covariation describes the behaviorbetweenvariables.Covariationis 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. How you do that should again depend on the type of variables involved.

A categorical and continuous variable

It's common to want to explore the distribution of a continuous variable broken down by a categorical variable, as in the previous frequency polygon. The default appearance ofgeom_freqpoly()is not that useful for that sort of comparison because the height is given by the count. That means if one of the groups is much smaller than the others, it's hard to see the differences in shape. For example, let's explore how the income of an individual varies depending on their educational level:

ggplot(data = india04_2, mapping = aes(x = incwage_adj)) + geom_freqpoly(mapping = aes(colour = edattain), binwidth = 500) 

It's hard to see the difference in distribution because the overall counts differ so much:

ggplot(india04_2) + geom_bar(mapping = aes(x = edattain)) 

To make the comparison easier we need to swap what is displayed on the y-axis. Instead of displaying count, we'll displaydensity, which is the count standardised so that the area under each frequency polygon is one.

ggplot(data = india04_2, mapping = aes(x = incwage_adj, y = ..density..)) + geom_freqpoly(mapping = aes(colour = edattain), binwidth = 500) 

There's something rather surprising about this plot - it appears that individuals that could not complete primary school and those that did have a similar distribution of income. This might mean that primary education does not have an effect on earnings! But maybe that's because frequency polygons are a little hard to interpret - there's a lot going on in this plot.

Another alternative to display the distribution of a continuous variable broken down by a categorical variable is the boxplot. Aboxplotis 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.

Let's take a look at the distribution of price by cut usinggeom_boxplot():

ggplot(data = india04_2, mapping = aes(x = edattain, y = incwage_adj)) + geom_boxplot() 

We see much less information about the distribution, but the boxplots are much more compact so we can more easily compare them (and fit more on one plot).

Section 3 Exercises

Use theindia04_new3dataframe that you created in Section 2 Exercise 2 to answer this question.

1. (10 pts)Using two different plot types among the ones seen above, explore the relationship between educational attainment, log income and sex.
2. (5 pts)Summarize your findings on the differences in the distribution of log income by sex and educational level.
3. (5 pts)It can be sometimes be difficult to interpret categorical variables likesexandurban. Dummy variables are often used to keep track of categorical variable using numbers. A dummy variable with a value of 1 is typically interpreted as being "TRUE", while a value of 0 is interpreted as "FALSE". Add two new variables to theindia04_new3dataframe and call the resulting dataframeindia04_new4:

• A new dummy variable calledfemale_dummythat has value of 1 if an individual's sex is categorized as "Female".
• A new dummy variable calledurban_dummythat has value of 1 if an individual's urban status is categorized as "Urban".

1. (5 pts)Create a new object calledindia04_summary, by selecting the log adjusted income variable, the age variable and the two dummy variables created in the previous question from theindia04_new4dataframe. Notice that these are all numerical variables.
2. (15 pts)Install the packagestargazerand import this package into your library (remember that install statements should not be in your final script). Use this package to save a summary statistics table ofindia04_summary. Some tips:

• Helpful instructions on how to use this package can be foundhere.
• Stargazer normally requires dataframes as an input; you can use the commandas.data.frame()to make sure your input data has been read correctly.
• Use the title option to give your table an appropriate name.
• Use thetype = htmloption to save the output file in your project folder. Open the html file (using a web browser like Chrome or Safari) to make sure this worked. Submit the html file with the rest of the files in your submission folder.

Two categorical variables

To visualise the covariation between categorical variables, you'll need to count the number of observations for each combination. One way to do that is to rely on the built-ingeom_count():

ggplot(data = india04_2) + geom_count(mapping = aes(x = sex, y = edattain)) 

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.

Another approach is to compute the count with dplyr:

india04_count <- india04_2 %>% count(sex, edattain) 

Then visualise withgeom_tile()and the fill aesthetic:

india04_count %>% ggplot(mapping = aes(x = sex, y = edattain)) + geom_tile(mapping = aes(fill = n)) 

If the categorical variables are unordered, you might want to use the seriation package to simultaneously reorder the rows and columns in order to more clearly reveal interesting patterns. For larger plots, you might want to try the d3heatmap or heatmaply packages, which create interactive plots.

Exercises Section 4 Exercises

Use theindia04_new3dataframe that you created in Section 2 Exercise 2 to answer this question.

1. (5 pts)Usegroup_byto display the proportion of each educational attainment. Use stargazer to report your results in a table which your code should produce and save as html and include in your submission folder. (Hint: You can use the optionssummary = FALSEandrownames = FALSEtogether with the options from Section 3 Exercise 5).
2. (10 pts)Usegroup_byto display the proportion of each education attainment within sex. Use stargazer to report your results in a table which you will also save as html and include in your submission folder.

Two continuous variables

You've already seen one great way to visualise the covariation between two continuous variables: draw a scatterplot withgeom_point(). You can see covariation as a pattern in the points. For example, if we look at working age individuals (aged 15 to 60) that work as "Professionals", we can see the positive relationship between the income and the age of an individual.

india04_3 <- filter(india04_2, between(age, 15, 59), occisco == "Professionals") ggplot(data = india04_3) + geom_point(mapping = aes(y = incwage_adj, x = age)) + geom_smooth(mapping = aes(y = incwage_adj, x = age)) 

As illustrated above, scatterplots become less useful as the size of your dataset grows, because points begin to overplot, and pile up into areas of uniform black (as above). You've already seen one way to fix the problem: using thealphaaesthetic to add transparency.

ggplot(data = india04_3) + geom_point(mapping = aes(y = incwage_adj, x = age), alpha = 1 / 10) 

But using transparency can be challenging for very large datasets. Another solution is to use bin. Previously you usedgeom_histogram()andgeom_freqpoly()to bin in one dimension. Now you'll learn how to usegeom_bin2d()andgeom_hex()to bin in two dimensions.

geom_bin2d()andgeom_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. You will need to install the hexbin package to usegeom_hex().

Some examples using thediamondsdata:

smaller <- filter(diamonds, carat < 3) ggplot(data = smaller) + geom_bin2d(mapping = aes(x = carat, y = price)) 

# install.packages("hexbin") library(hexbin) ggplot(data = smaller) + geom_hex(mapping = aes(x = carat, y = price)) 

Another option is to bin one continuous variable so it acts like a categorical variable. Then you can use one of the techniques for visualising the combination of a categorical and a continuous variable that you learned about. For example, you could bincaratand then for each group, display a boxplot:

ggplot(data = smaller, mapping = aes(x = carat, y = price)) + geom_boxplot(mapping = aes(group = cut_width(carat, 0.1))) 

Exercises Section 5 Exercises

Use theindia04_new3dataframe that you created in Section 2 Exercise 2 to answer this question.

1. (5 pts)Combine two of the techniques you've learned to visualise the combined distribution of log income, age, and sex.

ggplot2 calls

As we move on from these introductory chapters, we'll transition to a more concise expression of ggplot2 code. So far we've been very explicit, which is helpful when you are learning:

ggplot(data = faithful, mapping = aes(x = eruptions)) + geom_freqpoly(binwidth = 0.25) 

Typically, the first one or two arguments to a function are so important that you should know them by heart. The first two arguments toggplot()aredataandmapping, and the first two arguments toaes()arexandy. In the remainder of the book, we won't supply those names. That saves typing, and, by reducing the amount of boilerplate, makes it easier to see what's different between plots. That's a really important programming concern.

Rewriting the previous plot more concisely yields:

ggplot(faithful, aes(eruptions)) + geom_freqpoly(binwidth = 0.25) 

Sometimes we'll turn the end of a pipeline of data transformation into a plot. Watch for the transition from%>%to+. I wish this transition wasn't necessary but unfortunately ggplot2 was created before the pipe was discovered.

diamonds %>% count(cut, clarity) %>% ggplot(aes(clarity, cut, fill = n)) + geom_tile()