Descriptive Statistics in R: A Practical Approach

Introduction

In our previous post, we covered the theory behind descriptive statistics — the formulas, notation, and intuition. Here, we put that theory into practice, computing each measure in R, using base R functions along with a few custom functions for statistics.

The analysis separates variables into numeric and categorical types, since each requires different summary techniques:

For numeric variables, the following are calculated:

Central tendency — mean, median, mode
Spread/dispersion — range, variance (sample & population), standard deviation (sample & population), interquartile range (IQR), coefficient of variation
Position & shape — percentiles, skewness, kurtosis (raw and excess), and IQR-based outlier detection
Bivariate relationships — covariance, Pearson correlation, and Spearman rank correlation

For categorical variables, the following are calculated:

Frequency analysis — absolute frequency, relative frequency, percentage frequency, and cumulative frequency
Contingency tables — cross-tabulations between categorical pairs, tested for association using the Chi-square test and quantified using Cramér’s V

R code:

Variable & Parameter declaration

# Descriptive Statistics # CONFIG DATA_PATH <- "/dataset.csv" DATE_COL <- "date_col" NUMERIC_COLS <- "num_col" CATEGORICAL_COLS <- "cat_col" # NULL = auto-generate every pair from CATEGORICAL_COLS or set like: list(c("season","weathersit"), c("yr","workingday")) CATEGORICAL_PAIRS <- NULL PERCENTILE_PROBS <- c(0.10, 0.25, 0.50, 0.75) # percentile positions IQR_OUTLIER_MULTIPLIER <- 1.5 SAMPLE_DDOF <- 1 # sample : divide by (n-1) POPULATION_DDOF <- 0 # population: divide by n ROUND_DECIMALS <- 2

# Descriptive Statistics 
# CONFIG 

DATA_PATH            <- "/dataset.csv"
DATE_COL             <- "date_col"
NUMERIC_COLS         <- "num_col"
CATEGORICAL_COLS     <- "cat_col"

# NULL  = auto-generate every pair from CATEGORICAL_COLS or set like: list(c("season","weathersit"), c("yr","workingday"))
CATEGORICAL_PAIRS    <- NULL

PERCENTILE_PROBS     <- c(0.10, 0.25, 0.50, 0.75) # percentile positions

IQR_OUTLIER_MULTIPLIER <- 1.5

SAMPLE_DDOF          <- 1    # sample  : divide by (n-1)
POPULATION_DDOF      <- 0    # population: divide by n

ROUND_DECIMALS       <- 2

Code block

Show code

# ============================================================================= # 0. LOAD DATA + METADATA # ============================================================================= df <- read.csv(DATA_PATH, stringsAsFactors = FALSE) df[[DATE_COL]] <- as.Date(df[[DATE_COL]]) cat("0. METADATA\n") cat(sprintf("Shape (rows, columns): %d x %d\n", nrow(df), ncol(df))) cat(sprintf("Date range: %s to %s\n", min(df[[DATE_COL]]), max(df[[DATE_COL]]))) meta <- data.frame( dtype = sapply(df, class), non_null = sapply(df, function(x) sum(!is.na(x))), missing_count = sapply(df, function(x) sum(is.na(x))), missing_pct = round(sapply(df, function(x) mean(is.na(x)) * 100), ROUND_DECIMALS), n_unique = sapply(df, function(x) length(unique(x))) ) cat("\nColumn metadata:\n") print(meta) cat("\nNumeric columns :", paste(NUMERIC_COLS, collapse = ", "), "\n") cat("Categorical cols :", paste(CATEGORICAL_COLS, collapse = ", "), "\n") # ============================================================================= # 1. NUMERIC DATA # ============================================================================= # Helper: population std (R's sd() uses n-1 by default) pop_sd <- function(x) sqrt(sum((x - mean(x))^2) / length(x)) pop_var <- function(x) sum((x - mean(x))^2) / length(x) cv <- function(x) (sd(x) / mean(x)) * 100 # sd() = sample std # ── Central Tendency ──────────────────────────────────────────────────── cat(" CENTRAL TENDENCY (Mean, Median, Mode)\n") first_mode <- function(x) as.numeric(names(sort(table(x), decreasing = TRUE)[1])) central <- data.frame( mean = round(sapply(df[NUMERIC_COLS], mean), ROUND_DECIMALS), median = round(sapply(df[NUMERIC_COLS], median), ROUND_DECIMALS), mode = round(sapply(df[NUMERIC_COLS], first_mode), ROUND_DECIMALS) ) print(central) # ── Spread ─────────────────────────────────────────────────────────────── cat(" SPREAD / DISPERSION\n") spread <- data.frame( minimum = round(sapply(df[NUMERIC_COLS], min), ROUND_DECIMALS), maximum = round(sapply(df[NUMERIC_COLS], max), ROUND_DECIMALS), range = round(sapply(df[NUMERIC_COLS], function(x) diff(range(x))), ROUND_DECIMALS), sample_variance = round(sapply(df[NUMERIC_COLS], var), ROUND_DECIMALS), population_variance= round(sapply(df[NUMERIC_COLS], pop_var), ROUND_DECIMALS), sample_std = round(sapply(df[NUMERIC_COLS], sd), ROUND_DECIMALS), population_std = round(sapply(df[NUMERIC_COLS], pop_sd), ROUND_DECIMALS), IQR = round(sapply(df[NUMERIC_COLS], function(x) IQR(x)), ROUND_DECIMALS), CV_pct = round(sapply(df[NUMERIC_COLS], cv), ROUND_DECIMALS) ) print(spread) # ── Position & Shape ──────────────────────────────────────────────────── cat(" POSITION & SHAPE (Percentiles, Skewness, Kurtosis)\n") # Percentile table (rows = columns, cols = probabilities) pct_mat <- t(sapply(df[NUMERIC_COLS], function(x) quantile(x, probs = PERCENTILE_PROBS))) colnames(pct_mat) <- paste0("P", as.integer(PERCENTILE_PROBS * 100)) # Skewness (bias = TRUE ) skew_b <- function(x) { n <- length(x) m3 <- mean((x - mean(x))^3) m2 <- mean((x - mean(x))^2) m3 / m2^1.5 # biased g1 } # Kurtosis: raw (Pearson) and excess (Fisher) kurt_raw <- function(x) { m4 <- mean((x - mean(x))^4) m2 <- mean((x - mean(x))^2) m4 / m2^2 } kurt_excess <- function(x) kurt_raw(x) - 3 shape <- data.frame( skewness_g1 = round(sapply(df[NUMERIC_COLS], skew_b), ROUND_DECIMALS), kurtosis_raw = round(sapply(df[NUMERIC_COLS], kurt_raw), ROUND_DECIMALS), excess_kurtosis= round(sapply(df[NUMERIC_COLS], kurt_excess), ROUND_DECIMALS) ) position <- cbind(round(pct_mat, ROUND_DECIMALS), shape) print(position) # Outlier fences cat(sprintf("\nOutlier fences (Q1 - %.1f*IQR, Q3 + %.1f*IQR) and outlier counts:\n", IQR_OUTLIER_MULTIPLIER, IQR_OUTLIER_MULTIPLIER)) for (col in NUMERIC_COLS) { q1 <- quantile(df[[col]], 0.25) q3 <- quantile(df[[col]], 0.75) iqr <- q3 - q1 lower <- q1 - IQR_OUTLIER_MULTIPLIER * iqr upper <- q3 + IQR_OUTLIER_MULTIPLIER * iqr n_out <- sum(df[[col]] < lower | df[[col]] > upper) cat(sprintf(” %-12s: lower=%9.3f upper=%9.3f outliers=%d\n”, col, lower, upper, n_out)) } # ── Counts ─────────────────────────────────────────────────────────────── cat(” COUNTS (Total observations)\n”) cat(sprintf(“Total observations (N): %d\n”, nrow(df))) non_missing <- sapply(df[NUMERIC_COLS], function(x) sum(!is.na(x))) print(data.frame(non_missing_count = non_missing)) # ── Bivariate ──────────────────────────────────────────────────────────── cat(" BIVARIATE (Covariance, Pearson r, Spearman rho)\n") cat("\nSample covariance matrix:\n") print(round(cov(df[NUMERIC_COLS]), ROUND_DECIMALS)) cat("\nPearson correlation matrix:\n") print(round(cor(df[NUMERIC_COLS], method = "pearson"), ROUND_DECIMALS)) cat("\nSpearman rank correlation matrix:\n") print(round(cor(df[NUMERIC_COLS], method = "spearman"), ROUND_DECIMALS)) # ============================================================================= # 2. CATEGORICAL DATA # ============================================================================= # ── Frequency Analysis ─────────────────────────────────────────────────── cat(" FREQUENCY ANALYSIS (categorical / coded columns)\n") for (col in CATEGORICAL_COLS) { freq <- sort(table(df[[col]])) # sorted by value rel <- freq / nrow(df) pct <- rel * 100 cum_freq <- cumsum(freq) cum_pct <- cumsum(pct) tbl <- data.frame( abs_freq_f = as.integer(freq), rel_freq_p = round(as.numeric(rel), ROUND_DECIMALS), pct_freq_pct = round(as.numeric(pct), ROUND_DECIMALS), cum_freq_F = as.integer(cum_freq), cum_rel_freq_pct = round(as.numeric(cum_pct), ROUND_DECIMALS), row.names = names(freq) ) cat(sprintf("\n--- %s ---\n", col)) print(tbl) } # ── Contingency Tables ─────────────────────────────────────────────────── cat(" CONTINGENCY TABLES (categorical x categorical)\n") # Cramér's V cramers_v <- function(ct) { ch <- chisq.test(ct, correct = FALSE) n <- sum(ct) r <- nrow(ct) k <- ncol(ct) sqrt((ch$statistic / n) / (min(r - 1, k - 1))) } # Build pairs list if (is.null(CATEGORICAL_PAIRS)) { pairs <- combn(CATEGORICAL_COLS, 2, simplify = FALSE) } else { pairs <- CATEGORICAL_PAIRS } for (pair in pairs) { col1 <- pair[1]; col2 <- pair[2] cat(sprintf("\n--- %s x %s ---\n", col1, col2)) ct <- table(df[[col1]], df[[col2]]) dimnames(ct) <- list(col1 = rownames(ct), col2 = colnames(ct)) cat("\nObserved frequencies:\n") print(ct) # Row % row_pct <- round(prop.table(ct, margin = 1) * 100, ROUND_DECIMALS) cat("\nRow % :\n") print(row_pct) # Chi-square + Cramér's V ch <- chisq.test(ct, correct = FALSE) v <- cramers_v(ct) sig <- ifelse(ch$p.value < 0.05, "significant (p < 0.05)", "not significant (p >= 0.05)”) cat(sprintf(“\nChi-square = %.2f, dof = %d, p-value = %.4f, Cramer’s V = %.2f\n”, ch$statistic, ch$parameter, ch$p.value, v)) cat(sprintf(“Association is %s\n”, sig)) } cat(“\nDone.\n”)

# ===========================================================================
# 0. LOAD DATA + METADATA
# ===========================================================================
df             <- read.csv(DATA_PATH, stringsAsFactors = FALSE)
df[[DATE_COL]] <- as.Date(df[[DATE_COL]])


cat("0. METADATA\n")

cat(sprintf("Shape (rows, columns): %d x %d\n", nrow(df), ncol(df)))
cat(sprintf("Date range: %s  to  %s\n",
            min(df[[DATE_COL]]), max(df[[DATE_COL]])))

meta <- data.frame(
  dtype         = sapply(df, class),
  non_null      = sapply(df, function(x) sum(!is.na(x))),
  missing_count = sapply(df, function(x) sum(is.na(x))),
  missing_pct   = round(sapply(df, function(x) mean(is.na(x)) * 100), ROUND_DECIMALS),
  n_unique      = sapply(df, function(x) length(unique(x)))
)
cat("\nColumn metadata:\n")
print(meta)

cat("\nNumeric columns  :", paste(NUMERIC_COLS,    collapse = ", "), "\n")
cat("Categorical cols :", paste(CATEGORICAL_COLS, collapse = ", "), "\n")

# ===========================================================================
# 1. NUMERIC DATA
# ===========================================================================

# Helper: population std (R's sd() uses n-1 by default)
pop_sd  <- function(x) sqrt(sum((x - mean(x))^2) / length(x))
pop_var <- function(x) sum((x - mean(x))^2) / length(x)
cv      <- function(x) (sd(x) / mean(x)) * 100   # sd() = sample std

# ──  Central Tendency ────────────────────────────────────────────────────

cat(" CENTRAL TENDENCY  (Mean, Median, Mode)\n")

first_mode <- function(x) as.numeric(names(sort(table(x), decreasing = TRUE)[1]))

central <- data.frame(
  mean   = round(sapply(df[NUMERIC_COLS], mean),       ROUND_DECIMALS),
  median = round(sapply(df[NUMERIC_COLS], median),     ROUND_DECIMALS),
  mode   = round(sapply(df[NUMERIC_COLS], first_mode), ROUND_DECIMALS)
)
print(central)

# ──  Spread ───────────────────────────────────────────────────────────────

cat(" SPREAD / DISPERSION\n")


spread <- data.frame(
  minimum            = round(sapply(df[NUMERIC_COLS], min),     ROUND_DECIMALS),
  maximum            = round(sapply(df[NUMERIC_COLS], max),     ROUND_DECIMALS),
  range              = round(sapply(df[NUMERIC_COLS], function(x) diff(range(x))), ROUND_DECIMALS),
  sample_variance    = round(sapply(df[NUMERIC_COLS], var),     ROUND_DECIMALS),
  population_variance= round(sapply(df[NUMERIC_COLS], pop_var), ROUND_DECIMALS),
  sample_std         = round(sapply(df[NUMERIC_COLS], sd),      ROUND_DECIMALS),
  population_std     = round(sapply(df[NUMERIC_COLS], pop_sd),  ROUND_DECIMALS),
  IQR                = round(sapply(df[NUMERIC_COLS], function(x) IQR(x)), ROUND_DECIMALS),
  CV_pct             = round(sapply(df[NUMERIC_COLS], cv),      ROUND_DECIMALS)
)
print(spread)

# ──  Position & Shape ────────────────────────────────────────────────────

cat(" POSITION & SHAPE  (Percentiles, Skewness, Kurtosis)\n")

# Percentile table (rows = columns, cols = probabilities)
pct_mat <- t(sapply(df[NUMERIC_COLS],
                    function(x) quantile(x, probs = PERCENTILE_PROBS)))
colnames(pct_mat) <- paste0("P", as.integer(PERCENTILE_PROBS * 100))

# Skewness (bias = TRUE )
skew_b <- function(x) {
  n  <- length(x)
  m3 <- mean((x - mean(x))^3)
  m2 <- mean((x - mean(x))^2)
  m3 / m2^1.5          # biased g1
}

# Kurtosis: raw (Pearson) and excess (Fisher)
kurt_raw <- function(x) {
  m4 <- mean((x - mean(x))^4)
  m2 <- mean((x - mean(x))^2)
  m4 / m2^2
}
kurt_excess <- function(x) kurt_raw(x) - 3

shape <- data.frame(
  skewness_g1    = round(sapply(df[NUMERIC_COLS], skew_b),      ROUND_DECIMALS),
  kurtosis_raw   = round(sapply(df[NUMERIC_COLS], kurt_raw),    ROUND_DECIMALS),
  excess_kurtosis= round(sapply(df[NUMERIC_COLS], kurt_excess), ROUND_DECIMALS)
)

position <- cbind(round(pct_mat, ROUND_DECIMALS), shape)
print(position)

# Outlier fences
cat(sprintf("\nOutlier fences (Q1 - %.1f*IQR,  Q3 + %.1f*IQR) and outlier counts:\n",
            IQR_OUTLIER_MULTIPLIER, IQR_OUTLIER_MULTIPLIER))

for (col in NUMERIC_COLS) {
  q1    <- quantile(df[[col]], 0.25)
  q3    <- quantile(df[[col]], 0.75)
  iqr   <- q3 - q1
  lower <- q1 - IQR_OUTLIER_MULTIPLIER * iqr
  upper <- q3 + IQR_OUTLIER_MULTIPLIER * iqr
  n_out <- sum(df[[col]] < lower | df[[col]] > upper)
  cat(sprintf("  %-12s: lower=%9.3f  upper=%9.3f  outliers=%d\n",
              col, lower, upper, n_out))
}

# ──  Counts ───────────────────────────────────────────────────────────────

cat(" COUNTS  (Total observations)\n")


cat(sprintf("Total observations (N): %d\n", nrow(df)))
non_missing <- sapply(df[NUMERIC_COLS], function(x) sum(!is.na(x)))
print(data.frame(non_missing_count = non_missing))

# ── Bivariate ────────────────────────────────────────────────────────────

cat(" BIVARIATE  (Covariance, Pearson r, Spearman rho)\n")


cat("\nSample covariance matrix:\n")
print(round(cov(df[NUMERIC_COLS]), ROUND_DECIMALS))

cat("\nPearson correlation matrix:\n")
print(round(cor(df[NUMERIC_COLS], method = "pearson"), ROUND_DECIMALS))

cat("\nSpearman rank correlation matrix:\n")
print(round(cor(df[NUMERIC_COLS], method = "spearman"), ROUND_DECIMALS))

# ===========================================================================
# 2. CATEGORICAL DATA
# ===========================================================================

# ──  Frequency Analysis ───────────────────────────────────────────────────

cat(" FREQUENCY ANALYSIS (categorical / coded columns)\n")


for (col in CATEGORICAL_COLS) {
  freq     <- sort(table(df[[col]]))          # sorted by value 
  rel      <- freq / nrow(df)
  pct      <- rel * 100
  cum_freq <- cumsum(freq)
  cum_pct  <- cumsum(pct)
  
  tbl <- data.frame(
    abs_freq_f     = as.integer(freq),
    rel_freq_p     = round(as.numeric(rel),      ROUND_DECIMALS),
    pct_freq_pct   = round(as.numeric(pct),      ROUND_DECIMALS),
    cum_freq_F     = as.integer(cum_freq),
    cum_rel_freq_pct = round(as.numeric(cum_pct), ROUND_DECIMALS),
    row.names      = names(freq)
  )
  cat(sprintf("\n--- %s ---\n", col))
  print(tbl)
}

# ── Contingency Tables ───────────────────────────────────────────────────

cat(" CONTINGENCY TABLES (categorical x categorical)\n")


# Cramér's V
cramers_v <- function(ct) {
  ch  <- chisq.test(ct, correct = FALSE)
  n   <- sum(ct)
  r   <- nrow(ct)
  k   <- ncol(ct)
  sqrt((ch$statistic / n) / (min(r - 1, k - 1)))
}

# Build pairs list
if (is.null(CATEGORICAL_PAIRS)) {
  pairs <- combn(CATEGORICAL_COLS, 2, simplify = FALSE)
} else {
  pairs <- CATEGORICAL_PAIRS
}

for (pair in pairs) {
  col1 <- pair[1]; col2 <- pair[2]
  cat(sprintf("\n--- %s x %s ---\n", col1, col2))
  
  ct <- table(df[[col1]], df[[col2]])
  dimnames(ct) <- list(col1 = rownames(ct), col2 = colnames(ct))
  
  cat("\nObserved frequencies:\n")
  print(ct)
  
  # Row %
  row_pct <- round(prop.table(ct, margin = 1) * 100, ROUND_DECIMALS)
  cat("\nRow % :\n")
  print(row_pct)
  
  # Chi-square + Cramér's V
  ch  <- chisq.test(ct, correct = FALSE)
  v   <- cramers_v(ct)
  sig <- ifelse(ch$p.value < 0.05,
                "significant (p < 0.05)", "not significant (p >= 0.05)")
  
  cat(sprintf("\nChi-square = %.2f, dof = %d, p-value = %.4f, Cramer's V = %.2f\n",
              ch$statistic, ch$parameter, ch$p.value, v))
  cat(sprintf("Association is %s\n", sig))
}

cat("\nDone.\n")

Usage

For hands on experience refer the folder– it contains both the code and the dataset. run the code and see the output results for your understanding

A Book of Mindless Misc

Descriptive Statistics in R: A Practical Approach

Introduction

R code:

Variable & Parameter declaration

Code block

Usage

A Book of Mindless Misc

Subscribe Now

A Book of Mindless Misc

Quick Links

Quantitative Sciences

Stay Connected