Dianne Cook and D. F. Swayne, 2008: "Interactive and Dynamic Graphics for Data Analysis: With Examples Using R and GGobi". http://www.ggobi.org/book/

Hadley Wickham and G. Grolemund, 2016: "R for Data Science". http://r4ds.had.co.nz/

Exploratory data analysis (EDA) is an iterative cycle. 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.

Your goal during EDA is to develop an understanding of your data.

two types of questions will always be useful for making discoveries within your data. You can loosely word these questions as:

• What type of variation occurs within my variables?
• What type of covariation occurs between my variables?

You will need to clean, transform, and visualization:

• tidyr for data cleaning
• dplyr for transformations
• ggplot2 for visualization

We use a list of 1738 all-time woman's best performances on 10,000m (track) from http://www.alltime-athletics.com/w_10kok.htm. The following figure presents kernel density estimates and observed values ("ticks" on x-axis). Remember that the sample contains only observations from the left tail.

File: athletics_10000W.txt

athl %>% print(n = 5)

## # A tibble: 1,738 × 9
##    best     time             name country date_of_birth position
##   <int>    <dbl>            <chr>   <chr>         <dbl>    <chr>
## 1     1 29.29083      Almaz Ayana     ETH          1991        1
## 2     2 29.52967      Wang Junxia     CHN          1973        1
## 3     3 29.54217 Vivian Cheruiyot     KEN          1983        2
## 4     4 29.70933  Tirunesh Dibaba     ETH          1985        3
## 5     5 29.89183   Alice Nawowuna     KEN          1994        4
## # ... with 1,733 more rows, and 3 more variables: place <chr>, date <dbl>,
## #   age <dbl>

athl$time %>% summary

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   29.29   31.45   31.77   31.65   31.99   32.17

• Observed times has median of 31.77 and standard deviation of 0.46 minutes.
• Current world record (August 2016) is equal to median - 5.4*standard error – seems to be clear outlier.

(Just for comparison my trail PB on 10,000m is slightly above 50 minutes.)

If we see data like this we should consider whether we see a mistake (bug) or information.

Hawkins (1980) defined an outlier as "an observation which deviates so much from other observations as to arouse suspicions that it was generated by a different mechanism"

Barnett and Lewis (1994) defined an outlier as "an observation (or subset of observations) which appears to be inconsistent with the remainder of that set of data"

• Measurement/coding error – e.g. a race in Sevilla in 1999 was actually 9,600m long (one lap short)
• Coding values – e.g. -99 is favourite number for indicating missing value
• Changes in coding –
• Errors in table joining – "This dataset uses ISO/Polity IV/… coding." Hardly.
• …

• Do your scripts do what they are supposed to do?
• Does your solution really works in all settings?
• Do you use proper IDs for observations?
• Are joined data (e.g. from different survey waves) compatible?

Remember that this kind of bugs is very often silent! Careful reading of code books is the way to avoid them and EDA is the way to discover them.

On the other hand it is fairly possible that data are correct. In this case the outliers often show us some important information about data generating process.

If it is the case we come up with a new question to investigate (Was doping involved? How about wind/altitude/…?)

In order to demonstrate that we will use famous iris data set with petal and sepal sizes of iris setosa and iris versicolor.

File: iris.Rdata

print(iris, n=5)

## # A tibble: 100 × 5
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## # ... with 95 more rows

We will focus on relationship between sepal dimensions:

iris %>% 
  ggplot(
      aes(
    x = Sepal.Length,
    y = Sepal.Width
  )
  ) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  coord_fixed() +
  ggtitle("Sepal dimensions") +
  xlab("Length") + ylab("Width") + theme_bw() -> p

OLS spoke out clearly but the result is still a bit suspicious for someone at least touched by biology. To reveal what is really going on we need to dig a bit deeper.

Apart measurements errors one can suspect that there is a third variable which drives regression results. One is right in this case:

iris %>% 
  ggplot(aes(
    x = Sepal.Length,
    y = Sepal.Width
  )) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  stat_ellipse(aes(color = Species)) +
  geom_smooth(method = "lm", se = FALSE, aes(color = Species)) +
  coord_fixed() +
  ggtitle("Sepal dimensions") + xlab("Length") + ylab("Width") +
  theme_bw() -> p

Careful carrying out of EDA will helps you to better understand data generating process and subsequently prevent modeling mistakes.

In this exercise we will investigate whether there is a gender gap – a significant difference in wages – in Czech Republic.

Dataset covers sample of observations from Czech Republic gathered by WageIndicator Foundation (http://www.wageindicator.org/).

File: wi_sample_cr.dta

Data are available in dta (Stata) format. Data from dta can be load using foreign::read.dta() (only older versions) or haven::read_dta() which supports only limited number of data classes (logical, integer, numeric, character, and factor).

library(haven)
read_dta("data/wi_sample_cr.dta") -> wages

Other statistical packages use very often labelled vectors (class labelled in R). Our dataset is not different:

lapply(wages,class) %>% unlist %>% head

##     casenum      locale     surveyy     survemm     surveqy     country 
##   "numeric" "character"  "labelled"  "labelled"  "labelled"  "labelled"

Object of class labelled contains coded observed values and codes descriptions. Labels are stored as attributes.

See example:

wages$depmale %>% attributes() %>% names

## [1] "label"        "format.stata" "class"        "labels"

Our table contains 103 variables which is a bit too much for our exercise. Therefore we will restrict out variables to 6:

wages %<>% 
    select(
        surveyy,    #year of survey
        gender,     #gender
        wagegrhr,   #gross hourly wage in nat. currency
        tenuexpe,   #years of service
        satlife,    #matrix: satisfaction with life as-a-whole
        sathhdiv    #satisfaction with division of housework
    )

As you can see from the table automatic coersion to factor might not be a good idea. A lot of variables in our selection is clearly continous:

read_dta("data/wi_sample_cr.dta") %>% 
    transmute(
        surveyy = as.integer(surveyy),
        gender = as.factor(gender),
        wagegrhr = as.numeric(wagegrhr),
        tenuexpe = as.numeric(tenuexpe),
        satlife = as.integer(satlife),
        sathhdiv = as.integer(sathhdiv)
    ) -> wages

It is worth of consideration not to use such plots as far as the data can be presented in more concise way:

wages %$% table(surveyy) 

## surveyy
##  2009  2010  2011  2012  2013  2014  2015 
## 20502  4288  1039  2146  2908  2474  1155

In this case a bar plot does not bring any added value – it presents the same information using more ink. (Do you remember Tufte?)

As far as we are interested in gender gap we should be worried about number of males/females in the sample:

wages %>% 
  ggplot(aes(gender)) +
  geom_bar() +
  theme_bw() -> p

You shouldn't use pie charts because:

• it is difficult to read values (join data representation with the scale)
• ït is difficult to compare sizes of slices (different items are actually shown using different shapes) – Try to order shares from a pie chart above.

Have you ever seen a bar plot or a pie chart made in 3D? If so then you know that it makes everything 3Times worse.

Our primary variable of interest is hourly wage which is clearly a continous variables. To get a basic picture it is useful to take a look at decriptive statistics:

wages %$% summary(wagegrhr) # gross hourly wage in nat. currency

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 0.000e+00 9.800e+01 1.380e+02 1.846e+13 1.940e+02 5.132e+17      5755

It seems that highest hourly wage in the sample is 5.132153110^{17} CZK. It is quite a suspicious figure. It seems that there might be a bug in the data.

Let's take a look at histogram:

wages %>% 
  ggplot(
    aes(
      x = wagegrhr 
    )
  ) +
  geom_histogram(bins = 50) +
  theme_bw() -> p

You can see that observed wage distribution follows a common pattern:

• it is unimodal
• it is asymetric (close to log-normal)
• mean is higer than median (as a consequence)
• there are very distant observations on the right tail

Histogram depicts observed values. There are also geoms which allow us to plot density estimates:

wages %>% 
  ggplot(
    aes(
      x = wagegrhr 
    )
  ) +
  geom_histogram(binwidth = 50) +
  geom_density(
    aes(y = 50 * ..count..),
    color = "blue"
    ) +
  theme_bw() -> p

Both distributions are clearly visible now, but it is extremely hard to compare them.

Maybe it is time to call for something completely different…

There are several ways to deal with the issue:

The first one is to change properties of scatter plot – set an alpha or make "points" as small as possible using shape=".".

Theory shows that using rectangle bins might be misleading. It is better to use hexagon bins. (geom_hex has an extra dependency.) Practically horizontal orientation of "highest" bins suggests that there is actually no correlation between hourly wages (wagegrhr) and legth of practice (tenuexpe).

Relationship between income and life satisfaction is a subject of many papers. Let's see what can we see in our data. Happiness (satlife) is (as usuall) evaluated using likert scale from 1 to 10 (bigger is better):

wages %>% 
  ggplot(
    aes(x=wagegrhr,y=satlife)
  ) + geom_point() +
  theme_bw()  -> p

Or we can accept, that life satisfaction is discrete and use a boxplot which basically technique to display relationship discrete and continous variable:

wages %>% 
  ggplot(
    aes(y=wagegrhr,x=factor(satlife))
  ) + geom_boxplot() +
  coord_flip(ylim = c(0,1000)) +
  theme_bw() -> p

Happy life starts with happy home. And happy home starts with washed dishes. This bold so called wife hypothesis can be tested using our dataset which contains variable sathhdiv – satisfaction with housework division. But there is a problem – both satisfactions are evaluated using likert scale (from 1 to 5 in the case of sathhdiv).

wages %>% 
  ggplot(
    aes(y=satlife,x=sathhdiv)
  ) + geom_point() +
  theme_bw() -> p 

A very definition of the word useless.

Use S2 database of migartion flows (Abel & Sander, 2014) from the file HW_abel.Rdata. The table abel contains 4 columns:

print(abel, n=5)

## # A tibble: 153,664 × 4
##   origin destination  flow period
##    <chr>       <chr> <int>  <int>
## 1    ABW         ABW     0   1990
## 2    AFG         ABW     0   1990
## 3    AGO         ABW     0   1990
## 4    ALB         ABW     0   1990
## 5    ARE         ABW     1   1990
## # ... with 1.537e+05 more rows

• origin and destination contain ISO-3 country codes of countries of origin and destination
• flow contains the estimated size of the migration flow from origin to destination
• Estimated flows are available for 4 five-years long periods (1990-1995, 1995-2000, 2000-2005, 2005-2010). Column period contains the first year of each period.

You are supposed to:

1. Compute total outflows from each origin in each period.
2. Compute emigration rates – data on population in the first year of each period are available in HW_population_un.Rdata. Population data are in thousands of persons!
3. Identify suspcious total outflows – total outflows which are less or equal to 10 % of the value of total outflows from previous period.
4. Report only suspicious observations in the following structure:

solution %>% print(n=5)

## Source: local data frame [46 x 3]
## Groups: origin [35]
## 
##   origin period    emig_rate
##    <chr>  <int>        <dbl>
## 1    BIH   1995 0.000000e+00
## 2    CYP   1995 2.221212e-05
## 3    DJI   1995 0.000000e+00
## 4    ERI   1995 4.028640e-03
## 5    ESH   1995 0.000000e+00
## # ... with 41 more rows

Names of columns are mandatory. Do not perform any rounding etc.

Assign resulting table into variable abel_outliers.