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

Case study: Olympics

  • Note: Warning message: Removed 1346 rows containing missing values (geom_point)
  • Features:
    • linear relationship (expected, more than?)
    • outliers
    • discretization
  • Substantial overplotting, >10000 athletes.
  • What is interesting? Are there some sport(s) where you would expect specific relationships?

Individual data for the Summer 2012 Olympic Games. Downloaded from http://www.guardian.co.uk/sport/series/london-2012-olympics-data in 2013-03; more recently it has changed to https://www.theguardian.com/sport/series/london-2012-olympics-data.

skimr::skim(oly12)
Data summary
Name oly12
Number of rows 10384
Number of columns 14
_______________________
Column type frequency:
Date 1
factor 6
numeric 7
________________________
Group variables None

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
DOB 6192 0.4 1947-06-01 1997-07-09 1986-09-11 2149

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Name 0 1 FALSE 10366 Lei: 3, Lin: 3, Ale: 2, Hao: 2
Country 0 1 FALSE 205 Gre: 523, Uni: 518, Rus: 414, Aus: 399
Sex 0 1 FALSE 2 M: 5756, F: 4628
PlaceOB 0 1 FALSE 4108 emp: 2690, Seo: 57, Bud: 54, Mos: 50
Sport 0 1 FALSE 42 Ath: 2119, Swi: 907, Foo: 596, Row: 524
Event 0 1 FALSE 763 Men: 336, Wom: 260, Wom: 210, Men: 206

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Age 0 1.00 26.07 5.44 13.0 22.0 25.0 29.0 71.0 ▆▇▁▁▁
Height 561 0.95 1.77 0.11 1.3 1.7 1.8 1.9 2.2 ▁▃▇▃▁
Weight 1280 0.88 72.85 16.07 36.0 61.0 70.0 81.0 218.0 ▇▆▁▁▁
Gold 0 1.00 0.02 0.14 0.0 0.0 0.0 0.0 2.0 ▇▁▁▁▁
Silver 0 1.00 0.02 0.13 0.0 0.0 0.0 0.0 2.0 ▇▁▁▁▁
Bronze 0 1.00 0.02 0.14 0.0 0.0 0.0 0.0 2.0 ▇▁▁▁▁
Total 0 1.00 0.05 0.25 0.0 0.0 0.0 0.0 5.0 ▇▁▁▁▁ 
data(oly12, package = "VGAMdata")
ggplot(oly12, aes(x = Height, y = Weight, label = Sport)) +
  geom_point()

🧩 Try this

Interactivity can be a useful tool for exploring relationships.

Cut and paste the code into your R console, and the resulting plot to examine the sport of the athlete.

library(tidyverse) 
library(plotly) 
data(oly12, package = "VGAMdata") 
p <- ggplot(oly12, aes(x = Height, y = Weight, label = Sport)) + 
  geom_point() 
ggplotly(p)

How many athletes in the different sports?

Categories need re-working:

  • so many different events grouped into athletics
  • cycling split among many categories

Consolidate factor levels

There are several cycling events that are reasonable to combine into one category. Similarly for gymnastics and athletics.


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 the by sport

What is interesting?

  • 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
  • There’s a lot to digest

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
ggplot(oly12, aes(x = Height, y = Weight)) +
  geom_point(alpha = 0.5) + 
  facet_wrap(~Sport, ncol = 8) +
  theme(aspect.ratio = 1)



Note: alpha transparency, and aspect ratio

Remove missings, explore difference by sex

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.
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))

Common ways to augment scatterplots

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

Comparing association



  • Weightlifters are much heavier relative to height
  • Swimmers are leaner relative to height
  • Tennis players are a bit mixed, shorter tend to be heavier, taller tend to be lighter
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"
  )

Comparing spread



  • Modern pentathlon athletes are uniformly height and weight related
  • Shooters are quite varied in body type
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"
  )

Case study: Olympics

We learned that association between height and weight is different strata, defined by categorical variables: sport, gender, and possibly country and age, too.

Some of the association may be due to unmeasured variables, for example, “Athletics” is masking different body types in throwing vs running. This is a lurking variable.



If you were just given the Height and Weight in this data could you have detected the presence of conditional relationships?

It may appear as multimodality.

Can you see conditional dependencies?



There is a barely hint of multimodality.

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.

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)

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}\]

cor(d1$x, d1$y)
[1] 0.52
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)

cor(d2$x, d2$y)
[1] -0.05
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.

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
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
cor(df$x, df$y)
[1] 0.79
cor(df$xr, df$yr)
[1] -0.036
cor(df$x, df$y, method = "spearman")
[1] -0.036
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
cor(df$x, df$y)
[1] 0.79
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

Scatterplot case studies

Case study: Soils (1/4)

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.

Case study: Soils (2/4)



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")

Case study: Soils (3/4)


Lurking variable?



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")

Case study: Soils (4/4)

Colour high calcium (>5200ppm) calcium values

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"
  )

If calcium levels in the soil are high, yield is consistently high. If calcium levels are low, then there is a positive relationship between yield and iron, with higher iron leading to higher yields.

Case study: COVID-19




Bubble plots, size of point is mapped to another variable.

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.

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")

Scales matter





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



Stay tuned!

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

Is what you see really association?

Checking association with visual inference

ggplot(
  lineup(null_permute("Corn97BU"), baker, n = 12),
  aes(x = B, y = Corn97BU)
) +
  geom_point() +
  facet_wrap(~.sample, ncol = 4)

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

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 = 12),
  aes(x = Height, y = Weight, colour = Sport)
) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_colour_brewer("", palette = "Dark2") +
  facet_wrap(~.sample, ncol = 4) +
  theme(legend.position = "none")

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

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?

Resources