The world is full of obvious things which nobody by any chance observes



Quote from Arthur Conan Doyle, The Hound of the Baskervilles

Outline

  • The humble but powerful scatterplot
  • Additions and variations
  • Transformations to linearity
  • (Robust) numerical measures of association
  • Simpson’s paradox
  • Making null samples to test for association
  • Imputing missings
  • Multiple associations, exploring network data

The scatterplot

Scatterplots are the natural plot to make to explore association between two continuous (quantitative or numeric) variables.

They are not just for linear relationships but are useful for examining nonlinear patterns, clustering and outliers.

We also can think about scatterplots in terms of statistical distributions: if a histogram shows a marginal distribution, a scatterplot allows us to examine the bivariate distribution of a sample.

History

  • Descartes provided the Cartesian coordinate system in the 17th century, with perpendicular lines indicating two axes.
  • It wasn’t until 1832 that the scatterplot appeared, when John Frederick Herschel plotted position and time of double stars.
  • Kopf argues that The scatter plot, by contrast, proved more useful for scientists, but it clearly is useful for economics today.

http://www.datavis.ca/milestones/

Language and terminology

Are the words “correlation” and “association” interchangeable?

In the broadest sense correlation is any statistical association, though it commonly refers to the degree to which a pair of variables are linearly related. Wikipedia



If the relationship is not linear, call it association, and avoid correlated.

Features of a pair of continuous variables (1/3)

Feature Example Description
positive trend Low value corresponds to low value, and high to high.
negative trend Low value corresponds to high value, and high to low.
no trend No relationship
strong Very little variation around the trend
moderate Variation around the trend is almost as much as the trend
weak A lot of variation making it hard to see any trend

Features of a pair of continuous variables (2/3)

Feature Example Description
linear form The shape is linear
nonlinear form The shape is more of a curve
nonlinear form The shape is more of a curve
outliers There are one or more points that do not fit the pattern on the others
clusters The observations group into multiple clumps
gaps There is a gap, or gaps, but its not clumped

Features of a pair of continuous variables (3/3)

Feature Example Description
barrier There is combination of the variables which appears impossible
l-shape When one variable changes the other is approximately constant
discreteness Relationship between two variables is different from the overall, and observations are in a striped pattern
heteroskedastic Variation is different in different areas, maybe depends on value of x variable
weighted If observations have an associated weight, reflect in scatterplot, e.g. bubble chart

Additional considerations (Unwin, 2015):

  • causation: one variable has a direct influence on the other variable, in some way. For example, people who are taller tend to weigh more. The dependent variable is conventionally on the y axis. It’s not generally possible to tell from the plot that the relationship is causal, which typically needs to be argued from other sources of information.
  • association: variables may be related to one another, but through a different variable, eg ice cream sales are positively correlated with beach drownings, is most likely a temperature relationship.
  • conditional relationships: the relationship between variables is conditionally dependent on another, such as income against age likely has a different relationship depending on retired or not.

Famous data examples

Famous scatterplot examples

Anscombe’s quartet

All four sets of Anscombe has same means, standard deviations and correlations, \(\bar{x}\) = 9, \(\bar{y}\) = 7.5, \(s_x\) = 3.3, \(s_y\) = 2, \(r\) = 0.82.


Numerical statistics are the same, for very different association.

Datasaurus dozen

And similarly all 13 sets of the datasaurus dozen have same means, standard deviations and correlations, \(\bar{x}\) = 54, \(\bar{y}\) = 48, \(s_x\) = 17, \(s_y\) = 27, \(r\) = -0.06.

Aspect ratio

Notice that the scatterplots of the historical examples were made using equal aspect ratio, that is the plot space is a SQUARE not a rectangle.

Most software and most scatterplots you see use unequal aspect ratio. This is an ERROR that is ubiquitous!

When showing two variables, to read association effectively REQUIRES the SQUARE (1:1) aspect ratio.

These slides have more than 50 uses of aspect.ratio=1.

Make sure you don’t fall into the trap of failing to control aspect ratio.

Scatterplot case studies

Activity: Olympics

Your turn

  1. Load the individual athlete data for the Summer 2012 Olympic Games (downloaded from the Guardian’s London 2012 Olympics data)
data(oly12, package = "VGAMdata")

  1. Make a scatterplot of Height vs Weight. You’ll see a warning about missing values being removed – what features do you notice? Look for a linear relationship, outliers, discretization, and overplotting (there are more than 10,000 athletes)

  2. Are there sport(s) where you would expect a different relationship between height and weight?

  3. Count how many athletes compete in each Sport. Are the categories well organised, or would you group some of them together?

Code 
ggplot(oly12, aes(x = Height, y = Weight, label = Sport)) +
  geom_point()

skimr::skim(oly12)

oly12 |>
  count(Sport, sort = TRUE)

# Consolidate factor levels
oly12 <- oly12 |>
  mutate(Sport = as.character(Sport)) |>
  mutate(Sport = ifelse(grepl("Cycling", Sport), 
    "Cycling", Sport
  )) |> 
  mutate(Sport = ifelse(grepl("Gymnastics", Sport),
    "Gymnastics", Sport
  )) |>
  mutate(Sport = ifelse(grepl("Athletics", Sport),
    "Athletics", Sport
  )) |>
  mutate(Sport = as.factor(Sport))

# Drill down by sport
ggplot(oly12, aes(x = Height, y = Weight)) +
  geom_point(alpha = 0.5) + 
  facet_wrap(~Sport, ncol = 8) +
  theme(aspect.ratio = 1) 

# Handle missings, and examine by demographics
oly12 |>
  filter(!(Sport %in% c("Boxing", "Gymnastics", "Synchronised Swimming", "Taekwondo", "Trampoline"))) |>
  mutate(Sport = fct_drop(Sport)) |>
  ggplot(aes(x = Height, y = Weight, colour = Sex)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~Sport, ncol = 7, scales = "free") +
  scale_colour_brewer("", palette = "Dark2") +
  theme(aspect.ratio = 1, axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Did you notice these interesting patterns?

  • Some sports have no data for height, weight
  • The positive association between height and weight is visible across sports
  • Nonlinear in wrestling?
  • An outlier in judo, and football, and archery
  • Maybe flatter among swimmers
  • Taller in basketball, volleyball and handball
  • Shorter in athletics, weightlifting and wrestling
  • Little variance in tennis players

Refine plots to make comparisons easier:

  • Remove sports with missings
  • Overlay model to assess linear relationships
  • Drill down into other strata, eg male/female athletes
  • Compare in small chunks, one group against the rest

Note: Because the focus is now on males vs females association shape within sport, make plots scale separately.

  • Athletics category should have been broken into several more categories like track, field: a shot-putter has a very different physique to a sprinter.
  • Generally, clustering of male/female athletes
  • Outliers: a tall skinny male archer, a medium height very light female athletics athlete, tall light female weightlifter, tall light male volleyballer
  • Canoe slalom athletes, divers, cyclists are tiny.

Common ways to augment scatterplots (1/2)

Modification Example Purpose
alpha-blend alleviate overplotting to examine density at centre
model overlay focus on the trend
model + data trend plus variation
density overall distribution, variation and clustering
filled density high density locations in distribution (modes), variation and clustering
colour relationship with conditioning and lurking variables

Common ways to augment scatterplots (2/2)

Hexagon binning

Colour maps to the count of points falling in each hexagonal bin, alleviating overplotting.

Bubble plot

Size of each point maps to an extra (weight) variable.

High density regions

The gghdr package colours points by highest density regions (HDRs), revealing modes and labelling outliers beyond the largest region.

Code 
ggplot(df, aes(x = x2, y = y3)) +
  geom_hex(bins = 15) +
  scale_fill_continuous_sequential(palette = "YlGnBu") +
  theme(aspect.ratio = 1, legend.position = "none")
Code 
library(survey)
data(api)
ggplot(apiclus2, 
  aes(x=hsg, 
      y=pct.resp, 
      size=pw)) + 
  geom_point() +
  ylab("Response %") +
  xlab("Parent high school grad %") +
  theme(aspect.ratio = 1, legend.position = "none")
Code 
# install.packages("remotes")
# remotes::install_github("Sayani07/gghdr")
library(gghdr)
ggplot(data = faithful, 
    aes(x = waiting, 
        y = eruptions)) +
  geom_point(aes(colour = hdr_bin(x = waiting, y = eruptions))) +
  #geom_hdr_rug() +
  scale_colour_viridis_d(direction = -1) +
  theme(aspect.ratio = 1, legend.position = "none")

Activity: Comparing association and spread

Your turn

  1. Using oly12, filter to female athletes in Swimming, Archery, Basketball, Handball, Hockey, Tennis, Weightlifting, and Wrestling. Fit a linear trend of Weight on Height for each Sport and plot them together. How does the association between height and weight differ across sports?

  2. Now filter to female athletes in Shooting, Modern Pentathlon, and Basketball, and overlay density contours of Height and Weight by Sport. How does the spread of the relationship differ across sports?

  3. Some of this association may be due to unmeasured (lurking) variables. For example, “Athletics” combines throwing and running events, which have quite different body types. If you were only given Height and Weight, without Sport, could you have detected this hidden structure? What pattern might reveal it?

  4. Make a scatterplot and overlaid density contours (try a couple of different bandwidths) of Height vs Weight for Sport == "Athletics" only. Can you detect any multimodality that would hint at the hidden sport subgroup structure?

Code 
oly12 |>
  filter(Sport %in% c(
    "Swimming", "Archery", "Basketball",
    "Handball", "Hockey", "Tennis",
    "Weightlifting", "Wrestling"
  )) |>
  filter(Sex == "F") |>
  mutate(Sport = fct_drop(Sport), Sex = fct_drop(Sex)) |>
  ggplot(aes(x = Height, y = Weight, colour = Sport)) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_discrete_divergingx(palette = "Zissou 1") +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.direction = "horizontal"
  )

oly12 |>
  filter(Sport %in% c("Shooting", "Modern Pentathlon", "Basketball")) |>
  filter(Sex == "F") |>
  mutate(Sport = fct_drop(Sport), Sex = fct_drop(Sex)) |>
  ggplot(aes(x = Height, y = Weight, colour = Sport)) +
  geom_density2d() +
  scale_color_discrete_divergingx(palette = "Zissou 1") +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.direction = "horizontal"
  )

p1 <- oly12 |>
  filter(Sport == "Athletics") |>
  ggplot(aes(x = Height, y = Weight)) +
  geom_point(alpha = 0.2, size = 4)
p2 <- oly12 |>
  filter(Sport == "Athletics") |>
  ggplot(aes(x = Height, y = Weight)) +
  geom_density2d_filled() +
  theme(legend.position = "none")
p3 <- oly12 |>
  filter(Sport == "Athletics") |>
  ggplot(aes(x = Height, y = Weight)) +
  geom_density2d(binwidth = 0.01)
p4 <- oly12 |>
  filter(Sport == "Athletics") |>
  ggplot(aes(x = Height, y = Weight)) +
  geom_density2d(binwidth = 0.001, color = "white", size = 0.2) +
  geom_density2d_filled(binwidth = 0.001) +
  theme(legend.position = "none")
grid.arrange(p1, p3, p2, p4, ncol = 2)
  • Weightlifters are much heavier relative to height, swimmers are leaner, and tennis players are mixed (shorter athletes tend to be heavier, taller athletes tend to be lighter): the association between height and weight differs by sport.
  • Modern pentathlon athletes are uniformly related in height and weight, while shooters are quite varied in body type: the spread of the relationship also differs by sport.
  • Some of the association is due to unmeasured variables. For example, “Athletics” masks different body types in throwing vs running events — this is a lurking variable.
  • With only Height and Weight, the hidden structure is hard to detect directly, but it may appear as multimodality in the joint distribution — as seen (faintly) in the density contours for Athletics. It’s not easy to detect the presence of the additional variable, and thus accurately describe the relationship between height and weight among Olympic athletes.

Numerical measures of association

Correlation

  • Correlation between variables \(x_1\) and \(x_2\), with \(n\) observations in each.

\[r = \frac{\sum_{i=1}^n (x_{i1}-\bar{x}_1)(x_{i2}-\bar{x}_2)}{\sqrt{\sum_{i=1}^n(x_{i1}-\bar{x}_1)^2\sum_{i=1}^n(x_{i2}-\bar{x}_2)^2}} = \frac{\mbox{covariance}(x_1, x_2)}{(n-1)s_{x_1}s_{x_2}}\]

  • Test for statistical significance, whether population correlation could be 0 based on observed \(r\), using a \(t_{n-2}\) distribution:

\[t=\frac{r}{\sqrt{1-r^2}}\sqrt{n-2}\]

Code 
cor(d1$x, d1$y)
[1] 0.52
Code 
cor.test(d1$x, d1$y)

    Pearson's product-moment correlation

data:  d1$x and d1$y
t = 9, df = 198, p-value = 2e-15
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.41 0.62
sample estimates:
 cor 
0.52

Problems with correlation (1/2)

Code 
cor(d2$x, d2$y)
[1] -0.05
Code 
cor.test(d2$x, d2$y)

    Pearson's product-moment correlation

data:  d2$x and d2$y
t = -0.7, df = 198, p-value = 0.5
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.187  0.089
sample estimates:
  cor 
-0.05


It does not summarise non-linear associations.

Problems with correlation (2/2)

All observations

$estimate
cor 
0.3 

$statistic
  t 
4.4 

$p.value
[1] 1.6e-05

Without outlier

$estimate
   cor 
-0.012 

$statistic
    t 
-0.17 

$p.value
[1] 0.87


It is affected by extreme values.

Perceiving correlation

Let’s play a game: Guess the correlation!


Generally, people don’t do very well at this task. Typically people under-estimate \(r\) from scatterplots, particularly when it is around 0.4-0.7. The variation in a scatterplot perceptually doesn’t vary is not linearly with \(r\).

When someone says correlation is 0.5 it sounds impressive. BUT when someone shows you a scatterplot of data that has correlation 0.5, you will say that’s a weak relationship.

Code 
set.seed(7777)
vc <- matrix(c(1, 0, 0, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p1 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, 0.4, 0.4, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p2 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, 0.6, 0.6, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p3 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, 0.8, 0.8, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p4 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, -0.2, -0.2, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p5 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, -0.5, -0.5, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p6 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, -0.7, -0.7, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p7 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
vc <- matrix(c(1, -0.9, -0.9, 1), ncol = 2, byrow = T)
d <- as_tibble(rmvnorm(500, sigma = vc))
p8 <- ggplot(d, aes(x = V1, y = V2)) +
  geom_point() +
  theme_void() +
  theme(
    aspect.ratio = 1,
    plot.background = element_rect(fill = "gray90")
  )
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, ncol = 4)

Robust correlation measures (1/2)

  • Spearman (based on ranks)
    • Sort each variable, and return rank (of actual value)
    • Compute correlation between ranks of each variable
Code 
set.seed(60)
df <- tibble(
  x = c(round(rnorm(5), 1), 3, 10),
  y = c(round(rnorm(5), 1), -3, 10)
) |>
  mutate(xr = rank(x), yr = rank(y))
df
# A tibble: 7 × 4
      x     y    xr    yr
  <dbl> <dbl> <dbl> <dbl>
1   0.7  -1.7     5     2
2   0.5   1.1     4     6
3  -0.6   0.3     2     4
4  -0.2  -0.9     3     3
5  -1.7   0.4     1     5
6   3    -3       6     1
7  10    10       7     7
Code 
cor(df$x, df$y)
[1] 0.79
Code 
cor(df$xr, df$yr)
[1] -0.036
Code 
cor(df$x, df$y, method = "spearman")
[1] -0.036
Code 
ggplot(df, aes(x, y)) +
  geom_point() +
  xlim(c(-3.5, 10.5)) + ylim(c(-3.5, 10.5)) +
  theme(aspect.ratio=1)

Robust correlation measures (2/2)

  • Kendall \(\tau\) (based on comparing pairs of observations)
    • Sort each variable, and return rank (of actual value)
    • For all pairs of observations \((x_i, y_i), (x_j, y_j)\) compare ranks to determine if concordant, relatively \(x_i < x_j, y_i < y_j\) or \(x_i > x_j, y_i > y_j\), or discordant, \(x_i < x_j, y_i > y_j\) or \(x_i > x_j, y_i < y_j\).

\[\tau = \frac{n_c-n_d}{\frac12 n(n-1)}\]

# A tibble: 7 × 4
      x     y    xr    yr
  <dbl> <dbl> <dbl> <dbl>
1   0.7  -1.7     5     2
2   0.5   1.1     4     6
3  -0.6   0.3     2     4
4  -0.2  -0.9     3     3
5  -1.7   0.4     1     5
6   3    -3       6     1
7  10    10       7     7
Code 
cor(df$x, df$y)
[1] 0.79
Code 
cor(df$x, df$y, method = "kendall")
[1] -0.14

Comparison of correlation measures

sample corr spearman kendall
0.52 0.512 0.355
-0.05 -0.087 -0.073
0.30 -0.023 -0.014



Robust calculation corrects outlier problems, but nothing measures the non-linear association.

Transformations

for skewness, heteroskedasticity and linearising relationships, and to emphasize association

Circle of transformations for linearising

Remember the power ladder:

-1, 0, 1/3, 1/2, 1, 2, 3, 4


  1. Look at the shape of the relationship.
  2. Imagine this to be a number plane, and depending on which quadrant the shape falls in, you either transform \(x\) or \(y\), up or down the ladder: +,+ both up; +,- x up, y down; -,- both down; -,+ x down, y up


If there is heteroskedasticity, try transforming \(y\), may or may not help

Activity: Fisherman’s Reach crabs

Your turn

  1. Load the Fisherman’s Reach crabs data, collected by measuring 62 mud crabs.

  2. Make a scatterplot of Wgt by Length.NSW. Describe the relationship. (Which variable should be the dependent variable, on the \(y\) axis?)

  3. Fit a linear model of Wgt on Length.NSW and plot the residuals. Is there evidence of nonlinearity or heteroskedasticity?

  4. Using the circle of transformations, try transforming Wgt (e.g. cube root, square root, log10) to linearise the relationship. Which works best, and why?

  5. Is there possibly a lurking variable, that can explain the nonlinearity?

Code 
library(broom)

fr_crabs <- read_csv("https://ddde.numbat.space/data/fr-crab.csv") |>
  mutate(Sex = factor(Sex, levels = c(1, 2), labels = c("m", "f")))

ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt)) +
  geom_point()

# Fit a linear model, and check the residuals
cr_lm <- lm(Wgt ~ Length.NSW, data = fr_crabs)
fr_crabs <- augment(cr_lm, fr_crabs)
ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
ggplot(fr_crabs, aes(x = Length.NSW, y = .resid)) +
  geom_point()

# Try transformations of Wgt to linearise the relationship
ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt^(1 / 3))) +
  geom_point() +
  ylab("Cube root Wgt")
ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt)) +
  geom_point() +
  scale_y_sqrt() +
  ylab("Square root Wgt")
ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt)) +
  geom_point() +
  scale_y_log10() +
  ylab("Log10 Wgt")

# Check for a lurking variable
ggplot(fr_crabs, aes(x = Length.NSW, y = Wgt, colour = Sex)) +
  geom_point() +
  scale_colour_brewer(palette = "Dark2") +
  theme(legend.position = "bottom")

# Fit the cube root model to males only
fr_crabs_m <- fr_crabs |>
  filter(Sex == "m")
cr_m_lm <- lm(Wgt^(1 / 3) ~ Length.NSW, data = fr_crabs_m)
fr_crabs_m <- augment(cr_m_lm, fr_crabs_m)
coefs <- tidy(cr_m_lm)
ggplot(fr_crabs_m, aes(x = Length.NSW, y = Wgt^(1 / 3))) +
  geom_point() +
  geom_abline(intercept = coefs$estimate[1], slope = coefs$estimate[2])
  • Weight is the dependent variable: it should be on the \(y\) axis, because weight is related to the volume of the crab, which scales with length\(^p\), not the reverse.
  • The raw scatterplot is positive, strong, and a little nonlinear. Fitting a linear model and examining the residuals reveals a U-shape and clear heteroskedasticity, confirming the nonlinearity.
  • Cube root, square root, and log10 transformations of Wgt all reduce the curvature to varying degrees, consistent with a physical model where weight scales roughly with length-cubed.
  • Colouring by Sex reveals a lurking variable: males and females follow somewhat different Length-Weight relationships, and restricting to males and using a cube root transformation gives the best-fitting, most linear relationship — matching the physical intuition that weight is proportional to length\(^3\).

Interference between skewness and association

Code 
baker <- read_csv(here("data/baker.csv"))
p <- ggplot(baker, aes(x = B, y = Corn97BU)) +
  geom_point() +
  xlab("Boron (ppm)") +
  ylab("Corn Yield (bushells)")
ggMarginal(p, type = "density")

Interplay between skewness and association

Data is from a soil chemical analysis of a farm field in Iowa. Is there a relationship between Yield and Boron?


You can get a marginal plot of each variable added to the scatterplot using ggMarginal. This is useful for assessing the skewness in each variable.


Boron is right-skewed Yield is left-skewed. With skewed distributions in marginal variables it is hard to assess the relationship between the two. Make a transformation to fix, first.

Code 
p <- ggplot(
  baker,
  aes(x = B, y = Corn97BU^2)
) + 
  geom_point() +
  xlab("log Boron (ppm)") +
  ylab("Corn Yield^2 (bushells)") +
  scale_x_log10() 
ggMarginal(p, type = "density")

Checking it with visual inference lineup

Code 
set.seed(932)
ggplot(
  lineup(null_permute("Corn97BU"), baker, n = 24),
  aes(x = B, y = Corn97BU)
) +
  geom_point(size=1) +
  facet_wrap(~.sample, ncol = 8)

All but one of the panels have had the association broken by permuting one variable. There is no association in these data sets, and hence plots. Does the data plot stand out as being different from the null (no association) plots?

Code 
set.seed(933)
data(oly12, package = "VGAMdata")
oly12_sub <- oly12 |>
  filter(Sport %in% c(
    "Swimming", "Archery",
    "Hockey", "Tennis"
  )) |>
  filter(Sex == "F") |>
  mutate(Sport = fct_drop(Sport), Sex = fct_drop(Sex))

ggplot(
  lineup(null_permute("Sport"), oly12_sub, n = 24),
  aes(x = Height, y = Weight, colour = Sport)
) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_colour_brewer("", palette = "Dark2") +
  facet_wrap(~.sample, ncol = 8) +
  theme(legend.position = "none")

All but one of the panels have had the association broken by permuting the Sport label. There is no difference in the association between weight and height across sports in these data sets, and hence plots. Does the data plot stand out as being different from the null (no association difference between sports) plots?

Activity: Case study: Soils

Your turn

  1. Using the baker soils data, make a scatterplot (with marginal densities, and a 2D density overlay) of Corn Yield\(^2\) (Corn97BU^2) against Iron (Fe). Is there a relationship? Could there be a lurking variable?

  2. Colour the points by whether Calcium (Ca) is high (> 5200 ppm) or low. How does this change your interpretation of the relationship between Iron and Yield?

Code 
p <- ggplot(
  baker,
  aes(x = Fe, y = Corn97BU^2)
) +
  geom_density2d(colour = "orange") +
  geom_point() +
  xlab("Iron (ppm)") +
  ylab("Corn Yield^2 (bushells)")
ggMarginal(p, type = "density")

ggplot(baker, aes(
  x = Fe, y = Corn97BU^2,
  colour = ifelse(Ca > 5200,
    "high", "low"
  )
)) +
  geom_point() +
  xlab("Iron (ppm)") +
  ylab("Corn Yield^2 (bushells)") +
  scale_colour_brewer("", palette = "Dark2") +
  theme(
    aspect.ratio = 1,
    legend.position = "bottom",
    legend.direction = "horizontal"
  )
  • The raw scatterplot of Yield\(^2\) against Iron appears to have two different patterns — this hints at a lurking variable.
  • Colouring by high vs low Calcium (> 5200 ppm) reveals the structure: if calcium levels in the soil are high, yield is consistently high regardless of iron. If calcium levels are low, there is a positive relationship between yield and iron, with higher iron leading to higher yields.

Scales, benchmarks, calibration matter




This bubble plot here shows total count of COVID-19 incidence (as of Aug 30, 2020) for every county in the USA, inspired by the New York Times coverage.

Code 
load(here("data/nyt_covid.rda"))
usa <- map_data("state")
ggplot() +
  geom_polygon(
    data = usa,
    aes(x = long, y = lat, group = group),
    fill = "grey90", colour = "white"
  ) +
  geom_point(
    data = nyt_county_total,
    aes(x = lon, y = lat, size = cases),
    colour = "red", shape = 1
  ) +
  geom_point(
    data = nyt_county_total,
    aes(x = lon, y = lat, size = cases),
    colour = "red", fill = "red", alpha = 0.1, shape = 16
  ) +
  scale_size("", range = c(1, 30)) +
  theme_map() +
  theme(legend.position = "none")

Where has COVID-19 hit the hardest?

Where are there more people?

This plot tells you NOTHING except where the population centres are in the USA.


To understand relative incidence/risk, report COVID numbers relative the population. For example, number of cases per 100,000 people.

Beyond quantitative variables

When variables are not quantitative

What do you do if the variables are not continuous/quantitative?

Type of variable determines the appropriate mapping.

  • Continuous and categorical: side-by-side boxplots, side-by-side density plots
  • Both categorical: faceted bar charts, stacked bar charts, mosaic plots, double decker plots

Mosaic plots (1/2)

  • Built from a contingency table of two (or more) categorical variables.
  • The plot area is recursively split into tiles:
    • Split horizontally by the first variable, with widths proportional to its marginal counts.
    • Split each column vertically by the second variable, with heights proportional to the conditional counts.
  • Tile area \(\propto\) cell count \(\Rightarrow\) area encodes joint frequency.
  • Column widths show mother’s education level 4 (finished upper secondary) is the most common in this sample; level 1 is rare.
  • Within each column, the split shows the conditional distribution of books-at-home category, given mother’s education.
  • Students whose mothers have more education more often report >100 books at home, but with only 47 students, tiles are unshaded (grey) — the residuals from independence aren’t large enough to be notable at this sample size.

Read off the proportions of one variable relative to the other.

Number of books at home by mother’s education (1 = less than ISCED 1, …, 4 = ISCED 3A), for 2022 PISA students in Turkey (learningtower package).

  • ISCED 1: Elementary/primary school completed
  • ISCED 2: Middle school / junior high completed
  • ISCED 3A: Academic high school diploma
Code 
library(vcd)
library(ggmosaic)
load(here::here("data/student_turkey.rda"))

tur_book <- student_tur |>
  filter(!is.na(book), 
         !is.na(mother_educ),
         !is.na(gender),
         !is.na(computer),
         !is.na(room)) |>
  mutate(book = as.character(book)) |>
  mutate(book = ifelse(book %in% c("0", "1-10"), "<=10", book)) |>
  mutate(
    book = factor(book, levels = c("<=10", "11-25",
      "26-100", "101-200", "201-500", "More than 500"),
      labels = c("<=10", "11-25",
      "26-100", "101-200", "201-500", ">500")),
    mother_educ = factor(mother_educ,
      levels = c("less than ISCED1", "ISCED 1", "ISCED 2", "ISCED 3A"),
      labels = c("1", "2", "3", "4")
    )
  )

ggplot(tur_book) +
  geom_mosaic(aes(x = product(book, mother_educ), 
    fill = book)) +
   scale_fill_discrete("books", palette = "Greens") +
  theme_mosaic()

Mosaic plots (2/2)

The third variable would be interleaved on the horizontal axis.

A fourth variable would be interleaved on the vertical axis.

But it gets increasingly difficult to digest the different proportions.

Code 
ggplot(tur_book) +
  geom_mosaic(aes(x = product(gender, book, 
         mother_educ), 
    fill = gender), 
    ) +
   scale_fill_discrete("Gender", palette = "Greens") +
  theme_mosaic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Code 
ggplot(tur_book) +
  geom_mosaic(aes(x = product(gender, room, book, 
         mother_educ), 
    fill = gender), 
    ) +
   scale_fill_discrete("Gender", palette = "Greens") +
  theme_mosaic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Double decker plots

  • A double decker plot is a mosaic plot where one variable is treated as the response, plotted last.
  • Explanatory (predictor) categorical variables are combined and shown side by side, with widths proportional to their combined counts.
  • Within each combination of explanatory variables, a single bar is split (and coloured) by the response, showing the conditional probability of each response category.

Use when: you have a couple of categorical predictors and one categorical response, and want to compare how the response’s conditional distribution changes across combinations of predictors.

  • Column width encodes how many students are in that gender / mother’s-education combination; relatively few students had mothers who had not finished Year 12.
  • The response is heavily imbalanced: almost all students report having their own room, so most of each bar is orange regardless of group.
  • With only 45 students and just 3 without their own room, the small grey slivers (e.g. 1 male student with a mother who didn’t finish Year 12) reflect a handful of individuals — not a reliable estimate of any group difference.

Compared to a mosaic plot, a double decker plot makes it easier to focus on one response variable at a time, at the cost of not showing area for the response split as clearly.

Whether a student has their own room (grey = no) by gender and mother’s education, for 2022 PISA students in Australia.

Code 
tur_book2 <- student_tur |>
  filter(!is.na(book), 
         !is.na(mother_educ),
         !is.na(gender)) |>
  mutate(
    book = factor(
      case_when(
        book %in% c("0-10", "11-25") ~ "<=25",
        book == "26-100" ~ "26-100",
        book == "101-200" ~ "101-200",
        TRUE ~ ">200"
      ),
      levels = c("<=25", "26-100", "101-200", ">200")
    ),
    mother_educ = factor(mother_educ,
      levels = c("less than ISCED1", "ISCED 1", "ISCED 2", "ISCED 3A"),
      labels = c("1", "2", "3", "4")
    )
  )

tur_dd_tbl <- xtabs(~ gender + mother_educ + room, data = tur_book2)
vcd::doubledecker(tur_dd_tbl,
  gp = grid::gpar(fill = c("grey80", "orangered")),
  labeling_args = list(gp_labels = grid::gpar(fontsize = 10)),
  margins = c(3, 7, 3, 3)
)

Code 
tur_dd_tbl3 <- xtabs(~ mother_educ + gender + room, data = tur_book2)
vcd::doubledecker(tur_dd_tbl3,
  gp = grid::gpar(fill = c("grey80", "orangered")),
  labeling_args = list(gp_labels = grid::gpar(fontsize = 10)),
  margins = c(3, 7, 3, 3)
)

Code 
tur_dd_tbl2 <- xtabs(~ gender + mother_educ + room + computer, data = tur_book2)
vcd::doubledecker(tur_dd_tbl2,
  gp = grid::gpar(fill = c("grey80", "orangered")),
  labeling_args = list(gp_labels = grid::gpar(fontsize = 10)),
  margins = c(3, 7, 3, 3)
)

Code 
tur_dd_tbl4 <- xtabs(~ mother_educ + computer + gender + room, data = tur_book2)
vcd::doubledecker(tur_dd_tbl4,
  gp = grid::gpar(fill = c("grey80", "orangered")),
  labeling_args = list(gp_labels = grid::gpar(fontsize = 10)),
  margins = c(3, 7, 3, 3)
)

Paradoxes

Simpsons paradox

There is an additional variable, which if used for conditioning, changes the association between the variables, you have a paradox.

Simpson’s paradox: famous example (1/2)


Did Berkeley discriminate against female applicants?

Example from Unwin (2015)

Simpson’s paradox: famous example (2/2)

Based on separately examining each department, there is no evidence of discrimination against female applicants.

Example from Unwin (2015)

Always examine the associations in each strata

Handling and imputing missings (1/2)

Check if missings on one variable are related to distribution of the other variable.

Code 
ggplot(oceanbuoys,
       aes(x = air_temp_c,
           y = humidity)) +
     geom_miss_point()

  • Missings plotted in the margins.
  • Missings on humidity only occur for lower values of air tempoerature.

Imputing missings, at least for humidity requires using air temperature values.

But the clustering is due to year

Code 
ggplot(oceanbuoys,
       aes(x = air_temp_c,
           y = humidity)) +
     geom_miss_point() +
     facet_wrap(~year, ncol=2) +
     theme(legend.position = "none")

Handling and imputing missings (2/2)

Use the mean of complete cases to impute the missings

Code 
ocean_imp_yr_mean <- oceanbuoys |>
  select(air_temp_c, humidity, year) |>
  bind_shadow() |>
  group_by(year) |>
  impute_mean_at(vars(air_temp_c, humidity)) |>
  ungroup() |>
  add_label_shadow()
  
ggplot(ocean_imp_yr_mean,
       aes(x = air_temp_c,
           y = humidity,
           colour = any_missing)) + 
  geom_miss_point() +
  scale_color_discrete_divergingx(palette = "Zissou 1") +
  theme(legend.title = element_blank(),
        legend.position = "none")

Use simulation from a bivariate normal distribution, for each year.

Code 
ocean_imp_yr_sim <- nabular(oceanbuoys) |>
  select(air_temp_c, humidity, year) |>
  bind_shadow() |>
  add_label_shadow()

# Need to operate on each subset
ocean_imp_yr_sim_93 <- ocean_imp_yr_sim |>
  filter(year == 1993)
ocean_imp_yr_sim_97 <- ocean_imp_yr_sim |>
  filter(year == 1997)

ocean_imp_yr_sim_93 <- VIM::hotdeck(ocean_imp_yr_sim_93) 
ocean_imp_yr_sim_97 <- VIM::hotdeck(ocean_imp_yr_sim_97) 
  
ocean_imp_yr_sim <- bind_rows(ocean_imp_yr_sim_93, ocean_imp_yr_sim_97)  

ggplot(ocean_imp_yr_sim,
       aes(x = air_temp_c,
           y = humidity,
           colour = any_missing)) + 
  geom_miss_point() +
  scale_color_discrete_divergingx(palette = "Zissou 1") +
  theme(legend.title = element_blank(),
        legend.position = "none")

Networks

Definition

What do you do if you only receive the association information?

For example, in the Harry Potter series of books, some characters show support and comfort to other characters at different times. (This example is from Sam Tyner’s geomnet which was originally taken from The Siena Program.)

Code 
load("../data/hpchars.rda")
glimpse(hp.chars)
Rows: 64
Columns: 4
$ name       <chr> "Adrian Pucey", "Alicia Spinnet", "Ange…
$ schoolyear <int> 1989, 1989, 1989, 1991, 1991, 1989, 198…
$ gender     <fct> M, F, F, M, M, M, M, F, M, M, M, F, M, …
$ house      <fct> Slytherin, Gryffindor, Gryffindor, Rave…
Code 
load("../data/hpedges.rda")
glimpse(hp.edges)
Rows: 434
Columns: 3
$ name1 <chr> "Dean Thomas", "Dean Thomas", "Dean Thomas",…
$ name2 <chr> "Harry James Potter", "Hermione Granger", "N…
$ book  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…



Any data where association between multiple pairs of variables can computed can be examined using networks, and explored using exploratory network analysis.

Step 1: explore the node data

Code 
ggplot(hp.chars, aes(x=schoolyear)) + 
  geom_bar() +
  scale_x_continuous("", breaks = 1986:1995)

What is school year in HP??

Code 
ggplot(hp.chars, aes(x=house)) + 
  geom_bar() +
  xlab("")

Most of the characters are associated with Gryffindor!!

Code 
ggplot(hp.chars, aes(x=gender)) + 
  geom_bar() +
  xlab("")

Is the book lacking in diversity?

Code 
ggplot(hp.chars, aes(x=house, fill=gender)) + 
  geom_bar(position="fill") +
  xlab("") + 
  scale_fill_discrete_divergingx(palette="Zissou 1")

Just Slytherin is lacking diversity.

Step 2: explore the edge data

Code 
p1 <- hp.edges |>
  count(name1) |>
  ggplot(aes(x=fct_reorder(name1, n), y=n)) + 
  geom_col() +
  xlab("") + ylab("") +
  coord_flip()

p2 <- hp.edges |>
  count(name2) |>
  ggplot(aes(x=fct_reorder(name2, n), y=n)) + 
  geom_col() +
  xlab("") + ylab("") +
  coord_flip()
#p1 + p2 + plot_layout(ncol=2)

# Construct data to make plot comparison better
tmp1 <- hp.edges |>
  count(name1) |>
  rename(n1 = n,
         name = name1)
tmp2 <- hp.edges |>
  count(name2) |>
  rename(n2 = n,
         name = name2)
hp.edges_count <- full_join(tmp1, tmp2) |>
  mutate(n1 = if_else(is.na(n1), 0, n1), 
         n2 = if_else(is.na(n2), 0, n2)) |>
  mutate(n = n1 + n2) 
  
p1 <- hp.edges_count |>
  ggplot(aes(x=fct_reorder(name, n), y=n1)) + 
  geom_col() +
  xlab("") + ylab("") + ggtitle("Name 1") +
  coord_flip() +
  theme(aspect.ratio = 1.5)

p2 <- hp.edges_count |>
  ggplot(aes(x=fct_reorder(name, n), y=n2)) + 
  geom_col() +
  xlab("") + ylab("") + ggtitle("Name 2") +
  coord_flip() +
  theme(aspect.ratio = 1.5)

p3 <- hp.edges_count |>
  ggplot(aes(x=fct_reorder(name, n), y=n)) + 
  geom_col() +
  xlab("") + ylab("") + ggtitle("Total") +
  coord_flip() +
  theme(aspect.ratio = 1.5)

p1 + p2 + p3 + plot_layout(ncol=3)

Similar distribution for who is providing support and who is receiving support. Small differences.

Step 3: Examine association

  • Typically only get the connected nodes, so you need to expand the data to include all possible connections, along with the association measure (here 1 or 0).
  • Rows/columns ordered by name. Heatmaps need to have rows and columns ordered in a sensible way.
Code 
library(tidygraph)
library(ggraph)
# Expand edge data
hp.edges_full <- hp.edges |>
  mutate(n = 1) |>
  complete(name1, name2, book, fill = list(n = 0))
# Check
hp.edges_full |> count(name1, name2) -> x
hp.edges_full |> count(name1) -> y
hp.edges_full |> count(name2) -> z
dplyr::full_join(y, z, by=c("name1"="name2")) -> w
# It's missing some items!
name2_miss <- w |>
  filter(is.na(n.y)) |>
  pull(name1)
name1_miss <- w |>
  filter(is.na(n.x)) |>
  pull(name1)
# Ugly manual fix!
n1_d <- hp.edges_full |>
  filter(name2 %in% name1_miss) |>
  rename(name1 = name2,
         name2 = name1) |>
  mutate(n = 0) |>
  select(name1, name2, book, n)
hp.edges_full <- bind_rows(hp.edges_full, n1_d)
n2_d <- hp.edges_full |>
  filter(name1 %in% name2_miss) |>
  rename(name1 = name2,
         name2 = name1) |>
  mutate(n = 0) |>
  select(name1, name2, book, n)
hp.edges_full <- bind_rows(hp.edges_full, n2_d)
  
he_allbooks_long <- hp.edges_full |>
  group_by(name1, name2) |>
  summarise(n = sum(n)) 

he_order <- he_allbooks_long |>
  group_by(name1) |>
  summarise(n = sum(n)) |>
  arrange(desc(n)) 

he_allbooks_long <- he_allbooks_long |>
  mutate(name1 = factor(name1, levels=he_order$name1),
         name2 = factor(name2, levels=he_order$name1))
Code 
p <- ggplot(he_allbooks_long, 
       aes(x=name1,
           y=name2, fill=n)) +
  geom_tile() +
  scale_fill_continuous_sequential(palette="YlGnBu") +
  theme(aspect.ratio=1,
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none")
ggplotly(p, width=600, height=600)

Step 4: Make network summaries

  • Lay out the nodes based on connections over all books
  • Set up the data to allow facetting by book, that shows different connections in different books
Code 
he_connect <- he_allbooks_long |> 
  filter(n > 0) |>
  rename(from = name1, to = name2) 
he_graph <- he_connect |> as_tbl_graph() 

he_g_plt <- ggraph(he_graph, layout = 'tree', circular=TRUE) + 
  geom_edge_link() +
  geom_node_point() 
  
he_nodes <- tibble(x = he_g_plt$data$x,
                   y = he_g_plt$data$y,
                   name = levels(he_allbooks_long$name1))
he_connect <- he_connect |>
  mutate(from_num = as.numeric(from),
         to_num = as.numeric(to))
he_edges <- tibble(
  x = he_g_plt$data$x[he_connect$from_num],
  y = he_g_plt$data$y[he_connect$from_num],
  xend = he_g_plt$data$x[he_connect$to_num],
  yend = he_g_plt$data$y[he_connect$to_num],
  n = he_connect$n
)
  
# Now set up to facet by book, keeping layout set by 
# the overall count
book1 <- hp.edges_full |> 
  filter(book == 1, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book1_nodes <- he_nodes |>
  mutate(book = 1)
book1_edges <- tibble(
  x = he_g_plt$data$x[book1$from_num],
  y = he_g_plt$data$y[book1$from_num],
  xend = he_g_plt$data$x[book1$to_num],
  yend = he_g_plt$data$y[book1$to_num],
  book = 1
)

book2 <- hp.edges_full |> 
  filter(book == 2, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book2_nodes <- he_nodes |>
  mutate(book = 2)
book2_edges <- tibble(
  x = he_g_plt$data$x[book2$from_num],
  y = he_g_plt$data$y[book2$from_num],
  xend = he_g_plt$data$x[book2$to_num],
  yend = he_g_plt$data$y[book2$to_num],
  book = 2
)

book3 <- hp.edges_full |> 
  filter(book == 3, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book3_nodes <- he_nodes |>
  mutate(book = 3)
book3_edges <- tibble(
  x = he_g_plt$data$x[book3$from_num],
  y = he_g_plt$data$y[book3$from_num],
  xend = he_g_plt$data$x[book3$to_num],
  yend = he_g_plt$data$y[book3$to_num],
  book = 3
)

book4 <- hp.edges_full |> 
  filter(book == 4, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book4_nodes <- he_nodes |>
  mutate(book = 4)
book4_edges <- tibble(
  x = he_g_plt$data$x[book4$from_num],
  y = he_g_plt$data$y[book4$from_num],
  xend = he_g_plt$data$x[book4$to_num],
  yend = he_g_plt$data$y[book4$to_num],
  book = 4
)

book5 <- hp.edges_full |> 
  filter(book == 5, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book5_nodes <- he_nodes |>
  mutate(book = 5)
book5_edges <- tibble(
  x = he_g_plt$data$x[book5$from_num],
  y = he_g_plt$data$y[book5$from_num],
  xend = he_g_plt$data$x[book5$to_num],
  yend = he_g_plt$data$y[book5$to_num],
  book = 5
)

book6 <- hp.edges_full |> 
  filter(book == 6, n > 0) |>
  mutate(name1 = factor(name1, levels = he_order$name1),
         name2 = factor(name2, levels = he_order$name1)) |>
  mutate(from_num = as.numeric(name1),
         to_num = as.numeric(name2))

book6_nodes <- he_nodes |>
  mutate(book = 6)
book6_edges <- tibble(
  x = he_g_plt$data$x[book6$from_num],
  y = he_g_plt$data$y[book6$from_num],
  xend = he_g_plt$data$x[book6$to_num],
  yend = he_g_plt$data$y[book6$to_num],
  book = 6
)

all_books_nodes <- bind_rows(book1_nodes, book2_nodes, 
                             book3_nodes, book4_nodes, 
                             book5_nodes, book6_nodes)
all_books_edges <- bind_rows(book1_edges, book2_edges, 
                           book3_edges, book4_edges, 
                           book5_edges, book6_edges)
Code 
p2 <- ggplot() +
  geom_point(data=all_books_nodes, aes(x=x, y=y, label=name)) +
  geom_segment(data=all_books_edges, 
               aes(x=x, y=y, xend=xend, yend=yend),
               linewidth = 0.3) +
  scale_color_viridis_c() +
  facet_wrap(~book, ncol=3) +
  theme(aspect.ratio=1,
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_rect(fill=NA, 
                                        colour="black"),
        legend.position = "none")
ggplotly(p2, tooltip = "label", width=900, height=600)

Key points from today?

  • Scatterplots reveal form, strength, and outliers in the relationship between two continuous variables
  • Correlation is a single number summary — always look at the plot too, and consider robust measures when the relationship is non-linear or has outliers
  • Transformations can linearise a relationship and simplify its description
  • Mosaic and double decker plots extend these ideas to categorical variables
  • Watch for Simpson’s paradox: check associations within subgroups, not just overall
  • Missing values need careful handling — understand why they’re missing before imputing
  • Networks let us examine association between entities connected by relationships

Resources