library(dplyr)
library(tidyr)
data(iris)
# numeric variables only
num <- iris |> dplyr::select(where(is.numeric))
# per-variable summary (pooled)
num |>
  summarise(across(everything(), list(
    N        = ~sum(!is.na(.)),
    Mean     = ~mean(., na.rm = TRUE),
    Median   = ~median(., na.rm = TRUE),
    SD       = ~sd(., na.rm = TRUE),
    Variance = ~var(., na.rm = TRUE),
    Min      = ~min(., na.rm = TRUE),
    Max      = ~max(., na.rm = TRUE),
    Range    = ~diff(range(., na.rm = TRUE)),
    Q1       = ~quantile(., 0.25, na.rm = TRUE),
    Q3       = ~quantile(., 0.75, na.rm = TRUE),
    IQR      = ~IQR(., na.rm = TRUE),
    Skewness = ~moments::skewness(., na.rm = TRUE),
    Kurtosis = ~moments::kurtosis(., na.rm = TRUE)
  ), .names = "{.col}__{.fn}")) |>
  pivot_longer(everything(),
               names_to  = c("Variable", "Statistic"),
               names_sep = "__",
               values_to = "Value") |>
  pivot_wider(names_from = Statistic, values_from = Value)

# group-wise summary (e.g., by Species) uses the same block with:
# iris |> group_by(Species) |> group_modify(~ summarise_stats(.x))

library(nortest)
data(iris)
num <- iris[sapply(iris, is.numeric)]
# Shapiro-Wilk (3 <= N <= 5000)
lapply(num, shapiro.test)
# Anderson-Darling (N > 7)
lapply(num, nortest::ad.test)

data(iris)
num <- iris[sapply(iris, is.numeric)]
cor(num,
    use    = "pairwise.complete.obs",
    method = "pearson")   # or "spearman", "kendall"

data(iris)
# first variable categorical (Species); second variable continuous ->
# discretised to Low / Medium / High at mean +/- 1 SD
classify_num <- function(z) {
  m <- mean(z, na.rm = TRUE)
  s <- sd(z, na.rm = TRUE)
  cut(z,
      breaks = c(-Inf, m - s, m + s, Inf),
      labels = c("Low", "Medium", "High"),
      right  = TRUE)
}
tbl <- table(iris$Species, classify_num(iris$Sepal.Length))
prop.table(tbl, 1) * 100   # row percentages
prop.table(tbl, 2) * 100   # column percentages
chisq.test(tbl)