Discover possibilities

When I get my hands on new data, I always spend time on descriptives. Often, this exploration is fairly unstructured. I explore each variable separately in univariate statistics. I then start to relate one to the other in bivariate and trivariate statistics. Some researchers claim that if the relationship they study is not visible in the bivariate, it is not a project worth pursuing! That’s a pretty extreme approach.

The point is to get a sense of what information is contained in the data: Variation is information. Variables with a lot of variation therefore contain more information.

My choice of model furthermore depends on how the dependent variable is measured and the structure of the data. Certainly, I can read the codebook (or design it), but I may also perceive new modeling opportunities by knowing what information is available to me. To begin with, I therefore explore the phenomenon I’m interested in (\(y\)).

For example, you may go into a project where you test some theoretical argument about when a government chooses to bomb their enemy. A possible operationalization of your dependent variable is the number of air strikes in a particular area. However, you might also realize that you have a time variable (date, time of day) of each strike. Now, you can also consider an event history model that analyzes each strike separately; it’s a more fine-grained analysis.

Flag potential problems

Data come with all kinds of problems and challenges that you have to handle. It is not always easy to know which problem you might have in the future, but by knowing the data, you will know the quirks of each variable and more easily spot problems (and solutions) when you run into them.

For example, there can also be too much variation. I could, for example, discover that my continuous variable on respondent income has a few very rich people in it. These are outliers. They may or may not skew your analysis later on. Thanks to the descriptive statistics, you will have the potential problem on your radar.

Error probe

Analyzing data usually involves a lot of data wrangling: You can create new variables from different information sources or change the measurement level of a variable. You can move up a notch by making indexes to get a more fine-grained variable or you can move down a measurement level because there is not enough/too much variation in the original data. After changing my data, I always check if my code resulted in what I wanted. The best way to do this is through some visualization of the univariate distribution of the new data.

Interprets and communicate

Knowledge of your descriptive statistics will come back once you’ve run your statistical model. Now you have to interpret it. It is time to figure out what scenarios are realistic and interesting to discuss. By that, I mean that you need to consider which values on the \(x\)-variables are common and what changes in these variables are likely. Now you can use your model to predict what the outcome is. This is at the very heart of your analytical venture. Descriptive statistics therefore helps you understand what you have found.

For example, it may be that your regression analysis reports a small coefficient for the effect of income on scepticism towards immigration. This does not automatically mean that the effect is not substantial. If your income is measured in DKK, then a one-unit increase is not a lot. However, if you see a salary bump of 100.000, then the shift – even if the regression coefficient is small – might also be substantial.

Getting started

We start by loading in the data packages we will use. If you haven’t already installed ggplot2, you will have to do that first.

#install.packages("ggplot2")
library(dplyr); library(ggplot2)

We will use the same data as in the previous session: It is available in my package or online on my website.

library(RiPraksis)
data(kap6)
df <- kap6

Main logic

The workflow in ggplot2 is a pipe. This pipe is read chronologically: I first establish a sheet, then I add elements.

1. specify the data you want to work with

2. establish a blank sheet

3. add elements to the sheet one by one

4. the visual adaptations (“theme”) can be based on a template (there is a grey default temple). The template already specifies all the esthetics for you. However, you can “tweak” what you want after having established the template.

Google and ChatGPT are your best friends, as usual.

Example Let’s make a scatterplot between income and scepticism.

The main function, ggplot(), establishes the plot, but no plot elements. If you are going to use the same data for all the graphical elements, you can specify the data object in this function. When you do that, you do not specify your dataset multiple times in the pipe, so you save time. From then on, you can refer directly to the variables in it.

You bind all elements together with a pipe using “+” when you use ggplot. However, you can mix tidyverse and ggplot dialects. Here, I start by a tidyverse pipe, then continue with ggplot to make an empty sheet.

df %>%
  ggplot()

The main elements of the graphics are included using a series of geom_...() functions. They specify the “geometry,” which is what information should be included (information from variables) and where it should be placed (coordinates).

Inside each of these functions, you will find the aesthetic function, aes(), where you specify the x/y coordinates that come from your variables.

Here, I am creating a scatterplot for the relationship between immigration skepticism and income. This means that I want to draw points (geom_point()).

df %>%
  ggplot +
  geom_point(aes(y = Skepsis,
                 x = Indtaegt))

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

A simple scatter plot gives an impression of the bivariate relationship between immigration skepticism and income.

Numeric summary: mean, quantiles…

Summarizing continuous variables usually means looking at their average, but also their minimum, maximum values and their spread.

summary(df$Skepsis)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   4.000   5.000   5.009   6.000  10.000       5

I can also check out their standard deviation to get a feel for how much variation I have in them.

sd(df$Skepsis, na.rm = TRUE)

## [1] 1.72577

Visual: histogram

base

df$Skepsis %>% hist

We can make it relative using a simple additional argument.

df$Skepsis %>% hist(., probability = T)

ggplot2 Histograms are easily implemented in ggplot too.

df %>%
  ggplot +
  geom_histogram(aes(Skepsis))

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

## Warning: Removed 5 rows containing non-finite values (`stat_bin()`).

You can tweak them any which way you want, but it requires some more coding skills. Here, I make the histogram relative by manipulating the.

df %>%
  ggplot(aes(x = Skepsis)) +
  geom_histogram(aes(y = after_stat(density)))

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

## Warning: Removed 5 rows containing non-finite values (`stat_bin()`).

Numeric summary: frequency table

Here, we are looking at what our respondents voted last election.

df$Parti %>%
  table

## .
##                       Andet            Dansk Folkeparti 
##                          17                         143 
## Det Konservative Folkeparti        Det Radikale Venstre 
##                          65                         134 
##                Enhedslisten         Kristendemokraterne 
##                          77                           8 
##            Liberal Alliance          Socialdemokraterne 
##                          48                         268 
##     Socialistisk Folkeparti                     Venstre 
##                         108                         311

Visual: barplot

df$Parti %>%
  table %>%
  barplot()

df$Parti %>%
  #Frequency table
  table %>%
  #Relative distribution
  prop.table %>%
  #Barplot
  barplot()

df %>%
  ggplot +
  geom_bar(aes(Parti))

df %>%
  #Filter out missing observations before plotting
  filter(!is.na(Parti)) %>%
  ggplot() +
  #Barplot
  geom_bar(aes(x = Parti))

Numeric summary: correlation

When we have two numeric variables, we can calculate their correlation using Pearsons R (\(r\)).

cor(df$Indtaegt, df$Skepsis,
    use = "pairwise.complete.obs")

## [1] -0.1526337

The negative correlation is apparent, but what does it mean? Can we put it into plain English?

If we take the square of this correlation, we get the “coefficient of determination” (\(R^2\)). We can interpret the result as a proportion of shared variation between the two variables. Sometimes we’d interpret it as a causal relationship.

cor(df$Indtaegt, df$Skepsis,
    use = "pairwise.complete.obs")^2

## [1] 0.02329706

We may say that 2 % of the variation in attitudes towards immigration is shared with (explained by?) income.

We can test the statistical significance between the two.

cor.test(df$Indtaegt, df$Skepsis)

## 
##  Pearson's product-moment correlation
## 
## data:  df$Indtaegt and df$Skepsis
## t = -5.6155, df = 1322, p-value = 2.387e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.20482310 -0.09957898
## sample estimates:
##        cor 
## -0.1526337

Considering numbers in this way is really abstract though. The most intuitive way of exploring a correlation is through visuals.

Visual: scatterplot

When both variables are numeric, I can create a scatterplot for the relationship between immigration skepticism and income.

df %>%
  select(Indtaegt, Skepsis) %>%
  plot

ggplot2

I can do the same thing in ggplot and take the advantage of its additional functions to make the plot more legible.

I want to draw points, so I add the elements in the geom_point() function to the blank sheet.

df %>%
  ggplot +
  geom_point(aes(y = Skepsis,
                 x = Indtaegt))

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

A simple scatter plot gives an impression of the bivariate relationship between immigration skepticism and income.

It is hard to see any relationship between the two variables. Also, they only shared some 2% of their variation. However, the income variable also consists only of integers. Therefore, the points overlap each other, and we get a poor idea of how many units exist within the different value combinations. We can’t see what value combinations are the most common, only that there are a lot of them.

We can solve this in two ways: by making the points partially transparent, or by moving them slightly apart.

Transparent color I can specify that I want transparent colors. I do this with the argument alpha=. The value I specify ranges from 0 (completely transparent/invisible) to 1 (not transparent).

df %>%
  ggplot +
  geom_point(aes(y = Skepsis,
                 x = Indtaegt),
             #all observations are transparent
             alpha = 0.2)

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

By using a transparent color, I can show where the majority of the observations are located in the distribution.

Note that when we want to adjust the points without the adjustment being dependent on information from the data, this is done outside of aes() but within the geom_…() function.

Jitter data Another alternative is to “jitter” the data points a little. This means that R adds a random variation to the coordinates of the points. The actual data points become more imprecise, but the point here is not to perform a precise analysis, but to present relationships in the data so that the human eye gets an intuitive understanding of what is going on.

df %>%
  ggplot +
  geom_jitter(aes(y = Skepsis,
                  x = Indtaegt),
              #Horizontal, but not vertical shake
              height = 0, width = 0.3)

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

I get a better grasp of where my observations are by ‘shaking’ the data points.

The additional arguments width= and height= specify how much the points should be jittered horizontally and vertically, respectively. Here, I say that I want precise data along the y-axis (“in height”), but I let the points randomly vary by 0.3 units along the x-axis (“in width”).

Trivariate graphics: Group the data.

I can easily add more groupings. Here I am coloring the points based on the value of a third variable: The variable can be continuous or categorical.

Numeric/continous grouping “Prekaritet” flags respondents that have a weak link to the labor market (i.e. they are unemployed or on a temporary work contract).

df %>%
  ggplot +
  geom_point(aes(y = Skepsis,
                 x = Indtaegt,
                 #Grouping
                 color = Prekaritet))

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

Now I get a “heat scale” for the respondent’s link to the labor market. ggplot also automatically provides a legend. Now, the variable only has two values, 0 or 1. We better use it as a categorical moderator.

Categorical grouping If I want a categorical grouping, I need to feed the function a categorical variable (“character” or “factor”). It’s best to do this by mutating the data frame using tidyverse, not in the ggplot() function. Right here, I use as.factor() to change the measurement level into categorical.

df %>%
  mutate(Prekaritet = as.factor(Prekaritet)) %>%
  ggplot +
  geom_point(aes(y = Skepsis,
                 x = Indtaegt,
                 #Grouping
                 color = Prekaritet))

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

Of course, I can combine all these solutions – jittering, transparency and recoding – until I get a plot that tells the story I want to convey.

#tidyverse to transform the data
df %>%
  mutate(Prekaritet = as.factor(Prekaritet)) %>%
  #ggplot2 to plot the data
  ggplot +
  geom_jitter(aes(y = Skepsis,
                 x = Indtaegt,
                 #Continuous moderating third variable
                 color = Prekaritet),
              #Shake data data horizontaly and vertically
              height = 0.1, width = 0.4,
              #Transparency
              alpha = 0.7)

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

It is hard to tell much from this graphic, so let’s move on to other visualizations.

Bivariate regression lines

I can also show the relationship between two variables using regression lines. ggplot does not report regression coefficients but instead illustrates the regression line and the uncertainty around it. This is particularly relevant if at least one of the variables (the x-variable) is continuous.

Now we can choose between local regression and ordinary linear regression. Which one we choose depends on what we want to achieve.

Local regression

Local regression is – in my eyes – the queen of all ggplot functions. I use it all the time when exploring my data. It draws a “smooth” line (locally, the sliding average over the different x-values). Perfectly suited for time trends or exploring non-linear relationships, for example.

The function that gives us bivariate regression lines is geom_smooth(). The default setting is local regression.

df %>%
ggplot +
  geom_smooth(aes(y = Skepsis,
                  x = Indtaegt))

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

## Warning: Removed 178 rows containing non-finite values (`stat_smooth()`).

Local regression: A simple way to present empirical relationships between two variables.

Here I illustrate the relationship between income and skepticism about immigration. ggplot adjusts the axis limits for us automatically so as to zoom in on the results. It means that they always come across as large effects.

Grouped: Is the relationship the same within different categories?

df %>%
  ggplot +
  geom_smooth(aes(y = Skepsis,
                  x = Indtaegt,
                  color = as.factor(Prekaritet)))

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

## Warning: Removed 178 rows containing non-finite values (`stat_smooth()`).

Local regression: The queen of trivariate plots?

Would you like to replace the “spaghetti” with straight lines? Ask for a (bivariate) linear model.

df %>%
  ggplot +
  geom_smooth(aes(y = Skepsis,
                  x = Indtaegt,
                  color = as.factor(Prekaritet)),
              method = "lm")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 178 rows containing non-finite values (`stat_smooth()`).

Here we clearly see that the association between immigration skepticism and income is primarily present for respondents who have a safe connection to the labor market.

Numeric summary: group averages

This is a good excuse to have a second look at tidyverse . dplyr has a very handy function that allows us to do operations on subset of the data depending on a moderating variable.

Let’s calculate the average scepticism among voters of different parties and create a summary statistic. To do so, we use reframe()

df %>%
  group_by(Parti) %>%
  reframe("Scepticism" = mean(Skepsis, na.rm = T))

Before moving on, I’m removing the NAs.

df %>%
  filter(!is.na(Parti)) %>%
  group_by(Parti) %>%
  reframe("Scepticism" = mean(Skepsis, na.rm = T))

The observations are ordered alphabetically. It is not intuitive. Let’s rearrange.

df %>%
  filter(!is.na(Parti)) %>%
  group_by(Parti) %>%
  reframe("Scepticism" = mean(Skepsis, na.rm = T)) %>%
  arrange(Scepticism)

No surprises here; Dansk Folkeparti has voters that are more sceptical than other parties. This being said, in general, there is not a massive spread here.

My results looked ok, so I’m storing them in tab.

tab <-
  df %>%
  filter(!is.na(Parti)) %>%
  group_by(Parti) %>%
  reframe("Scepticism" = mean(Skepsis, na.rm = T)) %>%
  arrange(Scepticism)

Visual: barplot

Let’s instead plot this.

barplot(tab$Scepticism,
        names.arg = tab$Parti,
        las = 2)

tab %>%
  ggplot +
  geom_col(aes(y = Scepticism,
               x = Parti))

tab %>% 
  mutate(Parti = reorder(Parti, Scepticism)) %>% 
  ggplot +
  geom_col(aes(y = Scepticism,
               x = Parti))

ggplot2

df %>%
  #Create a new variable with categorical names
  mutate(Gender = if_else(Kvinde == 1,
                         "Woman",
                         "Man")) %>%
  #Plot
  ggplot +
  geom_bar(aes(x = Parti,
               #Moderate by gender
               fill = Gender),
           position = "dodge")

Define colors

You can choose your own colors, and there are many ways to do so. If you have linked your data with color choices, it may be useful to define a “palette.”

The text (“Precarious” and “Safe”) is part of the data I provided, so it is located in the variable “dfp$Worksituation.” I can change it in ggplot if I want to. I do this when defining the colors.

p <- 
  p +  
  #Define colors
  scale_color_manual(
    #What colors?
    values = c("purple", "magenta"),
    #Color the entire bar
    aesthetics = "fill",
    #What's the title of the legend?
    name = "Link to labor market",
    #What are the category names?
    labels = c("Stable",
               "Unstable")
    ) 
p

Background

You can define different templates for the aesthetic aspects of the plot with theme_...().

p <- 
  p +
  # Define a theme
  theme_minimal()
p

Here, I am using ready-made functions to modify the visual expression of the plot.

You can alter the themes at will. You do that after you have defined the theme, using the generik theme() function.

p <-
  p +
  # Modify the minimal theme after defining it
  theme(
    # Move the legend down to below the plot
    legend.position = "bottom",
    # Bold font for the title in the legend
    legend.title = element_text(face = "bold"),
    # Italic for axis values
    axis.text = element_text(face = "italic"),
    # Remove the gridlines
    panel.grid = element_blank())
p

Here, I’ve changed the background and other thematical elements in my graphic.

The element_...() functions are used inside the theme() function. They specify everything except the data content/information for an element. Should it be a blank element element_blank()? Specific font? element_text()… etc.

For example, I can specify that the text used for the axis names should be displayed in italic element_text(face = "italic").