functional_programming_learning.Rmd

---
title: "Functional Programming Learning"
date: '`r format(Sys.Date(), "%Y-%m-%d")`'
output:
  html_document:
    code_folding: hide
    theme: flat
    toc: yes
    toc_depth: 5
    toc_float:
      collapsed: no
  pdf_document:
    toc: yes
    toc_depth: '5'
editor_options:
  chunk_output_type: console
---

# Library
```{r}
library(tidyverse)
install.packages("palmerpenguins")
library(palmerpenguins)


```

# attributes
```{r}
attributes(penguins[1:20, ])

attr(penguins, "class")
attr(penguins, "names")

attributes(penguins$year)
```

# dimensions
```{r}
m <- matrix(1:6, ncol = 2)
m
attributes(m)

# modify the attributes
attr(m, "dim") <- c(2, 3)

# create an atomic vector
v <- 1:6

# assign a dimension attribute
attr(v, "dim") <- c(3, 2)

dim_names <- list(
  c("the first", "second", "III"),
  c("index", "value")
)

attr(v, "dimnames") <- dim_names

attributes(v)
rownames(v) <- c("the first", "second", "III")
colnames(v) <- c("index", "value")

attr(v, "matrix_mean") <- mean(v)
```

# A function that uses attributes
```{r}
db <- list(
  data.frame(
    color = c("red", "orange", "green"),
    transparency = c(0.8, 0.65, 0.93)
  ),
  data.frame(
    color = c("blue", "pink", "cyan"),
    transparency = c(0.40, 0.35, 0.87)
  )
)

# write a function to grab table; if it's 1920, grab the first one, otherwise grab the second one
pull_color_data <- function(year){
  to_pull <- if(year == "1920"){
    out <- db[[1]]
  }else{
    out <- db[[2]]
  }
  out
}

pull_color_data(1921)
pull_color_data(1920)

# modify the function to store the year as an attribute
pull_color_data <- function(year){
  to_pull <- if(year == "1920"){
    out <- db[[1]]
  }else{
    out <- db[[2]]
  }
  attributes(out) <- c(
    attributes(out),
    db = year
  )
  out
}

pull_color_data(1920)
attributes(pull_color_data(1920))
attr(pull_color_data(2021), "db")

library(colorspace)
print_colors <- function(color_data){
  title <- paste0("Colors for ", attr(color_data, "db"))
  colorspace::swatchplot(color_data$color)
  mtext(title)
}

pull_color_data(2021) |> 
  print_colors()
```

# Another example of attributes
```{r}
library(lme4)
# fit a multilevel model and pull the variance-covariance matrix
m <- lmer(Reaction ~ 1 + Days + (1 + Days|Subject),
          data = sleepstudy)

VarCorr(m)$Subject
```

# Matrices
```{r}
set.seed(42)
m <- matrix(rnorm(100, 200, 10), ncol = 10)
sum(m)
m
rowSums(m)
colSums(m)

# standardize the matrix
z <- (m - mean(m))/sd(m)


```

# names
```{r}
names(v) <- letters[1:6]

v3a <- c(a = 5, b = 7, c = 12)
v3a
typeof(v3a)
class(v3a)
names(v3a)
attributes(v3a)

v["b"]
v["d"]
```

# factors
```{r}
fct <- factor(c("a", "a", "b", "c"))
typeof(fct)
attributes(fct)
str(fct)


# create integer vector
int <- c(1L, 1L, 2L, 3L, 1L, 3L)
# assign some levels
attr(int, "levels") <- c("red", "green", "blue")
# change the class to a factor
class(int) <- "factor"

age <- factor(sample(c("baby", 1:10), 100, replace = TRUE))
age
typeof(age)
class(age)
str(age)

data.frame(age) |> 
  count(age)
```  

# Classes
```{r}
summary(mtcars)
summary(gss_cat$marital)
summary(lm(mpg ~ cyl, mtcars))

summary.data.frame(mtcars)
```

# lists
```{r}
l <- list(
  c("a", "b", "c"),
  rnorm(5),
  c(7L, 2L),
  c(TRUE, TRUE, FALSE, TRUE)
)

typeof(l)
attributes(l)
str(l)

l_df <- list(
  a = c("red", "blue"),
  b = rnorm(2),
  c = c(7L, 2L),
  d = c(TRUE, FALSE)
)
data.frame(l_df)
```

```{r}
x <- c(a = 3, b = 5, c = 7)

l <- list(
  x = x,
  x2 = c(x, x),
  x3 = list(
    vect = x,
    squared = x^2,
    cubed = x^3)
)

# subsetting lists
l[1]
typeof(l[1])

l[[1]]
typeof(l[[1]])
l[[1]]["c"]

# named list
l$x2
l$x3$squared["a"]

l[c(TRUE, FALSE, TRUE)]
l[c(1, 3)]

mtcars[, 4, drop = FALSE]
mtcars_tbl <- as_tibble(mtcars)
mtcars_tbl[, 4]


```

# for loops
```{r}
a <- letters[1:26]

for (i in 1:5){
  print(a[i])
}

# simulate tossing a coin multiple times
# pre-allocate a vector
result <- rep(NA, 10)
result
# run the trial n times, storing the result in your pre-allocated vector
for (i in seq_along(result)){
  result[i] <- sample(c("Heads", "Tails"), 1)
}
result

# create a alphabet of upper/lower case
result <- rep(NA, length(letters))

for (i in seq_along(result)){
  result[i] <- paste0(LETTERS[i], letters[i])
}
result

# simulate 100 cases from random normal distribution,
# where we varied the standard deviation in increments of 0.2, ranging from 1 to 5
increments <- seq(1, 5, by = 0.2)
# allocate a vector; do it in a list
simulated <- vector("list", length(increments))

for (i in seq_along(simulated)){
  simulated[[i]] <- rnorm(100, 0, increments[i])
}

# provide meaningful names
names(simulated) <- paste0("sd_", increments)
str(simulated)

# convert to df
sim_d <- data.frame(simulated)

# plot
pd <- sim_d |> 
  pivot_longer(
    everything(),
    names_to = "sd",
    values_to = "sim",
    names_prefix = "sd_",
    names_ptypes = list(sd = factor())# change the type of sd
  )
pd |> 
  ggplot(aes(sim)) +
  geom_density(aes(color = sd)) +
  guides(color = "none")

# calculate all the densities
densities <- vector("list", length(sim_d))

for (i in seq_along(densities)){
  densities[[i]] <- density(sim_d[, 1])
}
str(densities)

# plot the first density
plot(densities[[1]])

# loop through all the other densities
plot(densities[[1]], xlim = c(-20, 20))

for (i in seq(2, length(densities))){
  lines(x = densities[[i]]$x,
        y = densities[[i]]$y)
}

# breaking loops
set.seed(1)

rand_unif <- vector("double", 10)

for (i in seq_along(rand_unif)) {
  rand_unif[i] <- runif(1, 0, 10)
  if(any(rand_unif > 5)){
    break
  }
}
```

# lapply 
```{r}
# the for loop version
increments <- seq(1, 5, by = 0.2)
# allocate a vector; do it in a list
simulated <- vector("list", length(increments))

for (i in seq_along(simulated)){
  simulated[[i]] <- rnorm(10, 0, increments[i])
}
simulated

# the lapply version
increments <- seq(1, 5, by = 0.2)
sim_l <- lapply(increments, function(sd) rnorm(10, 0, sd))
sim_l

# loop through a data frame because a data frame is a list
lapply(palmerpenguins::penguins, is.double)
lapply(mtcars, mean)

# add a condition
lapply(penguins, function(x){
  if(is.numeric(x)){
    mean(x, na.rm = TRUE)
  }
})

# add a second condition
lapply(penguins, function(x){
  if(is.numeric(x)){
    return(mean(x, na.rm = TRUE))
  }else if(is.character(x) | is.factor(x)){
    return(table(x))
  }
})

# passing arguments
lapply(airquality, mean, na.rm = TRUE)

lapply(seq(1, 5, 0.2), rnorm, n = 10, mean = 0)
```

# operations by group
```{r}

# replicate group_by
by_cyl <- split(mtcars, mtcars$cyl)
# by_cyl is a list of data frames
lapply(by_cyl, function(x) mean(x$mpg))

# another example
by_species <- split(penguins, penguins$species)
lapply(by_species, function(x) mean(x$bill_length_mm, na.rm = TRUE))
```

# produce separate plots
```{r}
lapply(by_cyl, function(x){
  ggplot(x, aes(disp, mpg)) +
    geom_point() +
    geom_smooth()
})

lapply(by_species, function(x){
  ggplot(x, aes(bill_length_mm, body_mass_g)) +
    geom_point() +
    geom_smooth()
})

# save plots to an object
penguin_plots <- lapply(by_species, function(x){
  ggplot(x, aes(bill_length_mm, body_mass_g)) +
    geom_point() +
    geom_smooth()
})

str(penguin_plots) # this is a list

# save file names/directory
install.packages("here")
library(here)

library(fs)
# create a plots folder and a penguins folder inside the plots folder
dir_create(here("plots", "penguins"))
names(penguin_plots)
files <- paste0(names(penguin_plots), ".pdf")
file_paths <- here("plots", "penguins", files)

# write a for loop to save these plots
for (plot in seq_along(file_paths)){
  ggsave(file_paths[plot], # single bracket
         penguin_plots[[plot]], # double bracket
         width = 6.5,
         height = 8)
}
```

```{r}
# create a plot of bill length vs. bill depth for each combination of species and island
split2 <- split(penguins, list(penguins$species, penguins$island))
str(split2)

plots2 <- lapply(split2, function(x){
  ggplot(x, aes(bill_length_mm, bill_depth_mm)) +
    geom_point()
})

# get rid of empty plots

# check if n row > 0
keep <- sapply(split2, function(x) nrow(x) > 0)
# use this to subset the list
split3 <- split2[keep]
plots3 <- lapply(split3, function(x){
  ggplot(x, aes(bill_length_mm, bill_depth_mm)) +
    geom_point()
})

# OR 
sapply(split2, nrow)
sapply(split2, function(x){
  if(nrow(x) != 0){
    ggplot(x, aes(bill_length_mm, bill_depth_mm)) +
    geom_point()
  }
})

```

# sapply
```{r}
!sapply(penguins, is.double)

penguins[, !sapply(penguins, is.double)]

```

# vapply
```{r}
vapply(mtcars, mean, FUN.VALUE = double(1))
vapply(mtcars, is.double, FUN.VALUE = logical(1))

# count missing data
vapply(airquality, function(col){
  sum(is.na(col))
  },
  FUN.VALUE = double(1))
```

# apply
apply(x, dimension, function, function_args)
x: the thing to loop over (usually a matrix or data frame)
dimension: 1 = loop over each row; 2 = loop over each column; n = apply to nth dimension of an array (rare)

```{r}
people <- data.frame(
  first = c("Frederick", "Anna", "Julia"),
  last = c("Douglass", "Murray", "Griffiths")
)
people

# want a new column that was the first and last name
apply(people, 1, paste, collapse = " ")

people |> 
  mutate(full_name = apply(people, 1, paste, collapse =))
```

# example
```{r}
# example: calculate top mpg manufactures
four_cyl <- filter(mpg, cyl == 4)

ninety <- four_cyl %>%
    filter(cty >= quantile(cty, probs = 0.9)) # get the 90%tile
ninety
count_mfr <- count(ninety, manufacturer)
count_model <- count(ninety, manufacturer, model)
count_class <- count(ninety, class)

# produce a plot for each
plot_mfr <- count_mfr |> 
  ggplot(aes(fct_reorder(manufacturer, -n), n)) +
  geom_col(aes(fill = manufacturer)) +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "Manufacturers", x = "", y = "") +
  guides(fill = "none")

plot_car <- count_model %>%
    unite(car, manufacturer, model, sep = " ") %>%
    ggplot(aes(fct_reorder(car, -n), n)) +
        geom_col(aes(fill = car)) +
        scale_fill_brewer(palette = "Set3") +
        labs(title = "Top 10% of city mpg",
             subtitle = "Car frequency",
             x = "",
             y = "") +
        guides(fill = "none")

plot_class <- 
  ggplot(count_class, aes(fct_reorder(class, -n), n)) +
        geom_col(aes(fill = class)) +
        scale_fill_brewer(palette = "Set3") +
        labs(title = "Car Class",
             x = "", 
             y = "") +
        guides(fill = "none")

# Assemble the plots
library(patchwork)
plot_car / (plot_mfr + plot_class)

# functional programming version

by_cyl <- split(mpg, mpg$cyl)
top_10 <- lapply(by_cyl, function(x){
  filter(x, cty > quantile(cty, probs = 0.9))
})
top_10
str(top_10)
# keep <- vapply(top_10, function(x) nrow(x) > 0, FUN.VALUE = logical(1))
# top_10_new <- top_10[keep]

# all counts
counts <- lapply(top_10, function(x){
  count_manufacturer <- count(x, manufacturer)
  count_class <- count(x, class)
  count_model <- count(x, manufacturer, model) |> 
    unite(car, manufacturer, model, sep = " ")
  
  return(list(mfr = count_manufacturer,
              car = count_model,
              class = count_class))
})
names(counts[[1]])[1]
# write the plot function
counts_plot <- function(counts_df){
  var <- names(counts_df)[1]
  
  p <- ggplot(counts_df, aes(fct_reorder(!!sym(var), -n), n)) +
    geom_col(aes(fill = !!sym(var))) +
    scale_fill_brewer(palette = "Set3") +
    labs(title = str_to_title(var),
         x = "",
         y = "") +
    guides(fill = "none")
  
  if (var == "car") {
    p <- p + labs(title = "Top 10% of city mpg", subtitle = var)
  }
  p
}

counts_plot(counts[["6"]]$car)

# compile plots function
full_plot <- function(l){
  counts_plot(l[["car"]]) / (
    counts_plot(l[["mfr"]]) +
      counts_plot(l[["class"]])
  )
}

full_plot(counts[[1]])

# finish
plots <- lapply(counts, full_plot)
```

# purrr
```{r}
map(1:3, rnorm)
lapply(1:3, rnorm)

map(mtcars, function(x) length(unique(x)))
lapply(mtcars, function(x) length(unique(x)))

map(mtcars, ~ length(unique(.x))) # ~ is used in place of function(.x)

# vary the n
map(1:3, ~ rnorm(n = .x))

# vary the mean
map(1:3, ~ rnorm(n = 1, mean = .x))

# vary the sd
map(1:3, ~ rnorm(n = 1, sd = .x))

# use map to extract elements from list
l <- list(
  list(-1, x = 1, y = 2, z = "a"),
  list(-2, x = 4, y = c(5, 6), z = "b"),
  list(-3, x = 8, y = c(9, 10, 11))
)

map(l, 2)
lapply(l, function(x) x[[2]])
lapply(l, `[[`, 2)

# extract by name
map(l, "y")
# extract the first element of y
map(l, list("y", 1))

# variants of map
map_dbl(l, "x")
map_chr(l, "z", .default = NA_character_)
```

```{r}
df_list <- list(
  data.frame(var1 = 1:5),
  data.frame(var1 = 1:3),
  data.frame(var1 = 1)
)

# calculate the mean
map(df_list, function(x) mean(x$var1))
map(df_list, ~ mean(.x$var1))
map(df_list, ~ mean(.x[["var1"]]))

# return a vector
map_dbl(df_list, ~ mean(.x$var1))
```

# another example
```{r}
# fit a simple model to each year
# the relationship between personal consumption expenditures (pce) and unemployment percentage over time?
econ <- economics |> 
  mutate(year = year(date))

# split by year
by_year <- split(econ, econ$year)

unemploy <- map(by_year, ~ .x |> mutate(unemploy_pct = (unemploy / pop) * 100))

# fit the model
lm_model <- map(unemploy, ~ lm(unemploy_pct ~ pce, data = .x))

# extract the coefficient

coef_df <- map(lm_model, ~ coef(.x))
coefs <- map(lm_model, coef)

# extract the slope
coefs[[1]][2]

slopes <- map_dbl(coefs, 2)

# put all together by piping
by_year |> 
  map(~ mutate(.x, percent = unemploy / pop)) |> 
  map(~ lm(percent ~ pce, data = .x)) |> 
  map(coef) |> 
  map_dbl(2)
```

# lab 1
```{r}
set.seed(123) # Set random number generation

m <- lm(mpg ~ hp + cyl + drat, mtcars)

l <- list(a = list(m = matrix(1:12, ncol = 3),
              	   v = 1:7),
          b = data.frame(student = 1:15,
                         score = rnorm(15, 100, 10)))
```

```{r}
# extract the qr tolerance
m$qr$tol
View(m)
# extract the term labels from the model
str(m)
attr(m$terms, "term.labels")

attr(attr(m$model, "terms"), "term.labels")

# from l
## Use at least two different methods to extract m.
l$a$m
l[["a"]][["m"]]
## Extract the third column of m. Maintain the dimensions (matrix structure).
l$a$m[, 3, drop = FALSE]
## Extract the score for student 7 in the data frame.
l$b[l$b$student == 7, ]
l$b[7, ]
```

```{r}
ins <- read.csv("insurance_coverage.csv")
head(ins)

# Split the data file by "age_bucket" and "sex_name"
by_age_sex <- split(ins, list(ins$age_bucket, ins$sex_name))
str(by_age_sex)

# Use a for loop to calculate the mean health care population (hc_pop) for each "age_bucket"/"sex_name" combination.
hc_pop_mean_list <- rep(NA, length(by_age_sex))
for (i in seq_along(by_age_sex)){
  hc_pop_mean_list[i] <- mean(by_age_sex[[i]]$hc_pop, na.rm = TRUE)
}

# Replicate this calculation with lapply, sapply, and vapply 
lapply(by_age_sex, function(x) mean(x$hc_pop, na.rm = TRUE))
sapply(by_age_sex, function(x) mean(x$hc_pop, na.rm = TRUE))
vapply(by_age_sex, function(x) mean(x$hc_pop, na.rm = TRUE), FUN.VALUE = double(1))

# Produce separate plots showing the change in hc_pop from 2014 to 2015 for each "age_bucket"/"sex_name" combination (with a single loop). Set the hc_pop axis so the limits equal the overall minimum and maximum values for hc_pop.
 

ins_plots <- lapply(by_age_sex, function(x){
  ggplot(x, aes(x = factor(year), y = hc_pop, group = 1)) +
    geom_point() +
    geom_line() +
    scale_y_continuous(limits = c(min(ins$hc_pop), max(ins$hc_pop)))
})

# Use a for loop to save the plots to a folder on your computer
# save file names/directory
install.packages("here")
library(here)

library(fs)
# create a plots folder and a penguins folder inside the plots folder
dir_create(here("plots", "insurance_coverate"))
names(ins_plots)
files <- paste0(names(ins_plots), ".pdf")
file_paths <- here("plots", "insurance_coverate", files)


# write a for loop to save these plots
for (plot in seq_along(file_paths)){
  ggsave(file_paths[plot], # single bracket
         ins_plots[[plot]], # double bracket
         width = 6.5,
         height = 8)
}
```

# map_dfr
```{r}
# create a function that simulates data
set.seed(8765309)
simulate <- function(n, mean = 0, sd = 1) {
  tibble(sample_id = seq_len(n),
         sample = rnorm(n, mean, sd))
}

simulate(10)
simulate(4, mean = -10, sd = 1.5)

# if we want to vary the sample size from 10 to 150 by increments of 5
# mean stays constant at 100, sd is constant at 10

sims <- map(seq(10, 150, 5), simulate, mean = 100, sd = 10)
sims_df <- map_dfr(seq(10, 150, 5), simulate, mean = 100, sd = 10,
                   .id = "iteration")

sims_df[1:15,]

sample_size <- seq(10, 150, 5)
install.packages("english")
library(english)

sample_size <- setNames(sample_size, english::english(seq(10, 150, 5)))
sample_size

sims_df3 <- map_dfr(sample_size, simulate, mean = 100, sd = 10,
                    .id = "n")

# broom
lm(tvhours ~ age, gss_cat) |> 
  broom::tidy()

split(gss_cat, gss_cat$year) |> 
  map_dfr(~lm(tvhours ~ age, .x) |> 
            tidy(),
          .id = "year")
```

# batch-loading data
useful when you have data separated by some categories
```{r}
library(fs) # file system package
dir_ls(here::here("data"))

# we only want .csv
dir_ls(here("data"), glob = "*.csv")

# loop through the directories and read in the csv files
files <- dir_ls(
  here("data"),
  glob = "*.csv"
)

batch <- map_dfr(files, read_csv, .id = fi)

batch2 <- map_dfr(files, read_csv, .id = "file")

# remove the here::here patch from string
batch2 <- batch2 |> 
  mutate(file = str_replace_all(
    file, here::here("data"), ""
  ))

# use regular expressions to pull out pieces you need
batch2 |> 
  count(file)

## pull grade; but some file names have 2 digits after g, some have 1 digit
## then pull year
# then pull content area
d <- batch2 |> 
  # use () to specify what thing to pull out
  mutate(
  grade = str_replace_all(
    file,
    "/g(\\d?\\d).+",
    "\\1"
  ),
  grade = as.integer(grade),
  
  year = str_replace_all(
    file,
    ".+files(\\d\\d)_sim.+",
    "\\1"
  ),
  year = as.integer(year),
  
  content = str_replace_all(
    file,
    "/g\\d?\\d(.+)pfiles.+",
    "\\1"
  )
  ) |> 
  select(-file) |> 
  select(ssid, grade, year, content, testeventid, 
      asmtprmrydsbltycd, asmtscndrydsbltycd, Entry:WMLE)

d |> 
  count(content)

```

```{r}
# what if we want only math files
# we can't use glob because that's only for extension
dir_ls(here::here("data"), regexp = "Math") # use regexp to pass a pattern you want to search for

# read in only the grade 5 data
file_g5 <- dir_ls(here::here("data"), regexp = "g5")

g5_df <- map_dfr(file_g5, read_csv, id = "file")

# base equivalents for ls
list.files(here("data"))
list.files(here("data"), full.names = TRUE, pattern = "*.csv")

# fs package plays nicer with purrr
```

# list column
```{r}
# want to compare three models
d <- d |> 
  mutate(primary = as.factor(asmtprmrydsbltycd),
         secondary = as.factor(asmtscndrydsbltycd))

# create a data frame with a list column
by_content <- d |> 
  group_by(content) |> 
  nest()

map_dbl(by_content$data, nrow)
map_dbl(by_content$data, ncol)
map_dbl(by_content$data, ~ mean(.x$Theta))

by_content |> 
  mutate(n = map_dbl(data, nrow))

# fit a model
a <- by_content |> 
  mutate(m1 = map(data, ~ lm(Theta ~ primary, data = .x)),
         coefs = map(m1, coef))

by_content |> 
  mutate(m1 = map(data, ~ lm(Theta ~ primary, data = .x)),
         coefs = map(m1, coef),
         tidy_coefs = map(coefs, tidy)) |> 
  unnest(tidy_coefs) |> 
  filter(names %in% c("(Intercept)", "primary74")) |> 
  select(content, names, x) |> 
  pivot_wider(
    values_from = x,
    names_from = names
  ) |> 
  mutate(TBI = primary74 + `(Intercept)`) |> 
  select(content, intercept = 2, TBI)

# OR
by_content |> 
  mutate(m1 = map(data, ~ lm(Theta ~ primary, data = .x)),
         coefs = map(m1, coef),
         no_disab = map_dbl(coefs, 1),
         tbi = no_disab + map_dbl(coefs, "primary74")) |> 
  select(content, no_disab, tbi)

# fit all three models
mods <- by_content |> 
  mutate(
    m1 = map(data, ~ lm(Theta ~ primary, data = .x)),
    m2 = map(data, ~ lm(Theta ~ primary + secondary, data = .x)),
    m3 = map(data, ~ lm(Theta ~ primary * secondary, data = .x)))

# pull the p value for the model comparison
mods |> 
  mutate(comp12 = map2(m1, m2, anova),
         p12 = map_dbl(comp12, list("Pr(>F)", 2)))

a <- mods |> 
  mutate(comp12 = map2(m1, m2, anova))
a$comp12[[1]]$`Pr(>F)`[2]
a$comp12[[1]][["Pr(>F)"]][2]
# write a function that pull the p-value from model comparison objects
extract_p <- function(anova_ob) {
  anova_ob[["Pr(>F)"]][2]
}

mods |> 
  mutate(comp12 = map2(m1, m2, anova),
         p12 = map_dbl(comp12, extract_p))

# create a function using compose()
center <- compose(~ .x - mean(.x, na.rm = TRUE))

penguins$bill_length_mm |> 
  center()

# compose a p-val extractor
p <- compose(~ .x[["Pr(>F)"]][2])
mods |> 
  mutate(comp12 = map2(m1, m2, anova),
         p12 = map_dbl(comp12, p))
```

```{r}
# operations by row
df <- tibble(
  name = c("Me", "You"),
  x = 1:2,
  y = 3:4,
  z = 5:6
)

df |> 
  mutate(m = mean(c(x, y, z)))

df |> 
  rowwise() |> 
  mutate(m = mean(c(x, y, z)))

# if you apply rowwise operation with a list column, you don't have to loop
df <- tibble(var = list(1, 2:3, 4:6))
df |> 
  mutate(lngth = map_int(var, length))
df |> 
  rowwise() |> 
  mutate(lngth = length(var))

d |> 
  nest_by(content) |> 
  mutate(m1 = list(lm(Theta ~ primary, data = data)),
         m2 = list(lm(Theta ~ primary + secondary, data = data)),
         m3 = list(lm(Theta ~ primary * secondary, data = data)))

```

```{r}
# extract all r-squared
r2 <- mods |> 
  pivot_longer(
    m1:m3,
    names_to = "model",
    values_to = "output"
  ) |> 
  mutate(r2 = map_dbl(output, ~ summary(.x)$r.squared))

# tidy the coefficients
tidied <- mods |> 
  pivot_longer(
    m1:m3,
    names_to = "model",
    values_to = "output"
  ) |> 
  mutate(tidied = map(output, tidy)) |> 
  select(content, model, tidied) |> 
  unnest(tidied)

to_plot <- names(coef(mods$m1[[1]]))  
  
```

# map2
Simultate data from a normal distribution
- vary n from 5 to 150 by increments of 5
- For each n, vary miu from -2 to 2 by increments of 0.25

expand.grid()
```{r}
conditions <- expand.grid(n = seq(5, 150, 5),
            mu = seq(-2, 2, 0.25))

# simulate
sim1 <- map2(conditions$n, conditions$mu, ~{
  rnorm(n = .x, mean = .y, sd = 10)
})

# add it as a list column
sim2 <- conditions |> 
  as_tibble() |> 
  mutate(sim = map2(n, mu, ~ rnorm(n = .x, mean = .y, sd = 10))) |> 
  unnest(sim)


# replicate with rowwise() approach
conditions |> 
  rowwise() |> 
  mutate(sim = list(rnorm(n = n, mean = mu, sd = 10))) |> 
  unnest(sim)
```

```{r}
install.packages("fivethirtyeight")
library(fivethirtyeight)
pulitzer

# prep data
pulitzer <- pulitzer |> 
  select(newspaper, starts_with("num")) |> 
  pivot_longer(
    -newspaper,
    names_to = "year_range",
    values_to = "n",
    names_prefix = "num_finals"
  ) |> 
  mutate(year_range = str_replace_all(year_range, "_", "-")) |> 
  filter(year_range != "1990-2014")

# one plot
wsj <- pulitzer |> 
  filter(newspaper == "Wall Street Journal")

ggplot(wsj, aes(n, year_range)) +
  geom_col(aes(fill = n)) +
  scale_fill_distiller(
    type = "seq",
    limits = c(0, max(pulitzer$n)),
    palette = "BuPu",
    direction = 1
  ) +
  scale_x_continuous(
    limits = c(0, max(pulitzer$n)),
    expand = c(0, 0)
  ) +
  guides(fill = "none") +
  labs(
    title = "Pulitzer Prize winners: Wall Street Journal",
    x = "Total number of winners",
    y = ""
  )

# nest data
by_newspaper <- pulitzer |> 
  group_by(newspaper) |> 
  nest()

library(glue)
p <- by_newspaper |> 
  mutate(plot = map2(
    data, newspaper, ~{
      ggplot(.x, aes(n, year_range)) +
        geom_col(aes(fill = n)) +
        scale_fill_distiller(
          type = "seq",
          limits = c(0, max(pulitzer$n)),
          palette = "BuPu",
          direction = 1
        ) +
        scale_x_continuous(
          limits = c(0, max(pulitzer$n)),
          expand = c(0, 0)
        ) +
        guides(fill = "none") +
        labs(
          title = glue("Pulitzer Prize winners: {.y}"),
          x = "Total number of winners",
          y = ""
        )
    }
  ))

## using rowwise approach
pulitzer %>%
nest_by(newspaper) %>%
    mutate(
      plot = list(
     ggplot(data, aes(n, year_range)) +
        geom_col(aes(fill = n)) +
        scale_fill_distiller(
          type = "seq",
          limits = c(0, max(pulitzer$n)),
          palette = "BuPu",
          direction = 1
        ) +
        scale_x_continuous(
          limits = c(0, max(pulitzer$n)), 
          expand = c(0, 0)
        ) +
        guides(fill = "none") +
        labs(
          title = glue("Pulitzer Prize winners: {newspaper}"),
          x = "Total number of winners",
          y = ""
        )
    )
  )

## add a caption to the plot
pulitzer <- pulitzer |> 
  group_by(newspaper) |> 
  mutate(tot = sum(n))

# create a column to represent exactly the label you want
pulitzer <- pulitzer |> 
  mutate(
    label = glue(
      "{str_to_title(as.english(tot))} Total Pulitzer Awards"
    )
  )

## produce one plot to test
wsj2 <- pulitzer |> 
  filter(newspaper == "Wall Street Journal")

ggplot(wsj2, aes(n, year_range)) +
  geom_col(aes(fill = n)) +
  scale_fill_distiller(
    type = "seq",
    limit = c(0, max(pulitzer$n)),
    palette = "BuPu",
    direction = 1
  ) +
  scale_x_continuous(
    limits = c(0, max(pulitzer$n)),
    expand = c(0, 0)
  ) +
  guides(fill = "none") +
  labs(
    title = glue("Pulitzer Prize winners: Wall Street Journal"),
    x = "Total number of winners",
    y = "",
    caption = unique(wsj2$label)
  )

# produce all plots
by_newspaper_label <- pulitzer |> 
  group_by(newspaper, label) |> 
  nest()

final_plots <- by_newspaper_label |> 
  mutate(plots = pmap(list(newspaper, label, data), ~{
    ggplot(..3, aes(n, year_range)) +
      geom_col(aes(fill = n)) +
      scale_fill_distiller(
        type = "seq",
        limits = c(0, max(pulitzer$n)),
        palette = "BuPu",
        direction = 1
      ) +
      scale_x_continuous(
        limits = c(0, max(pulitzer$n)),
        expand = c(0, 0)
      ) +
      guides(fill = "none") +
      labs(
        title = glue("Pulitzer Prize winners: {..1}"),
        x = "Total number of winners",
        y = "",
        caption = ..2
      )
  }))

## the rowwise approach
pulitzer %>%
  ungroup() |> 
  nest_by(newspaper, label) |> 
  mutate(plot = list(
    ggplot(data, aes(n, year_range)) +
      geom_col(aes(fill = n)) +
      scale_fill_distiller(
        type = "seq",
        limits = c(0, max(pulitzer$n)),
        palette = "BuPu",
        direction = 1
      ) +
      scale_x_continuous(
        limits = c(0, max(pulitzer$n)),
        expand = c(0, 0)
      ) +
      guides(fill = "none") +
      labs(title = glue("Pulitzer Prize winners :{newspaper}"),
           x = "Total number of winners",
           y = "",
           caption = label)
  ))

## save all plots
# create a directory
fs::dir_create(here::here("plots", "pulitzers"))

# create file paths
files <- str_replace_all(
  tolower(final_plots$newspaper),
  " ",
  "-"
)
paths <- here::here("plots", "pulitzers", glue("{files}.png"))

# add paths to data frame
final_plots |> 
  ungroup() |> 
  mutate(path = paths) |> 
  rowwise() |> 
  summarize(
    ggsave(path,
           plots,
           width = 9.5,
           height = 6.5,
           dpi = 500)
  )

```

#pmap
simulation data from a normal distribution
vary n from 5 to 150 by increments of 5
for each n, vary miu from -2 to 2 by increments of 0.25
for each sigma from 1 to 3 by increments of 0.1
```{r}
#simulation
full_conditions <- expand.grid(
  n = seq(5, 150, 5),
  mu = seq(-2, 2, 0.25),
  sd = seq(1, 3, 0.1)
)

fsim <- pmap(
  list(number = full_conditions$n,
       average = full_conditions$mu,
       stdev = full_conditions$sd
),
function(number, average, stdev) { # the argument names here need to match the list names above
  rnorm(n = number, mean = average, sd = stdev)
}
)
## alternative
fsim <- pmap(
  list(full_conditions$n,
       full_conditions$mu,
       full_conditions$sd),
  ~ rnorm(n = ..1, mean = ..2, sd = ..3)
)

## another alternative
fsim <- pmap(
  full_conditions,
  ~ rnorm(n = ..1, mean = ..2, sd = ..3)
)

## list column version
full_conditions |> 
  as_tibble() |> 
  mutate(sim = pmap(list(n, mu, sd), ~ rnorm(..1, ..2, ..3))) |> 
  unnest(sim)

## rowwise approach
full_conditions |> 
  rowwise() |> 
  mutate(sim = list(rnorm(n, mu, sd))) |> 
  unnest(sim)

```

# walk

```{r}
library(tidyverse)
library(english)
library(fs)
library(here)
library(glue)
library(broom)
library(fivethirtyeight)

# save plots

dir_create(here("plots", "pulitzers"))
newspapers <- str_replace_all(
  tolower(p$newspaper),
  " ",
  "-"
)
paths <- here(
  "plots",
  "pulitzers",
  glue("{newspapers}.png")
)

walk2(paths, p$plot, ggsave, width = 6.5, height = 9)
```

# modify()
modify() always returns the same type as the input type
```{r}
modify(mtcars, ~ as.numeric(scale(.x)))
map_dfr(mtcars, ~ as.numeric(scale(.x)))

modify2(LETTERS[1:3], letters[1:3], paste0)
map2_chr(LETTERS[1:3], letters[1:3], paste0)
```

# safely()
makes your iterations safer
```{r}
by_cyl <- mpg |> 
  group_by(cyl) |> 
  nest()
by_cyl |> 
  mutate(mod = map(data, ~ lm(hwy ~ displ + drv, data = .x)))

# define safe function
safe_lm <- safely(lm)

# loop the safe function
safe_models <- by_cyl |> 
  mutate(safe_mod = map(data, ~ safe_lm(hwy ~ displ + drv, data = .x)))

# create a new variable to filter for results with errors
safe_models |> 
  mutate(error = map_lgl(safe_mod, ~!is.null(.x$error))) |> 
  # inspect the data that have the error message
  filter(isTRUE(error)) |> 
  select(cyl, data) |> 
  unnest(data)

# pull results that worked
safe_models |> 
  mutate(results = map(safe_mod, "result"),
         tidied = map(results, tidy)) |> 
  select(cyl, tidied) |> 
  unnest(tidied)
```

Use safely in web scraping because URLs don't always work
```{r}
library(rvest)
links <- list(
  "https://en.wikipedia.org/wiki/FC_Barcelona",
  "https://nosuchpage",
  "https://en.wikipedia.org/wiki/Rome"
)

pages <- map(links, ~{
  Sys.sleep(0.1)
  read_html(.x)
})

safe_read_html <- safely(read_html)

pages <- map(links, ~{
  Sys.sleep(0.1)
  safe_read_html(.x)
})

# for non-results, we want to double check the pages where we got no results
errors <- map_lgl(pages, ~ !is.null(.x$error))
links[errors]
```

# reduce()
When you use reduce(), you will always have the result of length one
```{r}
l <- list(
  c(1, 3),
  c(1, 5, 7, 9),
  3,
  c(4, 8, 12, 2)
)

reduce(l, sum)
map(l, sum)

# the use case of reduce()
l_df <- list(
  tibble(id = 1:3, score = rnorm(3)),
  tibble(id = 1:5, treatment = rbinom(5, 1, .5)),
  tibble(id = c(1, 3, 5, 7), other_thing = rnorm(4))
)

# apply the join function on all three tibbles in the list
reduce(l_df, full_join)

# another example to use reduce and bind_rows
l_df2 <- list(
  tibble(id = 1:3, scid = 1, score = rnorm(3)),
  tibble(id = 1:5, scid = 2, score = rnorm(5)),
  tibble(id = c(1, 3, 5, 7), scid = 3, score = rnorm(4))
)

reduce(l_df2, bind_rows)
```

# function
Anything that carries out an operation in R is a function
Functions are also usually objects

```{r}
3 + 5
`+` (3, 5)

```

We can have one thing that we want to apply multiple functions to
```{r}
funs <- list(
  quarter = function(x) x/ 4,
  half = function(x) x/ 2,
  double = function(x) x * 2,
  quadruple = function(x) x * 4
)
funs$quarter(100)

# store functions in a list
smry <- list(
  n = length,
  n_miss = function(x) sum(is.na(x)),
  n_valid = function(x) sum(!is.na(x)),
  mean = function(x) mean(x, na.rm = TRUE),
  sd = function(x) sd(x, na.rm = TRUE)
)
# for one column
map_df(smry, ~.x(mtcars$mpg))

# for multiple columns
map_df(mtcars, function(col) map_df(smry, ~.x(col)), .id = "column")

# avoid nested loops
# write a function
summarize_col <- function(column) {
  map_df(smry, ~.x(column))
}

map_df(mtcars, summarize_col, .id = "col")

# wrap the whole thing in a function
summarize_df <- function(df) {
  map_df(df, summarize_col, .id = "column")
}

summarize_df(airquality)
```

```{r}
# pull out specific coefficients
mods <- mtcars |> 
  group_by(cyl) |> 
  nest() |> 
  mutate(
    model = map(
      data, ~ lm(mpg ~ disp + hp + drat, data = .x)
    )
  )
# pull just the first model
m <- mods$model[[1]]
coef(m)

coef(m)["disp"]

pull_coef <- function(model, coef_name){
  coef(model)[coef_name]
}

mods %>%
    mutate(intercept = map_dbl(model, pull_coef, "(Intercept)"),
           disp      = map_dbl(model, pull_coef, "disp"),
           hp        = map_dbl(model, pull_coef, "hp"),
           drat      = map_dbl(model, pull_coef, "drat"))

# make some default
pull_coef <- function(model, coef_name = "(Intercept)"){
  coef(model)[coef_name]
}

mods |> 
  mutate(intercept = map_dbl(model, pull_coef))

# return all coefficients
# check the coefficient names of one model
names(coef(m))
# pull them in a data frame
tibble(term = names(coef(m)),
       coefficient = coef(m))


pull_all_coefficient <- function(model){
  tibble(term = names(coef(model)),
         coefficient = coef(model))
}


mods |> 
  mutate(coefs = map(model, pull_all_coefficient)) |> 
  select(cyl, coefs) |> 
  unnest(coefs)

## another example
# rescale each column to 0/1 (min is 0 and max is 1)
set.seed(42)
df <- tibble::tibble(
  a = rnorm(10, 100, 150),
  b = rnorm(10, 100, 150),
  c = rnorm(10, 100, 150),
  d = rnorm(10, 100, 150)
)

# subtract the min value from each observation, dividing that
# by the difference between the min/max value

# for example,
tibble(v1 = c(3, 4, 5),
       numerator = v1 - 3,
       denominator = 5 - 3,
       scaled = numerator / denominator)

# for one column
df |> 
  mutate(
    a = (a - min(a, na.rm = TRUE)) /
      (max(a, na.rm = TRUE) - min(a, na.rm = TRUE))
  )

# if we don't want to write a function
map_df(df, ~ (.x - min(.x, na.rm = TRUE)) /
         (max(.x, na.rm = TRUE) - min(.x, na.rm = TRUE)))

# write a function
rescale_function <- function(col){
  (col - min(col, na.rm = TRUE)) /
    (max(col, na.rm = TRUE) - min(col, na.rm = TRUE))
}

# make it easy to read
rescale_function <- function(col){
  z <- na.omit(col)
  min_z <- min(z)
  max_z <- max(z)
  
  (z - min_z) / (max_z - min_z)
}

df |> 
  map_df(rescale_function)

# add conditions to let a function behave differently depending on the input

# calculate mean of all numeric columns
mean2 <- function(x, ...) {
  if(is.numeric(x)) {
    out <- mean(x, ...)
  } else {
    return()
  }
  out
}
# calculate all means
means_mpg <- map(mpg, mean2)
# drop nulls
is_null <- map_lgl(means_mpg, is.null)
means_mpg[!is_null]

as.data.frame(means_mpg[!is_null])

get_numeric_mean <- function(df, ...) {
  means_df <- map(df, mean2, ...)
  is_null <- map_lgl(means_df, is.null)
  as.data.frame(means_df[!is_null])
}

get_numeric_mean(airquality, na.rm = TRUE)
```

# purity
```{r}
# returns the total number of instances where the value is NA for both vectors
df <- tibble(
  a = c(1, NA, NA, 3, 3, 9, NA),
b = c(NA, 3, NA, 4, NA, NA, NA)
)

df |> 
  rowwise() |> 
  mutate(both_null = ifelse(is.na(a)&is.na(b), TRUE, FALSE)) |> 
  count(both_null) |> 
  filter(both_null == TRUE)

is.na(a)
is.na(b)

sum(is.na(a) & is.na(b))

both_na <- function(x, y) {
  sum(is.na(x) & is.na(y))
}

# check length, stop the function if the length condition has not been met
both_na <- function(x, y) {
  if(length(x) != length(y)) {
    stop("The lengths for the two vectors are not the same.")
  }
  sum(is.na(x) & is.na(y))
}

both_na(a, c(b, b))

# more meaningful error message?
# state the length of each
both_na <- function(x, y) {
  if(length(x) != length(y)) {
    v_lngths <- paste0("x = ", length(x),
                       ", y = ", length(y))
    stop("Vectors are of different lengths:", v_lngths)
  }
  sum(is.na(x) & is.na(y))
}

# we don't need echo the function; set call. = FALSE
# we can also state the lengths on new lines
both_na <- function(x, y) {
  if(length(x) != length(y)) {
    v_lngths <- paste0("x = ", length(x),
                       ", y = ", length(y))
    stop("Vectors are of different lengths:\n", 
         v_lngths,
         call. = FALSE)
  }
  sum(is.na(x) & is.na(y))
}

# for quick checks, use stopifnot
z_score <- function(x) {
  stopifnot(is.numeric(x))
  x <- x[!is.na(x)]
  (x - mean(x)) / sd(x)
}

z_score(c("a", "b", "c"))

# if you want to embed a warning, just swap out stop() for warning()
both_na <- function(x, y) {
  if(length(x) != length(y)) {
    v_lngths <- paste0("x = ", length(x),
                       ", y = ", length(y))
    warning("Vectors are of different lengths:", v_lngths)
  }
  sum(is.na(x) & is.na(y))
}


```

```{r}
# let's build out a warning if the vectors are different lengths, but they are recyclable
# %% modula operator

## one approach; note the double condition
both_na <- function(x, y) {
  if(length(x) != length(y)) {
    lx <- length(x)
    ly <- length(y)
    
    v_lngths <- paste0("x = ", ly, ", y = ", ly)
    
    if(lx %% ly == 0 | ly %% lx == 0){
      warning("Vectors were recycled (", v_lngths, ")")
    }
    else {
      stop("Vectors are of different lengths and are not recyclable:",
           v_lngths)
    }
  }
  sum(is.na(x) & is.na(y))
}

both_na(a, c(b, b))
```

# Write better functions
```{r}

# calculate if recyclable
recyclable <- function(x, y) {
  test1 <- length(x) %% length(x)
  test2 <- length(y) %% length(x)
  
  any(c(test1, test2) == 0)
}

# test the recyclable function
recyclable(a, b)

# make errors/warning a function
check_lengths <- function(x, y) {
  if(length(x) != length(y)) {
    if(recyclable(x, y)) {
      warning(
        "Vectors were recycled\n",
        "x = ", length(x), "\n",
        "y = ", length(y),
        call. = FALSE
      )
    }
    else {
      stop(
        "Vectors are of different lengths and are not recyclable:\n",
        "x = ", length(x), "\n",
        "y = ", length(y),
        call. = FALSE
      )
    }
  }
}

both_na <- function(x, y){
  check_lengths(x, y)
  sum(is.na(x) & is.na(y))
}

both_na(c(a, a), c(b, b, b))
```

# non-standard evaluation
```{r}
percentile_df <- function(x) {
  sorted <- sort(x)
  d <- data.frame(sorted,
                  percentile = ecdf(sorted)(sorted))
  names(d)[1] <- paste0(substitute(x), collapse = "_")
  d
}

percentile_df(rnorm(100, 5, 0.2)) |> 
  head()

quote(subset(df, select = var))

substitute(subset(df, select = var))

extract_var <- function(df, var) {
  subset(eval(substitute(df)),
         select = var)
}
extract_var(mtcars, "mpg")
```

```{r}
# use NSE for both arguments
extract_var <- function(df, var){
  eval(substitute(var), envir = df)
}
extract_var(mtcars, mpg)

extract_var <- function(df, var){
  df[, as.character(substitute(var)), drop = FALSE]
}
extract_var(mtcars, mpg)

# for n columns
extract_vars <- function(df, ...){
  vars <- substitute(alist(...))
  df[, as.character(vars)[-1]]
}

head(extract_vars(mtcars, mpg, cyl, disp))
```

```{r}
   
mtcars %>%
  group_by(cyl, gear) %>%
  summarize(mean = mean(mpg, na.rm = TRUE)) %>%
  pivot_wider(names_from = cyl, values_from = mean)

# {rlang} version
group_means <- function(data, outcome, group_1, group_2) {
  out <- enquo(outcome) # quote the inputs
  g1 <- enquo(group_1)
  g2 <- enquo(group_2)
  
  data |> 
    group_by(!!g1, !!g2) |> # evaluate 
    summarize(mean = mean(!!out, na.rm = TRUE)) |> 
    pivot_wider(names_from = !!g1, values_from = mean)
}

group_means(mtcars, mpg, cyl, gear)

# double curly approach -- quote and evaluate at the same spot
group_means3 <- function(data, outcome, group_1, group_2) {
  
  data |> 
    group_by({{group_1}}, {{group_2}}) |> # evaluate 
    summarize(mean = mean({{outcome}}, na.rm = TRUE)) |> 
    pivot_wider(names_from = {{group_1}}, values_from = mean)
}
group_means3(mtcars, mpg, cyl, gear)

```


```{r}
# write a function that summarize any numeric columns by returning the mean, standard deviation, min, and max 
diamonds |> 
  select(depth, table, price) |> 
  pivot_longer(everything()) |> 
  group_by(name) |> 
  summarize(mean = mean(value),
            sd = sd(value))

summarize_cols <- function(data, ...){
  data |> 
    select(...) |> 
    pivot_longer(everything(),
                 names_to = "var",
                 values_to = "val") |> 
    group_by(var) |> 
    summarize(mean = mean(val, na.rm = TRUE),
            sd = sd(val, na.rm = TRUE),
            min = min(val, na.rm = TRUE),
            max = max(val, na.rm = TRUE))
}

mtcars |> 
  summarize_cols(mpg, hp, drat)

library(palmerpenguins)

penguins |> 
  select_if(is.numeric) |> 
  summarize_cols(everything())

```

```{r}
# we want to check if an x/y relation is linear by showing the plots
check_linear <- function(data, x, y, points = FALSE, ...) {
  p <- ggplot(data, aes({{x}}, {{y}}))
  if (points) {
    p <- p + geom_point(color = "gray80")
  }
  p +
    geom_smooth(method = "lm", color = "magenta", ...) +
    geom_smooth(...)
}

check_linear(penguins, bill_length_mm, bill_depth_mm)

# create a function that calculates means and standard errors for every numeric column in a data frame
# then use this function to create a plotting function

# SE function
se <- function(x) {
  x <- x[!is.na(x)] # remove missing data
  sd(x) / sqrt(length(x))
}

# solve it on an individual dataset first
penguins |> 
  select_if(is.numeric) |> 
  pivot_longer(everything(),
               names_to = "var",
               values_to = "val") |> 
  group_by(var) |> 
  summarize(mean = mean(val, na.rm = TRUE),
            se = se(val))

# generalize into a function
estimate_means <- function(data){
  data |> 
    select_if(is.numeric) |> 
    pivot_longer(everything(),
                 names_to = "var",
                 values_to = "val") |> 
    group_by(var) |> 
    summarize(mean = mean(val, na.rm = TRUE),
              se = se(val))
}


gss_cat |> 
  estimate_means()

# plot
estimate_means(gss_cat) |> 
  ggplot(aes(mean, se)) +
  geom_point() +
  geom_linerange(
    aes(xmin = mean + (qnorm(0.025) * se),
        xmax = mean + (qnorm(0.975) *se)),
    color = "magenta"
  )

# create plot function
plot_means <- function(df) {
  means <- estimate_means(df) |> 
    mutate(var = reorder(factor(var), mean))
  
  ggplot(means, aes(mean, var)) +
    geom_point() +
    geom_linerange(
      aes(xmin = mean + (qnorm(0.025) * se),
          xmax = mean + (qnorm(0.975) * se)),
      color = "magenta"
    )
}

plot_means(penguins)
```