R Cheatsheet

# assignment operators (all equivalent for simple values)
name <- "Alice"      # preferred
age = 30L            # = also works (avoid in functions)
30L -> age2          # rightward assignment

# basic types (called "modes" in R)
class("text")        # "character"
class(42)            # "numeric" (double)
class(42L)           # "integer"
class(3.14)          # "numeric"
class(TRUE)          # "logical"
class(2 + 3i)        # "complex"

# check and convert types
is.numeric(age)      # TRUE
is.character(name)   # TRUE
as.integer(3.9)      # 3 (truncates, not rounds)
as.character(42)     # "42"
as.numeric("3.14")   # 3.14

# special values
NA                   # missing value (Not Available)
NULL                 # null object (no value)
NaN                  # Not a Number (0/0)
Inf                  # infinity (1/0)
is.na(NA)            # TRUE
is.null(NULL)        # TRUE

# create vectors with c() (combine)
nums <- c(1, 2, 3, 4, 5)
chars <- c("a", "b", "c")
logs <- c(TRUE, FALSE, TRUE)

# sequences
1:5                  # 1 2 3 4 5
seq(1, 10, by = 2)  # 1 3 5 7 9
seq(0, 1, length.out = 5)  # 0 0.25 0.5 0.75 1
rep(1:3, times = 2) # 1 2 3 1 2 3
rep(1:3, each = 2)  # 1 1 2 2 3 3

# indexing (1-indexed!)
nums[1]              # 1 (first element)
nums[length(nums)]   # 5 (last element)
nums[c(1, 3, 5)]     # 1 3 5 (multiple indices)
nums[-1]             # 2 3 4 5 (all EXCEPT first)
nums[-c(1, 2)]       # 3 4 5 (exclude first two)
nums[2:4]            # 2 3 4 (range)

# logical indexing
nums[nums > 2]       # 3 4 5
nums[nums %% 2 == 0] # 2 4 (even numbers)
nums[nums > 2] <- 0  # modify in place

# named vectors
ages <- c(alice = 30, bob = 25, carol = 28)
ages["alice"]        # 30
ages[c("alice", "bob")]  # 30 25

# arithmetic is element-wise (vectorized)
a <- c(1, 2, 3)
b <- c(4, 5, 6)
a + b                # 5 7 9
a * b                # 4 10 18
a ^ 2                # 1 4 9
a / b                # 0.25 0.4 0.5

# recycling: shorter vector repeats
c(1, 2, 3, 4) + c(10, 20)  # 11 22 13 24

# summary functions
sum(a)               # 6
mean(a)              # 2
median(a)            # 2
sd(a)                # 1
var(a)               # 1
min(a); max(a)       # 1; 3
range(a)             # 1 3
cumsum(a)            # 1 3 6
cumprod(a)           # 1 2 6

# sorting and ordering
sort(c(3, 1, 2))              # 1 2 3
sort(c(3, 1, 2), decreasing = TRUE)  # 3 2 1
order(c(3, 1, 2))             # 2 3 1 (indices that would sort)
v[order(v)]                   # sort using order

# useful functions
length(a)            # 3
unique(c(1, 1, 2, 2, 3))  # 1 2 3
rev(a)                # 3 2 1
head(a, 2)            # 1 2 (first n)
tail(a, 2)            # 2 3 (last n)

# paste and paste0 (concatenation)
paste("Hello", "World")           # "Hello World"
paste("Hello", "World", sep = "_") # "Hello_World"
paste0("a", "b", "c")             # "abc" (no separator)
paste(c("a", "b", "c"), collapse = "-")  # "a-b-c"
paste("file", 1:3, ".csv", sep = "")     # "file1.csv" "file2.csv" "file3.csv"

# case conversion
toupper("hello")    # "HELLO"
tolower("WORLD")    # "world"

# substring
substr("Hello World", 1, 5)  # "Hello"
nchar("Hello")               # 5 (character count)

# split and replace
strsplit("a,b,c", ",")[[1]]  # "a" "b" "c"
gsub("o", "0", "Hello World") # "Hell0 W0rld" (all matches)
sub("o", "0", "Hello World")  # "Hell0 World" (first match only)
gsub("[0-9]+", "N", "a1b22c333") # "aNbNcN" (regex)

# grep and grepl (pattern matching)
grep("^A", c("Alice", "Bob", "Anna"))  # 1 3 (indices)
grepl("^A", c("Alice", "Bob"))         # TRUE FALSE

# sprintf (C-style formatting)
sprintf("Pi = %.2f", pi)       # "Pi = 3.14"
sprintf("%s is %d", "Alice", 30) # "Alice is 30"

# NA propagation and handling
x <- c(1, 2, NA, 4, 5)
mean(x)              # NA (NA propagates!)
mean(x, na.rm = TRUE) # 3 (skip NA)
sum(x, na.rm = TRUE)  # 12
is.na(x)             # FALSE FALSE TRUE FALSE FALSE
sum(is.na(x))        # 1 (count of NAs)
x[!is.na(x)]         # 1 2 4 5 (remove NAs)

# NULL vs NA
length(c(1, NA, 3))  # 3 (NA is an element)
length(c(1, NULL, 3)) # 2 (NULL is dropped)

# type coercion hierarchy: logical < integer < double < character
c(TRUE, 1)           # 1 1 (logical -> numeric)
c(1L, 2.5)           # 1.0 2.5 (integer -> double)
c(1, "a")            # "1" "a" (numeric -> character)

# explicit coercion
as.numeric(c("1", "2", "abc"))  # 1 2 NA (with warning)
as.logical(c(0, 1, 2))          # FALSE TRUE TRUE
factor(c("low", "high", "low")) # factor with 2 levels

# ifelse vectorized conditional
ifelse(x > 2, "big", "small")  # "small" "small" NA "big" "big"

# lists can hold different types and sizes
user <- list(
    name = "Alice",
    age = 30,
    scores = c(90, 85, 88),
    active = TRUE
)

# access by name ($ or [[]])
user$name             # "Alice"
user[["age"]]         # 30
user[["scores"]][2]   # 85

# [] returns a sublist (list); [[]] returns the element
user["name"]          # list with one element
user[["name"]]        # "Alice" (the string itself)
user[1:2]             # sublist with first 2 elements

# modify and add
user$age <- 31
user$email <- "[email protected]"  # add new element
user[["scores"]] <- NULL           # remove element

# iterate over a list
for (key in names(user)) {
    cat(key, ":", user[[key]], "\n")
}

# lapply and sapply on lists
lapply(user$scores, sqrt)   # list of square roots
sapply(user$scores, sqrt)   # vector of square roots

# create a data frame (like a table/spreadsheet)
df <- data.frame(
    name = c("Alice", "Bob", "Carol"),
    age = c(30, 25, 28),
    score = c(90, 85, 88),
    stringsAsFactors = FALSE
)

# dimensions
nrow(df)              # 3
ncol(df)              # 3
dim(df)               # 3 3
names(df)             # "name" "age" "score"
str(df)               # structure summary
head(df, 2)           # first 2 rows

# access columns
df$name               # vector (by name)
df[["age"]]           # vector (by name, alternative)
df[, "age"]           # vector (column)
df[, 2]               # vector (by index)
df[2]                 # data frame with one column

# access rows
df[1, ]               # first row (as data frame)
df[1:2, ]             # first 2 rows
df[c(1, 3), ]         # rows 1 and 3

# filter rows (logical indexing)
df[df$age > 26, ]              # rows where age > 26
df[df$score >= 88, c("name", "score")]
subset(df, age > 26, select = c(name, score))

# add and modify columns
df$grade <- c("A", "B", "A")   # add column
df$age <- df$age + 1           # modify column
df <- cbind(df, pass = df$score > 60)  # column bind

# add rows
new_row <- data.frame(name = "Dan", age = 22, score = 75, grade = "C")
df <- rbind(df, new_row)

# create a matrix (2D, all same type)
m <- matrix(1:12, nrow = 3, ncol = 4)
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    7   10
# [2,]    2    5    8   11
# [3,]    3    6    9   12

# fill by row instead of column
m2 <- matrix(1:6, nrow = 2, byrow = TRUE)

# dimensions and attributes
dim(m)                # 3 4
nrow(m)               # 3
ncol(m)               # 4
rownames(m) <- c("r1", "r2", "r3")
colnames(m) <- c("c1", "c2", "c3", "c4")

# indexing
m[2, 3]               # 8 (single element)
m[1, ]                # 1 4 7 10 (first row)
m[, 2]                # 4 5 6 (second column)
m[1:2, 2:3]           # 2x2 submatrix
m["r1", "c2"]         # 4 (by name)

# matrix operations
t(m)                  # transpose
m * m                 # element-wise multiply
m %*% t(m)            # matrix multiplication
solve(m[, 1:3])       # inverse (square matrix)
det(matrix(1:4, 2))   # determinant
rowSums(m)            # sum of each row
colMeans(m)           # mean of each column
apply(m, 1, sum)      # row sums (general)
apply(m, 2, max)      # column maxima

# arrays (n-dimensional)
arr <- array(1:24, dim = c(2, 3, 4))  # 2x3x4 array
arr[1, 2, 3]          # single element

# factors represent categorical data (stored as integers with labels)
gender <- factor(c("male", "female", "male", "female"))
print(gender)
# [1] male   female male   female
# Levels: female male

levels(gender)        # "female" "male" (sorted alphabetically)
table(gender)         # frequency table
nlevels(gender)       # 2

# ordered factor (for ordinal data)
size <- factor(c("M", "S", "L", "XL"),
               levels = c("S", "M", "L", "XL"),
               ordered = TRUE)
size[1] > size[2]     # TRUE (M > S)

# convert to factor and back
nums <- factor(c(1, 2, 3, 2, 1))
as.numeric(nums)      # 1 2 3 2 1 (level codes, NOT original values!)
as.numeric(as.character(nums))  # 1 2 3 2 1 (correct way)

# relevel (change reference level)
gender <- relevel(gender, ref = "male")

# cut: convert numeric to factor (binning)
ages <- c(15, 25, 35, 45, 55, 65)
age_groups <- cut(ages, breaks = c(0, 18, 35, 50, 100),
                  labels = c("child", "young", "adult", "senior"))
table(age_groups)     # frequency per group

# interaction of factors
f1 <- factor(c("A", "A", "B", "B"))
f2 <- factor(c("X", "Y", "X", "Y"))
interaction(f1, f2)   # A.X A.Y B.X B.Y

# convert list to data frame
my_list <- list(
    a = 1:3,
    b = c("x", "y", "z")
)
df <- as.data.frame(my_list)

# stack multiple vectors
do.call(rbind, list(
    data.frame(name = "A", val = 1),
    data.frame(name = "B", val = 2)
))

# reshape: wide to long
library(tidyr)
wide_df <- data.frame(
    id = 1:2,
    q1 = c(10, 20),
    q2 = c(15, 25)
)
long_df <- pivot_longer(wide_df, cols = starts_with("q"),
                        names_to = "quarter", values_to = "value")

# long to wide
pivot_wider(long_df, names_from = quarter, values_from = value)

# base R reshape (without tidyr)
reshape(wide_df, direction = "long",
        varying = c("q1", "q2"), v.names = "value",
        timevar = "quarter", idvar = "id")

# split and combine
split_results <- split(df, df$group)  # list of data frames by group
combined <- do.call(rbind, split_results)  # recombine

# melt and cast (reshape2 package, legacy)
# library(reshape2)
# melt(wide_df, id.vars = "id")
# dcast(long_df, id ~ quarter, value.var = "value")

# if-else if-else
score <- 85
if (score >= 90) {
    grade <- "A"
} else if (score >= 80) {
    grade <- "B"
} else {
    grade <- "C"
}

# ifelse (vectorized ternary)
ages <- c(15, 25, 35, 45)
ifelse(ages >= 18, "adult", "minor")
# "minor" "adult" "adult" "adult"

# nested ifelse (avoid deep nesting!)
ifelse(ages < 18, "minor",
       ifelse(ages < 65, "adult", "senior"))

# dplyr::case_when is cleaner for multiple conditions
# case_when(
#     ages < 18 ~ "minor",
#     ages < 65 ~ "adult",
#     TRUE      ~ "senior"
# )

# switch (dispatch on a value)
day_type <- function(day) {
    switch(day,
        "Mon" = "weekday",
        "Tue" = "weekday",
        "Wed" = "weekday",
        "Thu" = "weekday",
        "Fri" = "weekday",
        "Sat" = "weekend",
        "Sun" = "weekend",
        "unknown"
    )
}

# for loop
for (i in 1:5) {
    print(i)
}

# iterate over a vector
fruits <- c("apple", "banana", "cherry")
for (fruit in fruits) {
    print(paste("Fruit:", fruit))
}

# iterate with index
for (i in seq_along(fruits)) {
    print(paste(i, fruits[i]))
}

# while loop
count <- 0
while (count < 5) {
    count <- count + 1
    if (count == 3) next    # skip (like continue)
    print(count)
}

# repeat loop (infinite, must break)
i <- 0
repeat {
    i <- i + 1
    if (i >= 3) break       # exit loop
    print(i)
}

# preallocate for speed (IMPORTANT!)
result <- numeric(1000)
for (i in 1:1000) {
    result[i] <- i^2
}

# nested loops
mat <- matrix(0, 3, 3)
for (i in 1:3) {
    for (j in 1:3) {
        mat[i, j] <- i * j
    }
}

# basic function (last expression is returned)
add <- function(a, b) {
    a + b
}
add(3, 4)            # 7

# explicit return
is_positive <- function(x) {
    if (x > 0) return(TRUE)
    FALSE
}

# default arguments
greet <- function(name, greeting = "Hello", punctuation = "!") {
    paste0(greeting, ", ", name, punctuation)
}
greet("Alice")                    # "Hello, Alice!"
greet("Bob", greeting = "Hi")     # "Hi, Bob!"
greet(name = "Carol", punctuation = "?")  # named args

# ... (ellipsis: pass arguments through)
my_plot <- function(x, y, ...) {
    plot(x, y, col = "blue", ...)
}
my_plot(1:10, 1:10, main = "My Plot", type = "l")

# return multiple values via list
stats <- function(x) {
    list(mean = mean(x), sd = sd(x), n = length(x))
}
result <- stats(1:10)
result$mean         # 5.5

# lazy evaluation (args evaluated only when used)
f <- function(a, b) {
    a * 2    # b is never used, so not evaluated
}
f(5)                # 10 (no error despite missing b)

# lapply: apply function to each element of a list, returns list
my_list <- list(a = 1:3, b = 4:6, c = 7:9)
lapply(my_list, mean)    # list: 2, 5, 8
lapply(my_list, sum)     # list: 6, 15, 24

# sapply: like lapply but simplifies to vector/matrix
sapply(my_list, mean)    # 2 5 8 (named vector)
sapply(my_list, range)   # matrix with 2 rows

# vapply: like sapply but with guaranteed output type (safer)
vapply(my_list, mean, numeric(1))  # always numeric vector

# apply: apply function over matrix/array margins
m <- matrix(1:12, nrow = 3)
apply(m, 1, sum)         # row sums: 22 26 30
apply(m, 2, mean)        # column means: 2 5 8 11
apply(m, c(1, 2), sqrt)  # element-wise sqrt

# mapply: multivariate apply (vectorized over multiple args)
mapply(rep, 1:3, 3:1)    # rep(1,3), rep(2,2), rep(3,1)

# tapply: apply function by groups
df <- data.frame(
    group = c("A", "A", "B", "B", "B"),
    value = c(10, 20, 30, 40, 50)
)
tapply(df$value, df$group, mean)  # A: 15, B: 40

# replicate: repeat an expression n times
replicate(3, mean(rnorm(10)))  # 3 random means

# lexical scoping: functions look up variables in defining environment
y <- 10
f <- function(x) {
    x + y    # y found in global env
}
f(5)                 # 15

# local variables shadow globals
g <- function(x) {
    y <- 100         # local y
    x + y
}
g(5)                 # 105
y                    # 10 (global unchanged)

# <<- assigns to parent environment (super-assignment)
counter <- function() {
    count <- 0
    function() {
        count <<- count + 1   # modifies count in enclosing scope
        count
    }
}
c1 <- counter()
c1()                 # 1
c1()                 # 2

# search path for variables
search()             # shows environments: globalenv, package namespaces
ls()                 # list objects in current environment
ls(envir = .GlobalEnv)  # explicitly global
exists("y")          # TRUE if variable exists

# new environment (isolated scope)
e <- new.env()
e$x <- 42
e$x                  # 42 (separate from global x)

library(dplyr)

df <- data.frame(
    name = c("Alice", "Bob", "Carol", "Dan"),
    age = c(30, 25, 28, 35),
    dept = c("Eng", "Sales", "Eng", "Sales"),
    salary = c(80000, 50000, 75000, 90000)
)

# filter rows (like WHERE in SQL)
filter(df, age > 26)
filter(df, dept == "Eng" & salary > 70000)
filter(df, dept %in% c("Eng", "Sales"))
filter(df, age > 26 | salary > 85000)

# select columns (like SELECT in SQL)
select(df, name, age)
select(df, name:dept)           # range of columns
select(df, -salary)             # exclude column
select(df, starts_with("s"))    # columns starting with 's'
select(df, ends_with("e"))      # columns ending with 'e'
select(df, contains("am"))      # columns containing 'am'
select(df, everything())        # all columns (useful for reordering)

# arrange (sort, like ORDER BY)
arrange(df, age)                # ascending
arrange(df, desc(age))          # descending
arrange(df, dept, desc(salary)) # sort by dept, then salary desc

# rename columns
rename(df, years = age, department = dept)

# distinct rows
distinct(df, dept)              # unique departments
distinct(df, dept, .keep_all = TRUE)  # keep all columns

library(dplyr)

# mutate: add/modify columns
df %>%
    mutate(
        bonus = salary * 0.1,
        total = salary + bonus,
        category = ifelse(age > 30, "senior", "junior")
    )

# transmute: like mutate but only keeps new columns
transmute(df, name, annual = salary, monthly = salary / 12)

# summarize (reduces to single row per group)
summarize(df,
    avg_salary = mean(salary),
    max_age = max(age),
    n = n()
)

# group_by + summarize (like GROUP BY in SQL)
df %>%
    group_by(dept) %>%
    summarize(
        count = n(),
        avg_salary = mean(salary),
        avg_age = mean(age)
    ) %>%
    arrange(desc(avg_salary))

# multiple summaries
df %>%
    group_by(dept) %>%
    summarize(across(everything(), list(mean, sd)))

# count and tally
count(df, dept)                  # count per dept
df %>% group_by(dept) %>% tally()

# window functions (within groups)
df %>%
    group_by(dept) %>%
    mutate(
        rank = rank(desc(salary)),
        cumsum = cumsum(salary)
    ) %>%
    arrange(dept, rank)

library(magrittr)  # or library(dplyr)

# without pipe: nested, hard to read
result <- arrange(
    filter(
        select(df, name, age, salary),
        age > 25
    ),
    desc(salary)
)

# with pipe: linear, readable
result <- df %>%
    select(name, age, salary) %>%
    filter(age > 25) %>%
    arrange(desc(salary))

# the dot (.) refers to the piped data
df %>% plot(.$age, .$salary)     # . = df

# %$% exposes column names (magrittr)
df %$% cor(age, salary)          # correlation

# %<>% assigns result back (compound assignment)
df %<>% filter(age > 20)

# T pipe: returns first argument (for side effects)
rnorm(100) %T>%
    hist() %>%                    # plot histogram (returns input)
    mean()                        # then compute mean

# native pipe (R 4.1+): |>
df |>
    subset(age > 25) |>
    colMeans()

library(dplyr)

employees <- data.frame(
    id = c(1, 2, 3, 4),
    name = c("Alice", "Bob", "Carol", "Dan")
)

salaries <- data.frame(
    id = c(1, 2, 3, 5),
    salary = c(80000, 50000, 75000, 60000)
)

# inner join: only matching rows
inner_join(employees, salaries, by = "id")
#   id  name salary
#   1  Alice 80000
#   2  Bob   50000
#   3  Carol 75000

# left join: all rows from left, NA for non-matches
left_join(employees, salaries, by = "id")
# Dan has NA salary

# right join: all rows from right
right_join(employees, salaries, by = "id")
# id 5 has NA name

# full join: all rows from both
full_join(employees, salaries, by = "id")

# different column names
left_join(employees, salaries, by = c("id" = "id"))

# semi join: rows in left that have a match (no right columns)
semi_join(employees, salaries, by = "id")

# anti join: rows in left with NO match
anti_join(employees, salaries, by = "id")
# Dan (id 4 has no salary)

# multiple keys
left_join(df1, df2, by = c("dept" = "department", "year" = "yr"))

library(tidyr)

# wide to long
wide_df <- data.frame(
    id = 1:2,
    q1_sales = c(100, 200),
    q2_sales = c(150, 250),
    q3_sales = c(120, 220)
)

long_df <- wide_df %>%
    pivot_longer(
        cols = starts_with("q"),
        names_to = "quarter",
        values_to = "sales"
    )

# long to wide
long_df %>%
    pivot_wider(
        names_from = quarter,
        values_from = sales
    )

# separate one column into multiple
df <- data.frame(id = 1, name = "Alice_Smith")
df %>% separate(name, into = c("first", "last"), sep = "_")

# unite multiple columns into one
df %>% unite("full_name", first, last, sep = " ")

# separate_rows: split delimited values into rows
df <- data.frame(id = 1:2, tags = c("a,b,c", "x,y"))
df %>% separate_rows(tags, sep = ",")

# drop NA values
df %>% drop_na()                    # drop rows with any NA
df %>% drop_na(salary)              # drop rows where salary is NA

# replace NA with a value
df %>% replace_na(list(salary = 0, name = "Unknown"))

# fill missing values (carry forward)
df %>% fill(salary, .direction = "down")

data <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

# central tendency
mean(data)           # 5.5
median(data)         # 5.5
# R has no built-in mode; use:
as.numeric(names(sort(table(data), decreasing = TRUE)[1]))

# spread
sd(data)             # 3.028 (sample standard deviation)
var(data)            # 9.167 (sample variance)
IQR(data)            # 5 (interquartile range)
range(data)          # 1 10
diff(range(data))    # 9

# quantiles
quantile(data)              # 0% 25% 50% 75% 100%
quantile(data, c(0.1, 0.9)) # 10th and 90th percentiles

# summary (all at once)
summary(data)
# Min. 1st Qu. Median  Mean  3rd Qu.  Max.
# 1.00   3.25    5.50   5.50   7.75    10.00

# correlation and covariance
x <- c(1, 2, 3, 4, 5)
y <- c(2, 4, 5, 4, 5)
cor(x, y)            # 0.77 (Pearson correlation)
cor(x, y, method = "spearman")  # rank correlation
cov(x, y)            # covariance

# handling NA
data_na <- c(1, 2, NA, 4, 5)
mean(data_na, na.rm = TRUE)  # 3
sd(data_na, na.rm = TRUE)    # 1.826

# R uses d/p/q/r prefix for distributions:
#   d = density (PDF/PMF)
#   p = cumulative (CDF: P(X <= x))
#   q = quantile (inverse CDF)
#   r = random samples

# Normal distribution
dnorm(0)                    # 0.399 (density at 0)
pnorm(1.96)                 # 0.975 (P(Z <= 1.96))
qnorm(0.975)                # 1.96 (97.5th percentile)
rnorm(5, mean = 100, sd = 15)  # 5 random samples

# Binomial distribution
dbinom(3, size = 10, prob = 0.5)  # P(X=3)
pbinom(3, size = 10, prob = 0.5)  # P(X<=3)
rbinom(10, size = 1, prob = 0.5)  # 10 coin flips

# Poisson distribution
dpois(2, lambda = 3)        # P(X=2) with rate 3
rpois(100, lambda = 5)      # 100 random samples

# Uniform distribution
runif(5, min = 0, max = 1)  # 5 uniform random numbers

# Exponential
rexp(10, rate = 0.5)        # 10 exponential samples

# set random seed for reproducibility
set.seed(42)
rnorm(3)                    # same results every time with same seed

# sample from a vector
sample(1:10, 5)             # 5 random numbers without replacement
sample(1:10, 5, replace = TRUE)  # with replacement
sample(c("A", "B", "C"), 2) # sample from categories

# one-sample t-test (is mean different from hypothesized?)
data <- c(5.1, 4.9, 5.0, 5.2, 4.8, 5.1, 5.0)
t.test(data, mu = 5.0)
#  t = 0.527, df = 6, p-value = 0.616
#  (fail to reject H0: mean = 5)

# two-sample t-test (are two means different?)
group1 <- c(5.1, 4.9, 5.0, 5.2)
group2 <- c(4.5, 4.7, 4.6, 4.4)
t.test(group1, group2)
t.test(group1, group2, var.equal = TRUE)  # assume equal variance

# paired t-test (before/after)
before <- c(70, 80, 65, 90, 75)
after  <- c(75, 85, 70, 92, 80)
t.test(before, after, paired = TRUE)

# chi-squared test (independence of categorical vars)
tbl <- table(c("A","A","B","B"), c("X","Y","X","Y"))
chisq.test(tbl)

# Wilcoxon (non-parametric alternative to t-test)
wilcox.test(group1, group2)

# ANOVA (compare means across multiple groups)
df <- data.frame(
    value = c(5, 6, 7, 8, 9, 10),
    group = factor(c("A","A","B","B","C","C"))
)
result <- aov(value ~ group, data = df)
summary(result)         # F-test p-value

# extract p-value from any test
test_result <- t.test(data, mu = 5)
test_result$p.value     # 0.616

# simple linear regression: y ~ x
x <- c(1, 2, 3, 4, 5)
y <- c(2, 4, 5, 4, 5)
model <- lm(y ~ x)

# view results
print(model)          # coefficients
summary(model)        # full summary with R², t-tests, F-test
coef(model)           # (Intercept) 2.2, x 0.6
fitted(model)         # predicted values
residuals(model)      # y - predicted
confint(model)        # 95% confidence intervals

# make predictions
predict(model, newdata = data.frame(x = c(6, 7)))

# multiple regression
df <- data.frame(
    y = c(10, 12, 15, 18, 20),
    x1 = c(1, 2, 3, 4, 5),
    x2 = c(5, 4, 3, 2, 1)
)
model2 <- lm(y ~ x1 + x2, data = df)
model3 <- lm(y ~ ., data = df)          # all predictors
model4 <- lm(y ~ x1 + I(x1^2), data = df)  # polynomial

# interactions
model5 <- lm(y ~ x1 * x2, data = df)    # x1 + x2 + x1:x2

# diagnostic plots
par(mfrow = c(2, 2))
plot(model)            # 4 diagnostic plots

# with factors (automatic dummy variables)
model6 <- lm(y ~ x1 + factor(group), data = df)

# logistic regression (binary outcome)
df <- data.frame(
    admit = c(0, 1, 0, 1, 1, 0),
    gre = c(380, 660, 800, 640, 520, 760),
    gpa = c(3.6, 3.7, 3.8, 3.9, 3.5, 3.6)
)
logit <- glm(admit ~ gre + gpa, data = df, family = binomial)
summary(logit)

# predict probabilities (type = "response")
predict(logit, newdata = df, type = "response")

# Poisson regression (count data)
counts <- c(2, 3, 1, 4, 2, 5)
pois <- glm(counts ~ x, family = poisson)

# stepwise model selection
full_model <- lm(y ~ x1 + x2 + x3, data = df)
step_model <- step(full_model)    # AIC-based selection

# anova to compare models
model1 <- lm(y ~ x1, data = df)
model2 <- lm(y ~ x1 + x2, data = df)
anova(model1, model2)             # is x2 significant?

# other models
# nls()  - nonlinear least squares
# glm.nb() - negative binomial (MASS package)
# lme()/lmer() - mixed effects (nlme/lme4)
# rpart()/randomForest() - tree-based methods
# coxph() - survival analysis (survival package)

# cross-validation
library(boot)
cv_result <- cv.glm(df, logit, K = 10)  # 10-fold CV
cv_result$delta    # cross-validation error

# scatter plot
x <- 1:10
y <- x^2
plot(x, y, type = "p",        # p=points, l=lines, b=both, o=overplotted
     main = "Quadratic",       # title
     xlab = "x", ylab = "y",   # axis labels
     col = "blue", pch = 16,   # color and point character
     xlim = c(0, 10), ylim = c(0, 100))

# add lines and points to existing plot
lines(x, x*10, col = "red", lty = 2)  # dashed line
points(x, y + 5, col = "green")
legend("topleft", legend = c("x^2", "10x"),
       col = c("blue", "red"), lty = c(NA, 2), pch = c(16, NA))

# histogram
data <- rnorm(1000)
hist(data, breaks = 30, col = "lightblue",
     main = "Normal Distribution", xlab = "Value", freq = FALSE)
lines(density(data), col = "red", lwd = 2)  # add density curve

# boxplot (compare groups)
boxplot(count ~ spray, data = InsectSprays,
        col = "lightgreen", main = "Insect Count by Spray")

# barplot
counts <- table(mtcars$cyl)
barplot(counts, col = c("red","green","blue"),
        main = "Cars by Cylinders", xlab = "Cylinders")

# multiple plots in one window
par(mfrow = c(2, 2))   # 2x2 grid
plot(x, y); hist(data); boxplot(data); plot(density(data))
par(mfrow = c(1, 1))   # reset

library(ggplot2)

# basic structure: data + aesthetic mapping + geometry
ggplot(data = mtcars, aes(x = wt, y = mpg)) +
    geom_point()

# add aesthetics (color, size, shape mapped to data)
ggplot(mtcars, aes(x = wt, y = mpg, color = cyl, size = hp)) +
    geom_point()

# add layers
ggplot(mtcars, aes(x = wt, y = mpg)) +
    geom_point(color = "blue", size = 3) +
    geom_smooth(method = "lm", se = TRUE) +  # regression line
    labs(title = "MPG vs Weight",
         x = "Weight (1000 lbs)", y = "MPG",
         color = "Cylinders") +
    theme_minimal()

# the + operator adds layers (like %>% but for ggplot)
p <- ggplot(mtcars, aes(x = wt, y = mpg)) +
    geom_point()
p + geom_smooth()           # add to saved plot
p + facet_wrap(~ cyl)       # facet by cylinders

# key components:
#   data       - the data frame
#   aes()      - aesthetic mappings (x, y, color, size, shape)
#   geom_*()   - geometry (point, line, bar, histogram, etc.)
#   scale_*()  - control axes, colors, legends
#   facet_*()  - small multiples (subplots)
#   theme_*()  - overall appearance
#   labs()     - titles and labels

library(ggplot2)

# bar chart (counts)
ggplot(mpg, aes(class)) +
    geom_bar(fill = "steelblue") +
    coord_flip()              # horizontal bars

# histogram
ggplot(mpg, aes(hwy)) +
    geom_histogram(binwidth = 2, fill = "orange", color = "black")

# density plot
ggplot(mpg, aes(hwy, fill = class)) +
    geom_density(alpha = 0.5)   # semi-transparent

# boxplot by group
ggplot(mpg, aes(class, hwy)) +
    geom_boxplot() +
    geom_jitter(width = 0.2, alpha = 0.5)  # overlay points

# line plot (time series)
ggplot(economics, aes(date, unemploy)) +
    geom_line(color = "red") +
    geom_area(fill = "pink", alpha = 0.3)

# scatter with smoothing
ggplot(mpg, aes(displ, hwy, color = drv)) +
    geom_point(size = 2) +
    geom_smooth(method = "loess", se = FALSE)

# error bars
ggplot(df, aes(group, value)) +
    stat_summary(fun = mean, geom = "bar") +
    stat_summary(fun.data = mean_se, geom = "errorbar")

# text labels
ggplot(mtcars, aes(wt, mpg, label = rownames(mtcars))) +
    geom_text(size = 3)

library(ggplot2)

# facets: small multiples (subplots by a variable)
ggplot(mpg, aes(displ, hwy)) +
    geom_point() +
    facet_wrap(~ class)            # one panel per class

# facet_grid: 2D grid of panels
ggplot(mpg, aes(displ, hwy)) +
    geom_point() +
    facet_grid(drv ~ cyl)          # rows ~ cols

# scales: control axes and colors
ggplot(mpg, aes(displ, hwy, color = class)) +
    geom_point() +
    scale_x_log10() +              # log scale
    scale_color_brewer(palette = "Set1") +  # color palette
    scale_size_continuous(range = c(2, 8))

# manual colors
ggplot(mpg, aes(class, fill = drv)) +
    geom_bar() +
    scale_fill_manual(values = c("red", "green", "blue"))

# themes: overall appearance
ggplot(mpg, aes(displ, hwy)) +
    geom_point() +
    theme_minimal() +              # or theme_classic, theme_bw
    theme(
        text = element_text(size = 14, family = "serif"),
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom",
        panel.grid.major = element_line(color = "gray90")
    )

# save plot
ggsave("plot.png", width = 8, height = 6, dpi = 300)
ggsave("plot.pdf", device = "pdf")

# save single object to RDS (binary, preserves type)
saveRDS(my_data, "data.rds")
restored <- readRDS("data.rds")

# save multiple objects to RData
save(df1, df2, model, file = "workspace.RData")
load("workspace.RData")    # loads all objects into workspace

# CSV (most portable)
write.csv(df, "data.csv", row.names = FALSE)
df <- read.csv("data.csv")
df <- read.csv("data.csv", stringsAsFactors = FALSE)
df <- read.csv("data.csv", na.strings = c("", "NA", "N/A"))

# readr (tidyverse): faster, smarter CSV I/O
library(readr)
write_csv(df, "data.csv")
df <- read_csv("data.csv")           # returns a tibble
df <- read_csv("data.csv", col_types = cols(
    age = col_integer(),
    name = col_character()
))

# Excel (requires readxl/writexl packages)
library(readxl)
df <- read_excel("data.xlsx", sheet = 1, range = "A1:D100")
library(writexl)
write_xlsx(df, "output.xlsx")

# RDS vs RData vs CSV
# RDS:   one object, binary, fast, preserves types
# RData: multiple objects, binary, fast
# CSV:   text, portable to other tools, loses type info

# save and load workspace
save.image("project.RData")  # save everything
load("project.RData")        # restore everything

# serialize to JSON (jsonlite package)
library(jsonlite)
write_json(df, "data.json")
df <- fromJSON("data.json")

# S3 is R's simplest OOP system: a class is just an attribute
# create an object: list + class attribute
account <- list(owner = "Alice", balance = 1000)
class(account) <- "BankAccount"

# generic function dispatches by class
print.BankAccount <- function(x, ...) {
    cat("Account owner:", x$owner, "\n")
    cat("Balance: $", x$balance, "\n", sep = "")
}
print(account)
# Account owner: Alice
# Balance: $1000

# define a custom method for an existing generic
deposit <- function(obj, amount) {
    UseMethod("deposit")   # dispatch based on class of obj
}
deposit.BankAccount <- function(obj, amount) {
    obj$balance <- obj$balance + amount
    obj
}
account <- deposit(account, 500)
account$balance   # 1500

# check available methods for a generic
methods(print)         # all print methods
methods(class = "lm")  # methods for lm objects

# default method (fallback)
deposit.default <- function(obj, amount) {
    stop("No deposit method for this class")
}

# S4 is formal: classes are defined with setClass and validated
setClass("Person",
    slots = list(
        name = "character",
        age  = "numeric"
    ),
    prototype = list(name = NA_character_, age = NA_real_)
)

# constructor: new(ClassName, ...)
alice <- new("Person", name = "Alice", age = 30)

# access slots with @ (not $)
alice@name       # "Alice"
slot(alice, "age")  # 30

# define a method
setMethod("show", "Person", function(object) {
    cat("Person:", object@name, "(", object@age, "y)\n")
})
show(alice)   # Person: Alice ( 30 y)

# validity check
setValidity("Person", function(object) {
    if (object@age < 0) return("age must be non-negative")
    TRUE
})
# new("Person", name="Bob", age=-5)  # error

# inheritance
setClass("Student", contains = "Person",
    slots = list(gpa = "numeric"))
stu <- new("Student", name="Bob", age=20, gpa=3.8)

# a closure is a function that remembers its enclosing environment
# function factory: creates functions with private state
make_counter <- function() {
    count <- 0
    function() {
        count <<- count + 1   # <<- modifies in enclosing env
        count
    }
}

c1 <- make_counter()
c1()   # 1
c1()   # 2
c1()   # 3
c2 <- make_counter()
c2()   # 1 (independent counter)

# memoization (cache expensive results)
memoize <- function(f) {
    cache <- new.env(parent = emptyenv())
    function(x) {
        key <- as.character(x)
        if (!exists(key, envir = cache)) {
            assign(key, f(x), envir = cache)
        }
        get(key, envir = cache)
    }
}
slow_sqrt <- function(x) { Sys.sleep(0.1); sqrt(x) }
fast_sqrt <- memoize(slow_sqrt)
fast_sqrt(16)  # slow first time
fast_sqrt(16)  # instant second time

# capturing state in a list (module pattern)
bank_account <- function(initial) {
    balance <- initial
    list(
        deposit = function(amt) { balance <<- balance + amt },
        withdraw = function(amt) { balance <<- balance - amt },
        get_balance = function() balance
    )
}
acc <- bank_account(100)
acc$deposit(50)
acc$get_balance()   # 150

# Map: apply a function to each element (returns list)
Map(function(x) x^2, 1:5)
# [[1]] 1 [[2]] 4 [[3]] 9 [[4]] 16 [[5]] 25

# Reduce: combine elements with a binary function
Reduce("+", 1:5)            # 15 (1+2+3+4+5)
Reduce("*", 1:5)            # 120 (factorial 5!)
Reduce(c, list(1:2, 3:4, 5:6))  # 1 2 3 4 5 6

# Filter: keep elements satisfying a predicate
Filter(function(x) x > 2, 1:5)   # 3 4 5
Filter(is.numeric, list(1, "a", 2, "b"))  # 1 2

# Negate: invert a predicate
Negate(is.null)
is_not_null <- Negate(is.null)

# Position: find first index satisfying a predicate
Position(function(x) x > 3, 1:10)  # 4

# Find: first element satisfying a predicate
Find(function(x) x > 3, 1:10)      # 4

# purrr (tidyverse): cleaner functional programming
library(purrr)
map_dbl(1:5, ~ .x^2)               # 1 4 9 16 25 (vector output)
map_chr(c("a","b"), toupper)       # "A" "B"
keep(1:5, ~ .x > 2)                # 3 4 5
reduce(1:5, `+`)                   # 15
walk(1:3, print)                   # side effects only

# browser(): pause execution and inspect
f <- function(x) {
    browser()         # pauses here; use n, c, Q, ls, print
    y <- x * 2
    if (y > 10) stop("y too big")
    y
}
# f(6)  # enters browser at the browser() line

# debug() / undebug(): debug a function
debug(lm)
# lm(y ~ x)  # steps through lm line by line
undebug(lm)

# traceback(): see call stack after an error
# log("abc")   # error
# traceback()  # shows the call stack

# tryCatch(): handle errors gracefully
safe_log <- function(x) {
    tryCatch(
        log(x),
        error = function(e) {
            message("Error: ", conditionMessage(e))
            NA_real_
        },
        warning = function(w) {
            message("Warning: ", conditionMessage(w))
            log(x)   # retry or handle
        }
    )
}
safe_log("abc")   # NA with message
safe_log(-1)      # NA with warning (NaN)

# withCallingHandlers(): handle warnings without stopping
withCallingHandlers(
    { warn("oops"); 42 },
    warning = function(w) { message("caught: ", w$message); invokeRestart("muffleWarning") }
)

# options for debugging
options(warn = 2)   # warnings become errors (for debugging)
options(error = browser)  # enter browser on error
options(warn = 0)   # reset

# ts(): create a regular time series
# frequency: 12=monthly, 4=quarterly, 52=weekly, 365.25=daily
monthly_sales <- ts(c(30, 35, 40, 38, 42, 45, 50, 48, 52, 55, 60, 58),
                    start = c(2020, 1), frequency = 12)
print(monthly_sales)
#      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
# 2020  30  35  40  38  42  45  50  48  52  55  60  58

# attributes
start(monthly_sales)    # 2020 1
end(monthly_sales)      # 2020 12
frequency(monthly_sales)  # 12
time(monthly_sales)     # time points

# window(): subset a time range
q1 <- window(monthly_sales, start = c(2020, 3), end = c(2020, 5))

# multiple series (multivariate)
ts1 <- ts(rnorm(12), start = 2020, frequency = 12)
ts2 <- ts(rnorm(12), start = 2020, frequency = 12)
mts <- cbind(ts1, ts2)   # multivariate ts

# lag and diff
lag(monthly_sales, k = 1)   # shift by 1 period
diff(monthly_sales)         # first differences (month-over-month change)
diff(monthly_sales, lag = 12)  # year-over-year change

# aggregate to lower frequency
aggregate(monthly_sales, nfrequency = 4, FUN = mean)  # monthly -> quarterly

# decompose(): split a series into trend/seasonal/random
decomp <- decompose(monthly_sales, type = "additive")
# type = "additive" (T+S+R) or "multiplicative" (T*S*R)
plot(decomp)   # shows observed, trend, seasonal, random panels
decomp$trend      # trend component
decomp$seasonal   # seasonal component
decomp$random     # remainder

# stl(): Seasonal-Trend-Loess decomposition (more flexible)
stl_fit <- stl(monthly_sales, s.window = "periodic")
plot(stl_fit)
stl_fit$time.series[, "trend"]     # extract trend

# moving average smoothing
ma3 <- filter(monthly_sales, filter = rep(1/3, 3), sides = 2)
ma12 <- filter(monthly_sales, filter = rep(1/12, 12), sides = 1)

# HoltWinters(): exponential smoothing
hw <- HoltWinters(monthly_sales, seasonal = "additive")
plot(hw)
hw$fitted       # fitted values
hw$coefficients # alpha, beta, gamma

# forecast with HoltWinters
library(forecast)
fc <- forecast(hw, h = 12)   # 12-step ahead forecast
plot(fc)                      # with prediction intervals
autoplot(fc)                  # ggplot2 version

# ACF (autocorrelation) and PACF (partial autocorrelation)
acf(monthly_sales, main = "ACF")     # correlation with lagged self
pacf(monthly_sales, main = "PACF")   # direct correlation at each lag

# interpret:
#   ACF tails off slowly  -> need differencing (non-stationary)
#   ACF cuts off at lag p -> AR(p) process
#   PACF cuts off at lag q -> MA(q) process

# stationarity test
library(tseries)
adf.test(monthly_sales)   # Augmented Dickey-Fuller
# p < 0.05 -> stationary

# auto.arima(): automatically select ARIMA order
library(forecast)
fit <- auto.arima(monthly_sales)
summary(fit)              # shows ARIMA(p,d,q)(P,D,Q)[m]
# ARIMA(0,1,1)(1,0,0)[12]  example output

# forecast
fc <- forecast(fit, h = 12)
plot(fc)
accuracy(fc)              # ME, RMSE, MAE, MAPE

# manual ARIMA
fit2 <- arima(monthly_sales, order = c(1, 1, 1),
              seasonal = list(order = c(1, 0, 0), period = 12))

# residuals should look like white noise
checkresiduals(fit)       # Ljung-Box test
tsdisplay(residuals(fit)) # ACF + PACF of residuals

library(xts)
library(zoo)

# create xts from a matrix and time index
dates <- as.Date(c("2024-01-01", "2024-01-03", "2024-01-10"))
prices <- xts(c(100, 102, 98), order.by = dates)
colnames(prices) <- "AAPL"

# subset by date range (very intuitive)
prices["2024-01-03"]              # specific date
prices["2024-01-01/2024-01-05"]   # date range
prices["2024-01"]                 # entire January
prices["/2024-01-05"]             # up to a date

# lag and diff (returns xts)
lag(prices, k = 1)                # previous day's price
diff(prices)                      # daily change
daily_returns <- diff(prices) / lag(prices, 1)

# rolling operations (zoo)
rollmean(prices, k = 3)           # 3-day rolling mean
rollapply(prices, 3, sd)          # 3-day rolling std dev
rollmax(prices, 3)                # 3-day rolling max

# period.apply: aggregate by period
ep <- endpoints(prices, on = "weeks")  # week endpoints
period.apply(prices, ep, mean)         # weekly averages

# merge multiple series (aligns by date)
aapl <- xts(c(100, 102, 98), as.Date(c("2024-01-01","2024-01-02","2024-01-03")))
msft <- xts(c(200, 201), as.Date(c("2024-01-01","2024-01-03")))
merged <- merge(aapl, msft)       # NA fills gaps
merged <- merge(aapl, msft, fill = na.locf)  # carry forward

# to.ts <- as.ts(prices)  # convert to base ts (loses dates)

library(forecast)
library(ggplot2)

# split into train/test
train <- window(monthly_sales, end = c(2020, 9))
test  <- window(monthly_sales, start = c(2020, 10))

# fit multiple models
fit_naive  <- naive(train, h = 3)        # naive: last value
fit_snaive <- snaive(train, h = 3)       # seasonal naive
fit_ets    <- ets(train) %>% forecast(h = 3)  # exponential smoothing
fit_arima  <- auto.arima(train) %>% forecast(h = 3)

# compare accuracy on test set
accuracy(fit_naive, test)
accuracy(fit_arima, test)
# metrics: ME, RMSE, MAE, MPE, MAPE, MASE

# time series cross-validation
ts_cv <- tsCV(train, forecastfunction = naive, h = 3)
sqrt(mean(ts_cv^2, na.rm = TRUE))   # CV RMSE

# visualize forecasts with ggplot
autoplot(monthly_sales, series = "Actual") +
    autolayer(fit_arima, series = "ARIMA", PI = FALSE) +
    autolayer(fit_ets, series = "ETS", PI = FALSE) +
    labs(title = "Forecast Comparison", x = "Time", y = "Sales") +
    theme_minimal()

# prediction intervals
fc <- forecast(auto.arima(train), h = 3, level = c(80, 95))
autoplot(fc) +
    autolayer(test, series = "Actual")

# accuracy measures explained:
#   MAE  = Mean Absolute Error (same units as data)
#   RMSE = Root Mean Squared Error (penalizes large errors)
#   MAPE = Mean Absolute Percentage Error (scale-free, %)
#   MASE = Mean Absolute Scaled Error (< 1 beats naive)

library(dplyr)
df <- tibble(
  name = c("Alice","Bob","Carol","Dave"),
  dept = c("Eng","Sales","Eng","Sales"),
  salary = c(90000, 70000, 95000, 72000),
  years = c(5, 3, 8, 4)
)

# filter rows
eng <- df %>% filter(dept == "Eng", salary > 85000)

# select columns (with helpers)
df %>% select(name, salary)
df %>% select(starts_with("s"))
df %>% select(-years)
df %>% select(name:dept)         # range
df %>% select(dept, everything())# reorder

# mutate: add/modify columns
df %>% mutate(
  bonus = salary * 0.1,
  total = salary + bonus,
  level = if_else(years >= 5, "Senior", "Junior")
)

# transmute: keep only new columns
df %>% transmute(name, annual_k = salary / 1000)

# summarize (with group_by)
df %>%
  group_by(dept) %>%
  summarize(
    avg_salary = mean(salary),
    max_years = max(years),
    n = n()
  )

employees <- tibble(id = 1:4, name = c("Al","Bo","Cy","Di"))
salaries  <- tibble(id = c(1,2,4,5), salary = c(50,60,80,90))

# mutating joins (combine columns)
inner_join(employees, salaries, by = "id")    # only matching ids
left_join(employees, salaries,  by = "id")    # all from left
right_join(employees, salaries, by = "id")    # all from right
full_join(employees, salaries,  by = "id")    # all rows

# filtering joins (filter rows, no new columns)
semi_join(employees, salaries, by = "id")     # rows in left with match
anti_join(employees, salaries, by = "id")     # rows in left WITHOUT match

# different column names
left_join(x, y, by = c("emp_id" = "id"))

# set operations (require identical columns)
union(a, b)          # unique rows in either
union_all(a, b)      # all rows (with dups)
intersect(a, b)      # rows in both
setdiff(a, b)        # rows in a but not b

# bind rows/columns
bind_rows(list(df1, df2))   # stack vertically
bind_cols(df1, df2)         # side by side (no key check!)

df <- tibble(
  dept = c("A","A","A","B","B","B"),
  name = c("x","y","z","p","q","r"),
  salary = c(50, 70, 60, 80, 90, 75)
)

# ranking within groups
df %>% group_by(dept) %>%
  mutate(
    rank = row_number(),            # 1, 2, 3 (ties get distinct)
    min_rank = min_rank(salary),    # ties get same rank, gaps
    dense_rank = dense_rank(salary),
    pct_rank = percent_rank(salary)
  )

# offsets (lead / lag)
df %>% group_by(dept) %>%
  mutate(
    prev = lag(salary),
    next = lead(salary),
    change = salary - lag(salary)
  )

# cumulative aggregates
df %>% group_by(dept) %>%
  mutate(
    cum_total = cumsum(salary),
    cum_max = cummax(salary),
    running_avg = cummean(salary)
  )

# top N per group
df %>% group_by(dept) %>%
  slice_max(salary, n = 2)        # or slice_min, slice_head, slice_sample

# row-wise operations
df %>% rowwise() %>%
  mutate(total = sum(c_across(where(is.numeric))))

df <- tibble(id = 1:3, a = c(1,2,3), b = c(4,5,6), c = c("x","y","z"))

# apply function to multiple columns
df %>% mutate(across(a:b, ~ .x * 10))
df %>% mutate(across(where(is.numeric), log))
df %>% mutate(across(everything(), as.character))

# summarize multiple columns
df %>% summarize(across(where(is.numeric), list(
  mean = ~mean(.x),
  sd   = ~sd(.x)
), .names = "{.col}_{.fn}"))

# rename multiple columns
df %>% rename_with(toupper, starts_with("a"))

# conditional transformation
df %>% mutate(across(where(is.numeric), ~ if_else(.x > 3, "high", "low")))

# use .data pronoun for programmatic access
col <- "a"
df %>% mutate(new = .data[[col]] * 2)

# pick columns by type in any verb
df %>% filter(if_all(a:b, ~ .x > 1))   # all must satisfy
df %>% filter(if_any(a:b, ~ .x > 4))   # any must satisfy

library(DBI)
library(dbplyr)

# connect to a database
con <- dbConnect(RSQLite::SQLite(), "my.db")
# or: con <- dbConnect(odbc::odbc(), "PostgreSQL")

# copy data to database
copy_to(con, mtcars, "mtcars_db", temporary = FALSE)

# reference a remote table
mtcars_remote <- tbl(con, "mtcars_db")

# dplyr verbs translate to SQL!
mtcars_remote %>%
  group_by(cyl) %>%
  summarize(avg_mpg = mean(mpg), n = n()) %>%
  arrange(desc(avg_mpg))

# see the generated SQL
mtcars_remote %>%
  filter(mpg > 20) %>%
  select(mpg, cyl, hp) %>%
  show_query()

# collect: pull results into local tibble
result <- mtcars_remote %>%
  filter(cyl == 4) %>%
  collect()

# write back to database
copy_to(con, result, "efficient_cars")

dbDisconnect(con)

library(ggplot2)
df <- data.frame(x = 1:10, y = (1:10)^2 + rnorm(10))

# basic structure: data + aesthetic + geom
ggplot(df, aes(x = x, y = y)) +
  geom_point()

# multiple layers
ggplot(df, aes(x, y)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "lm", se = TRUE)

# aesthetics map data to visual properties
ggplot(mpg, aes(displ, hwy, color = class)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE)

# facets: small multiples
ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  facet_wrap(~ class, ncol = 4)

ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  facet_grid(drv ~ cyl)   # rows ~ cols

# statistics (compute then plot)
ggplot(diamonds, aes(price)) +
  geom_histogram(bins = 50)              # stat_bin
ggplot(diamonds, aes(cut, price)) +
  stat_summary(fun = "mean", geom = "bar")

p <- ggplot(mpg, aes(displ, hwy, color = class)) +
  geom_point(size = 3)

# scales control how data maps to visuals
p + scale_x_log10() +
    scale_y_continuous(limits = c(0, 50), breaks = seq(0, 50, 10))

# manual colors
p + scale_color_manual(values = c("red","blue","green"))

# color brewer / viridis (colorblind-safe)
p + scale_color_brewer(palette = "Set1")
p + scale_color_viridis_d()

# labels
p + labs(
  title = "Fuel Efficiency",
  subtitle = "By vehicle class",
  x = "Engine Displacement (L)",
  y = "Highway MPG",
  color = "Class",
  caption = "Source: EPA"
)

# themes
p + theme_minimal()
p + theme_classic()
p + theme(
  plot.title = element_text(face = "bold", size = 16),
  panel.grid.minor = element_blank(),
  legend.position = "bottom"
)

# annotations
p + annotate("text", x = 5, y = 40, label = "Outlier", color = "red") +
    annotate("rect", xmin = 4, xmax = 6, ymin = 30, ymax = 45,
             alpha = 0.2, fill = "blue")

# boxplot by group
ggplot(mpg, aes(class, hwy)) +
  geom_boxplot() +
  coord_flip()

# violin + boxplot overlay
ggplot(mpg, aes(class, hwy, fill = class)) +
  geom_violin() +
  geom_boxplot(width = 0.1)

# density plot
ggplot(mpg, aes(hwy, fill = drv)) +
  geom_density(alpha = 0.5)

# 2D density / heatmap
ggplot(diamonds, aes(carat, price)) +
  geom_bin2d(bins = 50)            # or geom_hex()

# quantile-quantile plot
ggplot(mtcars, aes(sample = mpg)) +
  geom_qq() + geom_qq_line()

# error bars
df <- data.frame(group = c("A","B","C"),
                 mean = c(5, 7, 4), se = c(0.5, 0.8, 0.3))
ggplot(df, aes(group, mean)) +
  geom_col() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2)

# ridgeline plot (ggridges)
# library(ggridges)
# ggplot(diamonds, aes(price, cut, fill = cut)) + geom_density_ridges()

# free scales per facet
ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  facet_wrap(~ class, scales = "free")

# polar coordinates (pie chart from bar)
ggplot(mtcars, aes(x = factor(1), fill = factor(cyl))) +
  geom_bar(width = 1) +
  coord_polar(theta = "y")

# flipped coordinates
ggplot(mpg, aes(class, hwy)) +
  geom_boxplot() +
  coord_flip()

# map coordinates (sf)
# library(sf)
# ggplot(nc_sf) + geom_sf(aes(fill = AREA))

# patchwork: combine multiple plots
# library(patchwork)
# p1 <- ggplot(...)
# p2 <- ggplot(...)
# p1 + p2
# p1 / (p2 + p3)

# save plots
ggsave("plot.png", width = 8, height = 6, dpi = 300)
ggsave("plot.pdf", device = cairo_pdf)  # vector, supports fonts

# interactive (plotly)
# library(plotly)
# ggplotly(p)

# wrap plots in functions for reuse
make_scatter <- function(data, x, y, color = NULL) {
  ggplot(data, aes({{ x }}, {{ y }}, color = {{ color }})) +
    geom_point() +
    theme_minimal()
}
make_scatter(mpg, displ, hwy, class)

# using .data pronoun (string column names)
plot_col <- function(data, col) {
  ggplot(data, aes(.data[[col]])) +
    geom_bar() +
    theme_minimal()
}
plot_col(mpg, "class")

# save and load plot objects
p <- ggplot(mpg, aes(displ, hwy)) + geom_point()
saveRDS(p, "plot.rds")
p_loaded <- readRDS("plot.rds")

# modify saved plot
p_loaded + geom_smooth()

# themes as reusable objects
my_theme <- theme_minimal() +
  theme(text = element_text(family = "Helvetica"),
        plot.title = element_text(face = "bold"))
ggplot(mpg, aes(displ, hwy)) + geom_point() + my_theme

# programmatic plot generation
for (col in c("hwy","cty","cyl")) {
  p <- ggplot(mpg, aes(.data[[col]])) + geom_histogram()
  ggsave(paste0("hist_", col, ".png"), p)
}

library(tidyr)

# wide format (one row per observation, columns for each variable)
wide <- tibble(
  country = c("A","B","C"),
  `2020` = c(100, 150, 200),
  `2021` = c(110, 160, 210),
  `2022` = c(120, 170, 220)
)

# pivot_longer: wide -> long
long <- wide %>%
  pivot_longer(
    cols = -country,            # all columns except country
    names_to = "year",
    values_to = "gdp"
  )
# A tibble: 9 x 3

# pivot_wider: long -> wide
wide2 <- long %>%
  pivot_wider(names_from = year, values_from = gdp)

# multiple value columns
df <- tibble(id = 1:2, a_1 = c(1,2), a_2 = c(3,4), b_1 = c(5,6), b_2 = c(7,8))
df %>%
  pivot_longer(
    -id,
    names_to = c(".value", "time"),   # .value keeps a, b as columns
    names_sep = "_"
  )
# id  time   a     b
# 1   1      1     5
# 1   2      3     7

df <- tibble(
  name = c("John_Smith", "Jane_Doe", "Bob_Jones"),
  date_range = c("2020-01-01_2020-12-31", "2021-01-01_2021-12-31", "2022-01-01_2022-12-31")
)

# separate one column into many
df %>% separate(name, into = c("first", "last"), sep = "_")

# separate with conversion
df %>% separate(name, into = c("first","last"), sep = "_", convert = TRUE)

# extract with regex groups
df %>%
  extract(date_range,
          into = c("start", "end"),
          regex = "(.+)_(.+)")

# unite multiple columns into one
df2 <- tibble(first = c("John","Jane"), last = c("Smith","Doe"))
df2 %>% unite("full_name", first:last, sep = " ")

# separate_rows: split delimited values into rows
df3 <- tibble(id = 1:2, tags = c("a,b,c", "x,y"))
df3 %>% separate_rows(tags, sep = ",")
# id  tags
# 1   a
# 1   b
# 1   c
# 2   x
# 2   y

df <- tibble(
  x = c(1, 2, NA, 4, 5),
  y = c(NA, 2, 3, NA, 5),
  z = c("a", NA, "c", "d", NA)
)

# drop rows with any NA
df %>% drop_na()
df %>% drop_na(x)              # only column x

# replace NA with a value
df %>% replace_na(list(x = 0, y = -1, z = "unknown"))

# fill NA with previous/next value
df %>% fill(x, y, .direction = "down")   # carry forward
df %>% fill(x, y, .direction = "up")     # carry backward

# coalesce: first non-missing value
df %>% mutate(x_filled = coalesce(x, y, 0))

# na_if: replace specific value with NA
df %>% mutate(z = na_if(z, "a"))

# detect missing values
is.na(df$x)
sum(is.na(df))                 # total NAs
df %>% map_df(~ sum(is.na(.))) # NAs per column

# imputation (simple)
df %>% mutate(x = if_else(is.na(x), mean(x, na.rm = TRUE), x))

# nest: group rows into list-columns
nested <- mtcars %>%
  group_by(cyl) %>%
  nest()
# cyl  data
# 4    <tibble 11x10>
# 6    <tibble 7x10>
# 8    <tibble 14x10>

# access nested data
nested$data[[1]]               # first group's tibble

# fit models to each group
models <- nested %>%
  mutate(
    model = map(data, ~ lm(mpg ~ wt, data = .x)),
    glance = map(model, broom::glance),
    tidy = map(model, broom::tidy)
  )

# unnest results
models %>% unnest(glance)
models %>% unnest(tidy)

# unnest a list-column back to rows
df <- tibble(id = 1:2, vals = list(c(1,2,3), c(4,5)))
df %>% unnest(vals)
# id  vals
# 1   1
# 1   2
# 1   3
# 2   4
# 2   5

# chop/unchop (similar but keeps other columns as lists)
df %>% chop(vals)
df %>% unchop(vals)

library(jsonlite)
library(tidyr)

# parse JSON
json <- '[
  {"name":"Alice","age":30,"skills":["R","Python"]},
  {"name":"Bob","age":25,"skills":["SQL"]}
]'
parsed <- fromJSON(json, simplifyDataFrame = FALSE)

# convert list to tibble
tibble(person = parsed) %>%
  hoist(person,
    name = "name",
    age = "age",
    skills = "skills"
  )
# name   age  skills
# Alice  30   <chr [2]>
# Bob    25   <chr [1]>

# unnest longer for nested lists
tibble(person = parsed) %>%
  unnest_wider(person)              # spread list to columns

tibble(skills = list(c("R","Python"), c("SQL"))) %>%
  unnest_longer(skills)

# deeply nested: use purrr + tibble
flatten_df <- function(lst) {
  tibble(
    name = lst$name,
    age = lst$age,
    n_skills = length(lst$skills)
  )
}
map_dfr(parsed, flatten_df)

# rectangularize arrays
jsonlite::fromJSON(json, simplifyDataFrame = TRUE)  # auto-flatten

library(purrr)

# map: apply function to each element, return list
map(1:3, ~ .x ^ 2)             # list(1, 4, 9)

# type-specific variants (safer, faster)
map_dbl(1:3, ~ .x ^ 2)         # numeric vector: 1 4 9
map_chr(c(1,2,3), ~ paste("n=", .x))
map_int(c(1.5, 2.7), floor)
map_lgl(c(1, NA, 3), ~ !is.na(.x))

# map over two vectors in parallel
map2_dbl(c(1,2,3), c(10,20,30), ~ .x + .y)   # 11 22 33

# map over multiple vectors (pmap)
pmap_dbl(list(a = 1:3, b = 4:6, c = 7:9), sum)  # 12 15 18

# map over indices (imap)
imap_chr(c("a","b","c"), ~ paste0(.y, ":", .x))
# "1:a" "2:b" "3:c"

# walk: side effects (no return value)
walk(c("file1.csv","file2.csv"), ~ print(read.csv(.x)))

# map over columns of a data frame
map_dbl(mtcars, mean)          # mean of each column
map(mtcars, class)             # class of each column

# formula syntax (~ creates a function)
map(1:3, ~ .x + 10)              # .x is the argument
map2(1:3, 10:12, ~ .x + .y)      # .x, .y for two args
pmap(list(1:3, 4:6, 7:9), ~ ..1 + ..2 + ..3)  # ..1, ..2 for many

# full function syntax (for complex bodies)
map(1:3, function(x) {
  if (x > 2) return(x * 10)
  x
})

# named function
square <- function(x) x^2
map(1:3, square)

# extract by name/position (no function needed)
map(list(a = list(x = 1), b = list(x = 2)), "x")
map(list(list(1,2), list(3,4)), 1)   # first element of each

# pluck: safely extract deeply nested
deep <- list(a = list(b = list(c = 42)))
pluck(deep, "a", "b", "c")           # 42
pluck(deep, "a", "x", "y", .default = NA)  # NA (no error)

# chuck: same but errors if missing
# chuck(deep, "a", "x")  # error

# reduce: combine elements pairwise
reduce(c(1,2,3,4), `+`)            # 10 (sum)
reduce(c(1,2,3,4), `*`)            # 24 (product)
reduce(list(df1, df2, df3), full_join, by = "id")  # merge many

# accumulate: keep intermediate results
accumulate(c(1,2,3,4), `+`)        # 1 3 6 10 (running sum)
accumulate(c(2,3,4), `*`)          # 2 6 24 (running product)

# keep/discard: filter by predicate
keep(1:10, ~ .x %% 2 == 0)        # 2 4 6 8 10
discard(1:10, ~ .x %% 2 == 0)     # 1 3 5 7 9
keep(mtcars, is.numeric)          # keep numeric columns

# some/every/detect: test predicates
some(1:10, ~ .x > 5)              # TRUE
every(1:10, ~ is.numeric(.x))     # TRUE
detect(1:10, ~ .x > 5)            # 6 (first match)
detect_index(1:10, ~ .x > 5)      # 6 (position)

# head_while/tail_while
head_while(1:10, ~ .x < 5)        # 1 2 3 4
tail_while(10:1, ~ .x > 5)        # 10 9 8 7 6

# negate a predicate
negate(is.na)(5)                  # TRUE
keep(c(1, NA, 3), negate(is.na))  # 1 3

# safely: capture errors without stopping
safe_log <- safely(log)
safe_log(10)         # list(result = 2.3, error = NULL)
safe_log(-1)         # list(result = NULL,  error = <error>)

# process many, capturing failures
results <- map(c(10, -1, 0, "x"), safely(log))
successes <- map(results, "result") %>% discard(is.null)
failures  <- map(results, "error")  %>% discard(is.null)

# possibly: return default on error
safe_log2 <- possibly(log, otherwise = NA_real_)
map_dbl(c(10, -1, 0, "x"), safe_log2)   # 2.3 NA NA NA

# quietly: capture warnings/messages
quiet_log <- quietly(log)
quiet_log(10)   # list(result, warnings, messages)

# transpose: restructure list-of-lists
transposed <- transpose(results)
transposed$result   # all results
transposed$error    # all errors

# rate-limited / retried operations
# library(purrr)
# safely_slow <- slowly(safely(f), rate = rate_backoff())

# walk + safely for batch file processing
walk(files, ~ safely(read.csv)(.x))

# vectorized map (faster for simple operations)
# map_vec returns a vector, auto-detecting type
map_vec(1:3, ~ .x * 2)            # numeric vector
map_vec(c("a","b"), ~ toupper(.x)) # character vector

# modify: like map but preserves type
modify(mtcars, ~ .x * 2)         # still a data frame
modify_if(mtcars, is.numeric, ~ .x * 2)
modify_at(mtcars, c("mpg","cyl"), ~ .x * 2)

# list_modify / list_merge
l1 <- list(a = 1, b = 2)
list_modify(l1, b = 20, c = 30)  # a=1, b=20, c=30
list_merge(l1, list(b = 20, c = 30))

# parallel mapping with future + furrr
# library(furrr)
# plan(multisession)              # parallel backend
# future_map(1:100, slow_function, .options = furrr_options(seed = TRUE))

# progress bar
# walk(1:100, ~ Sys.sleep(0.1), .progress = TRUE)

# conditional execution in map
map(1:5, ~ if (.x > 3) .x * 10 else .x)
# 1 2 3 40 50

library(stringr)

text <- c("apple pie", "banana bread", "cherry tart")

# detect pattern
str_detect(text, "pie")           # TRUE FALSE FALSE
str_which(text, "pie")            # 1
str_count(text, "a")              # 1 3 1

# subset strings
str_subset(text, "pie")           # "apple pie"
str_extract(text, "[aeiou]")      # first vowel
str_extract_all(text, "[aeiou]")  # list of all vowels

# match groups
str_match(text, "(\\w+) (\\w+)")  # matrix with groups
str_match_all(text, "(\\w+)")

# locate pattern positions
str_locate(text, "a")             # start/end matrix
str_locate_all(text, "a")

# replace
str_replace(text, "a", "A")       # first match
str_replace_all(text, "a", "A")   # all matches
str_remove(text, "a")             # str_replace(text, "a", "")

# split
str_split("a,b,c", ",")           # list of vectors
str_split_fixed("a,b,c", ",", 3)  # matrix

library(stringr)

# case conversion
str_to_upper("hello")             # "HELLO"
str_to_lower("WORLD")             # "world"
str_to_title("the wind in the willows")  # "The Wind In The Willows"
str_to_sentence("hello world")    # "Hello world"

# trimming and padding
str_trim("  hello  ")             # "hello" (both sides)
str_trim("  hello  ", side = "left")
str_squish("  hello   world  ")   # "hello world" (collapse spaces)
str_pad("5", width = 3, pad = "0") # "005"
str_pad("5", width = 3, side = "left", pad = "0")  # "005"

# subsetting
str_sub("hello", 1, 3)            # "hel"
str_sub("hello", -3, -1)          # "llo" (negative from end)
str_sub("hello", 2)               # "ello" (to end)
str_sub("hello", 1, 1) <- "H"     # assignment: "Hello"

# length
str_length("hello")               # 5
str_length(c("a","bb","ccc"))     # 1 2 3

# combine and duplicate
str_c("a", "b", "c", sep = "-")   # "a-b-c"
str_c(c("x","y"), collapse = ",") # "x,y"
str_dup("ab", 3)                  # "ababab"

# truncate
str_trunc("Hello World", 8)       # "Hello..."

library(stringr)

# anchors
str_detect(c("cat","scat","catch"), "^cat")   # starts with: T F T
str_detect(c("cat","scat","catch"), "cat$")   # ends with:   T F F

# character classes
str_extract_all("a1 b2 c3", "[a-z]")          # letters
str_extract_all("a1 b2 c3", "[0-9]")          # digits
str_extract_all("a1 b2 c3", "[^0-9]")         # non-digits
str_extract_all("a1 b2 c3", "\\w")          # word chars
str_extract_all("a1 b2 c3", "\\s")          # whitespace

# quantifiers
str_detect("aaa", "a{2,3}")      # 2-3 a's
str_detect("color", "colou?r")   # u optional
str_detect("abbbbc", "ab+c")     # one or more
str_detect("ac", "ab*c")         # zero or more

# groups and alternation
str_match("2024-01-15", "(\\d{4})-(\\d{2})-(\\d{2})")
str_detect("cat", "cat|dog|bird")  # alternation

# lookahead/lookbehind
str_extract("abc123", "(?<=abc)\\d+")  # 123 (lookbehind)
str_extract("abc123", "\\d+(?=end)")   # nothing (no 'end')

# named capture (stringr 1.5+)
str_match("John 30", "(?<name>[A-Za-z]+) (?<age>\\d+)")

# use fixed() for literal matching
str_detect("a.b.c", fixed("."))  # TRUE (no regex)

library(lubridate)

# parse dates (auto-detect format)
ymd("2024-01-15")               # 2024-01-15
ymd("2024/01/15")
ymd("24/1/15")
mdy("01-15-2024")               # month-day-year
dmy("15/01/2024")               # day-month-year
ymd_hms("2024-01-15 14:30:00")
ymd_hm("2024-01-15 14:30")
hms("14:30:45")

# from components
make_date(2024, 1, 15)
make_datetime(2024, 1, 15, 14, 30, 0, tz = "UTC")

# current date/time
today()                          # 2024-01-15
now()                            # 2024-01-15 14:30:00 UTC

# extract components
d <- ymd("2024-01-15")
year(d)                          # 2024
month(d)                         # 1
month(d, label = TRUE)           # "Jan"
day(d)                           # 15
wday(d, label = TRUE)            # "Mon"
yday(d)                          # 15 (day of year)
quarter(d)                       # 1
semester(d)                      # 1

# set components
year(d) <- 2025
month(d) <- 6
d <- update(d, year = 2025, month = 6, day = 1)

library(lubridate)

# timezones
t <- ymd_hms("2024-01-15 14:30:00", tz = "UTC")
with_tz(t, "America/New_York")   # 09:30 EST
with_tz(t, "Asia/Shanghai")      # 22:30 CST
force_tz(t, "America/New_York")  # keeps clock, changes zone

# list timezones
OlsonNames()[1:5]

# durations (exact seconds)
d <- dseconds(60) + dminutes(2)  # 180 seconds
as.numeric(d, "seconds")
dhours(1) + ddays(1)

# periods (calendar-aware, variable length)
p <- minutes(5) + hours(2)
ymd("2024-01-15") + p            # 2024-01-15 02:05:00
ymd("2024-03-10") + days(1)      # handles DST
months(1)                        # variable length!

# intervals (start to end)
int <- interval(ymd("2024-01-01"), ymd("2024-02-01"))
int_length(int)                  # seconds
int %within% int2                # test overlap
as.period(int)                   # "1M 0S"

# arithmetic
ymd("2024-01-15") %--% ymd("2024-02-15")  # interval
ymd("2024-01-15") + weeks(2)               # period
ymd("2024-01-15") + dweeks(2)              # duration (exact)

# rounding dates
floor_date(t, "hour")            # round down
ceiling_date(t, "day")
round_date(t, "hour")

---
title: "Analysis Report"
author: "Data Team"
date: "`r format(Sys.Date(), '%B %d, %Y')`"
output: html_document
---

## Introduction

This is **Markdown** with embedded R.

Inline code: the mean is `r mean(mtcars$mpg)`.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(dplyr)
library(ggplot2)
```

```{r load-data}
df <- mtcars
head(df)
```

```{r plot, fig.width=8, fig.height=6, fig.cap="MPG by weight"}
ggplot(df, aes(wt, mpg)) + geom_point() + theme_minimal()
```

```{r table}
library(knitr)
kable(summary(df), caption = "Summary statistics")
```

---
title: "Quarterly Report"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    theme: flatly
    df_print: paged
  pdf_document:
    toc: true
    number_sections: true
  word_document: default
params:
  quarter: "Q1"
  region: "North"
---

## Report for `r params$quarter` - `r params$region`

```{r}
data <- filter(all_data, quarter == params$quarter, region == params$region)
```

# Render with parameters:
# rmarkdown::render("report.Rmd", params = list(quarter = "Q2"))

# parameterized batch rendering
quarters <- c("Q1","Q2","Q3","Q4")
for (q in quarters) {
  rmarkdown::render("report.Rmd",
    params = list(quarter = q),
    output_file = paste0("report_", q, ".html")
  )
}

library(knitr); library(kableExtra)

# basic kable
kable(head(mtcars), caption = "First 6 rows")

# styled kable
kable(head(mtcars), "html") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "lightblue") %>%
  column_spec(1, bold = TRUE)

# gt package (modern, flexible)
library(gt)
mtcars %>%
  head() %>%
  gt() %>%
  tab_header(title = "Car Specifications", subtitle = "Top 6") %>%
  fmt_number(columns = mpg, decimals = 1) %>%
  data_color(columns = mpg, colors = scales::col_numeric("Reds", NULL))

# DT: interactive tables
library(DT)
datatable(mtcars,
  filter = "top",
  options = list(pageLength = 10),
  caption = "Searchable table"
) %>%
  formatRound(columns = c("mpg","disp"), digits = 2)

# reactable for even more interactivity
# library(reactable)

---
title: "Modern Report"
format:
  html:
    toc: true
    code-fold: true
  pdf:
    documentclass: article
execute:
  echo: false
  warning: false
filters:
  - lightbox
---

## Cross-references

See Figure @fig-scatter and Table @tbl-summary.

```{r}
#| label: fig-scatter
#| fig-cap: "MPG vs Weight"
ggplot(mtcars, aes(wt, mpg)) + geom_point()
```

```{r}
#| label: tbl-summary
#| tbl-cap: "Summary by cylinder"
mtcars %>%
  group_by(cyl) %>%
  summarize(mean_mpg = mean(mpg)) %>%
  gt()
```

## Multiple languages

```{python}
import pandas as pd
df = pd.DataFrame({'a': [1,2,3]})
print(df)
```

```{sql}
SELECT * FROM mtcars WHERE mpg > 20
```

# Quarto supports Python, Julia, Observable, and more
# in the same document, sharing data via ojs.

# render a report programmatically
rmarkdown::render("daily_report.Rmd",
  output_dir = "reports",
  output_file = paste0("report_", Sys.Date(), ".html"),
  params = list(date = Sys.Date() - 1),
  quiet = TRUE
)

# email the report
# library(blastula)
# render_email("email_template.Rmd") %>%
#   smtp_send(
#     to = "[email protected]",
#     from = "[email protected]",
#     subject = paste("Daily Report -", Sys.Date())
#   )

# schedule with cron (Linux/Mac) or Task Scheduler (Windows)
# crontab -e:
# 0 8 * * * cd /path/to/project && Rscript -e 'rmarkdown::render("daily.Rmd")'

# or use the cronR package
# library(cronR)
# cmd <- cron_rscript("render_report.R")
# cron_add(cmd, frequency = "daily", at = "08:00")

# batch render multiple reports
purrr::walk(regions, function(r) {
  rmarkdown::render("regional.Rmd",
    params = list(region = r),
    output_file = paste0("report_", r, ".html")
  )
})

# version control: cache expensive computations
# knitr::opts_chunk$set(cache = TRUE)

library(shiny)

ui <- fluidPage(
  titlePanel("My First Shiny App"),
  sidebarLayout(
    sidebarPanel(
      sliderInput("n", "Number of points:", min = 1, max = 100, value = 50),
      selectInput("color", "Color:", choices = c("red","blue","green"))
    ),
    mainPanel(
      plotOutput("scatter"),
      verbatimTextOutput("summary")
    )
  )
)

server <- function(input, output, session) {
  data <- reactive({
    data.frame(x = rnorm(input$n), y = rnorm(input$n))
  })

  output$scatter <- renderPlot({
    plot(data()$x, data()$y, col = input$color, pch = 19,
         xlab = "X", ylab = "Y", main = paste("n =", input$n))
  })

  output$summary <- renderPrint({
    summary(data())
  })
}

shinyApp(ui, server)

# save as app.R in a folder, then runApp("folder")
# or run with shiny::runApp()

server <- function(input, output, session) {
  # reactive expression: lazy, cached
  filtered <- reactive({
    mtcars %>% filter(cyl == input$cyl)
  })

  # observe: eager, for side effects
  observe({
    updateSelectInput(session, "model",
      choices = unique(filtered()$model))
  })

  # observeEvent: triggered by specific input
  observeEvent(input$reset, {
    updateSliderInput(session, "cyl", value = 4)
  })

  # eventReactive: lazy, triggered by event
  result <- eventReactive(input$go, {
    run_analysis(input$param)
  })

  # reactiveVal: mutable value
  counter <- reactiveVal(0)
  observeEvent(input$add, counter(counter() + 1))

  # reactiveValues: multiple mutable values
  rv <- reactiveValues(data = NULL, status = "ready")
  observeEvent(input$load, {
    rv$data <- read.csv(input$file$datapath)
    rv$status <- "loaded"
  })

  # isolate: read without dependency
  output$plot <- renderPlot({
    plot(isolate(rv$data))
  })
}

# dynamic UI
ui <- fluidPage(
  uiOutput("dynamic_controls"),
  plotOutput("plot")
)

server <- function(input, output, session) {
  output$dynamic_controls <- renderUI({
    cols <- names(input$data)
    selectInput("xvar", "X variable:", choices = cols)
  })

  # insert/remove UI
  observeEvent(input$add, {
    insertUI("#placeholder", "beforeEnd",
      sliderInput(paste0("s", input$add), "Slider", 0, 100, 50))
  })

  # modules: reusable UI+server
  histogramUI <- function(id) {
    ns <- NS(id)
    tagList(
      selectInput(ns("var"), "Variable:", names(mtcars)),
      plotOutput(ns("hist"))
    )
  }
  histogramServer <- function(id) {
    moduleServer(id, function(input, output, session) {
      output$hist <- renderPlot(hist(mtcars[[input$var]]))
    })
  }

  # use module
  ui <- fluidPage(histogramUI("hist1"), histogramUI("hist2"))
  server <- function(input, output, session) {
    histogramServer("hist1")
    histogramServer("hist2")
  }
}

# file upload
fileInput("data", "Upload CSV:", accept = ".csv")
# in server: input$data$datapath (temp file path)

# date range
dateRangeInput("dates", "Period:", start = Sys.Date() - 30, end = Sys.Date())

# tabset with conditional panels
tabsetPanel(
  tabPanel("Plot", plotOutput("p")),
  tabPanel("Table", DT::dataTableOutput("t")),
  tabPanel("Summary", verbatimTextOutput("s"))
)

# render DataTable (interactive)
output$t <- DT::renderDataTable({
  datatable(filtered(), filter = "top", options = list(pageLength = 25))
})

# render UI from ggplot
output$p <- renderPlot({
  ggplot(filtered(), aes(.data[[input$x]], .data[[input$y]])) +
    geom_point() + theme_minimal()
})

# download handlers
output$download <- downloadHandler(
  filename = function() paste0("data_", Sys.Date(), ".csv"),
  content = function(file) write.csv(filtered(), file)
)

# render image (pre-generated)
output$img <- renderImage({
  list(src = "plot.png", contentType = "image/png", width = 400)
}, deleteFile = FALSE)

# performance: debounce/throttle rapid inputs
input_debounced <- debounce(reactive({ input$text }), 500)

# async with future
library(future)
plan(multisession)
result <- reactive({
  future_promise({ expensive_computation(input$params) })
})

# caching expensive computations
# in UI: add resourcePath or use bindCache
output$plot <- renderPlot({ ... }) %>% bindCache(input$dataset)

# deploy to shinyapps.io
# library(rsconnect)
# deployApp("myapp/")

# run multiple Shiny apps behind a load balancer
# (shinyapps.io / RStudio Connect handle this automatically)

# dockerize for self-hosting
# Dockerfile:
# FROM rocker/shiny
# COPY app.R /srv/shiny-server/
# EXPOSE 3838

# shiny server config (/etc/shiny-server/shiny-server.conf)
# run_as shiny;
# server { listen 3838; location / { site_dir /srv/shiny-server; } }

# monitor with shinylogs or shiny.telemetry

# one-sample t-test: is mean different from hypothesized value?
t.test(mtcars$mpg, mu = 20)
#  t = 0.085, df = 31, p-value = 0.933
#  95% CI: [17.92, 22.26]
#  mean of x: 20.09

# two-sample t-test
t.test(mpg ~ cyl, data = subset(mtcars, cyl %in% c(4, 6)))

# paired t-test
t.test(pre, post, paired = TRUE)

# Wilcoxon (non-parametric alternative)
wilcox.test(mpg ~ cyl, data = subset(mtcars, cyl %in% c(4,6)))

# interpret results:
#   p < 0.05: reject null hypothesis (significant)
#   p > 0.05: fail to reject (not enough evidence)
#   CI shows plausible range of effect
#   check assumptions: normality (shapiro.test), equal variance

# power analysis
library(pwr)
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.8)
# n = 64 per group needed to detect medium effect

# one-way ANOVA: compare means across 3+ groups
fit <- aov(mpg ~ factor(cyl), data = mtcars)
summary(fit)
#              Df Sum Sq Mean Sq F value  Pr(>F)
# factor(cyl)   2  824.8   412.4    39.7 4.98e-09

# check assumptions
plot(fit)                       # residuals vs fitted, QQ
shapiro.test(resid(fit))        # normality of residuals
library(car); leveneTest(fit)  # equal variance

# post-hoc tests (which groups differ?)
TukeyHSD(fit)                   # Tukey HSD
pairwise.t.test(mtcars$mpg, mtcars$cyl, p.adjust = "bonferroni")

# two-way ANOVA with interaction
fit2 <- aov(mpg ~ cyl * am, data = mtcars)
summary(fit2)

# Kruskal-Wallis (non-parametric ANOVA)
kruskal.test(mpg ~ factor(cyl), data = mtcars)

# repeated measures ANOVA
# library(afex)
# aov_ez(id = "subject", dv = "score", data = df, within = "condition")

# mixed-effects models (lme4)
library(lme4)
lmer(score ~ treatment + (1|subject), data = df)

# correlation
cor(mtcars$mpg, mtcars$wt)                    # -0.867
cor.test(mtcars$mpg, mtcars$wt)               # with p-value
cor(mtcars[, c("mpg","wt","hp","disp")])      # correlation matrix
library(corrplot); corrplot(cor(mtcars[,1:4]))# visualize

# linear regression
fit <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit)
confint(fit)                                   # CIs for coefficients

# diagnostics
par(mfrow = c(2,2)); plot(fit)
# 1. Residuals vs Fitted: check linearity
# 2. Normal Q-Q: check normality
# 3. Scale-Location: check homoscedasticity
# 4. Residuals vs Leverage: check influential points

# check for multicollinearity
library(car); vif(fit)                         # variance inflation factor

# influential observations
cooks.distance(fit) |> plot()
influence.measures(fit)

# predictions with intervals
predict(fit, newdata, interval = "confidence") # CI for mean
predict(fit, newdata, interval = "prediction") # PI for new obs

# compare nested models
anova(fit1, fit2)                              # F-test
AIC(fit1, fit2); BIC(fit1, fit2)              # information criteria

# contingency table
tbl <- table(mtcars$cyl, mtcars$am)
#    0  1
# 4  3  8
# 6  4  3
# 8 12  2

# chi-square test of independence
chisq.test(tbl)
# X-squared = 8.74, df = 2, p-value = 0.0126

# expected counts (should be > 5 for valid chi-square)
chisq.test(tbl)$expected

# Fisher's exact test (small samples)
fisher.test(tbl)

# goodness of fit (one variable vs expected proportions)
chisq.test(c(10, 20, 30), p = c(1/3, 1/3, 1/3))

# McNemar's test (paired nominal data)
mcnemar.test(table(pre, post))

# Cochran-Mantel-Haenszel (stratified)
mantelhaen.test(tbl3d)

# visualize
library(ggplot2); library(ggmosaic)
ggplot(as.data.frame(tbl)) +
  geom_mosaic(aes(weight = Freq, x = product(Var1), fill = Var2))

# effect size
library(vcd); assocstats(tbl)  # Cramer's V, phi

library(brms)

# Bayesian linear regression
fit <- brm(
  mpg ~ wt + hp,
  data = mtcars,
  family = gaussian(),
  prior = c(
    prior(normal(0, 10), class = "b"),       # weakly informative
    prior(cauchy(0, 2), class = "sigma")
  ),
  chains = 4, cores = 4, iter = 2000
)

# summary
summary(fit)
plot(fit)                       # posterior distributions
pp_check(fit)                   # posterior predictive check

# extract posterior samples
posterior <- as_draws_df(fit)
mean(posterior$b_wt > 0)        # P(coefficient > 0)

# predictions
posterior_predict(fit, newdata) |> as.data.frame() -> preds

# model comparison
fit_null <- brm(mpg ~ 1, data = mtcars, ...)
loo(fit, fit_null)              # leave-one-out CV
waic(fit, fit_null)

# hierarchical model
fit_h <- brm(
  score ~ treatment + (1 + treatment|subject),
  data = df, family = gaussian()
)

# Stan code behind the model
stancode(fit)

library(caret)

# split data
set.seed(42)
idx <- createDataPartition(iris$Species, p = 0.8, list = FALSE)
train <- iris[idx, ]; test <- iris[-idx, ]

# stratified split for classification
# createDataPartition preserves class proportions

# k-fold cross-validation
ctrl <- trainControl(
  method = "cv", number = 10,
  savePredictions = "final",
  classProbs = TRUE,           # for AUC
  summaryFunction = multiClassSummary
)

# preprocessing pipeline
preproc <- preProcess(train[, -5],
  method = c("center","scale","YeoJohnson","nzv"))
train_processed <- predict(preproc, train[, -5])
test_processed  <- predict(preproc, test[, -5])

# common preprocessing methods:
#   center, scale: standardize
#   BoxCox, YeoJohnson: transform to normality
#   pca: principal components
#   nzv: remove near-zero variance
#   knnImpute, bagImpute: impute missing values

# dummy variables for categorical
dummies <- dummyVars(Species ~ ., data = train)
predict(dummies, train)

library(caret)

# train a random forest
ctrl <- trainControl(method = "cv", number = 5)
rf_fit <- train(
  Species ~ ., data = train,
  method = "rf",
  trControl = ctrl,
  tuneGrid = expand.grid(mtry = 1:4),
  tuneLength = 5              # auto-tune grid
)
print(rf_fit)                  # shows accuracy by tuning param
plot(rf_fit)                   # tuning curve

# try multiple models
models <- c("rf","gbm","svmRadial","knn","rpart")
fits <- lapply(models, function(m) {
  train(Species ~ ., data = train, method = m, trControl = ctrl)
})
names(fits) <- models

# compare models
resamps <- resamples(fits)
summary(resamps)
dotplot(resamps, metric = "Accuracy")

# custom tuning grid
grid <- expand.grid(
  mtry = c(2, 4, 8),
  splitrule = c("gini","extratrees"),
  min.node.size = c(1, 5, 10)
)
fit <- train(Species ~ ., data = train, method = "ranger",
             trControl = ctrl, tuneGrid = grid)

# predictions
pred <- predict(rf_fit, test)
prob <- predict(rf_fit, test, type = "prob")

# confusion matrix
cm <- confusionMatrix(pred, test$Species)
print(cm)
#               Reference
# Prediction   setosa versicolor virginica
# setosa          10       0       0
# versicolor       0      10       1
# virginica        0       0       9
# Overall Accuracy: 0.9667
# byClass: Sensitivity, Specificity, etc.

# two-class metrics
cm$byClass[, c("Sensitivity","Specificity","Precision","Recall","F1")]

# ROC and AUC
library(pROC)
roc_curve <- roc(test$Species == "versicolor", prob$versicolor)
auc(roc_curve)
plot(roc_curve)

# multi-class AUC
library(pRoc)
multiclass.roc(test$Species, prob)

# precision-recall curve (better for imbalanced data)
library(PRROC)
pr <- pr.curve(scores.class0 = prob$versicolor,
               weights.class0 = test$Species == "versicolor",
               curve = TRUE)
plot(pr)

# custom metrics in trainControl
twoClassSummary  # built-in for 2-class
mnLogLoss         # multi-class log loss

library(caret)

# recursive feature elimination (RFE)
ctrl <- rfeControl(functions = rfFuncs, method = "cv", number = 10)
rfe_fit <- rfe(train[,-5], train$Species,
  sizes = c(1:4), rfeControl = ctrl)
print(rfe_fit)              # optimal subset
predictors(rfe_fit)         # selected features
plot(rfe_fit)

# variable importance from trained model
varImp(rf_fit)              # caret extracts importance
plot(varImp(rf_fit))

# filter methods (before training)
# remove highly correlated features
cor_matrix <- cor(train[,-5])
high_cor <- findCorrelation(cor_matrix, cutoff = 0.9)
train_filtered <- train[,-high_cor]

# remove linear dependencies
findLinearCombos(train[,-5])

# genetic algorithm / simulated annealing selection
# gafsFit <- gafs(x, y, iters = 20)
# safsFit <- safs(x, y, iters = 20)

# SHAP values (model-agnostic interpretation)
# library(fastshap)
# explain(rf_fit, X = test, nsim = 100)

library(caretEnsemble)

# train multiple models with same resampling
ctrl <- trainControl(method = "cv", number = 5,
                     savePredictions = "final",
                     classProbs = TRUE)

models <- caretList(
  Species ~ ., data = train,
  trControl = ctrl,
  methodList = c("rf","gbm","svmRadial","knn")
)

# simple ensemble (average predictions)
ens <- caretEnsemble(models)
summary(ens)
plot(ens)

# stack: train a meta-model on base predictions
stack <- caretStack(models, method = "glm",
  metric = "Accuracy", trControl = ctrl)
print(stack)

# predict with ensemble
pred_ens <- predict(ens, test)
pred_stack <- predict(stack, test)

# weighted ensemble (custom weights)
weights <- c(0.4, 0.3, 0.2, 0.1)
probs <- lapply(models, function(m) predict(m, test, type = "prob"))
final_prob <- Reduce('+', Map(function(p, w) p * w, probs, weights))

# bagging (bootstrap aggregating)
bag_fit <- train(Species ~ ., data = train, method = "treebag",
                 trControl = ctrl)

# apply: over array margins (rows or columns)
m <- matrix(1:12, 3, 4)
apply(m, 1, sum)          # row sums: 22 26 30
apply(m, 2, mean)         # column means
apply(m, c(1,2), sqrt)    # element-wise (silly but valid)

# built-in row/col functions are faster
rowSums(m); rowMeans(m)
colSums(m); colMeans(m)

# lapply: list in, list out
lapply(list(a=1:3, b=4:6), mean)   # list(a=2, b=5)
lapply(mtcars, class)               # class of each column

# sapply: list in, vector out (simplifies)
sapply(mtcars, mean)                # named numeric vector
sapply(mtcars, is.numeric)          # logical vector

# vapply: type-safe sapply (specify output template)
vapply(mtcars, mean, numeric(1))    # always numeric(1)
vapply(mtcars, class, character(1)) # always character(1)

# replicate: repeat expression n times
replicate(5, rnorm(3))              # 3x5 matrix
replicate(100, mean(rexp(30)))      # 100 sample means

# mapply: multivariate map
mapply(rep, 1:3, 3:1)               # list(1,1,1, 2,2, 3)
mapply(function(a,b) a+b, 1:3, 4:6) # 5 7 9

# BAD: element-wise loop
x <- 1:1e6
y <- numeric(length(x))
for (i in seq_along(x)) y[i] <- x[i]^2

# GOOD: vectorized
y <- x^2   # ~100x faster

# ifelse vs vectorized
# BAD
y <- numeric(length(x))
for (i in seq_along(x)) {
  y[i] <- if (x[i] > 5) x[i] else 0
}
# GOOD
y <- ifelse(x > 5, x, 0)
# BEST (fastest)
y <- x * (x > 5)

# row-wise operations (avoid if possible)
df <- data.frame(a = 1:1000, b = 1001:2000)
# BAD
df$sum <- apply(df, 1, sum)
# GOOD
df$sum <- df$a + df$b

# preallocate!
n <- 1000
result <- numeric(n)       # not result <- c()
for (i in 1:n) result[i] <- i^2

# use Rcpp for hot loops
# library(Rcpp)
# cppFunction('double sumc(NumericVector x) { return sum(x); }')

# benchmark
library(microbenchmark)
microbenchmark(
  loop = {y <- numeric(length(x)); for(i in seq_along(x)) y[i] <- x[i]^2},
  vectorized = x^2,
  times = 10
)

# data.table: faster, memory-efficient alternative to data.frame
library(data.table)
dt <- data.table(mtcars)

# syntax: dt[i, j, by]
dt[mpg > 20, .(mean_hp = mean(hp)), by = .(cyl, gear)]
#   cyl gear mean_hp
# 1:   4    4    76.0
# 2:   4    5   102.0

# reference semantics (modify in place, no copy)
dt[, mpg_dbl := mpg * 2]              # add column in place
dt[1:5, mpg := 0]                     # modify subset in place
dt[, c("a","b") := .(mpg*2, hp/10)]   # multiple columns

# setkey for fast lookups and joins
setkey(dt, cyl)
dt[.(4)]                               # all rows where cyl == 4 (binary search)
dt[.(4), mean(mpg)]                    # fast aggregation

# joins
dt1 <- data.table(id = 1:3, x = letters[1:3])
dt2 <- data.table(id = 2:4, y = LETTERS[2:4])
setkey(dt1, id); setkey(dt2, id)
dt1[dt2]                               # right join
dt2[dt1]                               # left join
merge(dt1, dt2, by = "id")             # inner join

# fread/fwrite: fast I/O
big <- fread("huge.csv")               # ~10x faster than read.csv
fwrite(big, "out.csv")

# parallel package (built-in)
library(parallel)
ncores <- detectCores() - 1

# parallel lapply
cl <- makeCluster(ncores)
results <- parLapply(cl, 1:100, function(i) slow_function(i))
stopCluster(cl)

# foreach with doParallel
library(foreach); library(doParallel)
registerDoParallel(ncores)
results <- foreach(i = 1:100, .combine = "c") %dopar% {
  slow_function(i)
}
stopImplicitCluster()

# future ecosystem (modern, flexible)
library(future); library(furrr)
plan(multisession, workers = ncores)   # or multicore (Linux/Mac)
results <- future_map(1:100, slow_function)

# parallel caret
library(caret)
ctrl <- trainControl(method = "cv", number = 10, allowParallel = TRUE)
fit <- train(y ~ ., data = train, method = "rf", trControl = ctrl)

# benchmark
library(microbenchmark)
microbenchmark(
  serial = lapply(1:100, slow_function),
  parallel = parLapply(makeCluster(4), 1:100, slow_function),
  times = 5
)

# export variables to workers
clusterExport(cl, c("my_data", "my_function"))
clusterEvalQ(cl, library(dplyr))

# profile to find bottlenecks
library(profvis)
profvis({
  result <- my_analysis(big_data)
})
# opens interactive flame graph

# Rprof (base R)
Rprof("profile.out")
my_analysis(big_data)
Rprof(NULL)
summaryRprof("profile.out")

# memory profiling
Rprofmem("mem.out")
my_analysis(big_data)
Rprofmem(NULL)

# object sizes
object.size(my_data)              # bytes
format(object.size(my_data), "MB")

# find memory hogs
sort(sapply(ls(), function(x) object.size(get(x))))

# garbage collection
gc()                              # force collection, show memory
gcinfo(TRUE)                      # report on auto GC

# common optimizations:
# 1. Vectorize (avoid loops)
# 2. Preallocate
# 3. Use data.table instead of data.frame
# 4. Use Rcpp for hot loops
# 5. Avoid growing objects (c(), rbind())
# 6. Use appropriate data types (integer vs double)
# 7. Cache expensive computations (memoise)
library(memoise)
slow_cached <- memoise(slow_function)

library(dplyr)

starwars %>%
  filter(species == "Human", height > 170) %>%
  select(name, height, mass, homeworld) %>%
  mutate(bmi = mass / (height/100)^2) %>%
  arrange(desc(height)) %>%
  slice_head(n = 5)

# group_by + summarize
starwars %>%
  group_by(species) %>%
  summarize(
    n = n(),
    avg_height = mean(height, na.rm = TRUE),
    avg_mass = mean(mass, na.rm = TRUE)
  ) %>%
  arrange(desc(n)) %>%
  filter(n >= 2)

library(dplyr)

# mutating joins (add columns from y to x)
left_join(x, y, by = "id")        # all rows in x
right_join(x, y, by = "id")       # all rows in y
inner_join(x, y, by = "id")       # only matching rows
full_join(x, y, by = "id")        # all rows from both

# different column names
left_join(x, y, by = c("id" = "user_id"))

# multiple keys
left_join(x, y, by = c("first", "last"))

# filtering joins (filter x, no columns added)
semi_join(x, y, by = "id")        # rows in x that match y
anti_join(x, y, by = "id")        # rows in x that DON'T match y

# check for duplicates
x %>% count(id) %>% filter(n > 1)

library(dplyr)

# row numbering and ranking
df %>% group_by(group) %>%
  mutate(
    row_num = row_number(),
    rank = rank(value),
    dense_rank = dense_rank(value),
    min_rank = min_rank(value)
  )

# offsets
df %>% group_by(group) %>%
  mutate(
    prev = lag(value),
    next = lead(value, 2),
    diff = value - lag(value)
  )

# cumulative
df %>% group_by(group) %>%
  mutate(
    cum_sum = cumsum(value),
    cum_mean = cummean(value),
    cum_max = cummax(value),
    cum_min = cummin(value)
  )

# top n per group
df %>% group_by(group) %>%
  slice_max(value, n = 3)         # top 3
df %>% group_by(group) %>%
  slice_min(value, n = 2)

library(dplyr)

# apply function to multiple columns
starwars %>%
  mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

# summarize multiple columns
starwars %>%
  summarize(across(where(is.numeric), mean, na.rm = TRUE))

# rename multiple columns
starwars %>%
  rename_with(tolower, everything())

# transform specific columns
df %>%
  mutate(across(c(height, mass), ~ .x * 0.453592))  # lb to kg

# multiple functions
df %>%
  summarize(across(value, list(mean = mean, sd = sd, n = ~ n()), na.rm = TRUE))

# conditional transformation
df %>%
  mutate(across(where(is.character), tolower))

library(dplyr)

# basic
df %>%
  group_by(category) %>%
  summarize(
    count = n(),
    distinct = n_distinct(id),
    total = sum(value, na.rm = TRUE),
    avg = mean(value, na.rm = TRUE),
    median = median(value, na.rm = TRUE),
    sd = sd(value, na.rm = TRUE),
    min = min(value, na.rm = TRUE),
    max = max(value, na.rm = TRUE),
    first = first(value),
    last = last(value),
    q25 = quantile(value, 0.25, na.rm = TRUE)
  )

# multiple summary stats at once
df %>%
  group_by(category) %>%
  summarize(
    across(value, list(mean = mean, sd = sd, n = ~ n()), na.rm = TRUE)
  )

# count and proportions
df %>%
  count(category) %>%
  mutate(pct = n / sum(n))

library(ggplot2)

ggplot(mpg, aes(displ, hwy, color = class)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE) +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Fuel efficiency by engine size",
    subtitle = "Source: EPA mpg dataset",
    x = "Engine displacement (L)",
    y = "Highway MPG",
    color = "Vehicle class"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

# facets
ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  facet_wrap(~ class, ncol = 4)

ggplot(mpg, aes(displ, hwy)) +
  geom_point() +
  facet_grid(drv ~ cyl)

ggplot(mpg, aes(displ, hwy, color = cty)) +
  geom_point() +
  scale_x_log10() +
  scale_y_continuous(limits = c(10, 50), breaks = seq(10, 50, 5)) +
  scale_color_viridis_c(option = "plasma") +
  coord_cartesian(xlim = c(1, 7))   # zoom without removing data

# coord flip
ggplot(mpg, aes(class, hwy)) +
  geom_boxplot() +
  coord_flip()

# coord polar (pie chart from bar)
ggplot(mtcars, aes(factor(1), fill = factor(cyl))) +
  geom_bar(width = 1) +
  coord_polar(theta = "y")

# date axis
ggplot(economics, aes(date, unemploy)) +
  geom_line() +
  scale_x_date(date_labels = "%Y", date_breaks = "5 years")

library(ggplot2)

p <- ggplot(mpg, aes(class, hwy)) + geom_boxplot()

# built-in themes
p + theme_minimal()
p + theme_classic()
p + theme_bw()

# custom theme
p + theme(
  panel.background = element_rect(fill = "white"),
  panel.grid.major = element_line(color = "gray90"),
  panel.grid.minor = element_blank(),
  axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
  axis.title = element_text(face = "bold"),
  plot.title = element_text(size = 16, hjust = 0.5),
  plot.subtitle = element_text(color = "gray40"),
  legend.position = "right",
  legend.key = element_blank(),
  strip.background = element_rect(fill = "lightblue"),
  strip.text = element_text(face = "bold")
)

# save
ggsave("plot.png", width = 8, height = 6, dpi = 300)
ggsave("plot.pdf", device = cairo_pdf)

library(ggplot2)

# stat_summary for custom summaries
ggplot(mtcars, aes(factor(cyl), mpg)) +
  stat_summary(fun = "mean", geom = "point", size = 3) +
  stat_summary(fun.data = "mean_se", geom = "errorbar")

# smooth
ggplot(mtcars, aes(wt, mpg)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  geom_smooth(method = "loess", se = FALSE, color = "red")

# density
ggplot(iris, aes(Sepal.Length, fill = Species)) +
  geom_density(alpha = 0.5)

# 2d density
ggplot(diamonds, aes(carat, price)) +
  geom_density_2d() +
  scale_x_log10() + scale_y_log10()

# histograms with binning
ggplot(diamonds, aes(price)) +
  geom_histogram(bins = 50, fill = "steelblue") +
  scale_x_log10()

# boxplot and violin
ggplot(iris, aes(Species, Sepal.Length)) +
  geom_boxplot() +
  geom_violin(alpha = 0.3, fill = "orange")

library(ggplot2)
library(patchwork)

p1 <- ggplot(mtcars, aes(wt, mpg)) + geom_point()
p2 <- ggplot(mtcars, aes(factor(cyl))) + geom_bar()
p3 <- ggplot(mtcars, aes(mpg)) + geom_histogram(bins = 10)
p4 <- ggplot(mtcars, aes(factor(cyl), mpg)) + geom_boxplot()

# combine plots
p1 + p2                          # side by side
p1 / p2                          # stacked
(p1 + p2) / p3                   # grouped
p1 + p2 + p3 + plot_layout(nrow = 2, byrow = FALSE)

# annotations
p1 + p2 + plot_annotation(
  title = "MT cars analysis",
  tag_levels = "A"
)

# other extensions:
#   ggrepel: non-overlapping labels
#   ggridges: joy plots
#   gganimate: animated plots
#   ggiraph: interactive
#   ggplot2::geom_sf: maps

library(ggrepel)
ggplot(mtcars, aes(wt, mpg, label = rownames(mtcars))) +
  geom_point() +
  geom_text_repel()

library(tidyr)

# wide to long
#   id  q1  q2  q3
#   A   10  20  30
wide_data %>%
  pivot_longer(
    cols = starts_with("q"),
    names_to = "quarter",
    values_to = "score"
  )
#   id  quarter  score
#   A   q1       10
#   A   q2       20
#   A   q3       30

# long to wide
long_data %>%
  pivot_wider(
    names_from = quarter,
    values_from = score
  )

# multiple value columns
df %>%
  pivot_longer(
    cols = c(starts_with("q"), starts_with("r")),
    names_to = c(".value", "quarter"),
    names_pattern = "([a-z])([0-9])"
  )

library(tidyr)

# split a column
df <- tibble(x = c("a_1", "b_2", "c_3"))
df %>% separate(x, c("letter", "number"), sep = "_")
#   letter number
#   a      1
#   b      2

# split into rows
df %>% separate_rows(x, sep = "_")
#   x
#   a
#   1
#   b
#   2

# unite columns
df %>%
  separate(x, c("letter", "number")) %>%
  unite("combined", letter, number, sep = "-")
#   combined
#   a-1
#   b-2

# extract with regex
df %>% extract(x, c("letter", "number"), "([a-z])_([0-9])")

library(tidyr)

# replace NA
df %>% replace_na(list(value = 0, name = "unknown"))

# fill NA with previous/next value
df %>% fill(value, .direction = "down")    # forward fill
df %>% fill(value, .direction = "up")      # backward fill
df %>% fill(value, .direction = "downup")  # down then up

# drop rows with NA
df %>% drop_na()
df %>% drop_na(value)              # only specific columns

# detect NA
df %>% filter(!is.na(value))
sum(is.na(df$value))
df %>% mutate_all(~ sum(is.na(.))) # NA count per column

# coalesce: first non-NA
df %>% mutate(value = coalesce(value, fallback, 0))

# complete: expand combinations
df %>% complete(group, year)       # all group-year combos, NA filled
df %>% complete(group, year, fill = list(value = 0))

library(tidyr)
library(dplyr)

# nest: pack rows into list columns
nested <- mtcars %>%
  nest(data = -cyl)
#   cyl  data
#   4    <tibble 11x10>
#   6    <tibble 7x10>
#   8    <tibble 14x10>

# operate on each group
nested %>%
  mutate(model = map(data, ~ lm(mpg ~ wt, data = .x)))

# unnest
df %>%
  unnest(data)

# unnest specific column
df %>%
  unnest_wider(data)              # list-col of named lists -> cols
df %>%
  unnest_longer(data)             # list-col of vectors -> rows

# chop (similar to nest)
df %>% chop(value)

# pack (opposite of unnest)
df %>% pack(x = c(a, b))

library(tibble)
library(readr)

# create tibble
tibble(
  x = 1:5,
  y = c("a", "b", "c", "d", "e"),
  z = runif(5)
)

# tribble: row-wise construction
tribble(
  ~name, ~age, ~score,
  "Alice", 30, 90,
  "Bob",   25, 80
)

# read csv
df <- read_csv("data.csv", col_names = TRUE, col_types = "cdd")
df <- read_csv("data.csv", na = c("", "NA", "N/A"))
df <- read_tsv("data.tsv")
df <- read_delim("data.txt", delim = "|")

# write
write_csv(df, "out.csv")
write_tsv(df, "out.tsv")

# read from clipboard (Excel)
df <- read_clipboard()

# column types
#   c = character, d = double, i = integer, l = logical, f = factor, D = date

# one-sample t-test
t.test(x, mu = 5)
# H0: mean of x equals 5

# two-sample t-test
t.test(x, y)
t.test(x, y, var.equal = TRUE)    # Student's (assume equal var)
t.test(x, y, alternative = "greater")

# paired t-test
t.test(before, after, paired = TRUE)

# output interpretation
result <- t.test(x, y)
result$p.value                    # < 0.05 -> reject H0
result$conf.int                   # 95% CI of difference
result$statistic                  # t value

# check assumptions
shapiro.test(x)                   # normality
var.test(x, y)                    # equal variances

# effect size (effsize package)
library(effsize)
cohen.d(x, y)

# one-way ANOVA
result <- aov(yield ~ fertilizer, data = df)
summary(result)

# two-way ANOVA with interaction
result <- aov(yield ~ fertilizer + density + fertilizer:density, data = df)
result <- aov(yield ~ fertilizer * density, data = df)  # shorthand

# Tukey post-hoc
TukeyHSD(result)

# check assumptions
plot(result)                      # diagnostic plots
library(car)
leveneTest(yield ~ fertilizer, data = df)  # equal variances

# Kruskal-Wallis (non-parametric alternative)
kruskal.test(yield ~ fertilizer, data = df)

# repeated measures
result <- aov(score ~ time + Error(subject/time), data = df)

# mixed effects (lme4)
library(lme4)
model <- lmer(score ~ time + (1 | subject), data = df)
anova(model)

# chi-square test of independence
tbl <- table(df$gender, df$vote)
chisq.test(tbl)

# goodness of fit
chisq.test(table(df$color), p = c(0.25, 0.25, 0.25, 0.25))

# Fisher's exact (small samples)
fisher.test(tbl)

# McNemar (paired binary)
mcnemar.test(tbl)

# expected frequencies
ct <- chisq.test(tbl)
ct$expected
ct$observed
ct$residuals                     # Pearson residuals

# association measures
library(vcd)
assocstats(tbl)

# visualize
library(ggplot2)
ggplot(df, aes(gender, fill = vote)) +
  geom_bar(position = "fill")

# Pearson (linear, normal)
cor(x, y, method = "pearson")
cor.test(x, y, method = "pearson")

# Spearman (rank, non-parametric)
cor(x, y, method = "spearman")

# Kendall (rank, smaller samples)
cor(x, y, method = "kendall")

# correlation matrix
cor(mtcars[, c("mpg", "wt", "hp", "disp")])

# with p-values
library(psych)
corr.test(mtcars[, 1:6])

# visualize
library(corrplot)
M <- cor(mtcars)
corrplot(M, method = "circle", type = "lower", order = "hclust")

# GGally for scatterplot matrix
library(GGally)
ggpairs(mtcars[, c("mpg", "wt", "hp", "cyl")])

# handle missing
cor(df, use = "complete.obs")     # listwise deletion
cor(df, use = "pairwise.complete.obs")

# Mann-Whitney U (alternative to two-sample t)
wilcox.test(x, y)
wilcox.test(x, y, paired = FALSE)

# Wilcoxon signed-rank (alternative to paired t)
wilcox.test(before, after, paired = TRUE)

# Kruskal-Wallis (alternative to one-way ANOVA)
kruskal.test(y ~ group, data = df)

# Friedman (alternative to repeated measures ANOVA)
friedman.test(y ~ group | subject, data = df)

# permutation test
library(coin)
oneway_test(y ~ group, data = df, distribution = approximate(nresample = 9999))

# bootstrap CI
library(boot)
boot_mean <- function(data, idx) mean(data[idx])
b <- boot(x, boot_mean, R = 10000)
boot.ci(b, type = "perc")

# sign test
library(BSDA)
SIGN.test(x, md = 5)              # median different from 5?

# simple linear
model <- lm(mpg ~ wt, data = mtcars)
summary(model)
coef(model)
confint(model)
fitted(model)
residuals(model)
predict(model, newdata = new_df, interval = "confidence")

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

# interactions
model <- lm(mpg ~ wt * hp, data = mtcars)
model <- lm(mpg ~ wt + hp + wt:hp, data = mtcars)

# transformations
model <- lm(log(mpg) ~ wt, data = mtcars)
model <- lm(mpg ~ poly(wt, 2), data = mtcars)  # quadratic
model <- lm(mpg ~ I(wt^2) + wt, data = mtcars)

# categorical predictors
model <- lm(mpg ~ factor(cyl), data = mtcars)

# diagnostics
plot(model)                       # 4 diagnostic plots
library(car)
vif(model)                        # variance inflation factor

# logistic regression (binary)
model <- glm(am ~ wt + hp, data = mtcars, family = binomial)
summary(model)

# predicted probabilities
predict(model, type = "response")

# odds ratios
exp(coef(model))

# Poisson regression (counts)
model <- glm(count ~ group, data = df, family = poisson)

# negative binomial (overdispersed counts)
library(MASS)
model <- glm.nb(count ~ group, data = df)

# compare models
model1 <- glm(y ~ x1, family = binomial, data = df)
model2 <- glm(y ~ x1 + x2, family = binomial, data = df)
anova(model1, model2, test = "Chisq")

# goodness of fit
1 - pchisq(model$deviance, model$df.residual)

# McFadden pseudo R²
1 - model$deviance / model$null.deviance

# stepwise selection
full <- lm(mpg ~ ., data = mtcars)
step(full, direction = "both")
step(full, direction = "backward")
step(full, direction = "forward", scope = list(lower = ~1, upper = full))

# AIC and BIC
AIC(model1, model2)
BIC(model1, model2)

# all subsets
library(leaps)
leaps <- regsubsets(mpg ~ ., data = mtcars, nvmax = 10)
plot(leaps, scale = "adjr2")
plot(leaps, scale = "Cp")

# cross-validation
library(caret)
train(mpg ~ ., data = mtcars, method = "lm",
      trControl = trainControl(method = "cv", number = 10))

# regularized (glmnet)
library(glmnet)
X <- model.matrix(mpg ~ ., mtcars)[, -1]
y <- mtcars$mpg
cv_model <- cv.glmnet(X, y, alpha = 0)  # ridge
cv_model <- cv.glmnet(X, y, alpha = 1)  # lasso

library(lme4)
library(lmerTest)

# random intercept
model <- lmer(score ~ treatment + (1 | subject), data = df)

# random intercept and slope
model <- lmer(score ~ time + (time | subject), data = df)

# nested random effects
model <- lmer(score ~ 1 + (1 | school/class), data = df)

# crossed random effects
model <- lmer(score ~ 1 + (1 | subject) + (1 | item), data = df)

# generalized mixed model
glmer(outcome ~ treatment + (1 | subject), data = df, family = binomial)

# summary and CIs
summary(model)
confint(model, method = "Wald")
confint(model, method = "boot")

# anova
anova(model)

# random effects
ranef(model)
VarCorr(model)

# predict
predict(model, newdata = df, allow.new.levels = TRUE)

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

# diagnostic plots
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))

# residuals vs leverage
plot(model, which = 5)

# influence measures
influence.measures(model)
hatvalues(model)                   # leverage
cooks.distance(model)              # Cook's distance
dfbeta(model)                      # change in coef per obs

# identify outliers
which(cooks.distance(model) > 4 / nrow(mtcars))

# predictions
new_data <- data.frame(wt = c(2, 3, 4), hp = c(100, 150, 200))
predict(model, newdata = new_data, interval = "confidence")
predict(model, newdata = new_data, interval = "prediction")

# marginal effects
library(ggeffects)
ggeffect(model, terms = "wt")
plot(ggeffect(model, terms = c("wt", "hp")))

# R squared variants
library(rsq)
rsq(model, type = "v")             # variance-based
rsq(model, type = "kl")            # Kullback-Leibler

---
title: "Quarterly Report"
author: "Data Team"
date: "2024-12-31"
output: html_document
---

## Section

Inline code: the mean is `r mean(x)`.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(dplyr)
library(ggplot2)
data <- read_csv("data.csv")
```

```{r plot, fig.cap="Sales over time"}
ggplot(data, aes(date, sales)) + geom_line()
```

```{r table}
library(knitr)
kable(head(data), caption = "First 6 rows")
```

```{r chunk-name, options}
# code here
```

# Common options:
#   echo=FALSE       hide code, show output
#   eval=FALSE       show code, don't run
#   include=FALSE    run, hide code and output
#   results="hide"   show code, hide text output
#   warning=FALSE    hide warnings
#   message=FALSE    hide messages
#   fig.width=6      figure width (inches)
#   fig.height=4     figure height
#   fig.cap="..."    figure caption
#   fig.path="figs/" where to save figures
#   cache=TRUE       cache results (faster re-knits)
#   dependson="..."  cache dependencies
#   dev="svg"        output device

# global options
knitr::opts_chunk$set(
  echo = TRUE,
  fig.width = 6,
  fig.height = 4,
  fig.path = "figures/",
  cache = TRUE
)

# language engines
```{python}
# Python code
```
```{bash}
# Shell commands
```
```{sql, connection=db}
# SQL queries
```

---
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    number_sections: true
    theme: flatly
    highlight: tango
    code_folding: hide
    df_print: paged
---

---
output:
  pdf_document:
    toc: true
    number_sections: true
    fig_caption: true
    latex_engine: xelatex
---

---
output:
  word_document:
    toc: true
    reference_docx: template.docx
---

---
output:
  ioslides_presentation:
    widescreen: true
---

---
output:
  revealjs::revealjs_presentation:
    theme: solarized
    transition: slide
---

# Multiple formats
---
output:
  html_document: default
  pdf_document: default
---

library(knitr)
library(kableExtra)
library(gt)

# basic kable
kable(head(mtcars))

# styled kable
kable(head(mtcars), format = "html", caption = "Top cars") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  add_header_above(c(" " = 1, "Specs" = 3, "Performance" = 7))

# conditional formatting
kable(df) %>%
  cell_spec(value, color = ifelse(value > 0, "green", "red")) %>%
  row_spec(1, bold = TRUE, background = "yellow")

# gt (modern alternative)
mtcars %>%
  head() %>%
  gt() %>%
  tab_header(title = "MT Cars", subtitle = "First 6 rows") %>%
  cols_label(mpg = "MPG", cyl = "Cylinders") %>%
  tab_options(table.width = pct(80))

# DT for interactive tables
library(DT)
datatable(mtcars, filter = "top", options = list(pageLength = 5))

---
title: "Regional Sales Report"
output: html_document
params:
  region: "West"
  year: 2024
  data_file: "sales.csv"
---

## Report for `r params$region` in `r params$year`

```{r}
data <- read_csv(params$data_file) %>%
  filter(region == params$region, year == params$year)
```

# Render from R
rmarkdown::render("report.Rmd",
                  params = list(region = "East", year = 2024),
                  output_file = "east_2024.html")

# Render from command line
# Rscript -e 'rmarkdown::render("report.Rmd", params = list(region = "East"))'

# Loop over parameters
for (r in c("West", "East", "North")) {
  rmarkdown::render("report.Rmd",
                    params = list(region = r),
                    output_file = paste0(r, "_report.html"))
}

library(purrr)

# map: list output
map(1:5, ~ .x^2)                  # list(1, 4, 9, 16, 25)
map_dbl(1:5, ~ .x^2)              # numeric vector
map_chr(1:5, ~ paste0("n", .x))   # character vector
map_int(1:5, ~ .x * 2L)           # integer vector
map_lgl(1:5, ~ .x > 3)            # logical vector
map_dfr(1:3, ~ data.frame(id = .x, val = .x^2))  # row-bound df
map_dfc(1:3, ~ data.frame(x = .x))               # col-bound df

# multiple arguments
map2(1:3, 4:6, ~ .x + .y)
pmap(list(1:3, 4:6, 7:9), ~ ..1 + ..2 + ..3)

# walk (for side effects)
walk(1:5, ~ print(.x))
walk(files, ~ process(.x))
walk2(plots, filenames, ~ ggsave(.y, .x))

# in pipelines
mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .x)) %>%
  map(summary) %>%
  map_dbl(~ .x$r.squared)

library(purrr)

# safely: capture errors without stopping
safe_log <- safely(log)
safe_log(10)                       # list(result = 2.3, error = NULL)
safe_log("a")                      # list(result = NULL, error = <error>)

# process many, keep going on error
results <- map(values, safely(my_function))
successes <- map(results, "result") %>% compact()
errors <- map(results, "error") %>% compact()

# quietly: capture messages/warnings
quiet_log <- quietly(log)
quiet_log(10)                      # list(result, output, warnings, messages)

# possibly: return default on error
safe_log <- possibly(log, otherwise = NA_real_)
map(values, safe_log)              # never errors

# transpose for cleaner output
results <- map(values, safely(fun)) %>% transpose()
results$result                     # all results
results$error                      # all errors

# insistently: retry on failure
slow_fn <- insistently(api_call, rate = rate_backoff())

library(purrr)

# reduce: combine pairwise
reduce(c(1, 2, 3, 4), )        # 10
reduce(c(1, 2, 3, 4), ~ .x * .y)  # 24
reduce(list(df1, df2, df3), full_join, by = "id")

# accumulate: keep intermediate results
accumulate(c(1, 2, 3, 4), )    # 1 3 6 10
accumulate(c(2, 3, 4), ~ .x * .y) # 2 6 24

# right reduce
reduce(c(1, 2, 3, 4), , .dir = "backward")

# with init
reduce(c(1, 2, 3), ~ c(.x, .y), .init = numeric(0))

# detect
detect(1:10, ~ .x > 5)            # 6 (first match)
detect_index(1:10, ~ .x > 5)      # 6 (index)
detect(1:10, ~ .x > 5, .right = TRUE)  # 10 (last match)

# keep/discard
keep(1:10, ~ .x %% 2 == 0)        # 2 4 6 8 10
discard(1:10, ~ .x %% 2 == 0)     # 1 3 5 7 9
every(1:10, ~ .x > 0)             # TRUE
some(1:10, ~ .x > 5)              # TRUE

library(purrr)

# partial: pre-fill arguments
add_one <- partial(, 1)
add_one(5)                         # 6

round2 <- partial(round, digits = 2)
round2(3.14159)                    # 3.14

# compose: chain functions
clean_string <- compose(
  str_trim,
  str_to_lower,
  ~ str_replace_all(.x, "[^a-z]", " ")
)
clean_string("  Hello, World!  ")  # "hello  world "

# %>% (magrittr pipe) vs compose
# pipe: data %>% f %>% g %>% h
# compose: g %>% f creates a new function

# negate
is_missing <- negate(is.na)
is_missing(5)                      # TRUE

# lift (vectorize)
sum2 <- lift()                  # takes a list of 2 args
sum2(list(1, 2))                   # 3

# walk with index
imap(letters[1:3], ~ paste0(.y, ": ", .x))
# list("1: a", "2: b", "3: c")

library(purrr)
library(dplyr)
library(tidyr)

# nest data, fit models, predict
mtcars %>%
  group_by(cyl) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(mpg ~ wt, data = .x)),
    rsq = map_dbl(model, ~ summary(.x)$r.squared),
    pred = map2(model, data, predict)
  )

# unnest predictions
mtcars %>%
  group_by(cyl) %>%
  nest() %>%
  mutate(model = map(data, ~ lm(mpg ~ wt, .x))) %>%
  mutate(pred = map2(model, data, predict)) %>%
  select(cyl, data, pred) %>%
  unnest(c(data, pred))

# apply function to each column
mtcars %>% map_dbl(mean)
mtcars %>% map_dbl(sd)

# apply to columns of specific type
iris %>%
  map_if(is.numeric, mean) %>%
  map_dbl(round, 2)

# pluck
list(a = list(b = list(c = 42))) %>% pluck("a", "b", "c")
# 42

# modify (returns same type)
modify_if(iris, is.numeric, ~ .x * 2)
modify_depth(nested_list, 2, ~ .x + 1)

library(shiny)

ui <- fluidPage(
  titlePanel("Hello Shiny"),
  sidebarLayout(
    sidebarPanel(
      sliderInput("n", "Number of points", 1, 100, 50)
    ),
    mainPanel(
      plotOutput("scatter")
    )
  )
)

server <- function(input, output, session) {
  output$scatter <- renderPlot({
    plot(1:input$n, rnorm(input$n))
  })
}

shinyApp(ui, server)

# Save as app.R in a folder, run with:
#   shiny::runApp("path/to/folder")
# Or click "Run App" in RStudio

library(shiny)

ui <- fluidPage(
  numericInput("n", "N", value = 10, min = 1, max = 100),
  textInput("label", "Label", value = "Data"),
  selectInput("color", "Color", choices = c("red", "blue", "green")),
  dateInput("date", "Date"),
  dateRangeInput("range", "Range"),
  checkboxInput("show", "Show?", value = TRUE),
  checkboxGroupInput("vars", "Variables", choices = names(mtcars)),
  radioButtons("dist", "Distribution",
               choices = c("Normal", "Uniform", "Exponential")),
  fileInput("file", "Upload CSV", accept = ".csv"),
  actionButton("go", "Go!")
)

server <- function(input, output, session) {
  # event-triggered
  observeEvent(input$go, {
    showNotification(paste("Clicked", input$go))
  })

  # reactive expression (cached)
  data <- reactive({
    rnorm(input$n)
  })

  # reactive value
  val <- reactiveVal(0)
  observeEvent(input$go, val(val() + 1))
}

library(shiny)
library(ggplot2)
library(DT)

ui <- fluidPage(
  plotOutput("plot", click = "plot_click"),
  tableOutput("table"),
  DT::dataTableOutput("dt"),
  verbatimTextOutput("summary"),
  textOutput("text"),
  uiOutput("dynamic")
)

server <- function(input, output, session) {
  output$plot <- renderPlot({
    ggplot(mtcars, aes(wt, mpg)) + geom_point()
  })

  output$table <- renderTable({
    head(mtcars)
  })

  output$dt <- DT::renderDataTable({
    datatable(mtcars, filter = "top")
  })

  output$summary <- renderPrint({
    summary(mtcars)
  })

  output$text <- renderText({
    paste("You clicked:", input$plot_click$x)
  })

  output$dynamic <- renderUI({
    if (input$n > 5) strong("Big!") else em("Small")
  })
}

# observe vs reactive
#   observe: side effects (output$ assignments, updates)
#   reactive: returns a value, used by other reactives

library(shiny)

server <- function(input, output, session) {
  # reactive: cached, lazy
  filtered <- reactive({
    mtcars %>% filter(cyl == input$cyl)
  })

  # observe: eager, side effects
  observe({
    updateSelectInput(session, "model",
                      choices = unique(filtered()$model))
  })

  # eventReactive: triggered by event
  result <- eventReactive(input$go, {
    run_analysis(input$params)
  })

  # reactiveVal: single mutable value
  counter <- reactiveVal(0)
  observeEvent(input$add, counter(counter() + 1))

  # reactiveValues: multiple values
  rv <- reactiveValues(data = NULL, status = "idle")
  observeEvent(input$load, {
    rv$data <- read.csv(input$file$datapath)
    rv$status <- "loaded"
  })

  # debounce (throttle rapid changes)
  search <- reactive({ input$search }) %>% debounce(500)

  # isolate (use value without dependency)
  output$plot <- renderPlot({
    plot(isolate(rv$data))  # not reactive on rv$data
  })
}

# Deploy to shinyapps.io
#   1. Create account at shinyapps.io
#   2. Install rsconnect
install.packages("rsconnect")
#   3. Authorize (paste token from shinyapps.io)
rsconnect::setAccountInfo(name = "<name>", token = "<token>", secret = "<secret>")
#   4. Deploy
rsconnect::deployApp("path/to/app")

# Run on local network
shiny::runApp("app.R", host = "0.0.0.0", port = 3838)

# Shiny Server (self-hosted on Linux)
# Install shiny-server, put apps in /srv/shiny-server/
# Config: /etc/shiny-server/shiny-server.conf
# Access at http://server:3838/appname

# Docker
#   FROM rocker/shiny
#   COPY app.R /srv/shiny-server/myapp/
#   EXPOSE 3838

# Performance tips:
#   - Use renderCachedPlot for expensive plots
#   - Use future + promises for async work
#   - Move heavy compute to reactive({...})
#   - Avoid re-reading files in render functions

# Debugging
options(shiny.fullstacktrace = TRUE)
options(shiny.reactlog = TRUE)  # press Ctrl+F3 in app

library(profvis)

# profile a block of code
profvis({
  data <- read.csv("large.csv")
  result <- lapply(data, function(x) {
    x[x > 0]
  })
  plot(result)
})

# profile a function call
profvis(my_function(arg1, arg2))

# save profile
p <- profvis({ ... })
htmlwidgets::saveWidget(p, "profile.html")

# Rprof (base R)
Rprof("profile.out")
# ... code to profile ...
Rprof(NULL)
summaryRprof("profile.out")

# tips:
#   - profile realistic inputs (not too small)
#   - run multiple times for stable results
#   - look for the widest bars in the flame graph
#   - focus on hotspots, not micro-optimizations

# BAD: loop with growing result
slow <- function(n) {
  result <- numeric(0)
  for (i in 1:n) {
    result <- c(result, i^2)
  }
  result
}

# BETTER: preallocate
better <- function(n) {
  result <- numeric(n)
  for (i in 1:n) {
    result[i] <- i^2
  }
  result
}

# BEST: vectorize
fast <- function(n) {
  (1:n)^2
}

# benchmarks
microbenchmark::microbenchmark(
  slow = slow(1000),
  better = better(1000),
  fast = fast(1000),
  times = 10
)

# apply family
#   sapply/lapply: list apply
#   vapply: typed sapply (safer, faster)
#   map_dbl: tidyverse, typed
#   rowMeans/colSums: faster than apply(x, 1, mean)

library(Rcpp)

# inline C++ in R
cppFunction('
  double sumC(NumericVector x) {
    double s = 0;
    for (int i = 0; i < x.size(); i++) {
      s += x[i];
    }
    return s;
  }
')
sumC(1:1e6)

# source from file
# file: code.cpp
# #include <Rcpp.h>
# using namespace Rcpp;
# // [[Rcpp::export]]
# double meanC(NumericVector x) {
#   return std::accumulate(x.begin(), x.end(), 0.0) / x.size();
# }

sourceCpp("code.cpp")

# use Rcpp sugar (vectorized C++)
cppFunction('
  NumericVector funC(NumericVector x) {
    return x * 2 + 1;  // vectorized via Rcpp sugar
  }
')

# data frames
cppFunction('
  DataFrame subsetC(DataFrame df, LogicalVector mask) {
    return df[mask];
  }
')

library(data.table)

# data.table is much faster than data.frame for large data
dt <- fread("huge.csv")            # fast CSV reader
fwrite(dt, "out.csv")              # fast CSV writer

# syntax: dt[i, j, by]
dt[, mean(value), by = group]      # group by + summarize
dt[order(-value)]                  # sort
dt[value > 100, .N, by = group]    # filter + count by group
dt[, .(avg = mean(value), n = .N), by = group]

# reference semantics (no copy)
dt[, new_col := value * 2]         # add column in place
dt[value < 0, value := 0]          # modify in place

# keys for fast lookups
setkey(dt, id)
dt["abc"]                          # fast lookup
dt["abc", mult = "first"]

# joins
dt1[dt2, on = "id"]                # left join
dt1[dt2, on = .(id), nomatch = 0]  # inner join

# memory tips
#   - gc() to force garbage collection
#   - object.size(x) to check size
#   - rm() large objects when done
#   - read data in chunks for huge files

library(parallel)
library(future)
library(furrr)

# parallel lapply
cl <- makeCluster(detectCores() - 1)
result <- parLapply(cl, items, function(x) {
  # work on x
})
stopCluster(cl)

# future: modern, simpler
library(future)
plan(multisession)                 # parallel backend
result <- future_lapply(items, fun)
# or
plan(multicore)                    # fork (Linux/Mac, faster)
plan(cluster, workers = 4)

# furrr: parallel purrr
library(furrr)
plan(multisession, workers = 4)
result <- future_map(items, fun)
result <- future_map_dbl(items, ~ .x^2)

# foreach
library(foreach)
library(doParallel)
registerDoParallel(4)
result <- foreach(i = 1:10, .combine = c) %dopar% {
  i^2
}

# always benchmark — parallel overhead can outweigh gains
# for small tasks
microbenchmark::microbenchmark(
  serial = lapply(1:100, slow_fn),
  parallel = future_lapply(1:100, slow_fn),
  times = 5
)