R for Data Analysis: A Practical Guide for Engineers

1
2
3
4
5
6
7
8
9

# magrittr pipe (package required)
library(magrittr)
mtcars %>% filter(cyl == 6) %>% summarise(mean_mpg = mean(mpg))

# native pipe (base R 4.1+)
mtcars |> subset(cyl == 6) |> (\(d) mean(d$mpg))()

# native pipe with placeholder (R 4.2+)
mtcars |> subset(cyl == 6) |> lm(mpg ~ wt, data = _)

1
2
3
4
5
6
7
8
9

# old way
sapply(1:5, function(x) x^2)

# new way
sapply(1:5, \(x) x^2)

# practical in pipes
list(a = 1:10, b = 11:20) |>
  lapply(\(v) c(mean = mean(v), sd = sd(v)))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# Installing packages
install.packages("tidyverse")  # from CRAN
install.packages("data.table") # from CRAN

# Bioconductor packages
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DESeq2")    # differential expression analysis
BiocManager::install("limma")     # linear models for microarray/RNA-seq

# Check what's installed
installed.packages()[, c("Package", "Version")] |> head(10)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# No scalars - these are all vectors of length 1
typeof(5L)       # "integer"
typeof(5.0)      # "double"
typeof(TRUE)     # "logical"
typeof("hello")  # "character"

length(42)       # 1, not a "scalar"

# Vector arithmetic is element-wise and automatic
x <- c(1, 2, 3, 4, 5)
y <- c(10, 20, 30, 40, 50)
x + y       # c(11, 22, 33, 44, 55) - no loop needed
x * 2       # c(2, 4, 6, 8, 10) - scalar recycled
x > 3       # c(FALSE, FALSE, FALSE, TRUE, TRUE)

# Recycling - shorter vector wraps around (often useful, sometimes a footgun)
c(1, 2, 3, 4, 5, 6) + c(10, 20)  # c(11, 22, 13, 24, 15, 26)
# R warns if recycling is not an exact multiple

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

v <- c("a", "b", "c", "d", "e")
v[1]      # "a"  (NOT "b")
v[5]      # "e"
v[6]      # NA  (out of bounds returns NA, not an error)
v[-1]     # c("b","c","d","e") - negative index means "exclude"
v[c(1,3)] # c("a","c") - vector indexing
v[2:4]    # c("b","c","d") - range indexing, both ends inclusive

# Matrix indexing: [row, col]
m <- matrix(1:9, nrow = 3)
m[2, 3]   # row 2, column 3
m[1, ]    # entire first row
m[, 2]    # entire second column

# Data frame: similar, but also supports $ and [[ for columns
df <- data.frame(x = 1:3, y = c("a","b","c"))
df[1, ]       # first row
df[, "x"]     # column x as vector
df$x          # column x (common idiom)
df[["x"]]     # column x, drops to vector

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# Convention: <- for assignment
x <- 42
name <- "R"
results <- list()

# = works but style convention reserves it for function args
mean(x = c(1, 2, 3))   # x is a named argument here

# The difference matters inside complex expressions
# This assigns x=3 in the global environment AND passes x to median:
median(x <- 3)   # x is now 3 in global env, result is 3

# This just passes x=3 to median without creating a global variable:
median(x = 3)

# <<- is the superassignment operator: assigns to parent environment
f <- function() {
  outer_var <<- "I escaped the function scope"
}
f()
outer_var  # "I escaped the function scope"

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# Creating factors
status <- factor(c("low", "high", "medium", "low", "high"))
levels(status)   # c("high", "low", "medium") - alphabetical default
nlevels(status)  # 3

# Ordered factors
severity <- factor(c("low", "high", "medium"),
                   levels = c("low", "medium", "high"),
                   ordered = TRUE)
severity[1] < severity[3]  # TRUE: "low" < "high"

# Factors in modeling: baseline is always the first level
# Changing the reference level matters for interpretation
status <- relevel(status, ref = "low")  # low is now baseline

# Factors can be a footgun when you don't want them
# read.csv historically created factors by default; readr does not
df <- read.csv("data.csv", stringsAsFactors = FALSE)  # safe option

# forcats package (tidyverse) handles factors sensibly
library(forcats)
fct_infreq(status)    # reorder by frequency
fct_rev(status)       # reverse level order
fct_lump_n(status, 2) # collapse infrequent levels to "Other"

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# NA propagates
x <- c(1, 2, NA, 4, 5)
mean(x)          # NA -- cannot compute mean with missing values
mean(x, na.rm = TRUE)  # 3.0 -- explicit removal

sum(x)           # NA
sum(x, na.rm = TRUE)   # 12

# Checking for NA (never use == NA)
is.na(x)         # c(FALSE, FALSE, TRUE, FALSE, FALSE)
x == NA          # c(NA, NA, NA, NA, NA) -- wrong way

# NA has typed variants
NA_integer_
NA_real_
NA_complex_
NA_character_  # important for maintaining column types in data frames

# Handling NAs in filtering
df <- data.frame(x = c(1, 2, NA, 4), y = c("a","b","c",NA))
df[!is.na(df$x), ]          # rows where x is not NA
na.omit(df)                  # drop any row with any NA
complete.cases(df)           # logical: which rows have no NA

# tidyverse approach
library(dplyr)
df |> filter(!is.na(x))
df |> tidyr::drop_na()           # drop rows with any NA
df |> tidyr::drop_na(x)          # drop rows where x is NA

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# Lexical scoping example
x <- 10
f <- function() {
  x <- 20    # local x
  g()        # g sees the DEFINING environment's x, not f's x
}
g <- function() x   # g was defined in global env where x = 10
f()  # returns 10, not 20

# Closures capture the enclosing environment
make_adder <- function(n) {
  function(x) x + n   # n is captured from make_adder's environment
}
add5 <- make_adder(5)
add5(3)   # 8

# This is the foundation of R's apply/functional patterns
# and why factory functions work cleanly in R

# Environments are first-class objects
e <- new.env()
e$value <- 42
ls(e)           # "value"
get("value", envir = e)  # 42

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# Basic formula: response ~ predictors
lm(mpg ~ wt, data = mtcars)          # simple regression
lm(mpg ~ wt + hp, data = mtcars)     # multiple regression
lm(mpg ~ wt * hp, data = mtcars)     # interaction (wt + hp + wt:hp)
lm(mpg ~ wt:hp, data = mtcars)       # interaction only, no main effects
lm(mpg ~ ., data = mtcars)           # all other columns as predictors
lm(mpg ~ . - carb, data = mtcars)    # all except carb

# One-sided formulas: just predictors (used in model.matrix, recipes, etc.)
~ wt + hp

# Formula objects
f <- mpg ~ wt + hp
class(f)       # "formula"
terms(f)       # terms object with attribute information
all.vars(f)    # c("mpg", "wt", "hp")

# Formulas in non-modeling contexts
library(ggplot2)
facet_wrap(~ cyl)   # ggplot2 uses formulas for faceting

1
2
3
4
5
6
7
8
9

# Install the whole collection
install.packages("tidyverse")

# Or individual packages
install.packages(c("dplyr", "tidyr", "ggplot2", "purrr",
                   "readr", "stringr", "forcats", "lubridate"))

# Load it
library(tidyverse)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

library(dplyr)

# Using nycflights13 as a realistic dataset
library(nycflights13)

# FILTER: keep rows matching a condition
flights |>
  filter(carrier == "AA", month == 1, dep_delay > 30)

# Multiple conditions: comma means AND
flights |>
  filter(month %in% c(6, 7, 8), dest %in% c("LAX", "SFO"))

# SELECT: choose and rename columns
flights |>
  select(year, month, day, dep_time, arr_time, carrier, flight)

# Select helpers
flights |>
  select(year:day, starts_with("dep"), ends_with("time"))

flights |>
  select(-year, -month, -day)  # drop columns with -

# Rename while selecting
flights |>
  select(departure = dep_time, arrival = arr_time)

# MUTATE: create or modify columns
flights |>
  mutate(
    gain = arr_delay - dep_delay,
    speed = distance / (air_time / 60),
    on_time = dep_delay <= 0
  )

# transmute() returns only the new/modified columns
flights |>
  transmute(
    gain = arr_delay - dep_delay,
    hours = air_time / 60
  )

# ARRANGE: sort rows
flights |>
  arrange(desc(dep_delay))  # most delayed first

flights |>
  arrange(year, month, day, dep_time)  # multi-column sort

# SUMMARISE: collapse rows to aggregates
flights |>
  summarise(
    n = n(),
    mean_delay = mean(dep_delay, na.rm = TRUE),
    median_delay = median(dep_delay, na.rm = TRUE),
    max_delay = max(dep_delay, na.rm = TRUE),
    pct_on_time = mean(dep_delay <= 0, na.rm = TRUE)
  )

# GROUP_BY: the power multiplier
flights |>
  group_by(carrier) |>
  summarise(
    n_flights = n(),
    mean_dep_delay = mean(dep_delay, na.rm = TRUE),
    mean_arr_delay = mean(arr_delay, na.rm = TRUE)
  ) |>
  arrange(desc(mean_dep_delay))

# group_by + mutate: window functions (rank within group, etc.)
flights |>
  group_by(carrier) |>
  mutate(
    rank_by_delay = rank(desc(dep_delay)),
    delay_vs_carrier_mean = dep_delay - mean(dep_delay, na.rm = TRUE)
  )

# Chaining all verbs: a realistic pipeline
flights |>
  filter(!is.na(dep_delay), !is.na(arr_delay)) |>
  group_by(carrier, origin) |>
  summarise(
    n = n(),
    mean_dep_delay = mean(dep_delay),
    mean_arr_delay = mean(arr_delay),
    pct_delayed = mean(dep_delay > 15),
    .groups = "drop"    # ungroup after summarise (good practice)
  ) |>
  filter(n >= 100) |>
  arrange(desc(pct_delayed))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

# Sample data
employees <- tibble(
  id = 1:5,
  name = c("Alice", "Bob", "Carol", "Dave", "Eve"),
  dept_id = c(1, 2, 1, 3, NA)
)

departments <- tibble(
  dept_id = 1:3,
  dept_name = c("Engineering", "Marketing", "Finance")
)

# left_join: keep all left rows, match right where possible
employees |>
  left_join(departments, by = "dept_id")
# Eve (dept_id = NA) still appears, dept_name will be NA

# inner_join: only matching rows
employees |>
  inner_join(departments, by = "dept_id")
# Eve excluded (no dept_id match)

# full_join: all rows from both
employees |>
  full_join(departments, by = "dept_id")

# anti_join: rows in left that have NO match in right (very useful)
employees |>
  anti_join(departments, by = "dept_id")
# Returns Eve: the employee with no department

# semi_join: rows in left that DO have a match, but don't add right columns
employees |>
  semi_join(departments, by = "dept_id")

# Joining on columns with different names
sales <- tibble(employee_id = c(1, 2, 3), revenue = c(100, 200, 150))
employees |>
  left_join(sales, by = c("id" = "employee_id"))

# Multiple key columns
flights |>
  left_join(planes, by = "tailnum")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

library(tidyr)

# Wide data: each time point is a column (common in spreadsheets)
wide_data <- tibble(
  patient = c("P001", "P002", "P003"),
  week_0  = c(10.2, 12.5, 9.8),
  week_4  = c(8.7, 11.2, 8.1),
  week_8  = c(7.3, 10.8, 7.9),
  week_12 = c(6.9, 10.1, 7.2)
)

# pivot_longer: wide -> long (tidy)
long_data <- wide_data |>
  pivot_longer(
    cols = starts_with("week"),
    names_to = "week",
    values_to = "score",
    names_prefix = "week_"   # strip the prefix
  ) |>
  mutate(week = as.integer(week))

# pivot_wider: long -> wide
long_data |>
  pivot_wider(
    names_from = week,
    values_from = score,
    names_prefix = "week_"
  )

# unnest: flatten list-columns
nested <- tibble(
  group = c("A", "B"),
  data = list(
    tibble(x = 1:3, y = c(2.1, 3.4, 1.8)),
    tibble(x = 4:6, y = c(5.2, 4.7, 6.1))
  )
)
nested |> unnest(data)

# nest: create list-columns (useful for grouped modeling)
flights |>
  group_by(carrier) |>
  nest()    # returns a tibble with one row per carrier and a 'data' list-column

# separate and unite for splitting/combining columns
df <- tibble(date_str = c("2024-01-15", "2024-02-20", "2024-03-05"))
df |>
  separate(date_str, into = c("year", "month", "day"), sep = "-", convert = TRUE)

df2 <- tibble(year = 2024, month = 1:3, day = c(15, 20, 5))
df2 |>
  unite("date_str", year, month, day, sep = "-")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

library(purrr)

# map returns a list
map(1:5, \(x) x^2)              # list of results
map_dbl(1:5, \(x) x^2)          # numeric vector
map_chr(1:5, \(x) paste0("item_", x))  # character vector
map_lgl(1:5, \(x) x > 3)        # logical vector
map_int(1:5, \(x) x * 2L)       # integer vector

# map_df / map_dfr: return a data frame (binding rows)
map_dfr(c("EWR", "JFK", "LGA"), \(airport) {
  flights |>
    filter(origin == airport) |>
    summarise(
      origin = airport,
      n = n(),
      mean_delay = mean(dep_delay, na.rm = TRUE)
    )
})

# map2: map over two inputs in parallel
x <- list(1:5, 6:10, 11:15)
y <- list(10, 20, 30)
map2_dbl(x, y, \(v, mult) mean(v) * mult)

# pmap: map over arbitrary number of inputs
params <- list(
  mean = c(0, 5, 10),
  sd   = c(1, 2, 3),
  n    = c(100, 200, 300)
)
pmap(params, rnorm)   # generates 3 vectors with different params

# walk: map for side effects (printing, writing files, plotting)
walk(c("EWR", "JFK", "LGA"), \(airport) {
  cat("Processing", airport, "\n")
  # write a file, create a plot, etc.
})

# reduce: fold a list into a single value
reduce(list(1:5, 6:10, 11:15), c)   # concatenate: 1:15
reduce(c(2, 3, 4), `*`, .init = 1)  # 24: factorial-style

# keep and discard: filter lists
x <- list(a = 1, b = "hello", c = TRUE, d = 3.14)
keep(x, is.numeric)     # list(a = 1, d = 3.14)
discard(x, is.character) # everything except b

# pluck: safe element extraction (no NULL errors)
nested_list <- list(a = list(b = list(c = 42)))
pluck(nested_list, "a", "b", "c")  # 42, safe even if path doesn't exist

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# Native pipe (no package required, R 4.1+)
mtcars |>
  filter(cyl == 6) |>
  select(mpg, wt, hp) |>
  arrange(desc(mpg))

# The key difference: placeholder for non-first-arg piping
# Native pipe (R 4.2+): uses _ placeholder in named argument position
mtcars |> lm(mpg ~ wt, data = _)

# Magrittr: uses . placeholder, more flexible
library(magrittr)
mtcars %>% lm(mpg ~ wt, data = .)

# For the vast majority of pipelines, they are identical
# Use |> for new code unless you need the . placeholder frequently

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

library(ggplot2)

# The minimum: data + aesthetic mapping + geom
ggplot(data = mtcars, aes(x = wt, y = mpg)) +
  geom_point()

# Add more layers
ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point(aes(color = factor(cyl), size = hp), alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "black", linetype = "dashed") +
  labs(
    title = "Fuel Efficiency vs Weight",
    subtitle = "32 automobile models from Motor Trend (1974)",
    x = "Weight (1000 lbs)",
    y = "Miles per Gallon",
    color = "Cylinders",
    size = "Horsepower",
    caption = "Source: mtcars dataset"
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# Global aes: shared across all layers
ggplot(mtcars, aes(x = wt, y = mpg, color = factor(cyl))) +
  geom_point() +                          # inherits color
  geom_smooth(method = "lm", se = FALSE)  # also inherits color (per group)

# Local aes: override or extend for one layer
ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point(aes(color = factor(cyl))) +  # color only for points
  geom_smooth(method = "lm", se = FALSE,
              color = "black")            # single black line overall

# Fixed vs mapped aesthetics
ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point(color = "steelblue", size = 3, alpha = 0.8)  # fixed values, outside aes()

# Aesthetic mappings available: x, y, color/colour, fill, size, shape,
# alpha, linetype, linewidth, label, group

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

library(nycflights13)

# Scatter: geom_point
ggplot(flights |> sample_n(5000), aes(x = dep_delay, y = arr_delay)) +
  geom_point(alpha = 0.1) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")

# Line: geom_line (time series)
flights |>
  filter(carrier == "AA") |>
  count(month) |>
  ggplot(aes(x = month, y = n)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 1:12)

# Bar: geom_bar (counts) vs geom_col (precomputed values)
flights |>
  count(carrier) |>
  ggplot(aes(x = reorder(carrier, -n), y = n)) +
  geom_col(fill = "steelblue") +
  labs(x = "Carrier", y = "Number of Flights")

# Histogram
ggplot(flights |> filter(!is.na(dep_delay)), aes(x = dep_delay)) +
  geom_histogram(bins = 60, fill = "steelblue", color = "white") +
  xlim(-30, 120)

# Density
ggplot(flights |> filter(!is.na(dep_delay), carrier %in% c("AA","UA","DL")),
       aes(x = dep_delay, fill = carrier)) +
  geom_density(alpha = 0.4) +
  xlim(-30, 120)

# Boxplot + jitter: show distribution AND points
ggplot(flights |> filter(!is.na(dep_delay), carrier %in% c("AA","UA","DL")),
       aes(x = carrier, y = dep_delay)) +
  geom_boxplot(outlier.shape = NA) +  # hide outlier points (shown by jitter)
  geom_jitter(width = 0.2, alpha = 0.05, color = "steelblue") +
  ylim(-30, 120)

# Violin plot
ggplot(flights |> filter(!is.na(dep_delay), carrier %in% c("AA","UA","DL")),
       aes(x = carrier, y = dep_delay)) +
  geom_violin(fill = "steelblue", alpha = 0.7) +
  geom_boxplot(width = 0.1, fill = "white") +
  ylim(-30, 120)

# Heatmap with geom_tile
flights |>
  filter(!is.na(dep_delay)) |>
  group_by(month, day_of_week = wday(time_hour, label = TRUE)) |>
  summarise(mean_delay = mean(dep_delay), .groups = "drop") |>
  ggplot(aes(x = month, y = day_of_week, fill = mean_delay)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                       midpoint = 0) +
  scale_x_continuous(breaks = 1:12)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# facet_wrap: wrap panels into rows and columns
ggplot(flights |> filter(!is.na(dep_delay), carrier %in% c("AA","UA","DL","B6")),
       aes(x = dep_delay)) +
  geom_histogram(bins = 40, fill = "steelblue") +
  facet_wrap(~ carrier, ncol = 2) +  # 2 columns
  xlim(-30, 120)

# facet_grid: explicit row and column facet variables
ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point() +
  facet_grid(am ~ cyl,   # rows = am, columns = cyl
             labeller = labeller(
               am = c("0" = "Automatic", "1" = "Manual"),
               cyl = c("4" = "4 cyl", "6" = "6 cyl", "8" = "8 cyl")
             ))

# scales = "free" allows axis scales to vary per facet
ggplot(flights |> filter(!is.na(dep_delay)), aes(x = dep_delay)) +
  geom_histogram(bins = 40) +
  facet_wrap(~ origin, scales = "free_y")  # y-axis varies, x-axis fixed

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

# Built-in themes
p <- ggplot(mtcars, aes(x = wt, y = mpg, color = factor(cyl))) +
  geom_point(size = 2)

p + theme_minimal()
p + theme_bw()
p + theme_classic()
p + theme_void()   # no axes, no background (good for maps)

# Custom theme elements
p + theme(
  plot.title = element_text(size = 16, face = "bold"),
  axis.text = element_text(size = 11),
  legend.position = "bottom",
  panel.grid.minor = element_blank(),
  plot.background = element_rect(fill = "#f5f5f5", color = NA)
)

# Scales: control how data maps to visual properties
# Color scales
p + scale_color_brewer(palette = "Set1")    # ColorBrewer qualitative
p + scale_color_viridis_d()                  # viridis (perceptually uniform)
p + scale_color_manual(values = c("4" = "#E41A1C", "6" = "#377EB8", "8" = "#4DAF4A"))

# Continuous color scale
ggplot(mtcars, aes(x = wt, y = mpg, color = hp)) +
  geom_point(size = 3) +
  scale_color_viridis_c(option = "plasma")  # plasma colormap

# Scale transformations
ggplot(flights |> filter(!is.na(dep_delay)), aes(x = dep_delay)) +
  geom_histogram(bins = 50) +
  scale_y_log10()  # log scale on y-axis

# Axis limits and breaks
p +
  scale_x_continuous(
    breaks = seq(1, 6, by = 0.5),
    limits = c(1.5, 5.5),
    labels = scales::comma
  ) +
  scale_y_continuous(
    breaks = seq(10, 35, by = 5),
    labels = \(x) paste0(x, " mpg")
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

library(patchwork)

p1 <- ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point() + theme_minimal()
p2 <- ggplot(mtcars, aes(x = hp, y = mpg)) + geom_point() + theme_minimal()
p3 <- ggplot(mtcars, aes(x = factor(cyl))) + geom_bar() + theme_minimal()

# Side by side
p1 + p2

# Stacked
p1 / p2

# Complex layouts
(p1 + p2) / p3

# With shared title and annotations
(p1 + p2 + p3) +
  plot_annotation(
    title = "Motor Trend Car Road Tests",
    tag_levels = "A"   # adds (A), (B), (C) labels
  )

# inset_element: embed one plot inside another
p1 + inset_element(p3, left = 0.6, bottom = 0.6, right = 1, top = 1)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# ggsave: the correct way to save
ggsave("my_plot.png", width = 10, height = 7, dpi = 300)

# Save a specific plot object
ggsave("my_plot.pdf", plot = p1, width = 8, height = 6)

# For publications: PDF or SVG (vector formats)
ggsave("figure1.svg", width = 180, height = 120, units = "mm")
ggsave("figure1.eps", width = 7, height = 5)  # EPS for LaTeX

# For web: PNG at appropriate resolution
ggsave("plot.png", width = 12, height = 8, dpi = 150)

# Ragg device for better font rendering (install ragg package)
ggsave("plot.png", device = ragg::agg_png, width = 12, height = 8, dpi = 300)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

library(readr)

# read_csv: the readr version (NOT base R read.csv)
# - reads strings as characters (not factors)
# - parses dates automatically
# - faster for large files
# - returns a tibble
df <- read_csv("data.csv")
df <- read_csv("data.csv",
               col_types = cols(
                 date   = col_date(format = "%Y-%m-%d"),
                 price  = col_double(),
                 status = col_factor(levels = c("open", "closed"))
               ),
               skip = 1,        # skip header rows
               n_max = 10000,   # read only first 10k rows
               na = c("", "NA", "NULL", "-"))

# read_delim for other delimiters
df <- read_delim("data.tsv", delim = "\t")
df <- read_delim("data.txt", delim = "|")

# Reading from URLs
covid <- read_csv("https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv")

# write_csv: fast output
write_csv(df, "output.csv")
write_tsv(df, "output.tsv")

# readxl for Excel files
library(readxl)
df_excel <- read_excel("report.xlsx",
                        sheet = "Sheet2",
                        range = "A1:F100",
                        col_types = c("text", "numeric", "date", "text", "numeric", "numeric"))
excel_sheets("report.xlsx")  # list sheet names

# haven for SPSS/Stata/SAS files
library(haven)
df_spss <- read_sav("survey.sav")
df_stata <- read_dta("analysis.dta")
df_sas <- read_sas("data.sas7bdat")
# haven preserves labelled values from SPSS/Stata

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

library(janitor)

# clean_names: lowercase, underscores, removes special characters
messy_df <- data.frame(
  `First Name` = c("Alice", "Bob"),
  `Last Name`  = c("Smith", "Jones"),
  `2023 Revenue ($)` = c(100000, 200000),
  check.names = FALSE
)
clean_df <- messy_df |> clean_names()
# columns now: first_name, last_name, x2023_revenue

# tabyl: better table() with proportions
flights |>
  tabyl(carrier, origin) |>
  adorn_totals(where = c("row", "col")) |>
  adorn_percentages("row") |>
  adorn_pct_formatting(digits = 1) |>
  adorn_ns()

# remove_empty: drop empty rows and columns
df_clean <- messy_df |> remove_empty(which = c("rows", "cols"))

# get_dupes: find duplicate rows
flights |>
  get_dupes(flight, tailnum, origin, dest)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

library(lubridate)

# Parsing dates
ymd("2024-01-15")          # "2024-01-15"
mdy("January 15, 2024")    # "2024-01-15"
dmy("15/01/2024")          # "2024-01-15"
ymd_hms("2024-01-15 14:30:00")  # datetime with time

# Extracting components
d <- ymd_hms("2024-06-15 09:30:00", tz = "America/New_York")
year(d)    # 2024
month(d)   # 6
day(d)     # 15
hour(d)    # 9
minute(d)  # 30
wday(d, label = TRUE)   # "Sat"
yday(d)    # 167 (day of year)
week(d)    # ISO week number

# Arithmetic
d + days(7)           # add 7 days
d + months(1)         # add 1 month (handles month length correctly)
d - hours(2)          # subtract 2 hours
interval(ymd("2024-01-01"), ymd("2024-12-31")) / days(1)  # days in interval

# Time zones
with_tz(d, "UTC")           # convert to UTC
force_tz(d, "Europe/London") # declare time zone (no conversion)

# Practical example: calculate business days
flights |>
  mutate(
    flight_date = as.Date(time_hour),
    day_of_week = wday(flight_date, label = TRUE),
    is_weekend  = wday(flight_date) %in% c(1, 7),
    month_name  = month(flight_date, label = TRUE, abbr = FALSE),
    quarter     = quarter(flight_date)
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

library(stringr)

# str_detect: logical match (use with filter)
str_detect(c("apple", "banana", "apricot"), "^a")  # TRUE FALSE TRUE

# str_extract and str_extract_all
emails <- c("alice@company.com", "bob@personal.org", "invalid")
str_extract(emails, "[a-z]+@[a-z]+\\.[a-z]+")

# str_replace and str_replace_all
str_replace("hello world", "world", "R")
str_replace_all("2024-01-15", "-", "/")

# str_split
str_split("a,b,c,d", ",")[[1]]          # c("a","b","c","d")
str_split_fixed("a,b,c,d", ",", n = 2)  # matrix: "a" | "b,c,d"

# str_c: concatenate (vectorized paste)
str_c("item", 1:5, sep = "_")           # "item_1" ... "item_5"
str_c(c("x", "y", "z"), collapse = "+") # "x+y+z"

# str_pad and str_trim
str_pad("42", width = 6, side = "left", pad = "0")  # "000042"
str_trim("  hello world  ")                           # "hello world"

# str_glue: the glue-based interpolation (preferred over paste)
name  <- "Alice"
score <- 95.5
str_glue("Participant {name} scored {score} (grade: {if(score >= 90) 'A' else 'B'})")

# In a dplyr pipeline
flights |>
  mutate(
    route = str_c(origin, dest, sep = " -> "),
    carrier_flight = str_glue("{carrier}{flight}"),
    is_jfk_origin = str_detect(origin, "JFK")
  )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

# Linear regression
model <- lm(mpg ~ wt + hp + factor(cyl), data = mtcars)

# Comprehensive output
summary(model)
# Coefficients, standard errors, t-values, p-values, R², F-statistic

# Model diagnostics
par(mfrow = c(2, 2))  # 2x2 plot grid
plot(model)            # 4 diagnostic plots: residuals, Q-Q, scale-location, leverage
par(mfrow = c(1, 1))

# Confidence intervals for coefficients
confint(model, level = 0.95)

# Predictions
predict(model, newdata = data.frame(wt = 3.0, hp = 100, cyl = 4))
predict(model, newdata = data.frame(wt = 3.0, hp = 100, cyl = 4),
        interval = "confidence")  # confidence interval for mean
predict(model, newdata = data.frame(wt = 3.0, hp = 100, cyl = 4),
        interval = "prediction")  # prediction interval for individual

# Logistic regression with glm()
# Using a binary outcome
mtcars_binary <- mtcars |>
  mutate(high_mpg = as.integer(mpg > 20))

logit_model <- glm(high_mpg ~ wt + hp,
                   data = mtcars_binary,
                   family = binomial(link = "logit"))
summary(logit_model)

# Predicted probabilities
predict(logit_model, type = "response")  # probabilities, not log-odds

# Other GLM families
# Poisson: count data
glm(count ~ x + y, data = count_data, family = poisson(link = "log"))
# Negative binomial: overdispersed counts
library(MASS)
glm.nb(count ~ x + y, data = count_data)
# Gamma: positive continuous with variance proportional to mean squared
glm(cost ~ age + severity, data = claims, family = Gamma(link = "log"))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

library(broom)

model <- lm(mpg ~ wt + hp + factor(cyl), data = mtcars)

# tidy: coefficient-level summary as a tibble
tidy(model)
#   term           estimate std.error statistic  p.value
#   (Intercept)    38.75    1.79      21.6       ...
#   wt             -3.17    0.74      -4.28      ...

tidy(model, conf.int = TRUE)  # adds conf.low, conf.high columns

# glance: model-level summary
glance(model)
#   r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance

# augment: observation-level: adds fitted values, residuals, hat values
augment(model)
#   mpg  wt  hp  .fitted  .resid  .hat  .cooksd  .std.resid

# Coefficient plot: easily visualizable
tidy(model, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  ggplot(aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high)) +
  geom_pointrange() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_minimal()

# Fitting many models with purrr + broom
library(purrr)

# Fit a regression for each carrier in nycflights13
carrier_models <- flights |>
  filter(!is.na(arr_delay), !is.na(dep_delay)) |>
  group_by(carrier) |>
  nest() |>
  mutate(
    model    = map(data, \(d) lm(arr_delay ~ dep_delay + distance, data = d)),
    tidied   = map(model, tidy),
    glanced  = map(model, glance)
  )

# Extract all coefficients in one tibble
carrier_models |>
  unnest(tidied) |>
  filter(term == "dep_delay") |>
  select(carrier, estimate, std.error, p.value)

# Extract model fit statistics
carrier_models |>
  unnest(glanced) |>
  select(carrier, r.squared, adj.r.squared, sigma, AIC)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

library(tidymodels)

# Example: predicting mpg from car features
set.seed(42)
data(mtcars)

# 1. Split data
car_split <- initial_split(mtcars, prop = 0.75)
car_train <- training(car_split)
car_test  <- testing(car_split)

# 2. Define a recipe: preprocessing steps
car_recipe <- recipe(mpg ~ ., data = car_train) |>
  step_center(all_numeric_predictors()) |>    # center to mean 0
  step_scale(all_numeric_predictors()) |>     # scale to sd 1
  step_dummy(all_nominal_predictors())        # one-hot encode factors

# Preview what the recipe does
prep(car_recipe) |> bake(new_data = car_train) |> head()

# 3. Specify a model
# Linear regression
lm_spec <- linear_reg() |>
  set_engine("lm")

# Random forest
rf_spec <- rand_forest(
  trees = 500,
  mtry  = tune(),    # hyperparameter to tune
  min_n = tune()     # hyperparameter to tune
) |>
  set_engine("ranger") |>
  set_mode("regression")

# 4. Build a workflow
car_wflow <- workflow() |>
  add_recipe(car_recipe) |>
  add_model(lm_spec)

# 5. Fit
car_fit <- car_wflow |> fit(data = car_train)

# 6. Evaluate
car_predictions <- car_fit |>
  predict(new_data = car_test) |>
  bind_cols(car_test)

metrics(car_predictions, truth = mpg, estimate = .pred)
#   metric  estimator  estimate
#   rmse    standard   2.89
#   rsq     standard   0.87
#   mae     standard   2.23

# 7. Cross-validation for unbiased evaluation
car_folds <- vfold_cv(car_train, v = 10)

lm_cv_results <- fit_resamples(
  car_wflow,
  resamples = car_folds,
  metrics = metric_set(rmse, rsq, mae)
)

collect_metrics(lm_cv_results)

# 8. Hyperparameter tuning with grid search
rf_wflow <- workflow() |>
  add_recipe(car_recipe) |>
  add_model(rf_spec)

rf_grid <- grid_regular(
  mtry(range = c(2, 7)),
  min_n(range = c(2, 10)),
  levels = 5
)

rf_tuned <- tune_grid(
  rf_wflow,
  resamples = car_folds,
  grid = rf_grid,
  metrics = metric_set(rmse, rsq)
)

show_best(rf_tuned, metric = "rmse")
best_rf <- select_best(rf_tuned, metric = "rmse")

# 9. Finalize and fit on full training data
final_rf_wflow <- rf_wflow |> finalize_workflow(best_rf)
final_rf_fit <- final_rf_wflow |> last_fit(car_split)

# Final test-set performance
collect_metrics(final_rf_fit)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

---
title: "Quarterly Sales Analysis"
author: "Data Team"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    theme: flatly
params:
  quarter: "Q4"
  year: 2025
---

## Overview

This report covers **`r params$quarter` `r params$year`** sales data.

```r
# setup, include=FALSE
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 10,
  fig.height = 6
)
library(tidyverse)
library(knitr)
```

```r
sales <- read_csv("sales_data.csv") |>
  filter(quarter == params$quarter, year == params$year)
```

## Summary Statistics

```r
sales |>
  group_by(region) |>
  summarise(
    revenue = sum(revenue),
    units   = sum(units),
    avg_deal = mean(deal_size)
  ) |>
  kable(
    format.args = list(big.mark = ","),
    digits = 2,
    caption = "Revenue by Region"
  )
```

```r
ggplot(sales, aes(x = month, y = revenue, color = region)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::dollar) +
  theme_minimal() +
  labs(title = paste(params$quarter, params$year, "Revenue by Region"))
```

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# Render a report for each quarter
quarters <- c("Q1", "Q2", "Q3", "Q4")

purrr::walk(quarters, \(q) {
  rmarkdown::render(
    "sales_report.Rmd",
    params   = list(quarter = q, year = 2025),
    output_file = paste0("reports/sales_", q, "_2025.html")
  )
})

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

---
title: "Production Data Analysis"
author: "Jeff Moon"
date: last-modified
format:
  html:
    toc: true
    code-fold: true
    theme: cosmo
  pdf:
    toc: true
    number-sections: true
execute:
  echo: true
  warning: false
  cache: true
---

## Data Loading

```r
#| label: load-data
#| message: false

library(tidyverse)
data <- read_csv("analysis_data.csv")
glimpse(data)

1
2
3
4
5
6
7

#| label: python-analysis

import pandas as pd
import numpy as np

# Access R objects via the r module (Quarto handles the bridge)
# r.data contains the R tibble as a pandas DataFrame

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

library(shiny)

ui <- fluidPage(
  titlePanel("Flights Explorer"),

  sidebarLayout(
    sidebarPanel(
      selectInput(
        inputId = "carrier",
        label   = "Select Carrier:",
        choices = c("AA", "UA", "DL", "B6", "WN"),
        selected = "AA"
      ),
      sliderInput(
        inputId = "delay_threshold",
        label   = "Delay Threshold (minutes):",
        min = 0, max = 120, value = 30, step = 5
      ),
      dateRangeInput(
        inputId = "date_range",
        label   = "Date Range:",
        start   = "2013-01-01",
        end     = "2013-12-31"
      ),
      actionButton("update", "Update Plot", class = "btn-primary")
    ),
    mainPanel(
      plotOutput("delay_plot", height = "400px"),
      hr(),
      tableOutput("summary_table")
    )
  )
)

server <- function(input, output, session) {
  library(nycflights13)
  library(dplyr)
  library(ggplot2)

  # reactive(): a cached computation that reacts to inputs
  # Re-runs only when its inputs change
  filtered_data <- reactive({
    flights |>
      filter(
        carrier == input$carrier,
        !is.na(dep_delay),
        as.Date(time_hour) >= input$date_range[1],
        as.Date(time_hour) <= input$date_range[2]
      )
  })

  # eventReactive(): only re-runs when a specific event fires
  # Here: only recalculate when the button is clicked
  plot_data <- eventReactive(input$update, {
    filtered_data() |>
      mutate(delayed = dep_delay > input$delay_threshold)
  }, ignoreNULL = FALSE)  # run on initial load too

  # renderPlot: reactive output for ggplot2 graphics
  output$delay_plot <- renderPlot({
    ggplot(plot_data(), aes(x = dep_delay, fill = delayed)) +
      geom_histogram(bins = 60, alpha = 0.8) +
      geom_vline(xintercept = input$delay_threshold,
                 color = "red", linetype = "dashed", linewidth = 1) +
      scale_fill_manual(values = c("FALSE" = "steelblue", "TRUE" = "firebrick")) +
      xlim(-30, 180) +
      labs(
        title = paste("Departure Delays —", input$carrier),
        x = "Departure Delay (minutes)",
        y = "Number of Flights",
        fill = paste0("Delayed >", input$delay_threshold, "m")
      ) +
      theme_minimal(base_size = 14)
  })

  # renderTable: reactive output for tabular data
  output$summary_table <- renderTable({
    plot_data() |>
      group_by(origin) |>
      summarise(
        Flights         = n(),
        `Mean Delay`    = round(mean(dep_delay), 1),
        `Pct Delayed`   = paste0(round(mean(delayed) * 100, 1), "%"),
        `Max Delay`     = max(dep_delay),
        .groups = "drop"
      )
  }, striped = TRUE, hover = TRUE)
}

shinyApp(ui, server)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

server <- function(input, output, session) {

  # reactive(): cached value, like a computed property
  # Only recalculates when its reactive dependencies change
  expensive_data <- reactive({
    Sys.sleep(0.5)  # simulate an expensive query
    read_from_db(input$filter_value)
  })

  # observe(): side effects that react to inputs
  # Use for: logging, updating UI elements, writing files
  observe({
    cat("User selected:", input$carrier, "\n")
    # Does not return a value; runs for its side effect
  })

  # observeEvent(): side effect on a specific trigger
  observeEvent(input$download_btn, {
    # Only runs when download_btn is clicked
    write_csv(expensive_data(), "/tmp/export.csv")
    showNotification("File saved!", type = "message")
  })

  # eventReactive(): reactive value triggered by an event
  confirmed_data <- eventReactive(input$go_btn, {
    expensive_data() |> filter(value > input$threshold)
  })

  # renderUI(): dynamically build UI from server side
  output$dynamic_selector <- renderUI({
    choices <- unique(expensive_data()$category)
    selectInput("category", "Category:", choices = choices)
  })

  # Isolate: read a reactive value without creating a dependency
  output$info_text <- renderText({
    # This re-runs when 'carrier' changes, but NOT when 'threshold' changes
    carrier <- input$carrier
    threshold <- isolate(input$threshold)  # read without dependency
    paste("Carrier:", carrier, "| Threshold:", threshold)
  })
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

# Define a module: carrier_stats_ui + carrier_stats_server
carrier_stats_ui <- function(id) {
  ns <- NS(id)  # namespace function: ensures unique IDs
  tagList(
    selectInput(ns("carrier"), "Carrier:", choices = c("AA","UA","DL")),
    plotOutput(ns("plot")),
    tableOutput(ns("table"))
  )
}

carrier_stats_server <- function(id, data) {
  moduleServer(id, function(input, output, session) {
    filtered <- reactive({
      data |> filter(carrier == input$carrier)
    })

    output$plot <- renderPlot({
      ggplot(filtered(), aes(x = dep_delay)) +
        geom_histogram(bins = 40) +
        theme_minimal()
    })

    output$table <- renderTable({
      filtered() |>
        summarise(n = n(), mean_delay = mean(dep_delay, na.rm = TRUE))
    })
  })
}

# Use the module in an app — can instantiate multiple times
ui <- fluidPage(
  carrier_stats_ui("east_coast"),
  carrier_stats_ui("west_coast")
)

server <- function(input, output, session) {
  carrier_stats_server("east_coast", flights)
  carrier_stats_server("west_coast", flights)
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

# shinyapps.io: easiest option, free tier available
library(rsconnect)
rsconnect::deployApp(
  appDir = "myapp/",
  appName = "flights-explorer"
)

# Self-hosted: Shiny Server (open source) via Docker
# docker-compose.yml excerpt:
# shiny:
#   image: rocker/shiny:4.3.1
#   ports: ["3838:3838"]
#   volumes:
#     - ./apps:/srv/shiny-server
#     - ./logs:/var/log/shiny-server

# Testing with shinytest2
library(shinytest2)

test_that("carrier filter works", {
  app <- AppDriver$new(appDir = "myapp/", name = "flights-app")
  app$set_inputs(carrier = "UA")
  app$click("update")
  app$expect_values()     # snapshot test: captures current state
  app$expect_screenshot() # visual regression test

  # Assert specific output values
  delay_data <- app$get_value(output = "summary_table")
  expect_equal(nrow(delay_data), 3)  # 3 origins

  app$stop()
})

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

library(data.table)

# data.table syntax: DT[i, j, by]
# i = filter rows (like WHERE)
# j = select/transform columns (like SELECT)
# by = group by (like GROUP BY)
# This maps directly to SQL semantics

DT <- as.data.table(flights)   # convert data.frame to data.table
DT <- fread("large_file.csv")  # fast CSV read: 10-100x faster than read.csv

# Basic operations
DT[carrier == "AA", .(mean_delay = mean(dep_delay, na.rm = TRUE)), by = origin]

# In-place modification with := (no copy, memory efficient)
DT[, gain := arr_delay - dep_delay]
DT[, c("speed_mph", "on_time") := .(distance / (air_time / 60), dep_delay <= 0)]

# Multiple aggregations
DT[!is.na(dep_delay),
   .(
     n           = .N,
     mean_delay  = mean(dep_delay),
     p95_delay   = quantile(dep_delay, 0.95),
     pct_on_time = mean(dep_delay <= 0)
   ),
   by = .(carrier, origin)][order(-mean_delay)]

# Keys: sorted indexing for fast joins and lookups
setkey(DT, carrier, flight)
DT["AA"]  # all American Airlines flights (fast binary search)

# Joins
planes_DT <- as.data.table(planes)
setkey(DT, tailnum)
setkey(planes_DT, tailnum)
merged <- DT[planes_DT, on = "tailnum"]  # right join
merged <- planes_DT[DT, on = "tailnum"]  # left join

# Rolling joins: useful for as-of lookups
# Match each flight to the most recent weather observation
weather_DT[flights_DT, on = .(station = origin, time_hour),
           roll = TRUE]  # rolling match on time_hour

# Performance comparison reference (100M rows, 8-core machine):
# data.table group-by aggregation: ~1-2 seconds
# dplyr equivalent:                ~8-15 seconds
# pandas equivalent:               ~5-10 seconds
# The gap is largest on complex multi-group aggregations

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

library(Rcpp)

# Inline C++ function
cppFunction('
double sum_of_squares(NumericVector x) {
  double total = 0;
  for (int i = 0; i < x.size(); i++) {
    total += x[i] * x[i];
  }
  return total;
}
')

sum_of_squares(1:1000)  # fast C++ version

# More complex: implementing a custom rolling window function
cppFunction('
NumericVector rolling_mean(NumericVector x, int window) {
  int n = x.size();
  NumericVector result(n, NumericVector::get_na());
  for (int i = window - 1; i < n; i++) {
    double sum = 0;
    for (int j = i - window + 1; j <= i; j++) {
      sum += x[j];
    }
    result[i] = sum / window;
  }
  return result;
}
')

# Source a .cpp file
sourceCpp("fast_functions.cpp")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

library(future)
library(furrr)

# Set up a parallel plan (multicore on Linux/Mac, multisession on Windows)
plan(multisession, workers = 4)

# furrr is purrr with parallel execution — same API, just add future_ prefix
results <- future_map_dfr(
  c("AA", "UA", "DL", "B6", "WN"),
  \(carrier) {
    flights |>
      filter(carrier == !!carrier) |>
      summarise(
        carrier     = carrier,
        n           = n(),
        mean_delay  = mean(dep_delay, na.rm = TRUE)
      )
  },
  .options = furrr_options(seed = TRUE)  # reproducible randomness
)

# For heavy computation (model fitting, bootstrap, permutation tests)
bootstrap_means <- future_map_dbl(
  1:1000,
  \(i) mean(sample(flights$dep_delay, replace = TRUE, na.rm = TRUE))
)

# Reset to sequential
plan(sequential)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

library(arrow)

# Write a parquet file (columnar, efficient, cross-language)
write_parquet(flights, "flights.parquet")
write_parquet(flights, "flights.parquet",
              compression = "snappy")  # or "zstd", "gzip"

# Open without loading into memory — returns an Arrow Dataset
arrow_ds <- open_dataset("flights.parquet")

# Query with dplyr syntax — executed by Arrow, not R
result <- arrow_ds |>
  filter(carrier == "AA", !is.na(dep_delay)) |>
  group_by(origin) |>
  summarise(mean_delay = mean(dep_delay)) |>
  collect()  # only materializes the result into R memory

# Partitioned datasets: crucial for large data
# Write partitioned by year/month (like Hive partitioning)
write_dataset(flights |> mutate(year = 2013, month = month),
              "flights_partitioned/",
              partitioning = c("year", "month"),
              format = "parquet")

# Query a partitioned dataset efficiently
open_dataset("flights_partitioned/") |>
  filter(month == 6) |>  # partition pruning: only reads June data
  collect()

# DuckDB: analytical SQL engine, fast on local data
library(duckdb)
library(dplyr)

con <- dbConnect(duckdb())

# Register a data frame or parquet file as a table
duckdb_register(con, "flights", flights)

# Query with SQL
dbGetQuery(con, "
  SELECT carrier, origin,
         COUNT(*) as n,
         AVG(dep_delay) as mean_delay,
         PERCENTILE_CONT(0.95) WITHIN GROUP (ORDER BY dep_delay) as p95_delay
  FROM flights
  WHERE dep_delay IS NOT NULL
  GROUP BY carrier, origin
  ORDER BY mean_delay DESC
")

# Or use dplyr syntax (executed by DuckDB, not R)
tbl(con, "flights") |>
  filter(!is.na(dep_delay)) |>
  group_by(carrier) |>
  summarise(mean_delay = mean(dep_delay), n = n()) |>
  collect()

# DuckDB can query parquet directly (no loading needed)
dbGetQuery(con, "SELECT * FROM 'flights.parquet' LIMIT 10")
dbGetQuery(con, "SELECT * FROM 'flights_partitioned/**/*.parquet' WHERE month = 6")

dbDisconnect(con)