flightsdataset from nycflights package

pacman::p_load("nycflights13")Copy Code

Columns: 19
$ year           <int> 2013,…
$ month          <int> 1, 1,…
$ day            <int> 1, 1,…
$ dep_time       <int> 517, …
$ sched_dep_time <int> 515, …
$ dep_delay      <dbl> 2, 4,…
$ arr_time       <int> 830, …
$ sched_arr_time <int> 819, …
$ arr_delay      <dbl> 11, 2…
$ carrier        <chr> "UA",…
$ flight         <int> 1545,…
$ tailnum        <chr> "N142…
$ origin         <chr> "EWR"…
$ dest           <chr> "IAH"…
$ air_time       <dbl> 227, …
$ distance       <dbl> 1400,…
$ hour           <dbl> 5, 5,…
$ minute         <dbl> 15, 2…
$ time_hour      <dttm> 2013…Copy Code

Definition(s)

A sample is a sequence of items of the same type

In , if the type matches a basetype, a sample is typically represented using a vector

A numerical (quantitative) sample is just a sequence of numbers (either integers or floats), it is represented by a numeric vector.

A categorical (qualitative) sample is a sequence of values from some predefined finite set

In , categorical samples are handled using a special class: factor

In this lesson, we are concerned with univariate statistics, that is either with

• numerical samples,

• categorical samples taken from some finite (or seldom countable) set

The boundary between qualitative and quantitative is not always straightforward

flights dataset

• Columns carrier, dest, origin are considered categorical

levels(factor(flights$origin))Copy Code

## [1] "EWR" "JFK" "LGA"Copy Code

• Columns arr_delay, dep_delay, distance are numerical

• Column time_hour, ... contains timestamps

We will learn about:

• numerical summaries

• graphical displays

for univariate samples

Numerical summaries for quantitative variables

• Mean:

$$\overline{X}_n = \frac{1}{n}\sum_{i=1}^n x_i \qquad \textbf{(location)}$$

• Standard deviation:

$$s^2 = \frac{1}{n-1} \sum_{i=1}^n \big(x_i - \overline{X}_n\big)^2 \quad \text{or}\quad s^2 = \frac{1}{n} \sum_{i=1}^n \big(x_i - \overline{X}_n\big)^2\qquad\textbf{(dispersion)}$$

Beyond location et dispersion

• Skewness

$$\frac{\frac{1}{n} \sum_{i=1}^n \big(x_i - \overline{X}_n\big)^3}{s^3}$$

• Kurtosis

$$\kappa = \frac{\frac{1}{n} \sum_{i=1}^n \big(x_i -\overline{X}_n\big)^4}{s^4} - 2$$

Numerical summaries at work

s <- rnorm(1000)   # a random sample from the standard normal distribution
mean(s) - sum(s)/length(s)Copy Code

## [1] 0Copy Code

sd(s) - sqrt(sum((s-mean(s))^2)/(length(s)-1))Copy Code

## [1] 0Copy Code

mean((s- mean(s))^4)/sd(s)^4 -3Copy Code

## [1] 0.06372208Copy Code

Summarizing arr_delay

mean(flights$arr_delay)Copy Code

## [1] NACopy Code

mean(flights$arr_delay, na.rm=T)Copy Code

## [1] 6.895377Copy Code

sd(flights$arr_delay)Copy Code

## [1] NACopy Code

sd(flights$arr_delay, na.rm=T)Copy Code

## [1] 44.63329Copy Code

summary(flights$arr_delay) %>%
  broom::tidy() %>% 
  knitr::kable(digits = 2,
  caption = "Summary statistics for nycflights13 arrival delay")Copy Code

sum_funs <- list(mean, 
  sd, median, IQR)
purrr::map_dbl(sum_funs,
  ~ .x( flights$arr_delay, na.rm=T))Copy Code

## [1]  6.895377 44.633292 -5.000000 31.000000Copy Code

With dplyr::summarise

flights %>% 
  select(- year, -month, -day, -hour, -minute , -time_hour, -ends_with("time"), -flight) %>% 
  summarise(across(where(is.numeric),
                   .fns=list("mean"=mean, "median"=median , "IQR"=IQR),
                   na.rm=T,
                   .names="{.col}_{.fn}")) %>%
  pivot_longer(cols = everything()) %>% 
  knitr::kable(digits=2)Copy Code

Missing values and trimming

The code for mean.default goes beyond computing the weighted mean

function (x, trim = 0, na.rm = FALSE, ...)
{
    if (!is.numeric(x) && !is.complex(x) && !is.logical(x)) {
        warning("argument is not numeric or logical: returning NA")
        return(NA_real_)
    }
    if (na.rm)
        x <- x[!is.na(x)]
    if (!is.numeric(trim) || length(trim) != 1L)
        stop("'trim' must be numeric of length one")
    n <- length(x)
    if (trim > 0 && n) {
        if (is.complex(x))
            stop("trimmed means are not defined for complex data")
        if (anyNA(x))
            return(NA_real_)
        if (trim >= 0.5)
            return(stats::median(x, na.rm = FALSE))
        lo <- floor(n * trim) + 1
        hi <- n + 1 - lo
        x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
    }
    .Internal(mean(x))
}Copy Code

• What happens when na.rm is TRUE?

• What are the possible values of argument trim?

• If trim belongs to $(0, 1/2)$ what happens?

• Which lines actually compute the (uniformly) weighted mean?

• In a tidy sample, missing values are explicitely pointed as missing information using special objects NA, NA_integer_, NA_real_, ...

• To compute summary statistics when facing missing values, na.rm has to be TRUE.

• What is the class of NA, NA_integer_, NA_real_?

Why is var(sample) defined by

$$\sum_{i=1}^n \frac{1}{n-1} (x_i - \overline{X}_n)^2$$

when $\texttt{sample}=x_1, \ldots, x_n$?

• Order statistics: $x_{1:n} \leq x_{2:n} \leq \ldots \leq x_{n:n}$

• Median $x_{\lfloor n/2\rfloor+1:n}$ if $n$ is odd, $(x_{n/2:n}+ x_{n/2+1:n})/2$ if $n$ is even

• Quartiles $x_{\lfloor n/4\rfloor:n}$ and $x_{\lfloor 3 n/4\rfloor:n}$

• Quintiles $x_{\lfloor k n/5\rfloor:n}$ pour $k \in 1, \ldots, 4$

• Deciles $x_{\lfloor k n/10\rfloor:n}$ pour $k \in 1, \ldots, 9$

• Percentiles

• ...

Compute the derivative of

• the mean

• the median

• the trimmed mean

with respect to $x_{(n)}$

Tukey's boxplot assumes symmetry for the whiskers and normality for their length

What does this mean?

Assume that we boxplot the $\mathcal{N}(\mu, \sigma^2) \sim \mu + \sigma \mathcal{N}(0, 1)$ distribution instead of an empirical distribution.

Then:

• The median is equal to $\mu$ because of the symmetry of $\mathcal{N}(0,1)$,

• The quartiles are $\mu - \sigma \Phi^{\leftarrow}(3/4)$ and $\mu + \sigma \Phi^{\leftarrow}(3/4)$.

• The IQR is $2 \sigma \Phi^{\leftarrow}(3/4) \approx 1.35\times \sigma$

• The default length of the whiskers (under geom_boxplot) is $1.5 \times \text{IQR}$

• For a Gaussian population, $1.5 \text{IQR} = 3 \sigma \Phi^{\leftarrow}(3/4)$

• In a Gaussian population, the probability that a sample point lies outside the whiskers is $2 \times (1 - \Phi(4 \times \Phi^{\leftarrow}(3/4))) \approx 7/1000$

Boxplot of a Gaussian sample

If we boxplot a Gaussian sample, on average, $4\%$ of the points lie outside the whiskers. Such points are called extremes

Deciding which points should be considered as extremes is partly a matter of taste

Boxplots provide a visual assessment of the departure of a sample from normality/Gaussianity

geom_boxplot documentation

Boxplot of a non-Gaussian sample

The sampling distribution ( $t$-distribution with $3$ degrees of freedom) is heavy-tailed

Some sample points lie far away from the whiskers

title <- stringr::str_c("Student sample with 3 df, size : ", n)
ggplot(fake_data,
       aes(y=t, x=1L)) +
  geom_boxplot(outlier.shape=2, outlier.colour = 'red', outlier.alpha=.5, coef=1.5) +
  ylab("") +
  ggtitle(title) +
  ggplot2::theme(
  axis.text.x = element_blank(),
  axis.ticks.x = element_blank()
  )Copy Code

summary(fake_data$t)Copy Code

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -5.629301 -0.736732 -0.005991  0.038924  0.841745  5.314943Copy Code

Boxplot of a non-Gaussian sample

Plotting two boxplots side by side

n <- 1000
tibble(normal=rnorm(n), 
       t=rt(n, df=3)) %>%
  pivot_longer(cols=c(normal, t),
               names_to='Distribution') %>%
  ggplot() +
  aes(x=Distribution, y=value) +
  geom_boxplot(outlier.alpha = .5,
               outlier.colour = 'red',
               outlier.shape = 2) +
  ggtitle("Gaussian versus heavy-tailed samples")Copy Code

Boxplots for normal and heavy-tailed distributions

With usual whiskers length

• In Gaussian samples:

  • very few sample points outside the whiskers (few extremes/outliers),
  • extremes lie close to the whiskers endpoints

• In heavy-tailed distributions:

  • many sample points outside the whiskers, and
  • some extremes may lie far far away

Boxplots for delays

flights %>% 
  select(
    ends_with("_delay")
    ) %>% 
  pivot_longer(
    cols=ends_with("_delay"),
    names_pattern = "(.*)_delay", 
    values_to = "delay"
    ) %>% 
  filter(
    abs(delay) < 120
    ) %>% 
  ggplot() +
  aes(x=name, y=delay) +
  geom_boxplot()Copy Code

Construction

• Binning the sample values

  • Picking the right (?) number of bins
  • Placing the breaks

Histograms for (trimmed) delay distributions

flights %>% 
  select(
    ends_with("_delay")
    ) %>% 
  pivot_longer(
    cols=ends_with("_delay"),
    names_pattern = "(.*)_delay", 
    values_to = "delay"
    ) %>% 
  filter(
    abs(delay) < 120
    ) %>% 
  ggplot() +
  aes(x=delay,
      y=..density..,
      fill=name) +
  geom_histogram(alpha=.3,
                 bins=30)Copy Code

A real-valued sample defines a probability distribution, called the empirical distribution. Each sample point is assigned a probability equal to its relative frequency.

The cumulative distribution function associated with this empirical distribution is called the empirical cumulative distribution function (ECDF).

Let $x_1, \ldots, x_n$ denote the points from a real sample (they may not be pairwise distincts), let $P_n$ be the empirical distribution and $F_n$ be the ECDF.

$$\begin{array}{rl} P_n(A) & = \frac{1}{n} \sum_{i=1}^n \mathbb{I}_{x_i \in A} \\ F_n(z) & = \frac{1}{n} \sum_{i=1}^n \mathbb{I}_{x_i \leq z} \end{array}$$

Let $x_{1:n} \leq x_{2:n} \leq \ldots \leq x_{n:n}$ denote the order statistics of the sample (that is the non-decreasing rearrangement of sample, with possible repetitions).

The ECDF is easily computed from the order statistics:

$$F_n(z) = \begin{cases} 0 & \text{if } z < x_{1:n} \\ \frac{k}{n} & \text{if } x_{k:n} \leq z < x_{k+1:n}\\ 1 & \text{if } x_{n:n} \leq z \end{cases}$$

Recall that ECDFs are determined by order statistics (and conversely)

Example

Any numerical dataset may be used to illustrate ECDFs.

For example,

faithful: a data frame with 272 observations on 2 variables concerning waiting time between eruptions and the duration of the eruption for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA.

• eruptions: numeric, Eruption time in minutess

• wating: numeric Waiting time to next eruption (in minutes)

A closer look at faithful$eruptions reveals that these are heavily rounded times originally in seconds, where multiples of 5 are more frequent than expected under non-human measurement.

A glimpse at the faithfull data

data(faithful)
faithful %>%
  glimpse()Copy Code

## Rows: 272
## Columns: 2
## $ eruptions <dbl> 3.600, 1.800, 3.333, 2.283, 4.533, 2.883, 4.700, 3.600, 1.95…
## $ waiting   <dbl> 79, 54, 74, 62, 85, 55, 88, 85, 51, 85, 54, 84, 78, 47, 83, …Copy Code

Plotting an empirical CDF

faithful %>%
  ggplot() +
  aes(x=eruptions) +
  geom_step(stat="ecdf") +
  ylab("ECDF") +
  xlab("Eruptions duration in minutes") +
  ggtitle("Old Faithful Geyser Data")Copy Code

Computing the ECDF statistic

Fn <- ecdf(faithful$eruptions)
x_grid <- seq(min(faithful$eruptions)-1,
              max(faithful$eruptions)+1,
              0.01)
df <- data_frame(x=x_grid,
                 Fn=Fn(x_grid))Copy Code

## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Copy Code

df %>%
  ggplot() +
  aes(x=x, y=Fn) +
  geom_step(direction="hv")Copy Code

Comparing sample quantiles with distributional quantiles

When it comes to comparing the empirical distributions defined by two samples, quantile-quantile plots ( qqplot ) are favored.

• $F_n$ and $G_m$ : two ECDFs

• the quantile-quantile function is defined by

$$G_m^\leftarrow \circ F_n(z) = Y_{k+1:m}\qquad \text{if} \qquad \frac{k}{m}< F_n(z) \leq \frac{k+1}{m}$$

If $G_m = F_n$ (if the two samples are equal up to a permutation), then

$$F_n^\leftarrow \circ F_n(z) = X_{k+1:n} \text{if} \qquad X_{k:n} < z \leq X_{k+1:n}$$

The quantile-quantile plot lies above line $y=x$ and meets this line at every unique order statistics.

The departure of the quantile-quantile function from identity reflects the differences between the two empirical distribution functions.

We can for example compare the first and second half of the faithfull dataset.

Defining an empirical quantile function

equantile <- function(sample){
  sample <- sort(sample)
  m <- length(sample)
  function(p){
    ks <- pmax(1, ceiling(m * p))
    sample[ks]
  }
}Copy Code

]

QQ plot

not_cancelled <- flights %>% 
  filter(!is.na(arr_delay)) %>% 
  sample_n(100) %>% 
  arrange(dep_delay) %>% 
  filter(abs(arr_delay) < 120) 
nr <- nrow(not_cancelled)
not_cancelled %>% 
  ggplot() +
  aes(x=dep_delay, 
      y=quantile(arr_delay, 
                 seq(1, nr)/nr)
      ) +
  geom_point(size=.1) +
  coord_fixed() +
  geom_abline(slope=1, 
              intercept = 0, 
              alpha=.5) +
  ggtitle("QQ plot clipped arr_delay versus dep_delay") +
  ylab("arr_delay") +
  xlab("dep-delay")Copy Code