Chapter 9 Statistical Foundations infer

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, moderndive, infer)

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 |>
  dplyr::slice_sample(n = 25)
sf_25

Get the summary statistics for the sample taken.

sf_25 |>
  skim(arr_delay)

Data summary

Name	sf_25
Number of rows	25
Number of columns	19
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
arr_delay	0	1	-1.04	48.03	-44	-25	-15	3	182	▇▃▁▁▁

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

SF |>
  skim(arr_delay)

Data summary

Name	SF
Number of rows	13173
Number of columns	19
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
arr_delay	0	1	2.67	47.67	-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 flights 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 |>
  dplyr::slice_sample(n = n) |>
  summarize(mean_arr_delay = mean(arr_delay))

n <- 25
SF |>
  infer::rep_slice_sample(n = n, reps = 1) |>
  summarize(mean_arr_delay = mean(arr_delay))

Note that a different random sample produces a different value for the sample statistic.

SF |>
  infer::rep_slice_sample(n = n, reps = 1) |>
  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.

num_trials <- 500

n <- 25

sf_25_means <- SF |>
  infer::rep_slice_sample(n = n, reps = num_trials) |>
  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

Name	sf_25_means
Number of rows	500
Number of columns	3
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
mean_arr_delay	0	1	2.4	9.21	-17.68	-3.57	0.92	7.34	57.68	▅▇▂▁▁

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
  ) 

Alternatively, it can be calculated directly using a t-test.

sf_25_means |>
  infer::t_test(response = mean_arr_delay)

Using a larger sample size give a smaller Standard Error.

n <- 100

sf_100_means <- SF |> infer::rep_slice_sample(n = n, reps = num_trials) |> 
  group_by(replicate) |> 
  summarize(mean_arr_delay = mean(arr_delay)) |>
  mutate(n = n)

head(sf_100_means)

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 investigate the variability of these samples.

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

sf_200_bs <- orig_sample |> infer::rep_slice_sample(n = n, replace = TRUE, reps = num_trials) |>
  summarize(mean_arr_delay = mean(arr_delay)) |>
  mutate(n = n)

sf_200_bs |>
  skim(mean_arr_delay)

Data summary

Name	sf_200_bs
Number of rows	500
Number of columns	3
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
mean_arr_delay	0	1	1.64	3.19	-8.36	-0.58	1.68	3.63	13.01	▁▅▇▂▁

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

`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Compare with the entire population.

sf_200_pop <- SF |> infer::rep_slice_sample(n = n, replace = TRUE, reps = num_trials)  |>
  summarize(mean_arr_delay = mean(arr_delay)) |>
  mutate(n = n)

sf_200_pop |>
  skim(mean_arr_delay)

Data summary

Name	sf_200_pop
Number of rows	500
Number of columns	3
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
mean_arr_delay	0	1	2.44	3.35	-7.31	0.09	2.43	4.68	13.63	▁▆▇▃▁

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

`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Question: Do the histograms look somewhat similar?

Using the code from the Modern Dive book.

To learn about the infer package this is a good link to start with.

orig_sample |> 
  specify(response = arr_delay) |> 
  calculate(stat = "mean")

bootstrap_distribution <- orig_sample |> 
  specify(formula = arr_delay ~ NULL) |> 
  generate(reps = num_trials, type = "bootstrap") |> 
  calculate(stat = "mean")
bootstrap_distribution

bootstrap_distribution |> visualize()

See the Modern Dive book for the infer process.

percentile_ci <- bootstrap_distribution |> 
  get_confidence_interval(level = 0.95, type = "percentile")
percentile_ci

bootstrap_distribution |> visualize() + 
  shade_confidence_interval(endpoints = percentile_ci)

Example: Setting travel policy

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

n <- nrow(orig_sample)
sf_200_bs <- orig_sample |> infer::rep_slice_sample(n = n, replace = TRUE, reps = num_trials) |> 
  summarize(q98 = quantile(arr_delay, p = 0.98))

sf_200_bs |>
  skim(q98)

Data summary

Name	sf_200_bs
Number of rows	500
Number of columns	2
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
q98	0	1	156.48	50.67	58.08	101.36	183	188.26	267.12	▅▆▅▇▁

n_large <- 10000
sf_10000_bs <- SF |> slice_sample(n = n_large, replace = FALSE)

sf_200_bs <- sf_10000_bs |> infer::rep_slice_sample(n = n, replace = TRUE, reps = num_trials) |> 
  summarize(q98 = quantile(arr_delay, p = 0.98))
  

sf_200_bs |>
  skim(q98)

Data summary

Name	sf_200_bs
Number of rows	500
Number of columns	2
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
q98	0	1	142.96	34.04	70.06	114.42	140.64	166.61	267.04	▃▇▇▁▁

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

`geom_smooth()` using formula = 'y ~ x'

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