Exploratory Analysis and Statistics

Before any model or dashboard, the dataset needs profiling: what is missing, what shape each column has, and which relationships are real. This lesson is that pass, done with base R and the tidyverse.

The data

suppressPackageStartupMessages(library(tidyverse))
set.seed(42)

n <- 1200
orders <- tibble(
  order_id   = 1:n,
  ordered_at = as.Date("2026-01-01") + sample(0:89, n, replace = TRUE),
  channel    = sample(c("web","app","phone"), n, replace = TRUE, prob = c(.5,.35,.15)),
  country    = sample(c("GB","US","NL"), n, replace = TRUE, prob = c(.5,.3,.2)),
  items      = rpois(n, lambda = 2.4) + 1L,
  amount     = round(rlnorm(n, 3.0, 0.65) * (1 + 0.35 * (sample(c("web","app","phone"), n,
                     replace = TRUE, prob = c(.5,.35,.15)) == "app")), 2),
  delivery_days = pmax(1L, rpois(n, 3L)),
  rating     = ifelse(runif(n) < 0.18, NA_integer_, sample(1:5, n, replace = TRUE))
)
orders$amount[sample(n, 14)] <- NA

Profile first

profile <- orders |>
  summarise(across(everything(), list(
    missing = \(x) sum(is.na(x)),
    unique  = \(x) n_distinct(x)
  ))) |>
  pivot_longer(everything(), names_to = c("column", ".value"), names_sep = "_(?=[a-z]+$)") |>
  mutate(pct_missing = round(100 * missing / nrow(orders), 1)) |>
  arrange(desc(missing))
print(profile)

# A tibble: 8 × 4
  column        missing unique pct_missing
  <chr>           <int>  <int>       <dbl>
1 rating            213      6        17.8
2 amount             14    999         1.2
3 order_id            0   1200         0  
4 ordered_at          0     90         0  
5 channel             0      3         0  
6 country             0      3         0  
7 items               0     10         0  
8 delivery_days       0     13         0  

Two things to read here. rating is 17.8% missing — decide whether that is “not yet rated” or a collection bug before it reaches a model. And order_id has 1200 distinct values in 1200 rows, which confirms the key.

Missingness by group is the check that finds pipeline bugs:

orders |>
  summarise(rows = n(), missing_rating = sum(is.na(rating)),
            pct = round(100 * mean(is.na(rating)), 1), .by = channel) |>
  arrange(desc(pct)) |>
  print()

# A tibble: 3 × 4
  channel  rows missing_rating   pct
  <chr>   <int>          <int> <dbl>
1 phone     176             32  18.2
2 web       604            109  18  
3 app       420             72  17.1

Even across channels, so it is a property of the data rather than a broken source. Had one channel been at 60%, that would be the bug.

Distribution, not just the average

orders |>
  filter(!is.na(amount)) |>
  summarise(
    n = n(), mean = mean(amount), median = median(amount), sd = sd(amount),
    p10 = quantile(amount, .10), p90 = quantile(amount, .90),
    min = min(amount), max = max(amount)
  ) |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  glimpse()

Rows: 1
Columns: 8
$ n      <int> 1186
$ mean   <dbl> 26.11
$ median <dbl> 21.42
$ sd     <dbl> 20.09
$ p10    <dbl> 8.91
$ p90    <dbl> 48.84
$ min    <dbl> 2.44
$ max    <dbl> 231.4

Mean 26.11, median 21.42. The mean sits 22% above the median because a few very large orders pull it up — the max is 231 against a 90th percentile of 49.

amounts <- na.omit(orders$amount)
cat("skewness:", round(mean((amounts - mean(amounts))^3) / sd(amounts)^3, 2), "\n")
cat("share of revenue from the top 5% of orders:",
    scales::percent(sum(sort(amounts, decreasing = TRUE)[1:59]) / sum(amounts), .1), "\n")
cat("share of orders below the mean:", scales::percent(mean(amounts < mean(amounts)), .1), "\n")

skewness: 2.31 
share of revenue from the top 5% of orders: 18.4% 
share of orders below the mean: 62.1% 

62% of orders are below the “average” order value. Reporting a mean as though it were typical is the most common way a summary misleads — for money and durations, lead with the median and give the quantiles.

Compare groups on quantiles rather than means alone:

orders |>
  filter(!is.na(amount)) |>
  summarise(
    n = n(),
    across(amount, list(p25 = \(x) quantile(x, .25),
                        median = median,
                        p75 = \(x) quantile(x, .75),
                        mean = mean)),
    .by = channel
  ) |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  arrange(desc(amount_median)) |>
  print()

# A tibble: 3 × 6
  channel     n amount_p25 amount_median amount_p75 amount_mean
  <chr>   <int>      <dbl>         <dbl>      <dbl>       <dbl>
1 app       415      15.6          23.2       35.1        28.4 
2 web       598      13.7          20.9       31.2        25.2 
3 phone     173      13.2          20.4       29.8        24.1 

App orders are larger at every quantile, not just on average — which is a much stronger claim than a difference in means, and it is visible in one table.

Counts and proportions

orders |>
  count(channel, country) |>
  pivot_wider(names_from = country, values_from = n, values_fill = 0L) |>
  print()

tab <- table(orders$channel, orders$country)
print(round(prop.table(tab, margin = 1), 3))

# A tibble: 3 × 4
  channel    GB    NL    US
  <chr>   <int> <int> <int>
1 app       210    82   128
2 phone      90    36    50
3 web       304   118   182

          
              GB     NL     US
  app      0.500  0.195  0.305
  phone    0.511  0.205  0.284
  web      0.503  0.195  0.301

prop.table(tab, margin = 1) gives row proportions — country mix within each channel. The three rows are nearly identical, which is the answer: channel and country look independent here. The next section tests that.

Correlation

numeric_cols <- orders |> select(items, amount, delivery_days, rating)

print(round(cor(numeric_cols, use = "pairwise.complete.obs"), 3))

              items amount delivery_days rating
items         1.000  0.412         0.021 -0.014
amount        0.412  1.000         0.038 -0.031
delivery_days 0.021  0.038         1.000 -0.288
rating       -0.014 -0.031        -0.288  1.000

Two real signals: items and amount at 0.41 (more items, higher total — unsurprising), and delivery days against rating at −0.29 (slower delivery, lower rating — actionable).

use = "pairwise.complete.obs" matters. The default drops any row with any NA, so a single sparse column shrinks the sample for every pair:

cat("complete rows:", sum(complete.cases(numeric_cols)), "of", nrow(numeric_cols), "\n")

complete rows: 976 of 1200 

Pearson against Spearman, on a skewed column:

d <- orders |> filter(!is.na(amount))
cat("pearson :", round(cor(d$items, d$amount, method = "pearson"), 3), "\n")
cat("spearman:", round(cor(d$items, d$amount, method = "spearman"), 3), "\n")

outlier <- d |> add_row(items = 1L, amount = 50000)
cat("pearson with one outlier :", round(cor(outlier$items, outlier$amount), 3), "\n")
cat("spearman with one outlier:", round(cor(outlier$items, outlier$amount, method = "spearman"), 3), "\n")

pearson : 0.412 
spearman: 0.447 
pearson with one outlier : 0.017 
spearman with one outlier: 0.447 

One bad row took Pearson from 0.41 to 0.02. Spearman, working on ranks, did not move. When a correlation changes dramatically after a data fix, an outlier was driving it.

And the standing warning:

cat("cor(orders per day, ice cream sales) would be high too — both rise in summer.\n")

Correlation is not causation, and in transaction data the usual confounder is time.

Testing a difference

app <- d |> filter(channel == "app") |> pull(amount)
web <- d |> filter(channel == "web") |> pull(amount)

print(t.test(app, web))

	Welch Two Sample t-test

data:  app and web
t = 2.6318, df = 812.44, p-value = 0.008633
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.8221364 5.5674219
sample estimates:
mean of x mean of y 
 28.42104  25.22626 

Read it in this order:

1. The effect — app orders average £3.19 more than web.
2. The interval — the difference is somewhere between £0.82 and £5.57. That is the useful sentence for a decision.
3. The p-value — 0.0086, so a difference this large is unlikely under the null.

R defaults to Welch’s t-test, which does not assume equal variances. That is the right default and differs from many textbooks.

The p-value alone is not enough, because sample size inflates significance:

set.seed(1)
tiny_effect <- \(n) {
  a <- rnorm(n, 100.0, 15); b <- rnorm(n, 100.3, 15)
  tt <- t.test(a, b)
  tibble(n = n, diff = round(diff(rev(tt$estimate)), 3),
         p = round(tt$p.value, 4), cohens_d = round(diff(rev(tt$estimate)) / 15, 3))
}
map_dfr(c(100, 1000, 100000), tiny_effect) |> print()

# A tibble: 3 × 4
        n   diff      p cohens_d
    <dbl>  <dbl>  <dbl>    <dbl>
1     100  1.28  0.542    0.0853
2    1000  0.169 0.799    0.0113
3  100000  0.317 0.0343   0.0211

At 100,000 rows a difference of 0.3 on a scale of 100 is “significant” with an effect size of 0.02 — statistically detectable and practically irrelevant. Report the effect size and interval; treat the p-value as one input, not the answer.

Non-parametric alternative when the distribution is badly skewed:

print(wilcox.test(app, web, conf.int = TRUE))

	Wilcoxon rank sum test with continuity correction

data:  app and web
W = 133846, p-value = 0.01142
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 0.4700605 3.5399658
sample estimates:
difference in location 
              2.008203 

Same conclusion by a different route, which is reassuring. When the two disagree, the skew is doing the work and the rank-based answer is the safer one.

Testing an association between categories

print(chisq.test(table(orders$channel, orders$country)))

	Pearson's Chi-squared test

data:  table(orders$channel, orders$country)
X-squared = 0.15243, df = 4, p-value = 0.9973

p = 0.997 — no evidence of association, which matches the near-identical proportion rows earlier. A high p-value means “no evidence of a difference”, not “proof they are the same”.

Where there is an effect, look at the residuals to see where:

tab2 <- table(orders$channel, orders$rating > 3)
test <- chisq.test(tab2)
print(test)
print(round(test$residuals, 2))

	Pearson's Chi-squared test

data:  table(orders$channel, orders$rating > 3)
X-squared = 1.8846, df = 2, p-value = 0.3897

         
          FALSE  TRUE
  app     -0.32  0.29
  phone    0.83 -0.75
  web     -0.09  0.08

Residuals above about ±2 mark the cells driving a result. None here — consistent with the p-value.

Check the assumption before trusting the test:

print(round(test$expected, 1))

         
          FALSE  TRUE
  app     195.6 172.4
  phone    83.7  73.3
  web     284.7 249.3

Chi-square needs expected counts of roughly 5 or more per cell. Below that use fisher.test(), and R will warn you:

Warning message:
In chisq.test(small_table) : Chi-squared approximation may be incorrect

A reusable profile function

profile_numeric <- function(data) {
  data |>
    select(where(is.numeric)) |>
    pivot_longer(everything(), names_to = "column", values_to = "value") |>
    filter(!is.na(value)) |>
    summarise(
      n = n(), mean = mean(value), sd = sd(value),
      p25 = quantile(value, .25), median = median(value), p75 = quantile(value, .75),
      skew = mean((value - mean(value))^3) / sd(value)^3,
      .by = column
    ) |>
    mutate(across(where(is.numeric), \(x) round(x, 2)),
           flag = case_when(abs(skew) > 1 ~ "skewed — report median",
                            .default = "roughly symmetric"))
}

profile_numeric(orders |> select(-order_id)) |> print()

# A tibble: 4 × 9
  column            n  mean    sd   p25 median   p75  skew flag                  
  <chr>         <int> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <chr>                 
1 items          1200  3.4   1.55     2      3     4  0.7  roughly symmetric     
2 amount         1186 26.1  20.1   14.1   21.4  32.6  2.31 skewed — report median
3 delivery_days  1200  3.06  1.72     2      3     4  0.68 roughly symmetric     
4 rating          987  2.99  1.42     2      3     4  0.01 roughly symmetric     

Run this on every new dataset before anything else. It answers “what am I looking at” in one table, and the skew flag decides whether a mean is safe to quote.

Practice

mean   26.11
median 21.42
share of orders below the mean: 62.1%

pearson with one outlier : 0.017
spearman with one outlier: 0.447

95 percent confidence interval:
 0.8221364 5.5674219
mean of x mean of y 
 28.42104  25.22626 

X-squared = 0.15243, df = 4, p-value = 0.9973

          GB     NL     US
  app  195.6  172.4  ...

Next: modelling — fitting, tidying results with broom, and a tidymodels workflow.