Chapter 7 Statistical Foundations

Samples and Populations

The population considered in the book is contained in the nycflights13 dataset. In the nycflights13 dataset contains all of the departing flights from airports in the NYC Area.

library(pacman)
p_load(mdsr, nycflights13, tidyverse, skimr)

Recall the flights dataframe.

flights

Take a sample of flights from NYC to SF. Get all of the flights to SF and take a sample of 25 of them.

Note that I did not set a seed for the sample, so your answers will differ from what is in the book.

SF <- flights %>% 
  filter(dest == "SFO", !is.na(arr_delay))
SF

  
sf_25 <- SF %>%
  slice_sample(n = 25)

Get the summary statistics for the sample taken.

sf_25 %>%
  skim(arr_delay)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             25        
Number of columns          19        
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 arr_delay             0             1  47.4  94.5   -60    -8    16    70   317 ▇▅▁▁▁

Since the SF dataset contains all flights from NYC to SF in 2013, the statistics computed from the SF dataset are the calculated population paramters.

SF %>%
  skim(arr_delay)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             13173     
Number of columns          19        
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 arr_delay             0             1  2.67  47.7   -86   -23    -8    12  1007 ▇▁▁▁▁

Get the 98th percentile of the sample.

sf_25 %>%
  summarize(q98 = quantile(arr_delay, p = 0.98))

Get the estimated proportion of filights to SF with a delay less than 90 minutes.

SF %>%
  group_by(arr_delay < 90) %>%
  count() %>%
  mutate(pct = n / nrow(SF))

Compare with the 98th percentile. 90 minutes is the 95th percentile.

SF %>%
  summarize(q98 = quantile(arr_delay, p = 0.98))

Sample Statistics

The sampling distribution.

Usually sampling is done without replacement.

n <- 25
SF %>%
  slice_sample(n = n) %>%
  summarize(mean_arr_delay = mean(arr_delay))

Note that different random sample produce different values for the sample statistic.

SF %>%
  slice_sample(n = n) %>%
  summarize(mean_arr_delay = mean(arr_delay))

Using simulation we can approximate the sampling distribution of the sample statistics. Here we are computing the mean, but we could do this for any statistic we are interested in, the median, sample variance, sample standard deviation, etc.

Note that the do() function is from the mosaic package. We could just use a for loop.

num_trials <- 500
sf_25_means <- 1:num_trials %>%
  map_dfr(
    ~ SF %>%
      slice_sample(n = n) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)
head(sf_25_means)

Summarize the sample means.

sf_25_means %>%
  skim(mean_arr_delay)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             500       
Number of columns          2         
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable  n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 mean_arr_delay         0             1  2.62  9.61 -20.5 -4.59   1.3  8.75  57.3 ▃▇▃▁▁

Confidence intervals

sf_25_means %>%
  summarize(
    x_bar = mean(mean_arr_delay),
    se = sd(mean_arr_delay)
  ) %>%
  mutate(
    ci_lower = x_bar - 2 * se, # approximately 95% of observations 
    ci_upper = x_bar + 2 * se  # are within two standard errors
  ) 

Using a larger sample size give a smaller Standard Error.

n <- 100
sf_100_means <- 1:500 %>%
  map_dfr(
    ~ SF %>%
      slice_sample(n = n) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)

Plots to compare the sampling distributions with different sample sizes.

sf_25_means %>%
  bind_rows(sf_100_means) %>%
  ggplot(aes(x = mean_arr_delay)) + 
  geom_histogram(bins = 30) + 
  facet_grid( ~ n) + 
  xlab("Sample mean")

Boostrap

Usually for random sampling we sample without replacement.

three_flights <- SF %>%
  slice_sample(n = 3, replace = FALSE) %>%
  select(year, month, day, dep_time)
three_flights

With the Bootstrap we use sampling with replacement

three_flights %>% slice_sample(n = 3, replace = TRUE)

n <- 200
orig_sample <- SF %>% 
  slice_sample(n = n, replace = FALSE)
orig_sample

The reason for sampling with replacement is that we only have the Original Sample of size \(n\). We do not have the population. So if we sampled with out replacement we would only have one sample to look at. Sampling the Original Sample with replace allow us to general lots of samples and to investiage the variability of these samples.

orig_sample %>%
  slice_sample(n = n, replace = TRUE) %>%
  summarize(mean_arr_delay = mean(arr_delay))

sf_200_bs <- 1:num_trials %>%
  map_dfr(
    ~orig_sample %>%
      slice_sample(n = n, replace = TRUE) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)
sf_200_bs %>%
  skim(mean_arr_delay)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             500       
Number of columns          2         
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable  n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 mean_arr_delay         0             1  2.99  3.10 -5.32  1.06  2.84  4.82  14.2 ▁▇▇▂▁

sf_200_bs %>% ggplot(aes(x = mean_arr_delay)) +
  geom_histogram()

sf_200_pop <- 1:num_trials %>%
  map_dfr(
    ~SF %>%
      slice_sample(n = n, replace = TRUE) %>%
      summarize(mean_arr_delay = mean(arr_delay))
  ) %>%
  mutate(n = n)
sf_200_pop %>%
  skim(mean_arr_delay)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             500       
Number of columns          2         
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable  n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 mean_arr_delay         0             1  2.54  3.42 -7.16 0.232  2.52  4.77  11.5 ▁▅▇▆▁

sf_200_pop %>% ggplot(aes(x = mean_arr_delay)) +
  geom_histogram()

Question: Do the histograms look somewhat similar?

orig_sample %>%
  summarize(q98 = quantile(arr_delay, p = 0.98))

n <- nrow(orig_sample)
sf_200_bs <- 1:num_trials %>%
  map_dfr(
    ~orig_sample %>%
      slice_sample(n = n, replace = TRUE) %>%
      summarize(q98 = quantile(arr_delay, p = 0.98))
  )
sf_200_bs %>%
  skim(q98)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             500       
Number of columns          1         
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 q98                   0             1  148.  40.1  55.0  121.   140  162.   310 ▁▇▅▁▁

n_large <- 10000
sf_10000_bs <- SF %>% 
  slice_sample(n = n_large, replace = FALSE)
sf_200_bs <- 1:num_trials %>%
  map_dfr(~sf_10000_bs %>%
        slice_sample(n = n_large, replace = TRUE) %>%
        summarize(q98 = quantile(arr_delay, p = 0.98))
  )
sf_200_bs %>%
  skim(q98)

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             500       
Number of columns          1         
_______________________              
Column type frequency:               
  numeric                  1         
________________________             
Group variables                      

── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable n_missing complete_rate  mean    sd    p0   p25   p50   p75  p100 hist 
1 q98                   0             1  152.  4.91   137   151   152  154.  167. ▁▁▇▂▁

outliers

SF %>%
  filter(arr_delay >= 420) %>% 
  select(month, day, dep_delay, arr_delay, carrier)

SF %>% 
  filter(arr_delay < 420) %>%
  ggplot(aes(arr_delay)) + 
  geom_histogram(binwidth = 15) + 
  labs(x = "Arrival delay (in minutes)")

SF %>%
  group_by(hour) %>%
  count() %>%
  pivot_wider(names_from = hour, values_from = n) %>%
  data.frame()

SF %>%
  ggplot(aes(x = hour, y = arr_delay)) +
  geom_boxplot(alpha = 0.1, aes(group = hour)) +
  geom_smooth(method = "lm") + 
  xlab("Scheduled hour of departure") + 
  ylab("Arrival delay (minutes)") + 
  coord_cartesian(ylim = c(-30, 120)) 

mod1 <- lm(arr_delay ~ hour, data = SF)
broom::tidy(mod1)