Difference in means (simulation)

Complete code for class 10/7

Packages

Warning: package 'ggplot2' was built under R version 4.3.3
Warning: package 'tidyr' was built under R version 4.3.3
Warning: package 'readr' was built under R version 4.3.3
Warning: package 'dplyr' was built under R version 4.3.3
Warning: package 'stringr' was built under R version 4.3.3
Warning: package 'lubridate' was built under R version 4.3.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4          ✔ readr     2.1.5     
✔ forcats   1.0.0          ✔ stringr   1.5.1     
✔ ggplot2   3.5.1          ✔ tibble    3.2.1     
✔ lubridate 1.9.3          ✔ tidyr     1.3.1     
✔ purrr     1.0.2.9000     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Warning: package 'tidymodels' was built under R version 4.3.3
── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
✔ broom        1.0.5      ✔ rsample      1.2.1 
✔ dials        1.2.1      ✔ tune         1.2.0 
✔ infer        1.0.7      ✔ workflows    1.1.4 
✔ modeldata    1.3.0      ✔ workflowsets 1.1.0 
✔ parsnip      1.2.1      ✔ yardstick    1.3.1 
✔ recipes      1.0.10     
Warning: package 'dials' was built under R version 4.3.3
Warning: package 'scales' was built under R version 4.3.3
Warning: package 'infer' was built under R version 4.3.3
Warning: package 'modeldata' was built under R version 4.3.3
Warning: package 'parsnip' was built under R version 4.3.3
Warning: package 'recipes' was built under R version 4.3.3
Warning: package 'rsample' was built under R version 4.3.3
Warning: package 'tune' was built under R version 4.3.3
Warning: package 'workflows' was built under R version 4.3.3
Warning: package 'workflowsets' was built under R version 4.3.3
Warning: package 'yardstick' was built under R version 4.3.3
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()
• Dig deeper into tidy modeling with R at https://www.tmwr.org

Attaching package: 'palmerpenguins'

The following object is masked from 'package:modeldata':

    penguins

Hypothesis Testing

set.seed(12345)

null_dist <- penguins |>
  filter(species %in% c("Gentoo", "Chinstrap")) |>
  specify(response = bill_length_mm, explanatory = species) |>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  calculate(stat = "diff in means", order = c("Gentoo", "Chinstrap"))
Dropping unused factor levels Adelie from the supplied explanatory variable
'species'.
Warning: Removed 1 rows containing missing values.

p-value

null_dist |>
  ggplot(
    aes(x = stat)) +
  geom_histogram() +
  labs(title = "simulated distribution") +
  geom_vline(xintercept = -1.2 , color = "red", lwd = 3) +
  geom_vline(xintercept = 1.2 , color = "red", lwd = 3)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

get_p_value(null_dist, -1.2, "two-sided" )
# A tibble: 1 × 1
  p_value
    <dbl>
1   0.016

Confidence interval

set.seed(12345)

boot_df <- penguins |>
  filter(species %in% c("Gentoo", "Chinstrap")) |>
  specify(response = bill_length_mm, explanatory = species) |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "diff in means" , order = c("Gentoo", "Chinstrap"))
Dropping unused factor levels Adelie from the supplied explanatory variable
'species'.
Warning: Removed 1 rows containing missing values.
boot_df |>
  ggplot(
    aes(x = stat) 
  ) + 
  geom_histogram() +
  geom_vline(xintercept = -2.24, color = "red" , lwd = 3) +
  geom_vline(xintercept = -0.329, color = "red", lwd = 3)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

boot_df |>
  summarize(
    lower = quantile(stat, 0.025),
    upper = quantile(stat, 0.975)
  )
# A tibble: 1 × 2
  lower  upper
  <dbl>  <dbl>
1 -2.25 -0.338