Gaussian kernel

\[K(z) = \frac{1}{\sqrt{2 \pi}}e^{-\frac{1}{2} z^2}\]

Not to be confused with the (Gaussian) radial basis function kernel.

x <- rnorm(1000)

{
  qplot(x, geom = "blank") +
    stat_function(fun = dnorm) +
    labs(title = "Gaussian (normal) kernel",
         x = NULL,
         y = NULL)
} +
{
  ggplot(infant, aes(mortal)) +
    geom_density(kernel = "gaussian") +
    labs(title = "Infant mortality rate for 195 nations",
         x = "Infant mortality rate (per 1,000)",
         y = "Density")
}

Now we have a much smoother density function.

Rectangular (uniform) kernel

\[K(z) = \frac{1}{2} \mathbf{1}_{\{ |z| \leq 1 \} }\]

where \(\mathbf{1}_{\{ |z| \leq 1 \} }\) is an indicator function that takes on the value of 1 if the condition is true (\(|z| \leq 1\)) or 0 if the condition is false. This is the naive density estimator identified previously.

x <- runif(1000, -1.5, 1.5)
x_lines <- tribble(
  ~x, ~y, ~xend, ~yend,
  -1, 0, -1, .5,
  1, 0, 1, .5
)

{
  qplot(x, geom = "blank") +
    stat_function(fun = dunif, args = list(min = -1), geom = "step") +
    # geom_segment(data = x_lines, aes(x = x, y = y, xend = xend, yend = yend)) +
    labs(title = "Rectangular kernel",
         x = NULL,
         y = NULL)
} +
{
  ggplot(infant, aes(mortal)) +
    geom_density(kernel = "rectangular") +
    labs(title = "Infant mortality rate for 195 nations",
         x = "Infant mortality rate (per 1,000)",
         y = "Density")
}

Triangular kernel

\[K(z) = (1 - |z|) \mathbf{1}_{\{ |z| \leq 1 \} }\]

triangular <- function(x) {
  (1 - abs(x)) * ifelse(abs(x) <= 1, 1, 0)
}

{
  qplot(x, geom = "blank") +
    stat_function(fun = triangular) +
    labs(title = "Triangular kernel",
         x = NULL,
         y = NULL)
} +
{
  ggplot(infant, aes(mortal)) +
    geom_density(kernel = "triangular") +
    labs(title = "Infant mortality rate for 195 nations",
         x = "Infant mortality rate (per 1,000)",
         y = "Density")
}

Quartic (biweight) kernel

\[K(z) = \frac{15}{16} (1 - z^2)^2 \mathbf{1}_{\{ |z| \leq 1 \} }\]

biweight <- function(x) {
  (15 / 16) * (1 - x^2)^2 * ifelse(abs(x) <= 1, 1, 0)
}

{
  qplot(x, geom = "blank") +
    stat_function(fun = biweight) +
    labs(title = "Biweight kernel",
         x = NULL,
         y = NULL)
} +
{
  ggplot(infant, aes(mortal)) +
    geom_density(kernel = "biweight") +
    labs(title = "Infant mortality rate for 195 nations",
         x = "Infant mortality rate (per 1,000)",
         y = "Density")
}

Epanechnikov kernel

\[K(z) = \frac{3}{4} (1 - z^2) \mathbf{1}_{\{ |z| \leq 1 \} }\]

epanechnikov <- function(x) {
  (15 / 16) * (1 - x^2)^2 * ifelse(abs(x) <= 1, 1, 0)
}

{
  qplot(x, geom = "blank") +
    stat_function(fun = epanechnikov) +
    labs(title = "Epanechnikov kernel",
         x = NULL,
         y = NULL)
} +
{
  ggplot(infant, aes(mortal)) +
    geom_density(kernel = "epanechnikov") +
    labs(title = "Infant mortality rate for 195 nations",
         x = "Infant mortality rate (per 1,000)",
         y = "Density")
}

Comparison of kernels

qplot(x, geom = "blank") +
  stat_function(aes(color = "Gaussian"), fun = dnorm) +
  stat_function(aes(color = "Epanechnikov"), fun = epanechnikov) +
  stat_function(aes(color = "Rectangular"), fun = dunif,
                args = list(min = -1), geom = "step") +
  stat_function(aes(color = "Triangular"), fun = triangular) +
  stat_function(aes(color = "Biweight"), fun = biweight) +
  scale_color_brewer(type = "qual") +
  labs(x = NULL,
       y = NULL,
       color = NULL) +
  theme(legend.position = c(0.04, 1),
        legend.justification = c(0, 1),
        legend.background = element_rect(fill = "white"))

ggplot(infant, aes(mortal)) +
  geom_density(aes(color = "Gaussian"), kernel = "gaussian") +
  geom_density(aes(color = "Epanechnikov"), kernel = "epanechnikov") +
  geom_density(aes(color = "Rectangular"), kernel = "rectangular") +
  geom_density(aes(color = "Triangular"), kernel = "triangular") +
  geom_density(aes(color = "Biweight"), kernel = "biweight") +
  scale_color_brewer(type = "qual") +
  labs(title = "Density estimators of infant mortality rate for 195 nations",
       x = "Infant mortality rate (per 1,000)",
       y = "Density",
       color = "Kernel") +
  theme(legend.position = c(0.96, 1),
        legend.justification = c(1, 1),
        legend.background = element_rect(fill = "white"))

Selecting the bandwidth \(h\)

Even within the same kernel, different values for the bandwidth \(h\) will produce different density estimates because the moving window used to include observations in the local estimate will change.

ggplot(infant, aes(mortal)) +
  geom_density(kernel = "gaussian", adjust = 5) +
  geom_density(kernel = "gaussian", adjust = 1, linetype = 2) +
  geom_density(kernel = "gaussian", adjust = 1/5, linetype = 3) +
  labs(title = "Gaussian density estimators of infant mortality rate for 195 nations",
       subtitle = "Three different bandwidth parameters",
       x = "Infant mortality rate (per 1,000)",
       y = "Density")

If the underlying density of the sample is normal with standard deviation \(\sigma\), then for the Gaussian kernel estimation the most efficient bandwidth \(h\) will be:

\[h = 0.9 \sigma n^{-1 / 5}\]

As the sample size increases, the optimal window narrows and permits finer detail than a smaller sample. Of course we don’t actually know the population standard deviation \(\sigma\). Instead we know the sample standard deviation \(s\). If the underlying density is normal, we could just substitute \(s\) as an unbiased estimate for \(\sigma\). Of course the problem is that we assumed the underlying density is normal. If this is not true, then it’s possible that the sample standard deviation is inflated. In that case, we can adjust by using an “adaptive” estimator of spread:

\[A = \min \left( S, \frac{IQR}{1.349} \right)\]

where \(IQR\) is the interquartile range of the sample and 1.349 is the interquartile range of the standard normal distribution \(N(0,1)\).

This is the default method for calculating the bandwidth with a Gaussian kernel using the density() function in R. Other kernels use different functions for determining the optimal value for \(h\).

Violin plot

ggplot(infant, aes("", mortal)) +
  geom_violin() +
  labs(title = "Violin plot of infant mortality rate for 195 nations",
       x = NULL,
       y = "Infant mortality rate (per 1,000)")

A violin plot is a compact display of a continuous distribution. In essence, it is a mirrored density plot displayed in the same way as a boxplot. Because a boxplot only reports summary statistics, a violin plot is more useful if you want to compare or show the presence of different peaks, their position, and relative amplitude.

Galton’s heights

data(galton, package = "UsingR")

galton_tidy <-  as_tibble(galton) %>%
  gather(person, height)

ggplot(galton_tidy, aes(height)) +
  geom_histogram(binwidth = 1) +
  facet_wrap(~ person, scales = "free")

galton contains heights (in inches) for 928 children and 205 “midparents”. Each parent height is an average of the father’s height and 1.08 times the mother’s height. The daughters’ heights are aslo said to have been multiplied by 1.08. This dataset was famously used by Sir Francis Galton to develop theories of correlation and regression in the 19th century.

Above we see two histograms plotted, one for the children and one for the parents. Each has the same binwidth (1 inch). While they look symmetric with no outliers, it is difficult to directly compare the histograms because they use different scales for the \(x\) and \(y\) axes.

ggplot(galton_tidy, aes(height)) +
  geom_histogram(binwidth = .1) +
  facet_wrap(~ person, scales = "free")

ggplot(galton_tidy, aes(height, fill = person)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

galton_median <- galton_tidy %>%
  group_by(person) %>%
  summarize(median = median(height))

ggplot(galton_tidy) +
  geom_histogram(aes(height), binwidth = 1) +
  geom_vline(data = galton_median, aes(xintercept = median), color = "red") +
  facet_wrap(~ person, ncol = 1)

ggplot(galton_tidy, aes(height, color = person)) +
  geom_density()

ggplot(galton_tidy, aes(person, height)) +
  geom_boxplot()

Outlier detection

ggplot2movies contains a dataset called movies, with information on over 58,000 movies collected from IMDB. One of the variables it collects is length, or the runtime of the film in minutes. What is the distribution of this continuous variable?

data(movies, package = "ggplot2movies")

ggplot(movies, aes(length)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

This plot is not very useful, except that it tells us we have at least one outlier, or extremely long movie, in the data. This is why the data appears heavily skewed. With over 58,000 movies in the first two bins, one case will never be visible. This histogram is not very useful here.

A boxplot, on the other hand, does a nice job of summarizing and visualizing outliers.

ggplot(movies, aes("", length)) +
  geom_boxplot() +
  labs(x = NULL)

There are actually quite a few outliers in the data. We could now use this information to filter the dataset and find out which movies these are: are they errors in coding the observation, or just an extremely long film?

filter(movies, length > 1000)

## # A tibble: 3 x 24
##   title  year length budget rating votes    r1    r2    r3    r4    r5
##   <chr> <int>  <int>  <int>  <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cure…  1987   5220     NA    3.8    59  44.5   4.5   4.5   4.5     0
## 2 Four…  1967   1100     NA    3      12  24.5   0     4.5   0       0
## 3 Long…  1970   2880     NA    6.4    15  44.5   0     0     0       0
## # ... with 13 more variables: r6 <dbl>, r7 <dbl>, r8 <dbl>, r9 <dbl>,
## #   r10 <dbl>, mpaa <chr>, Action <int>, Animation <int>, Comedy <int>,
## #   Drama <int>, Documentary <int>, Romance <int>, Short <int>

In fact all of these movies exist, and their lengths are accurate.

If we want to investigate and explore this dataset further, we have to decide how to proceed. Since 99% of all the observations in the dataset have run times lower than three hours, let’s just look at those observations. Compare the following graphs: which appears the most useful?

movies_lite <- filter(movies, length <= 180)

ggplot(movies_lite, aes("", length)) +
  geom_boxplot() +
  coord_flip() +
  labs(x = NULL,
       y = "Movie length (in minutes)")

ggplot(movies_lite, aes(length)) +
  geom_density() +
  labs(x = "Movie length (in minutes)")

ggplot(movies_lite, aes(length)) +
  geom_histogram(binwidth = 1) +
  labs(x = "Movie length (in minutes)")

I think the third plot provides the most information:

1. There are few long films (i.e. over 2 and a half hours)
2. There is a distinct group of short films with a peak of 7 minutes
3. There is a sharp peak at a length of 90 minutes.
4. There is clear evidence of round numbers being favored.

Compare distributions of subgroups

Recall the earlier example of using a density estimate to compare the distribution of heights for children and parents:

ggplot(galton_tidy, aes(height, color = person)) +
  geom_density()

This worked because we were only comparing between two subgroups. What happens if we want to compare across more than that? At what point does the density estimate not prove effective anymore?

Let’s consider movies and compare the distribution of film length by genre using a density estimate:

genre <- movies_lite %>%
  gather(genre, value, Action:Short) %>%
  filter(value == 1) %>%
  select(-value)

ggplot(genre, aes(length, color = genre)) +
  geom_density() +
  scale_color_brewer(type = "qual")

This is really messy - can we learn much from it? There still appear to be two peaks for each genre of film, but beyond that I cannot really tell a story from this. What are our alternative options?

ggplot(genre, aes(length)) +
  geom_histogram(binwidth = 1) +
  labs(x = "Movie length (in minutes)") +
  facet_wrap(~ genre)

ggplot(genre, aes(length)) +
  geom_histogram(binwidth = 1) +
  labs(x = "Movie length (in minutes)") +
  facet_wrap(~ genre, scales = "free_y")

ggplot(genre, aes(genre, length)) +
  geom_boxplot() +
  labs(x = "Genre",
       y = "Movie length (in minutes)")

Bar chart

A bar chart is pretty simple: one bar for each value of the variable, and the height of the bar is equal to the number of obsrvations falling in that category.

Let’s consider a sample of data from the General Social Survey, which is a long-running US survey conducted by NORC at the University of Chicago. The survey has thousands of questions about social and political behavior, and is one of the most important surveys in social science.

library(forcats)
gss_cat

## # A tibble: 21,483 x 9
##     year marital     age race  rincome   partyid   relig   denom   tvhours
##    <int> <fct>     <int> <fct> <fct>     <fct>     <fct>   <fct>     <int>
##  1  2000 Never ma…    26 White $8000 to… Ind,near… Protes… Southe…      12
##  2  2000 Divorced     48 White $8000 to… Not str … Protes… Baptis…      NA
##  3  2000 Widowed      67 White Not appl… Independ… Protes… No den…       2
##  4  2000 Never ma…    39 White Not appl… Ind,near… Orthod… Not ap…       4
##  5  2000 Divorced     25 White Not appl… Not str … None    Not ap…       1
##  6  2000 Married      25 White $20000 -… Strong d… Protes… Southe…      NA
##  7  2000 Never ma…    36 White $25000 o… Not str … Christ… Not ap…       3
##  8  2000 Divorced     44 White $7000 to… Ind,near… Protes… Luther…      NA
##  9  2000 Married      44 White $25000 o… Not str … Protes… Other         0
## 10  2000 Married      47 White $25000 o… Strong r… Protes… Southe…       3
## # ... with 21,473 more rows

A simple bar chart tells us the racial makeup of the survey respondents:

ggplot(gss_cat, aes(race)) +
  geom_bar()

So we can see that there are approximately 16,000 whites, 3000 blacks, and 2000 “other”.

ggplot(gss_cat, aes(race)) +
  geom_bar() +
  scale_x_discrete(drop = FALSE)

relig_summary <- gss_cat %>%
  group_by(relig) %>%
  summarise(
    age = mean(age, na.rm = TRUE),
    tvhours = mean(tvhours, na.rm = TRUE),
    n = n()
  )

ggplot(relig_summary, aes(relig, tvhours)) +
  geom_col() +
  coord_flip()

ggplot(relig_summary, aes(factor(relig, levels = rev(sort(levels(relig)))), tvhours)) +
  geom_col() +
  coord_flip()

ggplot(relig_summary, aes(fct_reorder(relig, tvhours), tvhours)) +
  geom_col() +
  coord_flip()

rincome_summary <- gss_cat %>%
  group_by(rincome) %>%
  summarise(
    age = mean(age, na.rm = TRUE),
    tvhours = mean(tvhours, na.rm = TRUE),
    n = n()
  )

ggplot(rincome_summary, aes(age, fct_reorder(rincome, age))) +
  geom_point()

ggplot(rincome_summary, aes(age, fct_relevel(rincome, "Not applicable"))) +
  geom_point()

data(gun_deaths, package = "rcfss")

gun_deaths %>%
  count(intent) %>%
  na.omit() %>%
  ggplot(aes(fct_reorder(intent, n, .desc = TRUE), n)) +
  geom_col()

gun_deaths %>%
  count(intent, race) %>%
  na.omit() %>%
  mutate(intent = factor(intent, levels = c("Suicide", "Homicide",
                                            "Accidental", "Undetermined"))) %>%
  ggplot(aes(intent, n, fill = race)) +
  geom_col()

gun_deaths %>%
  count(intent, race) %>%
  na.omit() %>%
  mutate(intent = factor(intent, levels = c("Suicide", "Homicide",
                                            "Accidental", "Undetermined"))) %>%
  ggplot(aes(intent, n, fill = race)) +
  geom_col(position = "dodge")

gun_deaths %>%
  count(intent, race) %>%
  na.omit() %>%
  mutate(intent = factor(intent, levels = c("Suicide", "Homicide",
                                            "Accidental", "Undetermined"))) %>%
  ggplot(aes(intent, n, fill = race)) +
  geom_col(position = "fill")

Pie chart

gun_deaths %>%
  count(intent) %>%
  na.omit() %>%
  mutate(pct = n / sum(n)) %>%
  with(pie(pct, labels = intent, clockwise = TRUE))

Ahh yes, the dreaded pie chart. Pie charts are proportional graphs which use angles, arc lengths, and area to identify the frequency of observations in each group. This is different from bar graphs, which use height of the bars (or length) to identify the frequency of observations in each group. Pie charts have received a lot of criticisms from academics and visualization experts. The key question is which is easier to distinguish: line length or area, angle, and arc length?

Spence, Ian, and Stephan Lewandowsky. “Displaying proportions and percentages.” Applied Cognitive Psychology 5.1 (1991): 61-77.

Experimental evidence suggests angles and arc lengths are more difficult to perceive in their relative size than something like height (i.e. a bar chart).

\[\text{Subjective area} = \text{Area}^.86\]

• People misjudge/misestimate area
• This function is not derived from rigorous science
• Focus more on the subarea (area/size of individual slices rather than the entire circle)
• People also do not generally use pie charts to estimate precise magnitudes (if so, use a table instead)

pie <- data_frame(label = c("A", "B", "C", "D"),
                  value = c(10, 20, 40, 30))

# bar chart
pie_bar <- ggplot(pie, aes(label, value)) +
  geom_col(color = "black", fill = "white") +
  theme_void2()

# proportional area chart
pie_prop <- pie %>%
  mutate(ymin = 0,
         ymax = 1,
         xmax = cumsum(value),
         xmin = xmax - value,
         x = (xmax - xmin) / 2 + xmin)

pie_prop <- ggplot(pie_prop, aes(xmin = xmin, xmax = xmax,
             ymin = ymin, ymax = ymax)) +
  geom_rect(fill = "white", color = "black") +
  scale_x_continuous(breaks = pie_prop$x, labels = pie_prop$label) +
  theme_void2()

# pie chart
pie_pie <- ggplot(pie, aes(x = factor(1), y = value)) +
  geom_col(width = 1, color = "black", fill = "white") +
  geom_text(aes(label = label), position = position_stack(vjust = .5), size = 4) +
  coord_polar(theta = "y", direction = 1) +
  theme_void() +
  theme(legend.position = "none")

pie_bar +
  pie_prop +
  pie_pie +
  plot_layout(ncol = 1)

# table
pie %>%
  select(value) %>%
  t %>%
  knitr::kable(col.names = pie$label,
               row.names = FALSE)

movies example

ggplot(movies, aes(votes, rating)) +
  geom_point()

Insights that could be found include:

1. No films with lots of votes and a low average rating
2. For films with more than a small number of votes, the average rating increases with the number of votes
3. No film with lots of votes has an average rating close to the maximum (suggests a barrier)
4. A few films with high number of votes (over 50,000) look like outliers. See the points below the trend
5. Films with low number of votes have high variance in rating
6. The only films with very high ratings also have relatively few votes

Smoothing lines

To help identify if a pattern or trend really exists in the scatterplot, we can overlay the plot with a smoothing line, or a line derived from a statistical model. This model could be linear or non-linear, parametric or non-parametric - we will discuss specific methods next term.

ggplot(movies, aes(votes, rating)) +
  geom_point() +
  geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

• Adds a non-linear smoothing line to the graph
• Suggests a positive, non-linear association between votes and rating
• Less confidence in our estimates as votes increase - fewer observations from which to generalize

Adding jitter to the graph

ggplot(galton, aes(parent, child)) +
  geom_point()

Sometimes continuous variables don’t have many unique values due to rounding, such as Galton’s height dataset. In this situation, we can draw a scatterplot but we cannot see how many unique values actually exist in the dataset. One alternative is to jitter each value, or add a bit of random noise to the \(x\) and \(y\) coordinates for each data point. Because the noise is added at random, it should not unfairly distort or bias the graph.

ggplot(galton, aes(parent, child)) +
  geom_jitter()

We can now see that some combinations of parent-child heights are more common than others. Is there a positive trend?

ggplot(galton, aes(parent, child)) +
  geom_jitter(alpha = .2) +
  geom_abline(slope = 1, linetype = 2) +
  geom_smooth(method = "lm")

Here I added a linear (OLS) smoothing line, add a dashed reference line for a perfect 1-to-1 relationship between parent and child height, and fade the data points by adding some transparency, to make the graph less messy. We can see the relationship between parent and child height is positive and somewhat close to 1-to-1, though not precisely.

Note that smoothing lines are not definitive proof of correlation or causation. The relationship seen in the smoothing line could be driven by some third, ungraphed variable. Further exploration is certainly warranted, as well as theory to back up your claims.

Comparing groups within scatterplots

If you want to check for conditional relationships using scatterplots, you have a few options. If there is a small number of groups, you could add a color aesthetic to the graph to differentiate between each group.

ggplot(data = mpg, mapping = aes(x = displ, y = hwy, color = drv)) +
  geom_point()

Here we have a dataset of automobiles and want to understand the relationship between engine size and highway fuel efficiency, conditional on the drive-type. By changing the color of each data point based on its group membership, we can see a few things:

1. Front wheel drive cars do not have fuel efficiency below 20 MPG
2. Rear wheel drive cars tend to have larger engines
3. Within each group, the relationship still tend to be negative

However, as the number of groups increases, it becomes more difficult to visually distinguish between each group.

ggplot(data = mpg, mapping = aes(x = displ, y = hwy, color = class)) +
  geom_point()

With seven categories, this is becoming more tricky. What if we have even more categories? Consider a dataset of individual Olympians who competed in the 2012 Summer Olympics.

data(oly12, package = "VGAMdata")

ggplot(oly12, aes(Height, Weight, color = Sport)) +
  geom_point()

## Warning: Removed 1346 rows containing missing values (geom_point).

Too many categories. Instead, build a facet graph.

oly12 %>%
  mutate(Sport = abbreviate(Sport, 12)) %>%
  ggplot(aes(Height, Weight)) +
  geom_point(size = 1) +
  facet_wrap(~ Sport)

## Warning: Removed 1346 rows containing missing values (geom_point).

While still a bit cramped, it helps. Now we can see some different features of the data.

1. Some sports have no height/weight information for competitors
2. Some sports have far fewer competitors than others
3. Across sports, there still seems to be a positive relationship between height and weight

If you want to focus your comparison on just a handful of categories, filter the data first and then draw the plot.

oly12 %>%
  filter(Sport %in% c("Judo", "Weightlifting", "Wrestling")) %>%
  ggplot(aes(Height, Weight)) +
  geom_point(size = 1) +
  facet_wrap(~ Sport)

## Warning: Removed 69 rows containing missing values (geom_point).

Scatterplot matrix

To quickly visualize several variables in a dataset and their relation to one another, a scatterplot matrix is a quick and detailed tool for generating a series of scatterplots for each combination of variables. Consider credit.csv which contains a sample of individuals from a credit card company, identifying their total amount of credit card debt and other financial/demographic variables:

credit <- read_csv("data/Credit.csv") %>%
  # remove first ID column
  select(-X1)

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_integer(),
##   Income = col_double(),
##   Limit = col_integer(),
##   Rating = col_integer(),
##   Cards = col_integer(),
##   Age = col_integer(),
##   Education = col_integer(),
##   Gender = col_character(),
##   Student = col_character(),
##   Married = col_character(),
##   Ethnicity = col_character(),
##   Balance = col_integer()
## )

names(credit) <- stringr::str_to_lower(names(credit))   # convert column names to lowercase

If we want to quickly assess the relationship between all of the variables (in preparation for more advanced statistical learning techniques), we could generate a matrix of scatterplots using the base graphics package:

# only graph numeric variables
pairs(select_if(credit, is.numeric))

We can use this to look for trends and associations between each set of variables. This base graphics implementation is very minimal, as it is strictly for numeric variables and only shows tinyscatterplots.

We can also use GGally::ggpairs() to generate a scatterplot matrix. GGally is a package for R that extends ggplot2 by adding helper functions for common multivariate data structures. ggpairs() is a function that allows us to quickly generate a scatterplot matrix.

library(GGally)

## 
## Attaching package: 'GGally'

## The following object is masked from 'package:dplyr':
## 
##     nasa

ggpairs(select_if(credit, is.numeric))

When applied to strictly numeric variables, the lower triangle generates scatterplots, the upper triangle prints the correlation coefficient, and the diagonal panels are density plots of the variable.

ggpairs() also works on datasets with a mix of qualitative and quantitative variables, drawing appropriate graphs based on whether the variables are continuous or discrete:

ggpairs(select(rcfss::scorecard, type:debt))

Heatmap of correlation coefficients

Scatterplot matricies can provide lots of information, but can also be very densely packed. Perhaps instead we want to quickly visualize the correlation between each of the variables. GGally also contains the ggcorr() function which calculates Pearson’s correlation coefficient for each combination of continuous variables and visualizes them using a heatmap.

ggcorr(mpg)

## Warning in ggcorr(mpg): data in column(s) 'manufacturer', 'model', 'trans',
## 'drv', 'fl', 'class' are not numeric and were ignored

ggcorr(credit)

## Warning in ggcorr(credit): data in column(s) 'gender', 'student',
## 'married', 'ethnicity' are not numeric and were ignored

ggcorr(diamonds)

## Warning in ggcorr(diamonds): data in column(s) 'cut', 'color', 'clarity'
## are not numeric and were ignored

We quickly see for each dataset which variables are correlated most strongly with one another. ggcorr() also allows customization of the appearance of the heatmap. For example, we can display the actual correlation coefficient directly in the graph by setting label = TRUE:

ggcorr(credit, label = TRUE)

## Warning in ggcorr(credit, label = TRUE): data in column(s) 'gender',
## 'student', 'married', 'ethnicity' are not numeric and were ignored

ggcorr(credit, label = TRUE, label_round = 2, label_alpha = TRUE)

## Warning in ggcorr(credit, label = TRUE, label_round = 2, label_alpha =
## TRUE): data in column(s) 'gender', 'student', 'married', 'ethnicity' are
## not numeric and were ignored