Outline

  • Features to read from a plot of a single quantitative example
  • Variety of ways to display a single variable and the pros and cons
  • Methods to compute null samples
  • Numerical summaries
  • Transformations (to symmetrise)
  • Handling missing values
  • Plotting a single categorical variable
  • Ordering categories

Quantitative variables

Features of a single quantitative variable

Feature Example Description
Asymmetry The distribution is not symmetrical.
Outliers Some observations are that are far from the rest.
Multimodality There are more than one "peak" in the observations.
Gaps Some continuous interval that are contained within the range but no observations exists.
Heaping Some values occur unexpectedly often.
Discretized Only certain values are found, e.g. due to rounding.
Implausible Values outside of plausible or likely range.

Example: election (1/5)

First preference votes for the Greens from the 2019 Federal election

Code 
df1 <- read_csv(here::here("data/HouseFirstPrefsByCandidateByVoteTypeDownload-24310.csv"),
  skip = 1,
  col_types = cols(
    .default = col_character(),
    OrdinaryVotes = col_double(),
    AbsentVotes = col_double(),
    ProvisionalVotes = col_double(),
    PrePollVotes = col_double(),
    PostalVotes = col_double(),
    TotalVotes = col_double(),
    Swing = col_double()
  )
)
tdf3 <- df1 |>
  group_by(DivisionID) |>
  summarise(
    DivisionNm = unique(DivisionNm),
    State = unique(StateAb),
    votes_GRN = TotalVotes[which(PartyAb == "GRN")],
    votes_total = sum(TotalVotes)
  ) |>
  mutate(perc_GRN = votes_GRN / votes_total * 100)
tdf3 |>
  ggplot(aes(perc_GRN)) +
  geom_histogram(color = "white", 
    fill = "#00843D",
    breaks = seq(0, 50, 2)) +
  labs(
    x = "First preference votes %",
    y = "Count"
  )

  • What is the pattern?
  • What are the implications for an analysis?
  • Asymmetric: Skewed right
  • One outlier
  • Multimodal (?)
  • Greens won one seat in this election
  • To have a chance of winning an electorate a party needs to be one of the top first preferences %. Greens are not close on any other electorate.
  • Is the small bump a set of electorates that could be targeted by Greens in the next election?

Example: wordle (2/5)

Distribution of my daily wordle attempts

Code 
wordle <- tibble(rows = c(1:6), count = c(0, 13, 128, 204, 147, 49))
ggplot(wordle, aes(x=rows, y=count)) + 
  geom_col(fill = "#e78c45") +
  #scale_x_discrete(breaks=1:6, labels=1:6) +
  xlab("Number of rows to solve")

  • What is the pattern?
  • What are the implications for an analysis?
  • symmetric (?)
  • barrier (?)
  • Distribution of rows needed to solve puzzle, likely cannot be a binomial?

Example: olive oil (3/5)

% composition of eicosenoic acid in Italian olive oil samples

Code 
olive <- read_csv(here::here("data/olive.csv")) |>
  rename(id = `...1`)
olive |> 
  ggplot(aes(x=eicosenoic, y=1)) +
    geom_quasirandom( 
      colour = "#007BA5", 
      alpha = 0.5) +
    scale_x_continuous(breaks = seq(0, 100, 10)) +
    ylab("") +
    theme(axis.text.y = element_blank())

  • What is the pattern?
  • What are the implications for an analysis?
  • gap
  • Another variable might explain two groups

Example: air quality (4/5)

Concentration of SO2 at sensors across Melbourne throughout 2020.

Code 
melb_air <- read_csv(here::here("data/melb_air.csv"))
ggplot(melb_air, aes(x = 1, y = so2)) + 
  geom_quasirandom(colour = "#b9ca4a") +
  xlab("") +
  ylab("SO2 concentration (ppm)") +
  theme(axis.text.x = element_blank())

  • What is the pattern?
  • What are the implications for an analysis?
  • discrete: mix of two distributions
  • Resolution for analysis may need to be reduced to the lowest

Example: air travel (5/5)

Air traffic in and out of New York City airports

Code 
library(nycflights13)
ggplot(flights, aes(x = 1, y = dep_delay)) + 
  geom_boxplot(fill = "#c397d8", alpha = 0.5) +
  xlab("") +
  ylab("Departure delay (mins)") +
  theme(aspect.ratio = 2,
    axis.text.x = element_blank())

  • What is the pattern?
  • What are the implications for an analysis?
  • outliers
  • symmetric (?)
  • implausible (?)
  • 50% of observations are within a few minutes of scheduled departure time
  • Check the values of the extremes
  • Maybe bin up the data, or remove values close to 0, so that the distribution of delayed flights can be examined

Plot type can affect detectability

Code 
gp1 <- tdf3 |>
  ggplot(aes(perc_GRN)) +
  geom_histogram(color = "white", fill = "#00843D") +
  labs(
    x = "First preference votes %",
    y = "Count"
  )
gp2 <- tdf3 |>
  ggplot(aes(perc_GRN)) +
  geom_quasirandom(color = "#00843D") +
  labs(
    x = "First preference votes %",
    y = "Count"
  )
gp2 <- tdf3 |>
  ggplot(aes(x = perc_GRN, y = 1)) +
  geom_quasirandom(color = "#00843D") +
  labs(
    x = "First preference votes %",
    y = "Count"
  )
gp3 <- tdf3 |>
  ggplot(aes(x = perc_GRN, y = 1)) +
  geom_boxplot(fill = "#00843D", alpha=0.5) +
  labs(
    x = "First preference votes %",
    y = "Count"
  )
gp1 + gp2 + gp3 + plot_layout(ncol=3)

Code 
op1 <- olive |> 
  ggplot(aes(x=eicosenoic, y=1)) +
    geom_quasirandom( 
      colour = "#007BA5", 
      alpha = 0.5) +
    scale_x_continuous(breaks = seq(0, 100, 10)) +
    ylab("") +
    theme(axis.text.y = element_blank())
op2 <- olive |> 
  ggplot(aes(x=eicosenoic)) +
    geom_density( 
      colour = "#007BA5", 
      fill = "#007BA5", 
      alpha = 0.5) +
    scale_x_continuous(breaks = seq(0, 100, 10)) +
    ylab("") +
    theme(axis.text.y = element_blank())
op1 + op2 + plot_layout(ncol=2)

Code 
ap1 <- ggplot(melb_air, aes(x = 1, y = so2)) + 
  geom_quasirandom(colour = "#b9ca4a") +
  xlab("") +
  ylab("SO2 concentration (ppm)") +
  theme(axis.text.x = element_blank())
ap2 <- ggplot(melb_air, aes(x = 1, y = so2)) + 
  geom_violin(
    quantiles = c(0.25, 0.5, 0.75),
    quantile.colour = "white",
    colour = "#b9ca4a", fill = "#b9ca4a", alpha=0.5) +
  xlab("") +
  ylab("SO2 concentration (ppm)") +
  theme(axis.text.x = element_blank())
ap3 <- ggplot(melb_air, aes(x = so2, y = 1)) + 
  stat_halfeye(fill = "#b9ca4a", alpha=0.5) +  
  stat_dots(side = "bottom",
    colour = "#b9ca4a") +
  ylab("") +
  xlab("SO2 concentration (ppm)") +
  theme(axis.text.y = element_blank())
ap1 + ap2 + ap3 + plot_layout(ncol=3)

Code 
fp1 <- ggplot(flights, aes(x = 1, y = dep_delay)) + 
  geom_boxplot(fill = "#c397d8", alpha = 0.5) +
  xlab("") +
  ylab("Departure delay (mins)") +
  theme(aspect.ratio = 2,
    axis.text.x = element_blank())
fp2 <- flights |>
  filter(between(dep_delay, -15, 15)) |>
  ggplot(aes(x = 1, y = dep_delay)) + 
  geom_boxplot(fill = "#c397d8", alpha = 0.5) +
  xlab("") +
  ylab("Departure delay (mins)") +
  theme(aspect.ratio = 2,
    axis.text.x = element_blank())
fp3 <- flights |>
  filter(dep_delay > 15) |>
  ggplot(aes(x = dep_delay)) + 
  geom_density(colour = "#c397d8", 
    fill = "#c397d8", 
    alpha = 0.5) +
  scale_x_log10() +
  xlab("") +
  ylab("Departure delay (mins)") +
  theme(aspect.ratio = 0.8)
fp4 <- flights |>
  filter(dep_delay > 15) |>
  ggplot(aes(x = dep_delay, y = 1)) + 
  geom_quasirandom(colour = "#c397d8", 
    alpha = 0.5) +
  scale_x_log10() +
  xlab("") +
  ylab("Departure delay (mins)") +
  theme(aspect.ratio = 0.8, 
    axis.text.y = element_blank(),
    axis.title.y = element_blank())
fp1 + fp2 + fp3 + fp4 + plot_layout(ncol=4)

Activity: Movies

Your turn

  1. Grab a copy of the movies data
data(movies, package = "ggplot2movies")

  1. Make a plot of the length of the movies

  2. Describe the structure

  3. Play with the type of plot, and the scale, and describe some surprising aspects of movie length

Generating null samples

Simulation

  • Normal/Gaussian: rnorm(n, mean = 0, sd = 1)
  • Uniform: runif(n, min = 0, max = 1)
  • Exponential: rexp(n, rate = 1)
  • Log-normal: rlnorm(n, meanlog = 0, sdlog = 1)
  • Gamma: rgamma(n, shape, rate = 1)
  • Beta: rbeta(n, shape1, shape2)
  • Chi-square: rchisq(n, df)
  • Student’s t: rt(n, df)
  • Weibull: rweibull(n, shape, scale = 1)
  • Cauchy: rcauchy(n, location = 0, scale = 1)
  • Binomial: rbinom(n, size, prob)
  • Negative binomial: rnbinom(n, size, prob)
  • Geometric: rgeom(n, prob)

To generate null samples that might match your data, try fitdistrplus:

  • fit <- fitdist(your_data, "norm")
  • params <- fit$estimate
  • simulated <- rnorm(n, mean = params["mean"], sd = params["sd"])

Or compute summary statistics:

Normal

mu    <- mean(your_data)
sigma <- sd(your_data)

Exponential

rate <- 1 / mean(your_data)

Example: election (1/2)

Code 
set.seed(241)
l <- lineup(null_dist("perc_GRN", "exp"), tdf3, n=12, pos=1)
ggplot(l, 
       aes(x=perc_GRN)) +
  geom_histogram(color = "white", 
    fill = "#00843D", 
    breaks = seq(0, 50, 2)) +
  xlim(c(0, 50)) +
  facet_wrap(~.sample, ncol=6) +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid.major = element_blank())

Code 
library(fitdistrplus)
fit <- fitdist(tdf3$perc_GRN / 100, "beta")

set.seed(446)
l2 <- tdf3 |>
  select(perc_GRN) |>
  mutate(.sample = 1)
for (i in 1:11) {
  x <- rbeta(nrow(tdf3), fit$estimate[1], fit$estimate[2])*100
  d <- tibble(perc_GRN = x,
     .sample = i+1)
  l2 <- bind_rows(l2, d)
}
ggplot(l2, 
       aes(x=perc_GRN)) +
  geom_histogram(color = "white", 
    fill = "#00843D", 
    breaks = seq(0, 50, 2)) +
  xlim(c(0, 50)) +
  facet_wrap(~.sample, ncol=6, scale="free_y") +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid.major = element_blank())

  • The bump is may not be unusual for a sample from a standard distribution.
  • The single outlier is not typically seen in simulated samples.

Example: wordle (2/2)

Code 
wordle_uncount <- uncount(wordle, count)
p <- mean(wordle_uncount$rows)/6
set.seed(501)
wl <- wordle |>
  mutate(.sample = 1)
for (i in 1:11) {
  x <- rbinom(nrow(wordle_uncount), 6, p)
  d <- tibble(rows = x) |>
    count(rows) |>
    rename(count = n) |>
    mutate(.sample = i+1)
  wl <- bind_rows(wl, d)
}
ggplot(wl, 
       aes(x=rows, y=count)) +
  geom_col(fill = "#e78c45") +
  facet_wrap(~.sample, ncol=6) +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid.major = element_blank())

The distribution of actual wordle results is different from samples from a binomial

  • counts for 1-2 rows is too low
  • peak at 4 is too high

Numerical measures of a single quantitative variables

Central tendency and dispersion

  • A measure of central tendency, e.g. mean, median and mode

  • A measure of dispersion (also called variability or spread), e.g. variance, standard deviation and interquartile range

  • There are other measures, e.g. skewness and kurtosis that measures “tailedness”, but these are not as common as the measures of first two

  • The mean is also the first moment and variance, skewness and kurtosis are second, third, and fourth central moments

Robust measure of central tendency

  • Mean is a non-robust measure of location.
  • Median is the 50% quantile of the observations
  • Trimmed mean is the sample mean after discarding observations at the tails.
  • Winsorized mean is the sample mean after replacing observations at the tails with the minimum or maximum of the observations that remain.

Both trimmed and Winsorized mean trimmed 20% of the tails.

Plot Mean Median Trimmed Mean Winsorized Mean
1 0.109 0.114 0.120 0.103
2 0.054 -0.045 -0.016 -0.029
3 1.177 0.729 0.820 0.888
4 0.533 0.541 0.543 0.542
5 0.468 0.329 0.355 0.390
6 5.626 6.656 5.918 5.688

Robust measure of dispersion

  • Standard deviation or its square, variance, is a popular choice of measure of dispersion but is not robust to outliers
  • Standard deviation for sample \(x_1, ..., x_n\) is

\[\sqrt{\sum_{i=1}^n \frac{(x_i - \bar{x})^2}{n - 1}}\]

  • Interquartile range difference between 1st and 3rd quartile, more robust measure of spread
  • Median absolute deviance (MAD) is even more robust

\[\text{median}(|x_i - \text{median}(x_i)|)\]

Measure of dispersion
Plot SD IQR MAD Skewness Kurtosis
1 0.90 1.19 0.87 -0.072 3.0
2 0.99 1.41 1.08 0.358 2.2
3 1.33 1.18 0.79 1.944 7.2
4 0.29 0.45 0.34 -0.126 1.8
5 0.47 0.50 0.34 1.691 6.4
6 2.78 5.36 2.98 -0.351 1.7

Example: election (1/3)

% of first preference for the Greens
Mean Median SD MAD IQR Skewness Kurtosis
9.9 8.5 5.6 3.8 5 2.7 16

  • Why are the means and the medians different?

  • How are the standard deviations and the interquartile ranges similar or different?

  • Are there some other numerical statistics we should show?

Code 
tdf3 |>
 summarise(
    mean = mean(perc_GRN),
    median = median(perc_GRN),
    sd = sd(perc_GRN),
    mad = mad(perc_GRN),
    iqr = IQR(perc_GRN),
    skewness = moments::skewness(perc_GRN),
    kurtosis = moments::kurtosis(perc_GRN)
  ) |>
  knitr::kable(col.names = c("Mean", "Median", "SD", "MAD", "IQR", "Skewness", "Kurtosis"), digits = 3) |>
  kableExtra::kable_classic() |>
  kableExtra::add_header_above(c(" ", "% of first preference for the Greens" = 5, " "))

Example: olive oil (2/3)

% composition of eicosenoic acid
Mean Median SD MAD IQR Skewness Kurtosis
16 17 14 21 26 0.34 1.9

  • Why are the mean and the median so different?

  • How do the standard deviation and the interquartile range compare?

  • What do the skewness and kurtosis values suggest about the shape of the distribution?

Code 
olive |>
  summarise(
    mean = mean(eicosenoic),
    median = median(eicosenoic),
    sd = sd(eicosenoic),
    mad = mad(eicosenoic),
    iqr = IQR(eicosenoic),
    skewness = moments::skewness(eicosenoic),
    kurtosis = moments::kurtosis(eicosenoic)
  ) |>
  knitr::kable(col.names = c("Mean", "Median", "SD", "MAD", "IQR", "Skewness", "Kurtosis"), digits = 3) |>
  kableExtra::kable_classic() |>
  kableExtra::add_header_above(c("% composition of eicosenoic acid" = 7))

Example: air quality (3/3)

SO2 concentration (ppm)
Mean Median SD MAD IQR Skewness Kurtosis
0.00063 0.00041 0.00074 0.00037 0.00051 2.5 9.7

  • Why are the mean and the median so different?

  • How do the standard deviation and the interquartile range compare?

  • What do the skewness and kurtosis values suggest about the shape of the distribution?

Code 
melb_air |>
  summarise(
    mean = mean(so2, na.rm = TRUE),
    median = median(so2, na.rm = TRUE),
    sd = sd(so2, na.rm = TRUE),
    mad = mad(so2, na.rm = TRUE),
    iqr = IQR(so2, na.rm = TRUE),
    skewness = moments::skewness(so2, na.rm = TRUE),
    kurtosis = moments::kurtosis(so2, na.rm = TRUE)
  ) |>
  knitr::kable(col.names = c("Mean", "Median", "SD", "MAD", "IQR", "Skewness", "Kurtosis"), digits = 5) |>
  kableExtra::kable_classic() |>
  kableExtra::add_header_above(c("SO2 concentration (ppm)" = 7))

Inference for robust statistics

We have seen the re-sampling methods simulation and permutation used for generating null plots in a lineup. Re-sampling methods can be used with numeric statistics also.

Simulation from distribution, can be used to to check for outliers.

Code 
# Estimate the parameters without the outlier
est_r <- fitdist(tdf3$perc_GRN[tdf3$perc_GRN < 40]/100, "beta")
# Check fit
# ggplot(tdf3, aes(sample=perc_GRN)) + stat_qq(distribution = stats::qexp, dparams = est_r) + stat_qq_line(distribution = stats::qexp, dparams = est_r)
set.seed(912)
samp <- matrix(rbeta(n=151*100, 
  shape1=est_r$estimate[1], 
  shape2=est_r$estimate[2])*100, ncol=100, byrow=TRUE)
samp_max <- apply(samp, 2, max)
samp_max_df <- tibble(m = samp_max)
ggplot(samp_max_df, aes(x=m)) +
  geom_histogram(binwidth=2.5, fill="grey60", 
    colour="white") +
  xlim(c(0, 60)) +
  geom_vline(xintercept=
    tdf3$perc_GRN[tdf3$perc_GRN > 40], colour="#D93F00") +
  annotate("text", x=42, y=13, label="observed", colour="#D93F00") +
  xlab("Simulated maxima") +
  theme(aspect.ratio = 0.5)

We can also compute how many simulated values are more than the observed which gives a simulation \(p\)-value: 0.

For sample means, conventional tests provide a means for assessing what might be observed if different samples were taken.

Bootstrapping the current sample, can be used for robust statistics. If we have a sample of values:

[1] 2 2 3 6 7 7 8 8

to bootstrap sample with replacement:

Code 
sort(sample(x, replace=TRUE))
[1] 2 2 3 3 7 7 7 7
Code 
sort(sample(x, replace=TRUE))
[1] 2 3 6 6 6 6 8 8



Here’s an example of bootstrapping to get a confidence interval for a median.

[1] "Median: 6.34"
[1] "95% CI: ( 4.99 , 9.16 )"

Activity: Air quality

Your turn

  1. Look at the pm25 variable in the melb_air data (fine particulate matter, µg/m³)

  2. Compute both robust (median, MAD, IQR) and non-robust (mean, SD) summary statistics – how much do they disagree, and why?

  3. Bootstrap a 95% confidence interval for the median

  4. How might these quantities be used to determine reasonable levels of PM25, and indicate extreme values?

A comment about boxplots (1/2)

Code 
tdf3 |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(perc_GRN, State)) +
  geom_boxplot(varwidth = TRUE) +
  labs(
    x = "First preference votes %",
    y = "Count",
    title = "Greens party"
  )

Where are these electorates?

The width of the boxplot is proportional to the number of electoral districts in the corresponding state (which is roughly proportional to the population).

Code 
tdf3 |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(perc_GRN, State)) +
  geom_boxplot(varwidth = TRUE) +
  labs(
    x = "First preference votes %",
    y = "Count",
    title = "Greens party"
  )

A comment about boxplots (2/2)

Code 
tdf3 |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(perc_GRN, State)) +
  ggbeeswarm::geom_quasirandom(groupOnX = FALSE, varwidth = TRUE) +
  labs(
    x = "First preference votes %",
    y = "State",
    title = "Greens party"
  )

Now what do you notice from this graph that you didn’t notice before?

  • Only two electoral districts in NT.

  • And only 3 and 5 electoral districts in ACT and TAS, respectively!

  • Boxplots requires 5 points!

  • We should have summarised the number of electoral districts for each state with numerical statistics as a first step.

  • Also the outlier (yes, safe to call this an outlier!) and the cluster in the Victoria electorates.

Code 
tdf3 |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(perc_GRN, State)) +
  ggbeeswarm::geom_quasirandom(groupOnX = FALSE, varwidth = TRUE) +
  labs(
    x = "First preference votes %",
    y = "State",
    title = "Greens party"
  )

Transformations

Box-Cox

The Box-Cox transformation is a family of power transformations, indexed by a parameter \(\lambda\), used to make a variable more symmetric (and closer to normally distributed):

\[ y^{(\lambda)} = \begin{cases} \dfrac{y^\lambda - 1}{\lambda}, & \lambda \neq 0 \\[4pt] \log(y), & \lambda = 0 \end{cases} \]

  • Only defined for \(y > 0\)
  • \(\lambda\) can be chosen to maximise the log-likelihood of the transformed values being normal, e.g. MASS::boxcox()
  • The transformation is a guide: round \(\lambda\) to a nearby, interpretable value (e.g. \(0\) for \(\log\), \(0.5\) for \(\sqrt{\ }\))

The ladder of transformations

Tukey’s ladder of power transformations arranges the Box-Cox family by \(\lambda\). Moving down the ladder pulls in long right tails; moving up the ladder pulls in long left tails.

\(\lambda\) Transformation Typical use
2 \(y^2\) strongly left-skewed
1 \(y\) (no change) roughly symmetric
0.5 \(\sqrt{y}\) mild right-skew, count data
0 \(\log(y)\) moderate to strong right-skew
-0.5 \(1/\sqrt{y}\) strong right-skew
-1 \(1/y\) very strong right-skew, e.g. rates
  • Start close to \(\lambda=1\) (no transformation) and move up or down the ladder, checking symmetry with a plot after each step
  • The further \(\lambda\) moves from 1, the stronger the transformation

Example: air quality (PM2.5)

Original pm25 (µg/m³) – 34 zero values are shifted by \(+0.1\) before transforming, since Box-Cox requires \(y>0\):

Code 
pm25 <- melb_air$pm25[!is.na(melb_air$pm25)]
pm25_shift <- pm25 + 0.1
pm25p1 <- ggplot(tibble(pm25 = pm25), aes(x = pm25, y = 1)) +
  geom_quasirandom(colour = "#b9ca4a", alpha = 0.5) +
  ylab("") + xlab("PM2.5 (µg/m³)") +
  theme(axis.text.y = element_blank())

Box-Cox log-likelihood profile:

Code 
bc_pm25 <- MASS::boxcox(pm25_shift ~ 1, plotit = FALSE)
lambda_pm25 <- bc_pm25$x[which.max(bc_pm25$y)]
tibble(lambda = bc_pm25$x, loglik = bc_pm25$y) |>
  ggplot(aes(x = lambda, y = loglik)) +
  geom_line(colour = "#b9ca4a") +
  geom_vline(xintercept = lambda_pm25, linetype = "dashed", colour = "#D93F00") +
  labs(x = expression(lambda), y = "Log-likelihood") +
  theme(aspect.ratio = 0.4)

Optimal \(\lambda \approx\) 0.2

Transformed variable, using the optimal \(\lambda\) directly (not rounded to a ladder rung):

Code 
pm25_bc <- (pm25_shift^lambda_pm25 - 1) / lambda_pm25
pm25p2 <- ggplot(tibble(pm25_bc = pm25_bc), aes(x = pm25_bc, y = 1)) +
  geom_quasirandom(colour = "#b9ca4a", alpha = 0.5) +
  ylab("") + xlab("Box-Cox transformed PM2.5") +
  theme(axis.text.y = element_blank())
pm25p1  + pm25p2 + plot_layout(ncol=2)

  • Skewness before transforming: 11.04
  • Skewness after transforming: -0.17

\(\lambda \approx 0.2\) sits between \(\sqrt{y}\) (\(\lambda=0.5\)) and \(\log(y)\) (\(\lambda=0\)) on the ladder. Rounding to \(\log(y)\) actually overcorrects here (skewness becomes -1.33, i.e. left-skewed) – a reminder that the optimal \(\lambda\) is not always well-approximated by a “nice” rung when there are extreme outliers (here, from bushfire smoke).

Code 
pm25 <- melb_air$pm25[!is.na(melb_air$pm25)]
pm25_shift <- pm25 + 0.1
pm25p1 <- ggplot(tibble(pm25 = pm25), aes(x = pm25, y = 1)) +
  geom_quasirandom(colour = "#b9ca4a", alpha = 0.5) +
  ylab("") + xlab("PM2.5 (µg/m³)") +
  theme(axis.text.y = element_blank())
Code 
bc_pm25 <- MASS::boxcox(pm25_shift ~ 1, plotit = FALSE)
lambda_pm25 <- bc_pm25$x[which.max(bc_pm25$y)]
tibble(lambda = bc_pm25$x, loglik = bc_pm25$y) |>
  ggplot(aes(x = lambda, y = loglik)) +
  geom_line(colour = "#b9ca4a") +
  geom_vline(xintercept = lambda_pm25, linetype = "dashed", colour = "#D93F00") +
  labs(x = expression(lambda), y = "Log-likelihood") +
  theme(aspect.ratio = 0.4)
Code 
pm25_bc <- (pm25_shift^lambda_pm25 - 1) / lambda_pm25
pm25p2 <- ggplot(tibble(pm25_bc = pm25_bc), aes(x = pm25_bc, y = 1)) +
  geom_quasirandom(colour = "#b9ca4a", alpha = 0.5) +
  ylab("") + xlab("Box-Cox transformed PM2.5") +
  theme(axis.text.y = element_blank())
pm25p1  + pm25p2 + plot_layout(ncol=2)

Non box-cox example: Olive oil (1/2)

Code 
ggplot(olive, aes(x = eicosenoic, y = 1)) +
  geom_quasirandom(colour = "#007BA5", alpha = 0.5) +
  ylab("") +
  xlab("Eicosenoic acid %") +
  theme(axis.text.y = element_blank())

Two clusters

Check if there is another variable that explains the two groups

Code 
ggplot(olive, aes(x = eicosenoic, y = 1)) +
  geom_quasirandom(colour = "#007BA5", alpha = 0.5) +
  facet_wrap(~region, ncol=3) +
  ylab("") +
  xlab("Eicosenoic acid %") +
  theme(axis.text.y = element_blank())

Non box-cox example: Olive oil (2/2)

Mixture of discreteness and normal shape of continuous values. Why might this happen?

Check if there is a difference in the strata (here 1 thru 9), implying measurement policy differences.

Activity: Air quality

Your turn

When and why do you want to make a transformation of a variable?

Handling missing values

Melbourne Housing Prices (1/5)

Code 
df2 |>
  head(20) |>
  select(Suburb, Rooms, Type, Price, Date) |>
  mutate(
    Price = scales::comma(Price),
    Type = fct_recode(Type,
      "Home" = "h",
      "Townhouse" = "t",
      "Unit" = "u"
    )
  ) |>
  knitr::kable(
    col.names = c("Suburb", "Rooms", "Type", "Price ($)", "Date"),
    align = "lrlr"
  ) |>
  kableExtra::kable_classic() |>
  kableExtra::kable_styling(font_size = 24,
    full_width=FALSE)
Suburb Rooms Type Price ($) Date
Abbotsford 3 Home 1,490,000 2017-04-01
Abbotsford 3 Home 1,220,000 2017-04-01
Abbotsford 3 Home 1,420,000 2017-04-01
Aberfeldie 3 Home 1,515,000 2017-04-01
Airport West 2 Home 670,000 2017-04-01
Airport West 2 Townhouse 530,000 2017-04-01
Airport West 2 Unit 540,000 2017-04-01
Airport West 3 Home 715,000 2017-04-01
Albanvale 6 Home NA 2017-04-01
Albert Park 3 Home 1,925,000 2017-04-01
Albion 3 Unit 515,000 2017-04-01
Albion 4 Home 717,000 2017-04-01
Alphington 2 Home 1,675,000 2017-04-01
Alphington 4 Home 2,008,000 2017-04-01
Altona 2 Home 860,000 2017-04-01
Altona Meadows 4 Home NA 2017-04-01
Altona North 3 Home 720,000 2017-04-01
Armadale 2 Unit 836,000 2017-04-01
Armadale 2 Home 2,110,000 2017-04-01
Armadale 3 Home 1,386,000 2017-04-01
  • This data was scraped each week from domain.com.au from 2016-01-28 to 2018-10-13
  • In total there are 63,023 observations
  • All variables shown (there are more variables not shown here), except price, have complete records
  • The are 48,433 property prices across Melbourne (roughly 23% missing)

Data source: Tony Pio (2018) Melbourne Housing Market

How would you explore this data first?

Yes, with an overview plot.

Melbourne Housing Prices (2/5)

There are only missings on price, but a lot of them. Can we find some other collected variable which might help to explain the missing price?

Is missingness more likely for expensive houses?

Use a substitute like Rooms to examine the possible relationship. Show Rooms vs missing on Price and use a lineup before looking at the data.

Rooms is discrete, so a bar chart, possibly facaetted or side-by-side bars, would be a suitable plot.

  • To impute missings other variables will need to be used.

Note: Houses with more than 8 rooms removed. Why?

Code 
df2 <- read_csv(here::here("data/MELBOURNE_HOUSE_PRICES_LESS.csv"),
  col_types = cols(
    .default = col_character(),
    Rooms = col_double(),
    Price = col_double(),
    Date = col_date(format = "%d/%m/%Y"),
    Propertycount = col_double(),
    Distance = col_double()
  )
)
Code 
skimr::skim(df2)
── Data Summary ────────────────────────
                           Values
Name                       df2   
Number of rows             63023 
Number of columns          13    
_______________________          
Column type frequency:           
  character                8     
  Date                     1     
  numeric                  4     
________________________         
Group variables            None  

── Variable type: character ────────────────────────────────
  skim_variable n_missing complete_rate min max empty
1 Suburb                0             1   3  18     0
2 Address               0             1   7  27     0
3 Type                  0             1   1   1     0
4 Method                0             1   1   2     0
5 SellerG               0             1   1  27     0
6 Postcode              0             1   4   4     0
7 Regionname            0             1  16  26     0
8 CouncilArea           0             1  17  30     0
  n_unique whitespace
1      380          0
2    57754          0
3        3          0
4        9          0
5      476          0
6      225          0
7        8          0
8       34          0

── Variable type: Date ─────────────────────────────────────
  skim_variable n_missing complete_rate min       
1 Date                  0             1 2016-01-28
  max        median     n_unique
1 2018-10-13 2017-09-03      112

── Variable type: numeric ──────────────────────────────────
  skim_variable n_missing complete_rate      mean         sd
1 Rooms                 0         1          3.11      0.958
2 Price             14590         0.768 997898.   593499.   
3 Propertycount         0         1       7618.     4424.   
4 Distance              0         1         12.7       7.59 
     p0    p25      p50       p75       p100 hist 
1     1      3      3         4         31   ▇▁▁▁▁
2 85000 620000 830000   1220000   11200000   ▇▁▁▁▁
3    39   4380   6795     10412      21650   ▅▇▅▂▁
4     0      7     11.4      16.7       64.1 ▇▆▁▁▁
Code 
df2 |>
  select(Suburb, Rooms, Type, Price, Date) |>
  arrange(Suburb, Date) |>
  visdat::vis_miss()
Code 
df2 |>
  mutate(miss = ifelse(is.na(Price), 
    "Missing", "Recorded")) |>
  count(Rooms, miss) |>
  filter(Rooms < 8) |>
  group_by(miss) |>
  mutate(perc = n / sum(n) * 100) |>
  ggplot(aes(as.factor(Rooms), perc, fill = miss)) +
    geom_col(position = "dodge") +
    scale_fill_viridis_d(begin=0.3, end=0.7) +
    labs(x = "Rooms", y = "Percentage", fill = "Price") +
    theme(aspect.ratio = 0.8)

Is there a suspicious plot?

Break the association between Rooms and Missing/Not on Price, because the null hypothesis is that there is no difference in missing status for price based on the size of the house. Why?

Code 
library(nullabor)
df2_d <- df2 |>
  mutate(miss = ifelse(is.na(Price), "Missing", "Recorded")) |>
  select(Rooms, miss) |>
  filter(Rooms < 8)
df2_l <- lineup(null_permute("miss"), df2_d, n=10, pos=7) 
df2_l_agg <- df2_l |>
  group_by(.sample) |>
  count(Rooms, miss) |>
  ungroup() |>
  group_by(miss) |>
  mutate(perc = n / sum(n) * 100) |>
  mutate(Rooms = as.factor(Rooms))
ggplot(df2_l_agg, aes(x=Rooms, y=perc, fill = miss)) +
  geom_col(position = "dodge") +
  scale_fill_viridis_d(begin=0.3, end=0.7) +
  facet_wrap(~.sample, ncol=5) +
  theme(legend.position = "none", 
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid.major.x = element_blank())

Check the support of your data

If you have too few measurements in any region (extreme), summaries for these regions will be unreliable.

  • For quantitative variables, it may be necessary to remove extremes.
  • If the variable is categorical it might be best to combine levels.
  • It is important to script so decisions can be reversed or rare events are not ignored.


We removed houses with 8 or more rooms. What other way might we have handled these houses?

Melbourne Housing Prices (3/5)

What can we say from this plot?

  • The housing prices are right-skewed
  • There appears to be a lot of outlying housing prices (how can we tell?)

Note: We determined that it is likely that more higher price houses have not disclosed the sale price. The distribution of price will need to be checked again after imputation.

Code 
df2 <- read_csv(here::here("data/MELBOURNE_HOUSE_PRICES_LESS.csv"),
  col_types = cols(
    .default = col_character(),
    Rooms = col_double(),
    Price = col_double(),
    Date = col_date(format = "%d/%m/%Y"),
    Propertycount = col_double(),
    Distance = col_double()
  )
)
Code 
df2 |>
  ggplot(aes(Price / 1e6)) +
  geom_histogram(color = "white") +
  labs(
    x = "Price (mil)",
    y = "Count"
  )

Melbourne Housing Prices (4/5)

  • The x-axis has been \(\log_{10}\)-transformed in this plot
  • The plot appears more symmetrical now
  • What is a useful measure of central tendency here?
Code 
df2 <- read_csv(here::here("data/MELBOURNE_HOUSE_PRICES_LESS.csv"),
  col_types = cols(
    .default = col_character(),
    Rooms = col_double(),
    Price = col_double(),
    Date = col_date(format = "%d/%m/%Y"),
    Propertycount = col_double(),
    Distance = col_double()
  )
)
Code 
df2 |>
  ggplot(aes(Price / 1e6)) +
  geom_histogram(color = "white") +
  labs(
    x = "Price (mil)",
    y = "Count"
  ) +
  scale_x_log10()

Melbourne Housing Prices (5/5)

With no transformation:

Mean Median Trimmed Mean Winsorised Mean
$997,898 $830,000 $871,375 $903,823


With log transformation (and back-transformed to original scale):

Mean Median Trimmed Mean Winsorised Mean
$874,166 $830,000 $847,973 $859,325
Code 
df2 |>
  filter(!is.na(Price)) |>
  summarise(
    Mean = scales::dollar(mean(Price)),
    Median = scales::dollar(median(Price)),
    `Trimmed Mean` = scales::dollar(mean(Price, trim = 0.2)),
    `Winsorised Mean` = scales::dollar(psych::winsor.mean(Price))
  ) |>
  knitr::kable(align = "r") |>
  kableExtra::kable_classic() |>
  kableExtra::kable_styling(full_width=FALSE)
Code 
df2 |>
  filter(!is.na(Price)) |>
  mutate(lPrice = log10(Price)) |>
  summarise(
    Mean = scales::dollar(10^mean(lPrice)),
    Median = scales::dollar(10^median(lPrice)),
    `Trimmed Mean` = scales::dollar(10^mean(lPrice, trim = 0.2)),
    `Winsorised Mean` = scales::dollar(10^psych::winsor.mean(lPrice))
  ) |>
  knitr::kable(align = "r") |>
  kableExtra::kable_classic() |>
  kableExtra::kable_styling(full_width=FALSE)

Categorical variables

There are two types of categorical variables



Nominal where there is no intrinsic ordering to the categories
E.g. blue, grey, black, white.


Ordinal where there is a clear order to the categories.
E.g. Strongly disagree, disagree, neutral, agree, strongly agree.

Categorical variables in R

  • In R, categorical variables may be encoded as factors.
Code 
data <- c(2, 2, 1, 1, 3, 3, 3, 1)
factor(data)
[1] 2 2 1 1 3 3 3 1
Levels: 1 2 3
  • You can easily change the labels of the variables:
Code 
factor(data, labels = c("I", "II", "III"))
[1] II  II  I   I   III III III I  
Levels: I II III
  • Order of the factors are determined by the input:
Code 
# numerical input are ordered in increasing order 
factor(c(1, 3, 10))
[1] 1  3  10
Levels: 1 3 10
Code 
# character input are ordered by first char, alphabetically 
factor(c("1", "3", "10"))
[1] 1  3  10
Levels: 1 10 3
Code 
# you can specify order of levels explicitly 
factor(c("1", "3", "10"),
  levels = c("1", "3", "10")
)
[1] 1  3  10
Levels: 1 3 10

Order nominal variables meaningfully


Coding tip: use below functions to easily change the order of factor levels



stats::reorder(factor, value, mean)
forcats::fct_reorder(factor, value, median)
forcats::fct_reorder2(factor, value1, value2, func)

Numerical summaries

counts, proportions, percentages and odds

Tuberculosis counts in Australia

Code 
options(digits=2)
tb_oz |>
  filter(year >= 2000) |>
  mutate(p = count/sum(count),
         pct = p*100, 
         odds = count/count[year==2000]) |>
  print(n=100)
# A tibble: 22 × 7
   country   iso3   year count      p   pct  odds
   <chr>     <chr> <dbl> <dbl>  <dbl> <dbl> <dbl>
 1 Australia AUS    2000   982 0.0522  5.22 1    
 2 Australia AUS    2001   953 0.0507  5.07 0.970
 3 Australia AUS    2002  1008 0.0536  5.36 1.03 
 4 Australia AUS    2003   926 0.0493  4.93 0.943
 5 Australia AUS    2004  1036 0.0551  5.51 1.05 
 6 Australia AUS    2005  1030 0.0548  5.48 1.05 
 7 Australia AUS    2006  1127 0.0600  6.00 1.15 
 8 Australia AUS    2007  1081 0.0575  5.75 1.10 
 9 Australia AUS    2008  1182 0.0629  6.29 1.20 
10 Australia AUS    2009  1176 0.0626  6.26 1.20 
11 Australia AUS    2010  1146 0.0610  6.10 1.17 
12 Australia AUS    2011  1202 0.0640  6.40 1.22 
13 Australia AUS    2012  1259 0.0670  6.70 1.28 
14 Australia AUS    2013   512 0.0272  2.72 0.521
15 Australia AUS    2014   474 0.0252  2.52 0.483
16 Australia AUS    2015   438 0.0233  2.33 0.446
17 Australia AUS    2016   481 0.0256  2.56 0.490
18 Australia AUS    2017   524 0.0279  2.79 0.534
19 Australia AUS    2018   502 0.0267  2.67 0.511
20 Australia AUS    2019   554 0.0295  2.95 0.564
21 Australia AUS    2020   609 0.0324  3.24 0.620
22 Australia AUS    2021   593 0.0316  3.16 0.604

For qualitative data, compute

  • count/frequency,
  • proportion/percentage
  • and sometimes, an odds ratio. Here we have used ratio relative to the count in year 2000.



Note: For exploration, no rounding of digits was done, but to report you would need to make the numbers pretty.

Ways to plot a single categorical variable

Code 
wordle2 <- wordle |>
  mutate(attempts = factor(rows))

p_bar <- ggplot(wordle2, aes(x = attempts, y = count)) +
  geom_col(aes(fill = attempts)) +
  scale_fill_viridis_d() +
  labs(x = "Number of attempts", y = "Count") +
  theme(legend.position = "none") +
  ggtitle("Bar chart")

p_pie <- ggplot(wordle2, aes(x = "", y = count, fill = attempts)) +
  geom_col(width = 1, colour = "white") +
  coord_polar(theta = "y") +
  scale_fill_viridis_d() +
  labs(fill = "Attempts", x = NULL, y = NULL) +
  theme_void() +
  theme(legend.position = "none") +
  ggtitle("Pie chart")

p_rose <- ggplot(wordle2, aes(x = attempts, y = count, fill = attempts)) +
  geom_col(colour = "white") +
  coord_polar(theta = "x") +
  scale_fill_viridis_d() +
  labs(x = NULL, y = "Count", fill = "Attempts") +
  theme(axis.text.y = element_blank(), legend.position = "none") +
  ggtitle("Rose plot")

p_spine <- ggplot(wordle2, aes(x = "", y = count, fill = attempts)) +
  geom_col(position = "fill", colour = "white") +
  scale_fill_viridis_d() +
  coord_flip() +
  labs(x = NULL, y = "Proportion", fill = "Attempts") +
  theme(axis.text.y = element_blank()) +
  ggtitle("Spine plot")

p_bar + p_pie + p_rose + p_spine + plot_layout(ncol=2)

Number of attempts to solve Wordle:

  • Bar chart: height encodes count directly – easiest to compare categories accurately
  • Pie chart: angle (and area) encodes count – comparing categories by angle is harder than by length
  • Rose plot (nightingale/coxcomb): like a bar chart wrapped into a circle – radius encodes count, but the outer wedges look larger than their true count because area grows with the square of the radius
  • Spine plot: a single bar split into proportions – good for comparing proportions to the whole, or as a building block for two categorical variables (stacking within categories of a second variable)

Generally, bar charts (or spine plots) are the easiest to read accurately: humans compare lengths/positions along a common scale far better than angles or areas.

Activity: Melbourne housing

Your turn

  • Make the spine plot of Rooms and fill the bars based on missingness on Price.
  • If you can do this using a lineup, even better.

The code on slide 35 might help.

Imputing missings for univariate distributions

Quantitative variable: Simulate from a fitted distribution.

Code 
df2 <- df2 |>
  mutate(lPrice = log10(Price),
         price_miss = ifelse(is.na(Price), "yes", "no"))

df2 <- df2 |>
  filter(Rooms < 20) # remove one extreme
df2_fit <- lm(lPrice~Rooms, df2) 
coefs <- tidy(df2_fit)
fitstats <- glance(df2_fit)

set.seed(1003)  
df2_nomiss <- df2 |>
  filter(price_miss == "no")
df2_miss <- df2 |>
  filter(price_miss == "yes")
df2_miss <- df2_miss |>
  rowwise() |>
  mutate(lPrice = 
    coefs$estimate[1]+coefs$estimate[2]*Rooms +
    rnorm(1, 0, fitstats$sigma)) |>
  mutate(Price = ifelse(price_miss == "yes", 10^lPrice, Price))
df2_nomiss <- bind_rows(df2_nomiss, df2_miss)
df2_nomiss <- df2_nomiss |>
  mutate(price_miss = factor(price_miss, levels = c("yes", "no")))
Code 
mp1 <- ggplot(df2, aes(x=Rooms, y=lPrice)) + 
  geom_miss_point() 

mp2 <- ggplot(df2_nomiss, aes(x=Rooms, 
                               y=lPrice, 
                               colour=price_miss)) +
  geom_jitter(width=0.5, alpha=0.1)

mp3 <- df2_nomiss |>
  ggplot() +
  geom_histogram(aes(x=Price / 1e6, 
                     y=after_stat(density)), 
                 color = "white") +
  facet_wrap(~price_miss, ncol=2, 
              scales="free_y") +
  labs(
    x = "Price (mil)",
    y = "Count"
  ) +
  scale_x_log10()
mp1 + mp2 + plot_layout(ncol=2)

Code 
mp3

Categorical variable: Simulate from multinomial.

Code 
tb_oz_age <- tb |> 
  filter(iso3 == "AUS", year == 2012) |>
  select(contains("new_sp_f")) |>
  select(-new_sp_f04, -new_sp_f514, -new_sp_f014) |>
  pivot_longer(new_sp_f1524:new_sp_fu, 
    names_to="age", 
    values_to="count") |>
  mutate(age = str_remove(age, "new_sp_f"))
# Add some missing count
tb_oz_age$count[7] <- 12
tb_oz_age
# A tibble: 7 × 2
  age   count
  <chr> <dbl>
1 1524     27
2 2534     48
3 3544     15
4 4554     11
5 5564      9
6 65       15
7 u        12
Code 
tb_oz_age_long <- tb_oz_age |>
  uncount(count) |>
  mutate(age = ifelse(age == "u", NA, age))
set.seed(153)
fill_miss <- rbinom(tb_oz_age$count[7], size=5,
  prob=tb_oz_age$count[1:6]/sum(tb_oz_age$count[1:6]))+1
tb_oz_age_impute <- tb_oz_age 
for (i in 1:length(fill_miss)) 
  tb_oz_age_impute$count[fill_miss[i]] <-
    tb_oz_age_impute$count[fill_miss[i]] + 1
fill_miss
 [1] 5 3 1 1 2 1 1 2 1 1 1 1
Code 
tb_oz_age_impute
# A tibble: 7 × 2
  age   count
  <chr> <dbl>
1 1524     35
2 2534     50
3 3544     16
4 4554     11
5 5564     10
6 65       15
7 u        12

imputeMulti library can automate for multiple variables.

Key points

  • Make many plots of a single variable (histogram, density, boxplot, quasirandom) – the type of plot, and choices like bin width or bandwidth, change what patterns are visible
  • Use simulation and permutation to generate null samples and lineups, to judge whether a striking feature (outlier, cluster, shape) is more than what’s expected by chance
  • Compute both non-robust (mean, SD) and robust (median, MAD, IQR, trimmed/winsorised mean) statistics – skew and outliers can pull the non-robust ones a long way from the bulk of the data
  • Boxplots are a compact robust summary, but can hide multimodality – check against the raw data
  • Transformations (e.g. Box-Cox, the ladder of power transformations) can re-focus attention and symmetrise skewed data, but round the “optimal” \(\lambda\) with care, especially when extreme outliers are present
  • Check for missing values early: visualise the pattern of missingness (e.g. vis_miss(), a lineup) before deciding how, or whether, to impute – and script decisions so they can be reversed
  • For a categorical variable, compute counts, proportions/percentages and odds, and order categories meaningfully rather than relying on the default ordering

Resources

  • Unwin (2015) Graphical Data Analysis with R
  • Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0
  • Josse et al (2022) R-miss-tastic
  • Wilke Fundamentals of Data Visualization, chapters 6-11.