The tinyverse collection of packages makes things a lot more sane, IMO:

  penguins <- read_csv("penguins.csv") |>
    na.omit() |>
    select(species, island, bill_length_mm, body_mass_g) |>
    group_by(species, island) |>
    summarize(
      mean_bill_length = mean(bill_length_mm),
      mean_mass = mean(body_mass_g),
      n = n()
    ) |>
    arrange(species, desc(mean_bill_length))
  
  penguins |>
    ggplot(aes(x = species, y = mean_bill_length, fill = island)) +
    geom_col(position = "dodge") +
    labs(
      title = "Mean Bill Length by Species and Island",
      y = "Mean Bill Length (mm)"
    ) +
    theme_minimal()