ETC5521: Diving Deeply into Data Exploration

Making comparisons between groups and strata

Professor Di Cook

Department of Econometrics and Business Statistics

At the heart of quantitative reasoning is a single question: Compared to what?

-Edward Tufte

Making comparisons

Groups defined by strata labelled in categorical variables
Observations in strata, same or different?
Is there a baseline, or normal value?
What are the dependencies in the way the data was collected?
Are multiple samples recorded for the same individual, or recorded on different individuals?

Determining the appropriate comparison is not always easy, but important to document when decided.

How would you answer these questions?

Are housing prices increasing more in Sydney or Melbourne?

Is the rental price of the unit/apartment too high?

Is the Alfred or Epworth Hospital better for giving birth?

It’s hot and dry today, is the risk of bushfires too high to go hiking?

ED-ie

Outline

Comparing strata
Paired, matched or repeated measurements
Normalising, standardising, baseline, calibration
When there is no data
Including uncertainty
The problem of making thousands of comparisons
Generating comparison samples

Comparing strata

Case study: Melbourne’s daily maximum temperature (1/2)

📊
data
R

Melbourne’s daily maximum temperature from 1970 to 2020.

What are the strata in temporal data?

How are the temperatures different across months?
What about the temperature within a month?

Code

melb_temp <- read_csv(here::here("data", "melb_temp.csv")) |>
  janitor::clean_names() |>
  rename(temp = maximum_temperature_degree_c) |>
  filter(!is.na(temp)) |>
  dplyr::select(year, month, day, temp)
skimr::skim(melb_temp)

── Data Summary ────────────────────────
                           Values   
Name                       melb_temp
Number of rows             18310    
Number of columns          4        
_______________________             
Column type frequency:              
  character                2        
  numeric                  2        
________________________            
Group variables            None     

── Variable type: character ────────────────────────────────
  skim_variable n_missing complete_rate min max empty
1 month                 0             1   2   2     0
2 day                   0             1   2   2     0
  n_unique whitespace
1       12          0
2       31          0

── Variable type: numeric ──────────────────────────────────
  skim_variable n_missing complete_rate   mean    sd     p0
1 year                  0             1 1995.  14.5  1970  
2 temp                  0             1   19.9  6.48    5.7
     p25    p50    p75   p100 hist 
1 1983   1995   2008   2020   ▇▇▇▇▇
2   14.8   18.6   23.6   46.8 ▃▇▃▁▁

Code

ggplot(melb_temp, aes(x=month, y=temp)) +
  geom_violin(draw_quantiles=c(0.25, 0.5, 0.75), fill= "#56B4E9") +
  labs(x = "month", y = "max daily temp (°C)") +
  theme(aspect.ratio=0.5)

Case study: Melbourne’s daily maximum temperature (2/2)

Why can we make the comparison across months?

Because it is the same location, and same years, for each month subset.

Is some variation in temperature each month due to changing climate?

How would you check this?

Code

melb_temp |>
  group_by(year, month) |>
  summarise(temp = mean(temp)) |>
  ggplot(aes(x=year, y=temp)) +
  geom_point(alpha=0.5) +
  geom_smooth(se=F) + 
  facet_wrap(~month, ncol=4, scales="free_y") +
  scale_x_continuous("year", breaks=seq(1970, 2020, 20)) +
  ylab("max daily temp (°C)") +
  theme(aspect.ratio=0.7)

What is scales="free_y" for? ED-ie

Activity: Seasonality over years

Your turn

What does change in seasonality mean?

Using melb_temp, compute the average max daily temperature for each year and month.
Make a line plot of temperature by month (x-axis), with one line connecting the dots for each year. Is there any evidence that the shape of the seasonal cycle has changed over time?
Years differ in their overall average temperature, which can make it hard to compare the shape of the seasonal cycle rather than its level. Try calibrating each year by dividing its monthly temperatures by that year’s average temperature, then remake the plot. Does this change what you see?
What does this calibration achieve, and what does it hide?

Code
Learn

ED-ie

Code

melb_temp_month <- melb_temp |>
  group_by(year, month) |>
  summarise(temp = mean(temp), .groups = "drop")

# raw comparison
ggplot(melb_temp_month, 
    aes(x = month, 
    y = temp, 
    group = year,
    colour = year)) +
  geom_line(alpha = 0.4) +
  geom_point(alpha = 0.4) +
  ylab("max daily temp (°C)") +
  scale_x_discrete("", labels = c("J","F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

# calibrated by each year's average temperature
melb_temp_month |>
  group_by(year) |>
  mutate(temp_cal = temp / mean(temp)) |>
  ungroup() |>
  ggplot(aes(x = month, 
             y = temp_cal, 
             group = year, 
             colour = year)) +
  geom_line(alpha = 0.4) +
  geom_point(alpha = 0.4) +
  ylab("temp relative to yearly average") +   scale_x_discrete("", labels = c("J","F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))

The raw line plot is dominated by year-to-year differences in overall temperature: warmer years sit above cooler years across all months, which makes it hard to compare the shape of the seasonal cycle.
Dividing by each year’s average removes the level differences between years, so the lines overlap more and it becomes easier to compare the timing and shape of the seasonal peak and trough across years.
This calibration also removes information about whether some years were warmer overall than others — it isn’t appropriate if the question of interest is about overall warming rather than seasonal shape.

The olive oil data consists of the percentage composition of 8 fatty acids (palmitic, palmitoleic, stearic, oleic, linoleic, linolenic, arachidic, eicosenoic) found in the lipid fraction of 572 Italian olive oils.
There are 9 collection areas, 4 from southern Italy (North and South Apulia, Calabria, Sicily), two from Sardinia (Inland and Coastal) and 3 from northern Italy (Umbria, East and West Liguria).

Code

data(olives, package = "classifly")
df2 <- olives |>
  mutate(Region = factor(Region, labels = c("South", "Sardinia", "North")))

skimr::skim(df2)

── Data Summary ────────────────────────
                           Values
Name                       df2   
Number of rows             572   
Number of columns          10    
_______________________          
Column type frequency:           
  factor                   2     
  numeric                  8     
________________________         
Group variables            None  

── Variable type: factor ───────────────────────────────────
  skim_variable n_missing complete_rate ordered n_unique
1 Region                0             1 FALSE          3
2 Area                  0             1 FALSE          9
  top_counts                         
1 Sou: 323, Nor: 151, Sar: 98        
2 Sou: 206, Inl: 65, Cal: 56, Umb: 51

── Variable type: numeric ──────────────────────────────────
  skim_variable n_missing complete_rate   mean    sd   p0
1 palmitic              0             1 1232.  169.   610
2 palmitoleic           0             1  126.   52.5   15
3 stearic               0             1  229.   36.7  152
4 oleic                 0             1 7312.  406.  6300
5 linoleic              0             1  981.  243.   448
6 linolenic             0             1   31.9  13.0    0
7 arachidic             0             1   58.1  22.0    0
8 eicosenoic            0             1   16.3  14.1    1
     p25   p50    p75 p100 hist 
1 1095   1201  1360   1753 ▁▂▇▆▁
2   87.8  110   169.   280 ▂▇▅▃▁
3  205    223   249    375 ▂▇▃▁▁
4 7000   7302. 7680   8410 ▁▇▇▇▁
5  771.  1030  1181.  1470 ▃▅▃▇▃
6   26     33    40.2   74 ▂▅▇▂▁
7   50     61    70    105 ▂▁▇▇▂
8    2     17    28     58 ▇▃▅▂▁

Code

g1 <-
  df2 |>
  mutate(Area = fct_reorder(Area, palmitic)) |>
  ggplot(aes(Area, palmitic, color = Region)) +
  geom_boxplot() +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE, x = guide_axis(n.dodge = 2)) +
  theme(aspect.ratio=0.5)

g2 <- ggplot(df2, aes(Region, palmitic, color = Region)) +
  geom_boxplot() +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE) +
  theme(axis.text = element_blank())

g3 <- ggplot(df2, aes(palmitic, color = Region)) +
  geom_density() +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE) +
  theme(axis.text = element_blank())

g4 <- ggplot(df2, aes(palmitic, color = Region)) +
  stat_ecdf() +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE) +
  theme(axis.text = element_blank())

g5 <- g2 + g3 + g4 + plot_layout(ncol=3)
  
g1 + g5 + plot_layout(ncol=1, heights=c(2,1),
  guides = "collect")

Case study: olive oils (2/4)

📊
R

Colour is generally good to differentiate strata but if there are too many categories then it becomes hard to compare.

Code

ggplot(olives, aes(palmitoleic, palmitic, color = Area)) +
  geom_point() +
  scale_color_discrete_divergingx(palette="Zissou 1")

Case study: olive oils (3/4)

📊
R

It can be hard to compare across plots, because we need to remember what the previous pattern was when focusing on the new cell.

Code

ggplot(olives, aes(palmitoleic, palmitic, color = Area)) +
  geom_point() +
  facet_wrap(~Area) +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE)

Case study: olive oils (4/4)

📊
R

Comparison to all, by putting a shadow of all the data underneath the subset in each cell.

Code

ggplot(olives, aes(palmitoleic, palmitic)) +
  geom_point(data = dplyr::select(olives, -Area), color = "gray") +
  geom_point(aes(color = Area), size=2) +
  facet_wrap(~Area) +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  guides(color = FALSE)

Strata from quantitative variable

The coplot divides the numerical variable into chunks, and facets by these. The chunks traditionally we overlapping.

Becker, Cleveland and Shyu, (1996); Cleveland (1993)

Code

# Sizing of figure is difficult, so save it
library(ggcleveland)
olives_sard <- df2 |>
  filter(Region == "Sardinia")
p <- gg_coplot(olives_sard,
  x=arachidic, y=oleic, 
  faceting = linoleic,
  number_bins = 6, 
  overlap = 1/4) +
  theme(aspect.ratio=0.5)

Strata for categorical variables

Code

set.seed(42)

# Simulate a small employee-engagement survey
# 5 questions, 5-point Likert scale, 200 respondents

levels_5pt <- c(
  "Strongly disagree", "Disagree", "Neutral", "Agree", "Strongly agree"
)

survey_questions <- c(
  "My manager gives useful feedback",
  "I have the tools I need to do my job",
  "I understand how my work connects to company goals",
  "I feel comfortable raising concerns",
  "I would recommend this company as a place to work"
)

# Give each question a different underlying 
# sentiment so the charts actually show 
# a mix of positive/negative/mixed results
probs <- list(
  c(0.05, 0.10, 0.15, 0.40, 0.30),  # mostly positive
  c(0.20, 0.30, 0.20, 0.20, 0.10),  # skews negative
  c(0.05, 0.15, 0.50, 0.20, 0.10),  # mostly neutral
  c(0.15, 0.20, 0.10, 0.30, 0.25),  # mixed / bimodal-ish
  c(0.10, 0.15, 0.20, 0.30, 0.25)   # mildly positive
)

n_resp <- 200

survey_wide <- as_tibble(
  setNames(
    lapply(seq_along(survey_questions), function(i) {
      factor(
        sample(levels_5pt, n_resp, replace = TRUE, prob = probs[[i]]),
        levels = levels_5pt,
        ordered = TRUE
      )
    }),
    as.character(1:5)
  )
)

survey_long <- survey_wide |>
  pivot_longer(everything(), 
    names_to = "question_num", 
    values_to = "response") |>
  mutate(question = survey_questions[as.numeric(question_num)]) |>
  mutate(response = factor(response, 
    levels = levels_5pt, ordered = TRUE),
    question = factor(question, 
      levels = survey_questions))

stacked_data <- survey_long |>
  count(question_num, response) |>
  group_by(question_num) |>
  mutate(pct = n / sum(n)) |>
  ungroup()

Code

plain_stacked_bar <- ggplot(stacked_data, 
    aes(x = pct, 
        y = question_num, 
        fill = response)) +
  geom_col(position = position_stack(reverse = TRUE)) +
  scale_fill_brewer(palette = "RdBu", direction = 1) +
  scale_x_continuous(labels = scales::percent) +
  scale_y_discrete(limits = rev) +
  labs(
    x = "Percent of respondents", y = NULL, fill = NULL
  ) +
  theme_minimal()
plain_stacked_bar

The questions are: 1 My manager gives useful feedback, 2 I have the tools I need to do my job, 3 I understand how my work connects to company goals, 4 I feel comfortable raising concerns, 5 I would recommend this company as a place to work

Rank the 5 questions by “most disagreement”.

ED-ie

A Likert plot centres the bar on the neutral category.

gglikert
Sorted

Code

library(ggstats)
likert_plot <- gglikert(
  survey_wide,
  add_labels = FALSE,
  add_totals = FALSE
) +
  scale_fill_brewer(palette = "RdBu", direction = 1) 
likert_plot

NOTE: Also a good idea to sort questions (not done here).

Code

library(ggstats)
likert_plot <- gglikert(
  survey_wide,
  add_labels = FALSE,
  add_totals = FALSE,
  sort = "ascending"
) +
  scale_fill_brewer(palette = "RdBu", direction = 1) 
likert_plot

Simple rearranging comparisons: tuberculosis (1/3)

Code

tb_oz |>
  ggplot(aes(x=year, y=count, fill=sex)) +
    geom_col(position="fill") +
    scale_fill_discrete_divergingx(palette = "Zissou 1") +
    facet_wrap(~age, ncol=6) +
    xlab("") + ylab("proportion")

Primary comparison is sex, relative to yearly trend.

Simple rearranging comparisons: tuberculosis (2/3)

Code

tb_oz |>
  ggplot(aes(x=year, y=count, fill=age)) +
    geom_col(position="fill") +
    scale_fill_discrete_divergingx(palette = "Zissou 1") +
    facet_wrap(~sex, ncol=2) +
    xlab("") + ylab("proportion")

Primary comparison is age, relative to yearly trend.

Simple rearranging comparisons: tuberculosis (3/3)

Code

tb_oz |>
  ggplot(aes(x=year, y=count, fill=age)) +
    geom_col() +
    scale_fill_discrete_divergingx(palette = "Zissou 1") +
    facet_grid(sex~age, scales="free") +
    xlab("") + ylab("count") +
    theme(legend.position = "none")

Primary comparison is year trend, separately for age and sex.

Paired, matched or repeated measurements

Pairing adjusts for individual differences

If we were wanting to measure the effect of incorporating a data analytics competition on student learning which is the best design?

METHOD A

Divide students into two groups. Make sure that each group has similar types of students, so both groups are as similar as possible.
One group gets an extra traditional assignment, and the other group participates in a data competition.
Each student takes an exam on the content being taught.
The scores are compared using side-by-side boxplots and a two-sample permutation test

METHOD B

Each student takes an exam on the content being taught. We’ll call this their BEFORE score.
Divide students into two groups. Make sure that each group has similar types of students, so both groups are as similar as possible.
One group gets an extra traditional assignment, and the other group participates in a data competition.
Each student takes an exam on the content being taught. We’ll call this their AFTER score.
The difference between the before and after scores are compared using side-by-side boxplots and a two-sample permutation test

What other modifications to the design can you think of?

Case study: choropleth vs hexagon tile (1/3)

The goal is to demonstrate that the hexagon tile map is better than the choropleth for communicating disease incidence across Australia.

The choropleth fills geographic regions (LGAs, SA2s, …) with colour corresponding to the thyroid cancer relative difference from the overall mean. The hexagons, are also filled this way.

Kobakian et al

Case study: choropleth vs hexagon tile (2/3)

Each participant can only see the (same) data once.
Need to test for different types of spatial patterns.
Need to repeat measure each type of pattern, and each participant.

Pairing is done on the data set. Four different data sets used for each pattern.

Trial 1

Participant 1

Participant 2

Trial 2

Participant 1

Participant 2

Case study: choropleth vs hexagon tile (3/3)

Ignore the pairing

Code

hstudy <- read_csv("https://raw.githubusercontent.com/srkobakian/experiment/master/data/DAT_HexmapPilotData_V1_20191115.csv")
hstudy |>
  filter(trend == "three cities") |>
  ggplot(aes(x=detect)) + geom_bar() + facet_wrap(~type, ncol=2)

Looks like detection rate about 50-50 for hexagon tile map, which is better than almost zero for choropleth map.

Account for the pairing

Code

hstudy |>
  filter(trend == "three cities") |>
  select(type, replicate, detect) |>
  group_by(type, replicate) |>
  summarise(pdetect = length(detect[detect == 1])/length(detect)) |>
  ggplot(aes(x=type, y=pdetect)) +
    geom_point() +
    geom_line(aes(group=replicate)) +
    ylim(c(0,1)) +
    xlab("") +
    ylab("Proportion detected")

For each data set, the hexagon tile map performed better.

Normalising, standardising, baseline, calibration

Famous example: trade

📊
data
R

The export from England to the East Indies and the import to England from the East Indies in millions of pounds (A).
Import and export figures are easier to compare by plotting the difference like in (B).
Relative difference may be more of an interest: (C) plots the relative difference with respect to the average of export and import values.
The red area correspond to War of the Spanish Succession (1701-14), Seven Years’ War (1756-63) and the American Revolutionary War (1775-83).

Code

data(EastIndiesTrade, package = "GDAdata")
skimr::skim(EastIndiesTrade)

── Data Summary ────────────────────────
                           Values         
Name                       EastIndiesTrade
Number of rows             81             
Number of columns          3              
_______________________                   
Column type frequency:                    
  numeric                  3              
________________________                  
Group variables            None           

── Variable type: numeric ──────────────────────────────────
  skim_variable n_missing complete_rate  mean    sd   p0
1 Year                  0             1 1740   23.5 1700
2 Exports               0             1  518. 421.   100
3 Imports               0             1 1005. 320.   460
   p25  p50  p75 p100 hist 
1 1720 1740 1760 1780 ▇▇▇▇▇
2  145  370  840 1395 ▇▂▃▂▂
3  835  975 1000 1550 ▃▃▇▁▅

Code

g1 <- ggplot(EastIndiesTrade, aes(Year, Exports)) +
  annotate("rect",
    xmin = 1701, xmax = 1714,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1756, xmax = 1763,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1775, xmax = 1780,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  geom_line(color = "#339933", size = 2) +
  geom_line(aes(Year, Imports), color = "red", size = 2) +
  geom_ribbon(aes(ymin = Exports, ymax = Imports), fill = "gray") +
  labs(y = "<span style='color:#339933'>Export</span>/<span style='color:red'>Import</span>", tag = "(A)") +
  theme(aspect.ratio=0.7, axis.title.y = ggtext::element_markdown())

g2 <- ggplot(EastIndiesTrade, aes(Year, Imports - Exports)) +
  annotate("rect",
    xmin = 1701, xmax = 1714,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1756, xmax = 1763,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1775, xmax = 1780,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  geom_line(size = 2) +
  labs(tag = "(B)") +
  theme(aspect.ratio=0.7)

g3 <- ggplot(EastIndiesTrade, aes(Year, (Imports - Exports) / (Exports + Imports) * 2)) +
  annotate("rect",
    xmin = 1701, xmax = 1714,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1756, xmax = 1763,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  annotate("rect",
    xmin = 1775, xmax = 1780,
    ymin = -Inf, ymax = Inf,
    fill = "red", alpha = 0.3
  ) +
  geom_line(color = "#001a66", size = 2) +
  labs(y = "Relative difference", tag = "(C)") +
  theme(aspect.ratio=0.7)

g1 + g1 + g2 + g3 + plot_layout(ncol=2)

Am I short?

I am 165 cms tall.

Code

ggplot(df22, aes(x=height)) +
  geom_histogram(breaks = seq(132.5, 205, 5),
    colour="white") +
  geom_vline(xintercept = 165, colour="#D55E00",
    linewidth=2)

But, there are strata in humans, so compared to what would be better?

Code

ggplot(df22, aes(x=height)) +
  geom_histogram(breaks = seq(137.5, 205, 5),
    colour="white") +
  geom_vline(xintercept = 165, colour="#D55E00",
    linewidth=2) +
  facet_wrap(~sex, ncol=3, scales="free_y")

Nope, I’m average height.

Standardising

Within each strata convert values to a z-score.

\[ z = \frac{x-\bar{x}}{s} \]

Code

df22 <- df22 |>
  group_by(sex) |>
  mutate(zscore = (height -
    mean(height))/sd(height))

females: $\bar{x}=$ 164.43, $s=$ 6.41
males: $\bar{x}=$ 174.23, $s=$ 8.71
unknown: $\bar{x}=$ 165.74, $s=$ 18.15

My z-score is 0.09.

Rob’s height is 170 cms. His z-score is -0.49.

I am relatively TALLER than Rob.

Range of standardising methods

Put values on the same footing so different subgroups can be compared

Method	Formula	Use when
Z-score	$(x-\mu)/\sigma$	Roughly normal, no big outliers
Min-max	$(x-x_{min})/(x_{max}-x_{min})$	Need a bounded [0,1] range
Robust scaling	$(x-\text{median})/\text{IQR}$	Outliers / skew present
Unit vector norm	$x/\lVert x \rVert$	Care about direction, not magnitude
Rate / per-unit	$x/\text{exposure}$	Comparing groups of different size
Index / rebase	$100 \times x/x_{base}$	Comparing trends across units
Percentile rank	rank$(x)$	Compare relative standing

What happens when denominator is small with some of these metrics? ED-ie

Common ways to standardise observations

Put rows on the same footing so different individuals/units can be compared

Method	Formula	Use when
Ipsatization	$x_{ij}-\bar{x}_i$	Respondents differ in baseline scale use
Row z-score	$(x_{ij}-\bar{x}_i)/s_i$	Respondents differ in spread too
Compositional %	$x_{ij}/\sum_j x_{ij}$	Compare allocation pattern, not scale
Library-size norm.	adjust for exposure/depth	Unequal measurement effort per unit
Unit norm (row)	$x/\lVert x \rVert$	Compare pattern, not magnitude (e.g. text)
Rank-within-row	rank across own items	Avoid scale-use differences entirely

Comparing with a baseline

Code

data(anorexia, package="MASS")
ggplot(data=anorexia, 
 aes(x=Prewt, y=Postwt, 
    colour=Treat)) + 
 coord_equal() +
 xlim(c(70, 110)) + ylim(c(70, 110)) +
 xlab("Pre-treatment weight (lbs)") +  
 ylab("Post-treatment weight (lbs)") +
 geom_abline(intercept=0, slope=1,  
   colour="grey80", linewidth=1.25) + 
 geom_density2d() + 
 geom_point(size=3) +
 facet_grid(.~Treat) +
 theme(legend.position = "none")

Primary comparison is before treatment weight, and after treatment weight.
Three different treatments.

Unwin, Hofmann and Cook (2013)

Code

ggplot(data=anorexia, 
  aes(x=Prewt, colour=Treat,
    y=(Postwt-Prewt)/Prewt*100)) + 
  xlab("Pre-treatment weight (lbs)") +  
  ylab("Percent increase in weight") +
  geom_hline(yintercept=0, linewidth=1.25, 
    colour="grey80") + 
  geom_point(size=3) +   
  facet_grid(.~Treat) +
 theme(legend.position = "none")

Compute the difference
Compare difference relative to before weight
Before weight is used as the baseline

Famous calibration error: ovarian cancer proteomics (2002–2005)

A 2002 study claimed a mass-spectrometry blood test could detect ovarian cancer with near-perfect accuracy — using mass spectroscopy to display proteins in serum as a series of peaks, with a computer algorithm finding patterns unique to patients with the disease. It generated enormous excitement (a commercial test, “OvaCheck,” was being built on it).

Keith Baggerly and Kevin Coombes at MD Anderson reanalyzed the underlying data and could not reproduce the reported sensitivity and specificity, tracing the problem not to the biology but to how the data had been collected and processed.

See New York Times article.

Baggerly called this style of work forensic bioinformatics.

Discriminating peaks were spread across the entire spectrum, when true cancer-related protein changes should only affect a few specific peaks — a sign the “signal” wasn’t biological.
Some of the “important” peaks fell in regions of the spectrum where the machine couldn’t have been reliably sampling proteins at all — meaning the signal was probably just instrument noise.
Crucially, Baggerly pointed to incomplete randomization of samples as a likely source of these problems — cancer and control samples tended to be run in different batches, on different days, so any drift in the machine’s calibration between runs became perfectly confounded with disease status. The “biomarker” was picking up which day the sample went through the machine, not cancer.

When there is no data

Comparison when only have positives (1/2)

Code

library(ecotourism)
data(orchids)        # 302,123 rows: one row PER SIGHTING (presence-only)
data(weather)        # daily weather for every day, 2014-2024, per station
data(top_stations)

orchid_stations <- top_stations |>
  filter(organism == "orchids") |>
  pull(ws_id)

# Sightings by day and station
orchid_obs <- orchids |>
  filter(year == 2024) |>
  filter(ws_id %in% orchid_stations) |>
  summarise(count = n(), .by = c(date, ws_id))

# Weather by day and station
weather_2024 <- weather |>
  filter(year == 2024) |>
  filter(ws_id %in% orchid_stations)

Code

orchid_date <- orchid_obs |>
  left_join(weather_2024, by = c("ws_id", "date")) |>
  filter(!is.na(prcp)) |>
  mutate(zero_prcp = if_else(prcp < 0.001, "dry", "rain")) 

ggplot(orchid_date, aes(x=zero_prcp, y=count+0.1)) + 
  geom_lv(aes(fill = after_stat(LV))) +
  scale_fill_lv() +
  scale_y_log10() +
  xlab("")

Are there less sightings of orchids on rainy days?

Make sure all days are included.

Code

weather_2024_tsb <- as_tsibble(weather_2024, index = date, key = ws_id)

weather_2024_tsb |> has_gaps()

# A tibble: 3 × 2
  ws_id        .gaps
  <chr>        <lgl>
1 946300-99999 FALSE
2 956410-99999 FALSE
3 956470-99999 FALSE

Code

orchid_obs_tsb <- as_tsibble(orchid_obs, index = date, key = ws_id)

orchid_obs_tsb |> has_gaps()

# A tibble: 3 × 2
  ws_id        .gaps
  <chr>        <lgl>
1 946300-99999 TRUE 
2 956410-99999 TRUE 
3 956470-99999 TRUE

Code

ggplot(orchid_obs_tsb, aes(x=ws_id, y=date)) +
  geom_point(shape = "|") +
  xlab("") + ylab("") +
  coord_flip() +
  theme(aspect.ratio = 0.5)

Comparison when only have positives (2/2)

Join using the full time data as the primary set.

Code

weather_orchid <- weather_2024 |>
  left_join(orchid_obs, by = c("ws_id", "date")) |>
  filter(!is.na(prcp)) |>
  mutate(zero_prcp = if_else(prcp < 0.001, "dry", "rain")) |>
  mutate(count = replace_na(count, 0))

ggplot(weather_orchid, aes(x=zero_prcp, y=count+0.1)) + 
  geom_lv(aes(fill = after_stat(LV))) +
  scale_fill_lv() +
  scale_y_log10() +
  xlab("")

Orchids are a seasonal sighting, and there might be a difference in precipitation between seasons. Solution: Restrict to a month where orchids are typically visible.

Code

weather_orchid |>
  filter(month == 10) |>
  ggplot(aes(x=zero_prcp, y=count+0.1)) + 
  geom_lv(aes(fill = after_stat(LV))) +
  scale_fill_lv() +
  scale_y_log10() +
  xlab("")

Including uncertainty in plots

The usual approach: Error bars bolted on

Code

set.seed(1134)
county_name_sample <-
  toy_temp |> 
  select(county_name) |>
  distinct() |>
  sample_frac(0.5) |>
  pull(county_name)

toy_temp_eg <- toy_temp |>
  filter(county_name %in% county_name_sample) |>
  group_by(county_name) |>
  summarise(
    mean_temp = mean(recorded_temp),
    se_temp   = sd(recorded_temp) / sqrt(n())
  )

ggplot(toy_temp_eg, aes(x = fct_reorder(county_name, mean_temp), y = mean_temp)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean_temp - 2*se_temp,
                     ymax = mean_temp + 2*se_temp), width = 0.2) +
  labs(y = "Mean recorded temp (°C)", x = NULL) +
  coord_flip() +
  theme_minimal()

Uncertainty is a separate layer added after the fact — easy to skip, easy to ignore.

ggdibbler: Uncertainty is the data

Code

library(ggdibbler)
library(distributional)

# Same summary, but the estimate becomes a DISTRIBUTION, not a number
toy_temp_dist <- toy_temp |>
  filter(county_name %in% county_name_sample) |>
  group_by(county_name) |>
  summarise(
    temp_dist = 
      dist_normal(mu = mean(recorded_temp),
                             sigma = sd(recorded_temp) / sqrt(n()))
  ) |>
  mutate(county_name = fct_reorder(county_name, temp_dist, .fun = mean))

ggplot(toy_temp_dist, 
    aes(x = county_name, 
    y = temp_dist)) +
  geom_point_sample(times = 50, alpha = 0.4) +
  labs(y = "Mean recorded temp (°C)", x = NULL) +
  coord_flip() +
  theme_minimal()

Swap the number for a dist_normal(), use the _sample geom, and ggdibbler draws times plausible realisations — uncertainty shows up as visual spread/noise, the same plot you’d already make, just fed a random variable instead of a point estimate.

The problem of making thousands of comparisons

Case study: gene expression (1/2)

A typical experiment tests every gene for a difference between conditions: often 20,000+ tests at once, each with its own $p$-value.

At $\alpha = 0.05$, 5% of truly null genes will look “significant” by chance alone
With 20,000 genes and (say) 18,000 truly unaffected by the treatment:

$18{,}000 \times 0.05 \approx 900$ false positives, even if nothing real is happening

A single-test error rate ($\alpha$) controls the wrong thing here: it bounds the false-positive rate per test, not across the whole gene list you’ll actually report

The fix: control the False Discovery Rate, not $\alpha$ used in each test

Benjamini–Hochberg procedure: rank $p$-values, adjust each against its rank, relative to number of tests, keep everything below the corresponding threshold.

Case study: gene expression (1/2)

Gene expression studies rarely have more than a handful of replicates per group (often $n = 3$–$5$). A t-statistic depends on the variance estimate:

\[t = \frac{\bar{x}_1 - \bar{x}_2}{s_{\text{pooled}} \sqrt{2/n}}\]

By chance, some genes get an artificially small $s^2$ → an inflated $t$ → a “highly significant” result that’s really just an unlucky variance estimate, not a real effect
This is the same fragility problem as small-sample A/B tests — but multiplied across thousands of genes simultaneously

(Also can think of this as some genes have more impact than others with smaller expression.)

The fix: Empirical Bayes shrinkage of variance estimates, or moderated t-statistic

Each gene has too little data on its own, but thousands of other genes are measured in the same experiment
Borrow strength across genes: shrink each gene’s variance estimate toward the typical variance seen across all genes, weighted by how reliable that gene’s own estimate is.

The common thread: both are about not letting one noisy comparison speak for itself — FDR corrects for volume of tests, empirical Bayes corrects for scarcity of data per test.

Generating comparison samples

Bootstrap confidence intervals

📊
data
R

Confidence intervals show what might happen to estimates with different samples,
with the same dependence structure.
Sample the current sample, but don’t change anything else.
Reason for sampling with replacement, is to keep sample size the same - we know that variance decreases with smaller sample size.

For choropleth vs hexagon tiles, sample participants with replacement.

Code

hstudy |> filter(trend == "three cities") |> count(type, replicate)

# A tibble: 8 × 3
  type      replicate     n
  <chr>         <dbl> <int>
1 Geography         9    10
2 Geography        10    10
3 Geography        11    11
4 Geography        12    11
5 Hexagons          9    11
6 Hexagons         10    11
7 Hexagons         11    10
8 Hexagons         12    10

Code

hstudy_sub <- hstudy |>
  filter(trend == "three cities") |>
  select(id, type, replicate, detect) 

# Function to compute proportions
prop_func <- function(df) {
  df_smry <- df |>
    group_by(type, replicate) |>
    summarise(pdetect = length(detect[detect == 1])/length(detect)) |> 
    ungroup() |>
    pivot_wider(names_from = c(type, replicate),
               values_from = pdetect)
  df_smry
}

nboots <- 100
set.seed(1023)
bsamps <- tibble(samp="0", prop_func(hstudy_sub))
for (i in 1:nboots) {
  samp_id <- sort(sample(unique(hstudy_sub$id),
    replace=TRUE))
  hs_b <- NULL
  for (j in samp_id) {
    x <- hstudy_sub |>
      filter(id == j)
    hs_b <- bind_rows(hs_b, x)
  }
  bsamps <- bind_rows(bsamps,
    tibble(samp=as.character(i), prop_func(hs_b)))
}

bsamps_long <- bsamps |>
  pivot_longer(Geography_9:Hexagons_12, 
    names_to = "treatments", 
    values_to = "pdetect") |> 
  separate(treatments, into=c("type", "replicate"))

ggplot() +
    geom_line(data=filter(bsamps_long, samp != "0"),
      aes(x=type, 
          y=pdetect, 
          group=samp),
     linewidth=0.5, alpha=0.6, colour="grey70") +
    geom_line(data=filter(bsamps_long, 
                     samp == "0"),
      aes(x=type, 
          y=pdetect, 
          group=samp),
     linewidth=2) +
    facet_wrap(~replicate) +
    ylim(c(0,1)) +
    xlab("") +
    ylab("Proportion detected")

Lineups

Lineups show what might happen to estimates with null samples, where (by construction) there is no relationship.
Thus you need to break dependence structure.

For choropleth vs hexagon tiles, randomise the type of plot each participant received. This breaks any dependence between type and detection rate.

Code

n_nulls <- 11
set.seed(1110)
lsamps <- tibble(samp="0", prop_func(hstudy_sub))
for (i in 1:n_nulls) {

  hs_b <- hstudy_sub |>
    group_by(id) |>
    mutate(type = sample(type)) |>
    ungroup()
  lsamps <- bind_rows(lsamps,
    tibble(samp=as.character(i), prop_func(hs_b)))
}

lsamps_long <- lsamps |>
  pivot_longer(Geography_9:Hexagons_12, 
    names_to = "treatments", 
    values_to = "pdetect") |> 
  separate(treatments, into=c("type", "replicate"))

lsamps_long |> ggplot() +
    geom_line(aes(x=type, 
          y=pdetect, 
          group=replicate)) +
    facet_wrap(~samp, ncol=4) +
    ylim(c(0,1)) +
    xlab("") +
    ylab("Proportion detected")

Activity: Comparisons gone wrong, for real

On December 5, 2025, the University of Nebraska Board of Regents voted 7-1 to eliminate the Department of Statistics at UNL entirely, along with three other departments (Earth and Atmospheric Sciences, Educational Administration, and Textiles, Merchandising and Fashion Design), as part of a $6.74 million cut affecting 51.5 positions — Statistics alone accounted for $1.75 million and 12 eliminated positions. All degree programs (BS, MS, PhD) were closed, and tenured and tenure-track faculty were let go, with their last day of employment set for May 14, 2027.

The decision leaned heavily on a set of cross-department performance metrics — exactly the kind of comparison problems from today’s lecture.

Your turn

Look at the 7 slides under “Comparisons” in the UNL Statistics faculty’s rebuttal deck.

Can you see what comparison calculation was wrong?

?
Learn

ED-ie

Bundling every department’s SRI into one distribution ignores that fields differ enormously in citation norms, grant sizes, and publishing rates — comparable to the “fruit salad” mixture-distribution problem: a mixture is fine for describing the whole university, but not for ranking individual departments against each other.
A department’s position in a “shift” plot depends on which comparison group and time window it’s judged against; without a consistent, disclosed baseline for every department, claims of “equal treatment” can’t be verified — this is the same transparency issue as the OvaCheck calibration failure.
Taking a difference relative to a reference group (e.g. public AAU institutions) is exactly the baseline comparison idea from today: it controls for field-wide norms so the number reflects relative standing, not just which discipline you’re in.
Ranking departments on raw, un-standardised metrics compares apples to pumpkins — the fix is the same one we used all lecture: identify the correct strata, standardise within them, and always ask “compared to what?” before trusting a ranking.

Key points

Compared to what is especially important for exploring observational data — the strata, dependencies, and baseline all need to be identified before comparing.
Adjust for individual variation, by
- pairing (or repeated measures on the same unit)
- comparing relative to a baseline
- standardising onto a common scale (z-score, rate, index, …)
Calibration matters: comparisons across batches, instruments, or samples can produce spurious “signals” if not properly randomised or normalised — it’s the process being compared, not the effect of interest.

When data only records presence (not absence), join against the complete set of possibilities so missing combinations count as zero, not nothing.
Uncertainty should be part of the estimate, not an afterthought bolted on with error bars.
Making many comparisons at once inflates the chance of false discoveries — correct for multiplicity.
Bootstrap samples and permutation lineups let you simulate what the comparison would look like under resampling or under no real effect, to judge whether an observed difference is real.

Resources

Wilke (2019) Fundamentals of Data Visualization Chapters 9, 10, 11
Unwin (2015) Graphical Data Analysis with R Chapter 10
UNL Statistics faculty (2025) The Metrics: The Good, the Bad, and the Ugly — rebuttal seminar on the metrics used to justify the elimination of the UNL Department of Statistics

Method	Formula	Use when
Z-score	\((x-\mu)/\sigma\)	Roughly normal, no big outliers
Min-max	\((x-x_{min})/(x_{max}-x_{min})\)	Need a bounded [0,1] range
Robust scaling	\((x-\text{median})/\text{IQR}\)	Outliers / skew present
Unit vector norm	\(x/\lVert x \rVert\)	Care about direction, not magnitude
Rate / per-unit	\(x/\text{exposure}\)	Comparing groups of different size
Index / rebase	\(100 \times x/x_{base}\)	Comparing trends across units
Percentile rank	rank\((x)\)	Compare relative standing

Method	Formula	Use when
Ipsatization	\(x_{ij}-\bar{x}_i\)	Respondents differ in baseline scale use
Row z-score	\((x_{ij}-\bar{x}_i)/s_i\)	Respondents differ in spread too
Compositional %	\(x_{ij}/\sum_j x_{ij}\)	Compare allocation pattern, not scale
Library-size norm.	adjust for exposure/depth	Unequal measurement effort per unit
Unit norm (row)	\(x/\lVert x \rVert\)	Compare pattern, not magnitude (e.g. text)
Rank-within-row	rank across own items	Avoid scale-use differences entirely