## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)

## ----data---------------------------------------------------------------------
measurement_names <- c(
  "bill_length_mm",
  "bill_depth_mm",
  "flipper_length_mm",
  "body_mass_g"
)

penguins <- palmerpenguins::penguins
complete_rows <- stats::complete.cases(penguins[, measurement_names])

x <- as.matrix(penguins[complete_rows, measurement_names])

c(
  observations = nrow(x),
  measurements = ncol(x),
  omitted_rows = sum(!complete_rows)
)

## ----correlations-------------------------------------------------------------
round(stats::cor(x), 2)

## ----pca----------------------------------------------------------------------
penguin_pca <- stats::prcomp(
  x,
  center = TRUE,
  scale. = TRUE
)

eigenvalues <- penguin_pca$sdev^2
variance_explained <- 100 * eigenvalues / sum(eigenvalues)

pca_summary <- data.frame(
  component = paste0("PC", seq_along(eigenvalues)),
  eigenvalue = round(eigenvalues, 3),
  variance_percent = round(variance_explained, 1),
  cumulative_percent = round(cumsum(variance_explained), 1)
)

knitr::kable(pca_summary)

## ----fit----------------------------------------------------------------------
set.seed(2026)

penguin_sig <- sigPCA::sigPCA(
  x,
  method = "both",
  num_permutations = 999
)

penguin_sig$mp$mp_bounds
penguin_sig$combined

## ----results-table------------------------------------------------------------
component_results <- data.frame(
  component = paste0("PC", seq_along(eigenvalues)),
  eigenvalue = round(eigenvalues, 3),
  variance_percent = round(variance_explained, 1),
  beyond_mp_upper = seq_along(eigenvalues) %in%
    penguin_sig$mp$significant_components,
  permutation_p = round(penguin_sig$perm$pvalues, 4)
)

knitr::kable(component_results)

## ----spectrum, fig.cap="Observed eigenvalues and the Marchenko--Pastur bounds. Components above the upper bound are highlighted."----
sigPCA::plot_sigPCA(
  penguin_sig$mp$eigenvalues,
  penguin_sig$mp$mp_bounds
)

## ----spectrum-histogram, fig.cap="Empirical eigenvalue distribution with the Marchenko--Pastur bounds. Rug marks identify eigenvalues above the upper bound."----
sigPCA::plot_sigPCA_histogram(
  penguin_sig,
  bins = ncol(x)
)

## ----loadings, fig.cap="Absolute loadings for the retained first component. Fill indicates the loading sign."----
loading_data <- data.frame(
  measurement = rownames(penguin_pca$rotation),
  loading = penguin_pca$rotation[, 1]
)

ggplot2::ggplot(loading_data) +
  ggplot2::aes(
    x = abs(loading),
    y = stats::reorder(measurement, abs(loading)),
    fill = loading > 0
  ) +
  ggplot2::geom_col() +
  ggplot2::scale_fill_manual(
    values = c("#b6dfe2", "#0A537D"),
    labels = c("Negative", "Positive")
  ) +
  ggplot2::labs(
    x = "Absolute loading",
    y = NULL,
    fill = "Sign"
  ) +
  ggplot2::theme_minimal()

## ----biplot, fig.cap="Scores and loading vectors for the first two components. Species labels are used only for interpretation."----
score_data <- data.frame(
  PC1 = penguin_pca$x[, 1],
  PC2 = penguin_pca$x[, 2],
  species = penguins$species[complete_rows]
)

loading_vectors <- data.frame(
  measurement = rownames(penguin_pca$rotation),
  PC1 = penguin_pca$rotation[, 1],
  PC2 = penguin_pca$rotation[, 2]
)

score_span <- vapply(
  score_data[, c("PC1", "PC2")],
  function(z) diff(range(z)),
  numeric(1)
)

loading_span <- vapply(
  loading_vectors[, c("PC1", "PC2")],
  function(z) diff(range(z)),
  numeric(1)
)

arrow_multiplier <- 0.35 * min(score_span / loading_span)
loading_vectors$PC1 <- loading_vectors$PC1 * arrow_multiplier
loading_vectors$PC2 <- loading_vectors$PC2 * arrow_multiplier

penguin_colours <- c(
  Adelie = "darkorange",
  Chinstrap = "purple",
  Gentoo = "cyan4"
)

ggplot2::ggplot(score_data) +
  ggplot2::aes(x = PC1, y = PC2, colour = species, shape = species) +
  ggplot2::geom_point(alpha = 0.75, size = 2) +
  ggplot2::geom_segment(
    data = loading_vectors,
    ggplot2::aes(x = 0, y = 0, xend = PC1, yend = PC2),
    inherit.aes = FALSE,
    colour = "#0A537D",
    arrow = grid::arrow(length = grid::unit(0.08, "inches"))
  ) +
  ggplot2::geom_text(
    data = loading_vectors,
    ggplot2::aes(x = PC1, y = PC2, label = measurement),
    inherit.aes = FALSE,
    colour = "#0A537D",
    hjust = 0,
    vjust = 1,
    size = 3.5,
    check_overlap = TRUE
  ) +
  ggplot2::scale_colour_manual(values = penguin_colours) +
  ggplot2::labs(
    x = sprintf("PC1 (%.1f%% of variance)", variance_explained[1]),
    y = sprintf("PC2 (%.1f%% of variance)", variance_explained[2]),
    colour = "Species",
    shape = "Species"
  ) +
  ggplot2::theme_minimal()