library(ggplot2)           # ggplot2 for creating visualizations, 
library(dplyr)             # dplyr for data manipulation, 
library(palmerpenguins)    # palmerpenguins for the dataset


focal_penguins <- penguins %>%
  filter(species %in% c("Adelie", "Chinstrap")) %>%   # Keep only Adelie and Chinstrap penguins
  filter(!is.na(flipper_length_mm + body_mass_g)) %>% # Remove rows with missing values in flipper length or body mass
  select(species, flipper_length_mm, body_mass_g)     # Select only the species, flipper length, and body mass columns


Ztrans <- function(x) {
  # Z-transform function: standardizes the values by subtracting the mean and dividing by the standard deviation
  (x - mean(x)) / sd(x)
}

z_focal_penguins <- focal_penguins %>%
  group_by(species) %>%            # Group the data by species
  mutate_all(Ztrans)               # Apply Z transformation to all columns (flipper length and body mass)



gg_color_hue <- function(n) {
  # Function to generate a sequence of colors based on the number of groups (n)
  # This uses a color space that ensures distinct and attractive colors for each group
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}


ggplot(focal_penguins, aes(x = body_mass_g, y = flipper_length_mm, color = species)) +
  geom_point() +                                   # Create a scatter plot with body mass on x-axis and flipper length on y-axis, colored by species
  geom_smooth(method = 'lm') +                     # Add a linear regression line (method = 'lm') for each species
  scale_color_manual(values = gg_color_hue(3)[1:2], # Use custom colors for the species, reusing the color function defined earlier
                     guide = guide_legend(reverse = TRUE)) + # Reverse the order of species in the legend
  geom_label(data = . %>%
    group_by(species) %>%                          # For each species, compute the slope of the relationship between body mass and flipper length
    summarise(slope = cov(body_mass_g, flipper_length_mm) / var(body_mass_g),  # Calculate slope using covariance/variance
              slope_label = paste(unique(species), "\nslope:", round(slope, digits = 4)),  # Create label text with the species and slope value
              body_mass_g = 5200,                  # Set a fixed body mass for label placement (x-coordinate)
              flipper_length_mm = max(c(as.numeric(species == "Chinstrap") * 209, 197))), # Set y-coordinate for label placement (higher for Chinstrap)
             aes(label = slope_label), size = 6) + # Add the slope label to the plot with a specific font size
  coord_cartesian(xlim = c(2500, 5500)) +          # Set the x-axis limits to zoom in on a specific range of body mass values
  theme(legend.position = "none")                  # Remove the legend from the plot