中文译文

基础 R 编程技能

一份关于基础 R 编程的全面参考——涵盖数据结构、控制流、函数、I/O、统计计算和绘图。

快速参考

数据结构

# Vectors (atomic)
x <- c(1, 2, 3)              # numeric
y <- c("a", "b", "c")        # character
z <- c(TRUE, FALSE, TRUE)    # logical

# Factor
f <- factor(c("low", "med", "high"), levels = c("low", "med", "high"), ordered = TRUE)

# Matrix
m <- matrix(1:6, nrow = 2, ncol = 3)
m[1, ]       # first row
m[, 2]       # second column

# List
lst <- list(name = "ali", scores = c(90, 85), passed = TRUE)
lst$name      # access by name
lst[[2]]      # access by position

# Data frame
df <- data.frame(
  id = 1:3,
  name = c("a", "b", "c"),
  value = c(10.5, 20.3, 30.1),
  stringsAsFactors = FALSE
)
df[df$value > 15, ]    # filter rows
df$new_col <- df$value * 2  # add column

子集化

# Vectors
x[1:3]             # by position
x[c(TRUE, FALSE)]  # by logical
x[x > 5]           # by condition
x[-1]              # exclude first

# Data frames
df[1:5, ]                    # first 5 rows
df[, c("name", "value")]     # select columns
df[df$value > 10, "name"]    # filter + select
subset(df, value > 10, select = c(name, value))

# which() for index positions
idx <- which(df$value == max(df$value))

控制流

# if/else
if (x > 0) {
  "positive"
} else if (x == 0) {
  "zero"
} else {
  "negative"
}

# ifelse (vectorized)
ifelse(x > 0, "pos", "neg")

# for loop
for (i in seq_along(x)) {
  cat(i, x[i], "\n")
}

# while
while (condition) {
  # body
  if (stop_cond) break
}

# switch
switch(type,
  "a" = do_a(),
  "b" = do_b(),
  stop("Unknown type")
)

函数

# Define
my_func <- function(x, y = 1, ...) {
  result <- x + y
  return(result)  # or just: result
}

# Anonymous functions
sapply(1:5, function(x) x^2)
# R 4.1+ shorthand:
sapply(1:5, \(x) x^2)

# Useful: do.call for calling with a list of args
do.call(paste, list("a", "b", sep = "-"))

Apply 家族

# sapply — simplify result to vector/matrix
sapply(lst, length)

# lapply — always returns list
lapply(lst, function(x) x[1])

# vapply — like sapply but with type safety
vapply(lst, length, integer(1))

# apply — over matrix margins (1=rows, 2=cols)
apply(m, 2, sum)

# tapply — apply by groups
tapply(df$value, df$group, mean)

# mapply — multivariate
mapply(function(x, y) x + y, 1:3, 4:6)

# aggregate — like tapply for data frames
aggregate(value ~ group, data = df, FUN = mean)

字符串操作

paste("a", "b", sep = "-")    # "a-b"
paste0("x", 1:3)              # "x1" "x2" "x3"
sprintf("%.2f%%", 3.14159)    # "3.14%"
nchar("hello")                # 5
substr("hello", 1, 3)         # "hel"
gsub("old", "new", text)      # replace all
grep("pattern", x)            # indices of matches
grepl("pattern", x)           # logical vector
strsplit("a,b,c", ",")        # list("a","b","c")
trimws("  hi  ")              # "hi"
tolower("ABC")                # "abc"

数据输入/输出

# CSV
df <- read.csv("data.csv", stringsAsFactors = FALSE)
write.csv(df, "output.csv", row.names = FALSE)

# Tab-delimited
df <- read.delim("data.tsv")

# General
df <- read.table("data.txt", header = TRUE, sep = "\t")

# RDS (single R object, preserves types)
saveRDS(obj, "data.rds")
obj <- readRDS("data.rds")

# RData (multiple objects)
save(df1, df2, file = "data.RData")
load("data.RData")

# Connections
con <- file("big.csv", "r")
chunk <- readLines(con, n = 100)
close(con)

基础绘图

# Scatter
plot(x, y, main = "Title", xlab = "X", ylab = "Y",
     pch = 19, col = "steelblue", cex = 1.2)

# Line
plot(x, y, type = "l", lwd = 2, col = "red")
lines(x, y2, col = "blue", lty = 2)  # add line

# Bar
barplot(table(df$category), main = "Counts",
        col = "lightblue", las = 2)

# Histogram
hist(x, breaks = 30, col = "grey80",
     main = "Distribution", xlab = "Value")

# Box plot
boxplot(value ~ group, data = df,
        col = "lightyellow", main = "By Group")

# Multiple plots
par(mfrow = c(2, 2))  # 2x2 grid
# ... four plots ...
par(mfrow = c(1, 1))  # reset

# Save to file
png("plot.png", width = 800, height = 600)
plot(x, y)
dev.off()

# Add elements
legend("topright", legend = c("A", "B"),
       col = c("red", "blue"), lty = 1)
abline(h = 0, lty = 2, col = "grey")
text(x, y, labels = names, pos = 3, cex = 0.8)

统计信息

# Descriptive
mean(x); median(x); sd(x); var(x)
quantile(x, probs = c(0.25, 0.5, 0.75))
summary(df)
cor(x, y)
table(df$category)  # frequency table

# Linear model
fit <- lm(y ~ x1 + x2, data = df)
summary(fit)
coef(fit)
predict(fit, newdata = new_df)
confint(fit)

# t-test
t.test(x, y)                    # two-sample
t.test(x, mu = 0)               # one-sample
t.test(before, after, paired = TRUE)

# Chi-square
chisq.test(table(df$a, df$b))

# ANOVA
fit <- aov(value ~ group, data = df)
summary(fit)
TukeyHSD(fit)

# Correlation test
cor.test(x, y, method = "pearson")

数据操作

# Merge (join)
merged <- merge(df1, df2, by = "id")                  # inner
merged <- merge(df1, df2, by = "id", all = TRUE)      # full outer
merged <- merge(df1, df2, by = "id", all.x = TRUE)    # left

# Reshape
wide <- reshape(long, direction = "wide",
                idvar = "id", timevar = "time", v.names = "value")
long <- reshape(wide, direction = "long",
                varying = list(c("v1", "v2")), v.names = "value")

# Sort
df[order(df$value), ]              # ascending
df[order(-df$value), ]             # descending
df[order(df$group, -df$value), ]   # multi-column

# Remove duplicates
df[!duplicated(df), ]
df[!duplicated(df$id), ]

# Stack / combine
rbind(df1, df2)    # stack rows (same columns)
cbind(df1, df2)    # bind columns (same rows)

# Transform columns
df$log_val <- log(df$value)
df$category <- cut(df$value, breaks = c(0, 10, 20, Inf),
                   labels = c("low", "med", "high"))

环境与调试

ls()                  # list objects
rm(x)                 # remove object
rm(list = ls())       # clear all
str(obj)              # structure
class(obj)            # class
typeof(obj)           # internal type
is.na(x)              # check NA
complete.cases(df)    # rows without NA
traceback()           # after error
debug(my_func)        # step through
browser()             # breakpoint in code
system.time(expr)     # timing
Sys.time()            # current time

参考文件

如需更深入的了解，请阅读 references/ 中的参考文件：

函数陷阱与快速参考（摘自 R 4.5.3 参考手册）

不明显的行为、令人意外的默认值以及棘手的交互——只包含 Claude 尚不知道的内容：

data-wrangling.md — 阅读时机：子集提取返回错误的类型，对数据框使用 apply 时产生意外的强制转换，merge/split/cbind 行为异常，过滤后因子水平仍然保留，table/duplicated 的边界情况。
modeling.md — 阅读时机：公式语法令人困惑（I()、* 与 :、/），aov 给出错误的平方和类型，glm 静默地拟合 OLS，nls 无法收敛，predict 返回错误的尺度，optim/optimize 需要调参。
statistics.md — 阅读时机：假设检验给出令人意外的结果，需要选择正确的 p.adjust 方法，聚类参数似乎有误，分布函数命名令人困惑（d/p/q/r 前缀）。
visualization.md — 阅读时机：par 设置意外重置，layout/mfrow 的交互令人困惑，坐标轴标签被裁剪，颜色看起来不对，需要专业图形（contour、persp、mosaic、pairs）。
io-and-text.md — 阅读时机：read.table 静默丢弃数据或错误解析列，正则表达式行为与预期不同，sprintf 格式化很棘手，write.table 输出带有不需要的行名。
dates-and-system.md — 阅读时机：Date/POSIXct 转换给出错误的日期，时区导致差一错误，difftime 单位出乎意料，需要以编程方式查找/列出/测试文件。
misc-utilities.md — 阅读时机：do.call 的行为与直接调用不同，需要 Reduce/Filter/Map，tryCatch 处理程序未触发，all.equal 返回字符串而非逻辑值，时间序列函数需要设置。

编写优秀 R 代码的技巧

在生产代码中使用 vapply() 而非 sapply()——它会强制返回类型
优先使用 seq_along(x) 而非 1:length(x)——当 x 为空时后者会出错
在 read.csv() / data.frame() 中使用 stringsAsFactors = FALSE（默认值在 R 4.0 中已更改）
尽可能向量化操作，而不是编写循环
使用 stop()、warning()、message() 进行错误处理——而不是 print()
<<- 赋值给父环境——请谨慎且有意识地使用
with(df, expr) 避免到处重复 df$
使用 Sys.setenv() 和 .Renviron 设置环境变量 FILE:references/misc-utilities.md

杂项工具——快速参考

R 函数中不明显的行为、陷阱和棘手的默认值。只包含 Claude 还不知道的内容。

do.call

do.call(fun, args_list) — args 必须是列表，即使只有一个参数也是如此。
quote = TRUE 可防止在调用前对参数求值——在传递表达式/符号时需要。
do.call 内部 substitute 的行为与直接调用不同。这种情况的语义并未完全定义。
有用的模式：do.call(rbind, list_of_dfs) 用于合并数据框列表。

Reduce / Filter / Map / Find / Position

R 基础包中的函数式编程辅助函数——确实很不明显。

Reduce(f, x) 累积地应用二元函数 f：Reduce("+", 1:4) = ((1+2)+3)+4。对于非交换运算，方向很重要。
Reduce(f, x, accumulate = TRUE) 返回所有中间结果——等价于 Python 的 itertools.accumulate。
Reduce(f, x, right = TRUE) 从右侧折叠：f(x1, f(x2, f(x3, x4)))。
带 init 的 Reduce 会添加一个起始值：Reduce(f, x, init = v) = f(f(f(v, x1), x2), x3)。
Filter(f, x) 保留 f(elem) 为 TRUE 的元素。与 x[sapply(x, f)] 不同，它能正确处理 NULL/空值。
Map(f, ...) 是 mapply(f, ..., SIMPLIFY = FALSE) 的简单包装——始终返回列表。
Find(f, x) 返回 f(elem) 为 TRUE 的第一个元素。Find(f, x, right = TRUE) 用于最后一个。
Position(f, x) 返回第一个匹配项的索引（类似 Find，但返回位置而非值）。

lengths

lengths(x) 返回列表每个元素的长度。等价于 sapply(x, length) 但更快（用 C 实现）。
适用于任何类列表对象。返回整数向量。

conditions (tryCatch / withCallingHandlers)

tryCatch 展开调用栈——处理程序在调用环境中运行，而不是在错误发生处。无法恢复执行。
withCallingHandlers 不会展开——处理程序在条件被发出处运行。可以检查/记录，然后让条件继续传播。
tryCatch(expr, error = function(e) e) 返回错误条件对象。
tryCatch(expr, warning = function(w) {...}) 捕获第一个警告并退出。使用 withCallingHandlers + invokeRestart("muffleWarning") 来抑制警告但继续执行。
tryCatch 的 finally 子句始终运行（类似 Java 的 try/finally）。
globalCallingHandlers() 注册在会话期间持续存在的处理程序（对日志记录有用）。
自定义条件：stop(errorCondition("msg", class = "myError"))，然后用 tryCatch(..., myError = function(e) ...) 捕获。

all.equal

使用容差测试近似相等（默认 1.5e-8，即 sqrt(.Machine$double.eps)）。
返回 TRUE 或描述差异的字符串——而不是 FALSE。在条件判断中使用 isTRUE(all.equal(x, y))。
tolerance 参数控制数值容差。scale 用于绝对比较与相对比较。
检查属性、名称、维度——比 == 更彻底。

combn

combn(n, m) 或 combn(x, m)：生成从 x 中取 m 个项的所有组合。
返回一个矩阵，有 m 行；每一列是一个组合。
FUN 参数对每个组合应用一个函数：combn(5, 3, sum) 返回所有 3 元素子集的和。
simplify = FALSE 返回列表而不是矩阵。

modifyList

modifyList(x, val) 按名称用 val 中的元素替换列表 x 中的元素。
将某个值设为 NULL 会从列表中移除该元素。
会添加 x 中不存在的新名称——它在内部使用 x[names(val)] <- val，因此 val 中的任何名称都会被添加或替换。

relist

unlist 的逆操作：给定一个扁平向量和一个骨架列表，重建嵌套结构。
relist(flesh, skeleton)——flesh 是扁平数据，skeleton 提供形状。
适用于因子、矩阵和嵌套列表。

txtProgressBar

txtProgressBar(min, max, style = 3)——样式 3 显示百分比 + 进度条（最有用）。
用 setTxtProgressBar(pb, value) 更新。用 close(pb) 关闭。
样式 1：旋转的 |/-\，样式 2：简单进度。只有样式 3 显示百分比。

object.size

返回对象所用内存的估计值。对于共享引用并不总是精确。
format(object.size(x), units = "MB") 用于人类可读的输出。
不计算环境或外部指针的大小。

installed.packages / update.packages

installed.packages() 可能很慢（会扫描所有包）。若要检查某个特定包，请使用 find.package() 或 requireNamespace()。
update.packages(ask = FALSE) 会在不提示的情况下更新所有包。
lib.loc 指定要检查/更新哪个库。

vignette / demo

vignette() 列出所有 vignette；vignette("name", package = "pkg") 打开指定的一个。
demo() 列出所有演示；demo("topic") 以交互方式运行其中一个。
browseVignettes() 在 HTML 中打开 vignette 浏览器。

时间序列：acf / arima / ts / stl / decompose

ts(data, start, frequency)：frequency 是每单位时间内的观测数（月度用 12，季度用 4）。
acf 默认 type = "correlation"。PACF 请使用 type = "partial"。plot = FALSE 用于抑制自动绘图。
arima(x, order = c(p,d,q)) 用于 ARIMA 模型。seasonal = list(order = c(P,D,Q), period = S) 用于季节成分。
arima 能处理时间序列中的 NA 值（通过卡尔曼滤波）。
stl 需要 s.window（季节窗口）——必须指定，没有默认值。s.window = "periodic" 假定季节性固定。
decompose：比 stl 更简单，使用移动平均。type = "additive" 或 "multiplicative"。
stl 结果组件：$time.series 矩阵，其列为 seasonal、trend、remainder。 FILE:references/data-wrangling.md

数据整理 — 快速参考

R 函数中不明显的行文、陷阱和棘手的默认值。只包含 Claude 尚不知道的内容。

Extract / Extract.data.frame

base R 中的索引陷阱。

m[j = 2, i = 1] 是 m[2, 1] 而不是 m[1, 2]——在 [ 中参数名会被忽略，只按位置匹配。绝不要给索引参数命名。
因子索引：x[f] 使用因子 f 的整数编码，而不是其字符标签。若要基于标签索引，请使用 x[as.character(f)]。
x[[]] 不带索引总是会报错。x$name 默认会进行部分匹配；x[["name"]] 不会（默认精确匹配）。
通过 x[[i]] <- NULL 或 x$name <- NULL 赋值 NULL 会删除该列表元素。
数据框 [ 取单列时：df[, 1] 返回一个向量（列默认 drop=TRUE），但 df[1, ] 返回一个数据框（行默认 drop=FALSE）。请显式使用 drop = FALSE。
对数据框进行矩阵索引（df[cbind(i,j)]）会先强制转换为矩阵——避免这样做。

subset

仅用于交互式场景；用于编程不安全。

subset 参数使用非标准求值——列名在数据框中解析，在编程使用中可能悄无声息地选错变量。在函数中请使用带显式逻辑的 [。
逻辑条件中的 NA 会被当作 FALSE（行会被悄无声息地丢弃）。
因子在取子集后可能保留未使用的水平；请调用 droplevels()。

match / %in%

%in% 绝不返回 NA——这使其在 if() 条件中比 == 更安全。
match() 只返回第一个匹配的位置；table 中的重复项会被忽略。
因子、原始向量和列表在匹配前都会转换为字符。
NaN 匹配 NaN 但不匹配 NA；NA 只匹配 NA。

apply

作用于数据框时，apply 会先通过 as.matrix 强制转换为矩阵——混合类型会变成字符。
返回值的朝向是转置的：如果 FUN 返回长度为 n 的向量，结果的维度为 c(n, dim(X)[MARGIN])。行结果会变成列。
因子结果在输出数组中会被强制转换为字符。
... 参数不能与 X、MARGIN 或 FUN 同名（存在部分匹配风险）。

lapply / sapply / vapply

sapply 可能不可预测地返回向量、矩阵或列表——在非交互式代码中请使用 vapply 并显式提供 FUN.VALUE 模板。
在 lapply 中直接调用原语可能引发分派问题；请用 function(x) is.numeric(x) 包裹，而不是裸写 is.numeric。
sapply 配合 simplify = "array" 可能产生更高阶的数组（不只是矩阵）。

tapply

返回数组（不是数据框）。返回值上的类信息会被丢弃（例如，Date 对象会变成数值）。
传给 FUN 的 ... 参数不会被划分到各个单元格——它们全局生效，因此 FUN 不应期望有与 X 长度相同的额外参数。
default = NA 填充空单元格；对于求和类操作请设置 default = 0。在 R 3.4.0 之前，这被硬编码为 NA。
使用 array2DF() 将结果转换为数据框。

mapply

参数名是 SIMPLIFY（全大写），不是 simplify——与 sapply 不一致。
MoreArgs 必须是一个列表，其中的参数不会被向量化。
将较短的参数循环回收至共同长度；零长度参数会得到零长度结果。

merge

默认的 by 是 intersect(names(x), names(y))——如果数据框共享列名，可能会在非预期的列上悄悄进行合并。
by = 0 或 by = "row.names" 会按行名合并，并添加一个 "Row.names" 列。
by = NULL（或 by.x/by.y 长度均为 0）会产生笛卡尔积。
结果默认按 by 列排序（sort = TRUE）。若要得到未排序的输出，请使用 sort = FALSE。
重复的键匹配会产生所有组合（每个匹配对一行）。

split

如果 f 是一个因子列表，则会使用交互项；除非更改 sep，否则包含 "." 的水平可能导致意外的拆分。
drop = FALSE（默认）会将空的因子水平保留为空列表元素。
支持公式语法：split(df, ~ Month)。

cbind / rbind

对数据框调用 cbind 会调用 data.frame(...)，而不是 cbind.matrix。混合矩阵和数据框可能产生意外结果。
对数据框调用 rbind 时按名称匹配列，而不是按位置。缺失的列会得到 NA。
cbind(NULL) 返回 NULL（不是矩阵）。为保持一致，rbind(NULL) 也返回 NULL。

table

默认排除 NA（useNA = "no"）。使用 useNA = "ifany" 或 exclude = NULL 来统计 NA。
将 exclude 设为非空且非默认值会隐含 useNA = "ifany"。
结果始终是一个数组（即使是一维），类别为 "table"。用 as.data.frame(tbl) 转换为数据框。
两种 NA（因子水平 NA 与真正的 NA）会根据 useNA/exclude 被区别对待。

duplicated / unique

duplicated 将第二次及之后出现的值标记为 TRUE，而不是第一次。使用 fromLast = TRUE 可反转。
对于数据框，按整行操作。对于列表，则递归比较。
unique 保留每个值的第一次出现。

data.frame（易错点）

自 R 4.0.0 起，stringsAsFactors = FALSE 是默认值（之前是 TRUE）。
原子向量会循环以匹配最长的列，但仅当长度恰好为整数倍时。用 I() 保护以防止转换。
仅当 check.names = FALSE 时才允许重复的列名，但许多操作会静默地对它们去重。
矩阵参数会被展开为多个列，除非用 I() 保护。

factor（易错点）

as.numeric(f) 返回的是整数编码，而不是原始值。请使用 as.numeric(levels(f))[f] 或 as.numeric(as.character(f))。
因子之间只有 == 和 != 可用；因子必须具有完全相同的水平集合。有序因子支持 <、>。
对因子使用 c() 会合并水平集合（自 R 4.1.0 起），但更早的版本会转换为整数。
水平默认会被排序，但排序顺序在创建时依赖于区域设置。

aggregate

公式接口（aggregate(y ~ x, data, FUN)）默认会丢弃 NA 分组。
数据框方法要求 by 是一个列表（而不是向量）。
返回的列以分组变量命名，结果列保留原始名称。
如果 FUN 返回多个值，结果列是数据框内的一个矩阵列。

complete.cases

返回一个逻辑向量：对于所有列/参数中没有 NA 的行，为 TRUE。
可作用于多个参数（例如，complete.cases(x, y) 会同时检查两者）。

order

返回索引的排列向量，而不是排序后的值。请使用 x[order(x)] 来排序。
默认是升序；对数值使用 -x 实现降序，或使用 decreasing = TRUE。
对于字符排序，取决于区域设置。使用 method = "radix" 可实现与区域无关的快速排序。
对于大型整数/字符向量，使用 method = "radix" 的 sort.int() 要快得多。 FILE:references/dates-and-system.md

日期与系统 — 快速参考

R 函数中不明显的行、坑点与棘手的默认行为。只包含 Claude 尚不知道的内容。

日期（Date 类）

Date 对象存储为自 1970-01-01 以来的整数天数。算术以天为单位。
Sys.Date() 以 Date 对象返回当前日期。
seq.Date(from, to, by = "month") —— "month" 增量可能产生长度不等的区间。给 1 月 31 日加 1 个月得到 3 月 3 日（而非 2 月 28 日）。
diff(dates) 返回以天为单位的 difftime 对象。
format(date, "%Y") 取年份，"%m" 取月份，"%d" 取日，"%A" 取星期名（依赖区域设置）。
公元 1 年之前的年份可能无法正确处理。
若扩展，length(date_vector) <- n 会用 NA 填充。

DateTimeClasses（POSIXct / POSIXlt）

POSIXct：自 1970-01-01 UTC 以来的秒数（紧凑，一个数值向量）。
POSIXlt：列表，包含分量 $sec、$min、$hour、$mday、$mon（0-11！）、$year（自 1900 年起！）、$wday（0-6，周日=0）、$yday（0-365）。
在 POSIXct 与 Date 之间转换：as.Date(posixct_obj) 默认使用 tz = "UTC" —— 若原对象处于其他时区，可能得到与预期不同的日期。
Sys.time() 以当前时区返回 POSIXct。
strptime 返回 POSIXlt；用 as.POSIXct(strptime(...)) 得到 POSIXct。
difftime 算术：POSIXct 对象相减得到 difftime。单位自动选择（"secs"、"mins"、"hours"、"days"、"weeks"）。

difftime

difftime(time1, time2, units = "auto") —— 自动选择最小的合理单位。
显式单位："secs"、"mins"、"hours"、"days"、"weeks"。没有 "months" 或 "years"（长度可变）。
as.numeric(diff, units = "hours") 以指定单位提取数值。
units(diff_obj) <- "hours" 就地更改单位。

system.time / proc.time

system.time(expr) 返回 user、system 和 elapsed 时间。
gcFirst = TRUE（默认）：计时前运行垃圾回收，以获得更一致的结果。
proc.time() 返回自 R 启动以来的累计时间 —— 取差值以获得区间。
elapsed（挂钟时间）可能小于 user（多线程 BLAS）或大于它（I/O 等待）。

Sys.sleep

Sys.sleep(seconds) —— 允许小数秒。实际睡眠时间可能更长（操作系统调度）。
进程在睡眠期间会向操作系统让出（不会忙等）。

options（关键选项）

部分不太直观的选项：

options(scipen = n)：正值偏向固定记法，负值偏向科学记法。默认 0。作用于 print/format/cat，但不作用于 sprintf。
options(digits = n)：打印时的有效数字位数（1-22，默认 7）。仅为建议值。
options(digits.secs = n)：时间格式化中秒的最大小数位数（0-6，默认 0）。
options(warn = n)：-1 = 忽略警告，0 = 收集（默认），1 = 立即显示，2 = 转换为错误。
options(error = recover)：出错时进入调试器。options(error = NULL) 重置为默认值。
options(OutDec = ",")：更改输出中的小数分隔符（影响 format、print，不影响 sprintf）。
options(stringsAsFactors = FALSE)：data.frame 的全局默认值（自 R 4.0.0 起已无意义，因为那时它已经是 FALSE）。
options(expressions = 5000)：最大嵌套求值次数。深度递归时可增大。
options(max.print = 99999)：控制 print 输出中的截断。
options(na.action = "na.omit")：模型函数中默认的 NA 处理方式。
options(contrasts = c("contr.treatment", "contr.poly"))：无序/有序因子的默认对比。

file.path / basename / dirname

file.path("a", "b", "c.txt") → "a/b/c.txt"（平台适用的分隔符）。
basename("/a/b/c.txt") → "c.txt"。dirname("/a/b/c.txt") → "/a/b"。
file.path 不会规范化路径（不解析 ..）；如需此功能请使用 normalizePath()。

list.files

list.files(pattern = "*.csv") —— pattern 是正则表达式，不是通配符！请使用 glob2rx("*.csv") 或 "\\.csv$"。
full.names = FALSE（默认）仅返回基本文件名。如需完整路径请使用 full.names = TRUE。
recursive = TRUE 用于搜索子目录。
all.files = TRUE 用于包含隐藏文件（以 . 开头）。

file.info

返回包含 size、isdir、mode、mtime、ctime、atime、uid、gid 的数据框。
mtime：修改时间（POSIXct）。可用于 file.info(f)$mtime。
在某些文件系统上，ctime 是状态变更时间，而不是创建时间。

file_test

file_test("-f", path)：如果常规文件存在则为 TRUE。
file_test("-d", path)：如果目录存在则为 TRUE。
file_test("-nt", f1, f2)：如果 f1 比 f2 新则为 TRUE。
在区分文件和目录方面比 file.exists() 更可靠。 FILE:references/io-and-text.md

I/O 与文本处理 — 快速参考

R 函数中不明显的行为、陷阱和棘手的默认值。只包含 Claude 尚不知道的内容。

read.table（陷阱）

sep = ""（默认值）表示任意空白字符（空格、制表符、换行符）——而不是字面上的空字符串。
默认 comment.char = "#"——包含 # 的行会被截断。使用 comment.char = "" 禁用（也更快）。
header 自动检测：如果第一行的字段数比后续行少一个，则设为 TRUE（缺失的字段被假定为行名）。
colClasses = "NULL" 会完全跳过该列——对提升速度非常有用。
read.csv 的默认值与 read.table 不同：header = TRUE、sep = ","、fill = TRUE、comment.char = ""。
对于大文件：指定 colClasses 和 nrows 可大幅减少内存使用。read.table 对宽数据框（数百列）很慢；对于矩阵，使用 scan 或 data.table::fread。
自 R 4.0.0 起 stringsAsFactors = FALSE（之前为 TRUE）。

write.table（陷阱）

默认 row.names = TRUE——会生成一个未命名的第一列，导致重新读取时混淆。使用 row.names = FALSE 或 col.names = NA 以生成 Excel 兼容的 CSV。
write.csv 固定 sep = ","、dec = "."，并使用 qmethod = "double"——无法通过 ... 覆盖这些设置。
quote = TRUE（默认值）会为字符/因子列加引号。数值列从不加引号。
数据框中的矩阵式列会静默扩展为多列。
对于列数很多（数百列以上）的数据框很慢；每一列按类别分别处理。

read.fwf

读取固定宽度格式文件。widths 是字段宽度向量。
负宽度会跳过相应数量的字符（用于忽略字段很有用）。
buffersize 控制一次读取多少行；处理大文件时增大该值。
在拆分字段后，内部使用 read.table。

count.fields

统计文件中每行的字段数——用于诊断读取错误很有用。
sep 和 quote 参数与 read.table 的对应参数一致。

grep / grepl / sub / gsub（易错点）

三种正则模式：POSIX 扩展（默认）、perl = TRUE、fixed = TRUE。它们在边界情况下的行为不同。
显式命名参数——x/pattern 之后未命名的参数会按位置匹配到 ignore.case、perl 等。这是静默 bug 的常见来源。
sub 只替换第一个匹配；gsub 替换所有匹配。
反向引用：替换中用 "\\1"（R 字符串中需双反斜杠）。使用 perl = TRUE 时："\\U\\1" 用于转换为大写。
grep(value = TRUE) 返回匹配的元素；grep(value = FALSE)（默认）返回索引。
grepl 返回逻辑向量——过滤时首选。
regexpr 返回第一个匹配的位置 + 长度（作为属性）；gregexpr 以列表形式返回所有匹配。
regexec 返回匹配 + 捕获组位置；gregexec 对所有匹配执行此操作。
像 [:alpha:] 这样的字符类在 POSIX 模式下必须放在 [[:alpha:]] 内（双层方括号）。

strsplit

返回一个列表（每个输入字符串对应一个元素），即使只有一个字符串也是如此。
split = "" 或 split = character(0) 会拆分为单个字符。
在字符串开头匹配：结果的第一个元素是 ""。在末尾匹配：没有尾随的 ""。
fixed = TRUE 更快，且避免正则解释。
常见错误：未命名的参数会静默匹配到 fixed、perl 等。

substr / substring

substr(x, start, stop)：提取/替换子字符串。从 1 开始索引，两端均包含。
substring(x, first, last)：相同，但 last 默认为 1000000L（实际上相当于“到末尾”）。对 first/last 进行向量化。
赋值形式：substr(x, 1, 3) <- "abc" 就地替换（替换内容必须长度相同）。

trimws

which = "both"（默认）、"left" 或 "right"。
whitespace = "[ \\t\\r\\n]" —— 可自定义正则表达式，用于定义什么算作空白字符。

nchar

type = "bytes" 统计字节数；type = "chars"（默认）统计字符数；type = "width" 统计显示宽度。
nchar(NA) 返回 NA（不是 2）。nchar(factor) 作用于水平标签。
keepNA = TRUE（自 R 3.3.0 起的默认值）；设为 FALSE 可将 "NA" 计为 2 个字符。

format / formatC

format(x, digits, nsmall)：nsmall 强制指定最少小数位数。big.mark = "," 添加千位分隔符。
formatC(x, format = "f", digits = 2)：C 风格格式化。format = "e" 表示科学计数法，"g" 表示通用格式。
format 返回字符向量；默认始终右对齐（justify = "right"）。

type.convert

将字符向量转换为合适的类型（logical、integer、double、complex、character）。
as.is = TRUE（推荐）：保持字符为字符，而不是因子。
在数据框上按列应用。tryLogical = TRUE（R 4.3+）转换 "TRUE"/"FALSE" 列。

Rscript

commandArgs(trailingOnly = TRUE) 获取脚本参数（不包括 R/Rscript 标志）。
Unix 上的 #! 行：/usr/bin/env Rscript 或完整路径。
--vanilla 或 --no-init-file 跳过 .Rprofile 加载。
退出码：quit(status = 1) 表示错误退出。

capture.output

捕获来自 cat、print 或任何写入 stdout 的表达式的输出。
file = NULL（默认）返回字符向量。file = "out.txt" 直接写入文件。
type = "message" 改为捕获 stderr。

URLencode / URLdecode

URLencode(url, reserved = FALSE) 默认不编码保留字符（/、?、& 等）。
设置 reserved = TRUE 可编码 URL 组件（查询参数值）。

glob2rx

将 shell 通配符模式转换为正则表达式：glob2rx("*.csv") → "^.*\\.csv$"。
与 list.files(pattern = glob2rx("data_*.RDS")) 搭配使用很有用。 FILE:references/modeling.md

建模 — 快速参考

R 函数中不明显的行为、陷阱和棘手的默认值。只包含 Claude 尚不知道的内容。

formula

符号化模型规范的陷阱。

若要按字面意义使用算术运算符，必须使用 I()：y ~ x + I(x^2)。没有 I() 时，^ 表示交互交叉。
* = 主效应 + 交互：a*b 展开为 a + b + a:b。
(a+b+c)^2 = 所有主效应 + 所有二阶交互（不是平方）。
- 移除项：(a+b+c)^2 - a:b 只删除 a:b 交互。
/ 表示嵌套：a/b = a + b %in% a = a + a:b。
公式中的 . 表示“数据中的所有其他列”（在 terms.formula 上下文中）或“先前内容”（在 update.formula 中）。
公式对象携带一个用于变量查找的环境；as.formula("y ~ x") 使用 parent.frame()。

terms / model.matrix

model.matrix 创建设计矩阵，包括虚拟变量编码。默认对比：无序因子使用 contr.treatment，有序因子使用 contr.poly。
terms 对象属性：order（每个项的交互阶数）、intercept、factors 矩阵。
来自 model.matrix 的列名可能出人意料：例如 factorLevelName 拼接。

glm

默认 family = gaussian(link = "identity") — 不指定 family 的 glm() 会静默拟合 OLS（与 lm 相同，但更慢，且输出基于偏差）。
常见分布族：binomial(link = "logit")、poisson(link = "log")、Gamma(link = "inverse")、inverse.gaussian()。
binomial 接受响应为：0/1 向量、逻辑值、因子（第二水平 = 成功），或两列矩阵 cbind(success, failure)。
glm 中的 weights 表示先验权重（不是频数权重）—— 对于频数权重，请使用 cbind 技巧或 offset。
predict.glm(type = "response") 用于预测概率；默认 type = "link" 返回对数几率（对于逻辑回归）或对数率（对于泊松回归）。
anova(glm_obj, test = "Chisq") 用于基于偏差的检验；"F" 对于非高斯分布族无效。
拟分布族（quasibinomial、quasipoisson）允许过度离散 —— 不计算 AIC。
收敛：如果默认的 25 次迭代不够，使用 control = glm.control(maxit = 100)。

aov

aov 是 lm 的包装器，为平衡 ANOVA 存储额外信息。对于非平衡设计，计算 Type I SS（序贯型）——项的顺序很重要。
对于 Type III SS，使用 car::Anova() 或将对比设置为 contr.sum/contr.helmert。
重复测量的误差分层：aov(y ~ A*B + Error(Subject/B))。
summary.aov 给出 ANOVA 表；summary.lm(aov_obj) 给出回归风格的摘要。

nls

需要在 start = list(...) 中提供良好的起始值，否则会收敛失败。
自启动模型（SSlogis、SSasymp 等）会自动计算起始值。
算法 "port" 允许对参数设置边界（lower/upper）。
如果数据拟合得过于精确（没有残差噪声），收敛检查会失败——使用 control = list(scaleOffset = 1) 或对数据添加抖动。
weights 参数用于加权 NLS；na.action 用于缺失值处理。

step / add1

step 通过 AIC（默认）进行逐步模型选择。使用 k = log(n) 则为 BIC。
方向：direction = "both"（默认）、"forward" 或 "backward"。
add1/drop1 评估单项的添加/删除；step 迭代调用这些函数。
scope 参数定义搜索的上界/下界模型。
step 会就地修改模型对象——对于有许多候选项的大型模型可能会很慢。

predict.lm / predict.glm

predict.lm 使用 interval = "confidence" 给出均值响应的 CI；interval = "prediction" 给出新观测的 PI（更宽）。
newdata 必须包含与原始公式变量匹配的列——因子必须具有相同的水平。
predict.glm 使用 type = "response" 给出响应尺度上的预测（例如，logistic 的概率）；type = "link"（默认）给出链接尺度上的预测。
se.fit = TRUE 返回标准误；对于 predict.glm，无论 type 如何，这些标准误都在链接尺度上。
predict.lm 使用 type = "terms" 返回每一项的贡献。

loess

span 控制平滑程度（默认 0.75）。span < 1 时使用该比例的点；span > 1 时使用所有点并调整距离。
最多 4 个预测变量。内存占用大致与 n 成 二次关系（1000 个点约 10MB）。
允许 degree = 0（局部常数），但测试不充分——请谨慎使用。
与 S 的 loess 不完全相同；未实现条件化。
normalize = TRUE（默认）将预测变量标准化到统一尺度；对空间坐标设为 FALSE。

lowess vs loess

lowess 是较旧的函数；返回 list(x, y)——无法在新点上预测。
loess 是较新的公式接口，带有 predict 方法。
lowess 的参数是 f（span，默认 2/3）；loess 的参数是 span（默认 0.75）。
lowess 的 iter 默认为 3（稳健化迭代）；loess 默认 family = "gaussian"（无稳健性）。

smooth.spline

默认平滑参数由 GCV（广义交叉验证）选择。
cv = TRUE 改用普通留一交叉验证——不要在有重复 x 值时使用。
spar 和 lambda 控制平滑程度；df 可指定等价自由度。
返回带有 predict、print、plot 方法的对象。fit 分量包含节点和系数。

optim

默认 最小化。要最大化：设置 control = list(fnscale = -1)。
默认方法是 Nelder-Mead（无梯度，稳健但慢）。对 1D 效果差——请使用 "Brent" 或 optimize()。
"L-BFGS-B" 是唯一支持箱约束（lower/upper）的方法。设置边界会自动选择此方法并给出警告。
"SANN"（模拟退火）：收敛代码 始终为 0——它从不“失败”。maxit = 总函数求值次数（默认 10000），没有其他停止准则。
parscale：缩放参数，使每个参数的单位变化产生可比的目标函数变化。对混合尺度问题至关重要。
hessian = TRUE：返回 无约束 问题的数值 Hessian，即使箱约束处于活动状态。
fn 可以返回 NA/Inf（"L-BFGS-B" 除外，它始终要求有限值）。初始值必须是有限的。

optimize / uniroot

optimize：在有界区间上进行一维最小化。返回 minimum 和 objective。
uniroot：在 [lower, upper] 内寻找 f 的一个根。要求 f(lower) 和 f(upper) 符号相反。
使用 extendInt = "yes" 的 uniroot 可以自动扩展区间以寻找符号变化——但对于实际上不穿过零点的函数，可能会找到虚假的根。
nlm：牛顿型最小化器。梯度/海森矩阵作为 fn 返回值的属性（不寻常的接口）。

TukeyHSD

需要一个拟合好的 aov 对象（不是 lm）。
默认 conf.level = 0.95。返回所有两两比较的校正 p 值和置信区间。
仅对平衡或接近平衡的设计有意义；对于非常不平衡的数据可能过于宽松。

anova（用于 lm）

anova(model)：序贯（I 型）平方和——项的次序很重要。
anova(model1, model2)：比较嵌套模型的 F 检验。
如需 II 型或 III 型平方和，请使用 car::Anova()。 FILE:references/statistics.md

统计学——快速参考

R 函数中不明显的行文、陷阱和棘手的默认值。只包含 Claude 尚不知道的内容。

chisq.test

correct = TRUE（默认）仅对 2x2 表格应用 Yates 连续性校正。
simulate.p.value = TRUE：使用 B = 2000 次重复的蒙特卡洛模拟（最小 p 值约 0.0005）。模拟假设边际固定（Fisher 式抽样，而非卡方假设）。
对于拟合优度检验：传入向量，而非矩阵。p 必须求和为 1（或设置 rescale.p = TRUE）。
返回对象包含 $expected、$residuals（Pearson）和 $stdres（标准化）。

wilcox.test

对于无并列值的小样本，默认 exact = TRUE。存在并列值时，使用正态近似。
correct = TRUE 对正态近似应用连续性校正。
conf.int = TRUE 计算 Hodges-Lehmann 估计量和置信区间（不仅仅是 p 值）。
配对检验：paired = TRUE 使用符号秩检验（Wilcoxon），而非秩和检验（Mann-Whitney）。

fisher.test

对于大于 2x2 的表格，使用模拟（simulate.p.value = TRUE）或网络算法。
workspace 控制网络算法的内存；如果在大表格上遇到错误，请增大该值。
or 参数检验特定的优势比（默认 1）——仅适用于 2x2 表格。

ks.test

双样本检验或针对参考分布的单样本检验。
不能很好地处理并列值——会发出警告并使用渐近近似。
对于复合假设（参数从数据中估计），p 值是保守的（过大）。对于离散分布，请使用 dgof 或带 exact = NULL 的 ks.test。

p.adjust

方法："holm"（默认）、"BH"（Benjamini-Hochberg FDR）、"bonferroni"、"BY"、"hochberg"、"hommel"、"fdr"（BH 的别名）、"none"。
n 参数：假设的总数（如果某些 p 值被排除，可以大于 length(p)）。
处理 NA：当输入为 NA 时，调整后的 p 值为 NA。

pairwise.t.test / pairwise.wilcox.test

p.adjust.method 默认为 "holm"。改为 "BH" 以控制 FDR。
pool.sd = TRUE（t 检验的默认值）：使用所有组的合并标准差（假设方差相等）。
返回 p 值矩阵，而非检验统计量。

shapiro.test

样本量必须在 3 到 5000 之间。
检验正态性；p 值低 = 反对正态性的证据。

kmeans

建议 nstart > 1（例如 nstart = 25）：从多个随机起点运行算法，返回最佳结果。
默认 iter.max = 10——对于收敛可能太低。对于大型/复杂数据请增大该值。
默认算法是 "Hartigan-Wong"（通常最佳）。非常接近的点可能导致不收敛（警告并带有 ifault = 4）。
聚类编号是任意的；不同平台上的顺序可能不同。
当指定 k 时，总是返回 k 个聚类（除了 Lloyd-Forgy 可能返回更少）。

hclust

method = "ward.D2" 正确实现 Ward 准则（使用平方距离）。较旧的 "ward.D" 未对距离取平方（为向后兼容而保留）。
输入必须是 dist 对象。使用 as.dist() 转换对称矩阵。
plot() 中的 hang = -1 将所有标签对齐到底部。

dist

method = "euclidean"（默认）。其他选项："manhattan"、"maximum"、"canberra"、"binary"、"minkowski"。
返回一个 dist 对象（仅下三角）。使用 as.matrix() 获取完整矩阵。
"canberra"：分子和分母均为零的项会从求和中省略（不视为 0/0）。
Inf 值：涉及 Inf 的欧几里得距离为 Inf。同一观测中有多个 Inf 时，某些方法会得到 NaN。

prcomp 与 princomp

prcomp 使用 SVD（数值上更优）；princomp 对协方差使用 eigen（稳定性较差，N-1 与 N 的缩放差异）。
prcomp 中的 scale. = TRUE 会对变量进行标准化；当变量尺度差异很大时很重要。
princomp 的标准差与 prcomp 相差因子 sqrt((n-1)/n)。
两者都返回 $rotation（载荷）和 $x（得分）；不同运行之间成分的符号可能不同。

density

默认带宽：bw = "nrd0"（Silverman 经验法则）。对于多峰数据，考虑 "SJ" 或 "bcv"。
adjust：带宽的乘数因子。adjust = 0.5 将带宽减半（更不平滑）。
默认核："gaussian"。密度的范围会超出数据范围（由 cut 控制，默认为 3 个带宽）。
n = 512：评估点的数量。增大可获得更平滑的绘图。
from/to：显式限定评估范围。

quantile

九种 type 选项（1-9）。默认 type = 7（R 默认，线性插值）。Type 1 = 经验 CDF 的逆（SAS 默认）。Type 4-9 为连续型；1-3 为不连续型。
默认 na.rm = FALSE —— 若存在任何 NA 则返回 NA。
默认 names = TRUE，将 "0%"、"25%" 等添加为名称。

分布（所有分布的常见陷阱）

所有分布函数都遵循 d/p/q/r 模式。常见的非显而易见之处：

r*() 函数中的 n 参数：如果 length(n) > 1，则使用 length(n) 作为数量，而不是 n 本身。因此 rnorm(c(1,2,3)) 生成 3 个值，而不是 1+2+3。
log = TRUE / log.p = TRUE：在对数尺度上计算，以便在尾部保持数值稳定性。
lower.tail = FALSE 直接给出生存函数 P(X > x)（在尾部比 1 - pnorm() 更精确）。
Gamma：由 shape 和 rate（= 1/scale）参数化。默认 rate = 1。同时指定 rate 和 scale 会报错。
Beta：shape1（alpha）、shape2（beta）——没有 mean/sd 参数化形式。
Poisson dpois：x 可以是非整数（当 log = FALSE 时，对非整数值返回 0 并给出警告）。
Weibull：shape 和 scale（没有 rate）。R 的参数化形式：f(x) = (shape/scale)(x/scale)^(shape-1) exp(-(x/scale)^shape)。
Lognormal：meanlog 和 sdlog 是对数的均值/标准差，而不是分布本身的。

cor.test

默认方法："pearson"。还有 "kendall" 和 "spearman"。
返回 $estimate、$p.value、$conf.int（置信区间仅适用于 Pearson）。
公式接口：cor.test(~ x + y, data = df)——注意 ~ 左侧没有内容。

ecdf

返回一个函数（阶梯函数）。在新值上调用它：Fn <- ecdf(x); Fn(3.5)。
plot(ecdf(x)) 给出经验 CDF 图。
返回的函数是右连续且具有左极限的（cadlag）。

weighted.mean

处理权重中的 NA：如果权重为 NA，则丢弃该观测。
权重不需要求和为 1；它们会在内部被归一化。 FILE:references/visualization.md

可视化——快速参考

R 函数的非显而易见行为、陷阱和棘手的默认值。只包含 Claude 尚不知道的内容。

par（易错点）

par() 的设置是按设备生效的。打开新设备会重置一切。
设置 mfrow/mfcol 会把 cex 重置为 1，把 mex 重置为 1。在 2x2 布局下，基础 cex 会乘以 0.83；在 3 行/列及以上时，乘以 0.66。
mai（英寸）、mar（行数）、pin、plt、pty 之间会相互影响。在设备调整大小后恢复所有已保存的参数可能产生不一致的结果——按字母顺序最后的那个生效。
通过 par() 设置 bg 也会设置 new = FALSE。通过 par() 设置 fg 也会设置 col。
xpd = NA 裁剪到设备区域（允许在外边距中绘图）；xpd = TRUE 裁剪到图形区域；xpd = FALSE（默认）裁剪到绘图区域。
mgp = c(3, 1, 0)：控制标题行（mgp[1]）、标签行（mgp[2]）、轴线（mgp[3]）。全部以 mex 为单位。
las：0 = 与轴平行，1 = 水平，2 = 与轴垂直，3 = 竖直。不响应 srt。
tck = 1 会绘制贯穿整个绘图区的网格线。tcl = -0.5（默认）给出向外的刻度线。
对数刻度下的 usr：包含坐标范围的 log10 值，而不是原始值。
只读参数：cin、cra、csi、cxy、din、page。

layout

layout(mat)，其中 mat 是一个整数矩阵，用于指定图形排列方式。
widths/heights 接受 lcm()，以便将绝对尺寸与相对尺寸混合使用。
比 mfrow/mfcol 更灵活，但一旦设置后无法查询（不同于 par("mfrow")）。
layout.show(n) 可视化布局以便调试。

axis / mtext

axis(side, at, labels)：side 1=底部，2=左侧，3=顶部，4=右侧。
轴标签之间的默认间距由 par("mgp") 控制。如果不加管理，标签可能会重叠。
mtext：line 参数将文本定位在边距行中（0 = 紧邻绘图区，正值 = 向外）。adj 控制水平位置（0-1）。
带 outer = TRUE 的 mtext 会写入外边距（由 par(oma = ...) 设置）。

curve

第一个参数可以是 x 的表达式或函数：curve(sin, 0, 2*pi) 或 curve(x^2 + 1, 0, 10)。
add = TRUE 用于叠加到现有图形上。默认 n = 101 个求值点。
默认 xname = "x"；如果你的表达式使用不同的变量名，请更改。

pairs

panel 函数为每一对接收 (x, y, ...)。lower.panel、upper.panel、diag.panel 用于不同区域。
gap 控制面板之间的间距（默认 1）。
公式接口：pairs(~ var1 + var2 + var3, data = df)。

coplot

条件图：coplot(y ~ x | a) 或 coplot(y ~ x | a * b) 用于两个条件变量。
panel 函数可以自定义；rows/columns 控制布局。
默认面板绘制点；使用 panel = panel.smooth 叠加 loess 平滑曲线。

matplot / matlines / matpoints

将一个矩阵的列与另一个矩阵的列作图。跨列循环使用 col、lty、pch。
默认 type = "l"（不同于 plot 默认的 "p"）。
适用于同时绘制多个时间序列或拟合曲线。

contour / filled.contour / image

contour(x, y, z)：z 必须是 dim = c(length(x), length(y)) 的矩阵。
filled.contour 具有非标准布局——它为颜色键创建自己的绘图区域。不能与 par(mfrow) 一起使用。添加元素需要使用 plot.axes 参数。
image：将 z 值绘制为彩色矩形。默认配色方案可能具有误导性；请显式设置 col。
对于 image，x 和 y 根据上下文指定单元格边界或中点。

persp

persp(x, y, z, theta, phi)：theta = 方位角，phi = 余纬。
返回一个变换矩阵（不可见），用于将 3D 投影到 2D——使用 trans3d() 向透视图中添加点/线。
shade 和 col 控制曲面着色。border = NA 移除网格线。

segments / arrows / rect / polygon

全部接受向量化坐标；按需循环回收。
arrows：code = 1（箭头在起点），code = 2（箭头在终点，默认），code = 3（两端都有）。
polygon：最后一个点会自动连接到第一个点。用 col 填充；border 控制轮廓。
rect(xleft, ybottom, xright, ytop) —— 注意参数顺序与其他系统不同。

dev / dev.off / dev.copy

dev.new() 打开一个新设备。dev.off() 关闭当前设备（并刷新 pdf 等文件设备的输出）。
对最后一个打开的设备调用 dev.off() 会回退到空设备。
dev.copy(pdf, file = "plot.pdf") 后接 dev.off() 以保存当前绘图。
dev.list() 返回所有打开的设备；dev.cur() 返回当前活动设备。

pdf

必须调用 dev.off() 才能最终完成文件。否则，文件可能为空或损坏。
onefile = TRUE（默认）：多个页面放在一个 PDF 中。onefile = FALSE：每页一个文件（在文件名中使用 %d 进行编号）。
建议使用 useDingbats = FALSE，以避免某些 PDF 查看器和 pch 符号出现问题。
默认尺寸：7x7 英寸。family 控制字体族。

png / bitmap devices

res 控制 DPI（默认 72）。用于出版：res = 300，并配合以像素或英寸为单位的适当 width/height（使用 units = "in"）。
type = "cairo"（在支持 cairo 的系统上）比默认设置提供更好的抗锯齿效果。
bg = "transparent" 用于透明背景（PNG 支持 alpha）。

colors / rgb / hcl / col2rgb

colors() 返回全部 657 种命名颜色。col2rgb("color") 返回 RGB 矩阵。
rgb(r, g, b, alpha, maxColorValue = 255) —— 注意 maxColorValue 默认值是 1，而不是 255。
hcl(h, c, l)：感知均匀的颜色空间。推荐用于颜色标度。
adjustcolor(col, alpha.f = 0.5)：添加透明度的简便方法。

colorRamp / colorRampPalette

colorRamp 返回一个将 [0,1] 映射到 RGB 矩阵的函数。
colorRampPalette 返回一个接受 n 并返回 n 个插值颜色的函数。
space = "Lab" 比 "rgb" 提供感知上更均匀的插值。

palette / recordPlot

palette() 返回当前调色板（默认 8 种颜色）。palette("Set1") 设置一个内置调色板。
绘图中的整数颜色会索引到调色板（带循环）。索引 0 = 背景色。
recordPlot() / replayPlot()：保存和恢复完整图形——依赖设备，跨会话时很脆弱。 FILE:assets/analysis_template.R

============================================================

分析模板 — Base R

复制此文件，重命名，并填写你的详细信息。

============================================================

作者 :

日期 :

数据 :

目的 :

============================================================

── 0. 设置 ─────────────────────────────────────────────────

清空环境（可选——如果要加载到现有会话中，请注释掉）

rm(list = ls())

如需要，设置工作目录

setwd("/path/to/your/project")

可复现性

set.seed(42)

库——取消注释你需要的部分

library(haven) # 读取 .dta / .sav / .sas

library(readxl) # 读取 Excel 文件

library(openxlsx) # 写入 Excel 文件

library(foreign) # 较旧的 Stata / SPSS 格式

library(survey) # 调查加权分析

library(lmtest) # Breusch-Pagan、Durbin-Watson 等

library(sandwich) # 稳健标准误

library(car) # Type II/III ANOVA、VIF

── 1. 加载数据 ─────────────────────────────────────────────

df <- read.csv("your_data.csv", stringsAsFactors = FALSE)

df <- readRDS("your_data.rds")

df <- haven::read_dta("your_data.dta")

初次查看 — 始终运行这些

dim(df) str(df) head(df, 10) summary(df)

── 2. 数据质量检查 ────────────────────────────────────

缺失值

na_report <- data.frame( column = names(df), n_miss = colSums(is.na(df)), pct_miss = round(colMeans(is.na(df)) * 100, 1), row.names = NULL ) print(na_report[na_report$n_miss > 0, ])

重复值

n_dup <- sum(duplicated(df)) cat(sprintf("重复行：%d\n", n_dup))

分类列的唯一值

cat_cols <- names(df)[sapply(df, function(x) is.character(x) | is.factor(x))] for (col in cat_cols) { cat(sprintf("\n%s（%d 个唯一值）：\n", col, length(unique(df[[col]])))) print(table(df[[col]], useNA = "ifany")) }

── 3. 清洗与转换 ─────────────────────────────────────

重命名列（示例）

names(df)[names(df) == "old_name"] <- "new_name"

转换类型

df$group <- as.factor(df$group)

df$date <- as.Date(df$date, format = "%Y-%m-%d")

重新编码值（示例）

df$gender <- ifelse(df$gender == 1, "Male", "Female")

创建新变量（示例）

df$log_income <- log(df$income + 1)

df$age_group <- cut(df$age,

breaks = c(0, 25, 45, 65, Inf),

labels = c("18-25", "26-45", "46-65", "65+"))

筛选行（示例）

df <- df[df$year >= 2010, ]

df <- df[complete.cases(df[, c("outcome", "predictor")]), ]

删除未使用的因子水平

df <- droplevels(df)

── 4. 描述性统计 ────────────────────────────────

数值摘要

num_cols <- names(df)[sapply(df, is.numeric)] round(sapply(df[num_cols], function(x) c( n = sum(!is.na(x)), mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE), median = median(x, na.rm = TRUE), min = min(x, na.rm = TRUE), max = max(x, na.rm = TRUE) )), 3)

交叉表

table(df$group, df$category, useNA = "ifany")

prop.table(table(df$group, df$category), margin = 1) # 行比例

── 5. 可视化（探索性数据分析） ───────────────────────────────────

par(mfrow = c(2, 2))

主要结果变量的直方图

hist(df$outcome_var, main = "Distribution of Outcome", xlab = "Outcome", col = "steelblue", border = "white", breaks = 30)

按组别绘制箱线图

boxplot(outcome_var ~ group_var, data = df, main = "Outcome by Group", col = "lightyellow", las = 2)

散点图

plot(df$predictor, df$outcome_var, main = "Predictor vs Outcome", xlab = "Predictor", ylab = "Outcome", pch = 19, col = adjustcolor("steelblue", alpha.f = 0.5), cex = 0.8) abline(lm(outcome_var ~ predictor, data = df), col = "red", lwd = 2)

── 6. 分析 ───────────────────────────────────────────────

·· 6a. 均值比较 ··

t.test(outcome_var ~ group_var, data = df)

·· 6b. 线性回归 ··

fit <- lm(outcome_var ~ predictor1 + predictor2 + group_var, data = df) summary(fit) confint(fit)

检查多重共线性的 VIF（需要 car）

car::vif(fit)

稳健标准误（需要 lmtest + sandwich）

lmtest::coeftest(fit, vcov = sandwich::vcovHC(fit, type = "HC3"))

·· 6c. 方差分析 ··

fit_aov <- aov(outcome_var ~ group_var, data = df)

summary(fit_aov)

TukeyHSD(fit_aov)

·· 6d. 逻辑回归（二分类结果变量） ··

fit_logit <- glm(binary_outcome ~ x1 + x2,

data = df,

family = binomial(link = "logit"))

summary(fit_logit)

exp(coef(fit_logit)) # 优势比

exp(confint(fit_logit)) # 优势比置信区间

── 7. 模型诊断 ─────────────────────────────────────

par(mfrow = c(2, 2)) plot(fit) par(mfrow = c(1, 1))

残差正态性

shapiro.test(residuals(fit))

同方差性（需要 lmtest）

lmtest::bptest(fit)

── 8. 保存输出 ────────────────────────────────────────────

清洗后的数据

write.csv(df, "data_clean.csv", row.names = FALSE)

saveRDS(df, "data_clean.rds")

模型结果保存到文本文件

sink("results.txt")

cat("=== Linear Model ===\n")

print(summary(fit))

cat("\n=== Confidence Intervals ===\n")

print(confint(fit))

sink()

图形保存到文件

png("figure1_distributions.png", width = 1200, height = 900, res = 150)

par(mfrow = c(2, 2))

# ... 你的图形 ...

par(mfrow = c(1, 1))

dev.off()

============================================================

模板结束

============================================================

FILE:scripts/check_data.R

check_data.R — 任意 R 数据框的快速数据质量报告

用法：source("check_data.R") 然后调用 check_data(df)

或： source("check_data.R"); check_data(read.csv("yourfile.csv"))

check_data <- function(df, top_n_levels = 8) {

if (!is.data.frame(df)) stop("Input must be a data frame.")

n_row <- nrow(df) n_col <- ncol(df)

cat("══════════════════════════════════════════\n") cat(" DATA QUALITY REPORT\n") cat("══════════════════════════════════════════\n") cat(sprintf(" Rows: %d Columns: %d\n", n_row, n_col)) cat("══════════════════════════════════════════\n\n")

── 1. 列概览 ──────────────────────

cat("── COLUMN OVERVIEW ────────────────────────\n")

for (col in names(df)) { x <- df[[col]] cls <- class(x)[1] n_na <- sum(is.na(x)) pct <- round(n_na / n_row * 100, 1) n_uniq <- length(unique(x[!is.na(x)])) na_flag <- if (n_na == 0) "" else sprintf(" *** %d 个 NA (%.1f%%)", n_na, pct) cat(sprintf(" %-20s %-12s %d 个唯一值%s\n", col, cls, n_uniq, na_flag)) }

── 2. NA 摘要 ────────────────────────────

cat("\n── NA 摘要 ─────────────────────────────\n")

na_counts <- sapply(df, function(x) sum(is.na(x))) cols_with_na <- na_counts[na_counts > 0]

if (length(cols_with_na) == 0) { cat(" 没有缺失值。 \n") } else { cat(sprintf(" 含 NA 的列：%d / %d\n\n", length(cols_with_na), n_col)) for (col in names(cols_with_na)) { bar_len <- round(cols_with_na[col] / n_row * 20) bar <- paste0(rep("█", bar_len), collapse = "") pct_na <- round(cols_with_na[col] / n_row * 100, 1) cat(sprintf(" %-20s [%-20s] %d (%.1f%%)\n", col, bar, cols_with_na[col], pct_na)) } }

── 3. 数值列 ───────────────────────

num_cols <- names(df)[sapply(df, is.numeric)]

if (length(num_cols) > 0) { cat("\n── 数值列 ────────────────────────\n") cat(sprintf(" %-20s %8s %8s %8s %8s %8s\n", "列", "最小值", "均值", "中位数", "最大值", "标准差")) cat(sprintf(" %-20s %8s %8s %8s %8s %8s\n", "──────", "───", "────", "──────", "───", "──"))

for (col in num_cols) {
  x  <- df[[col]][!is.na(df[[col]])]
  if (length(x) == 0) next
  cat(sprintf("  %-20s  %8.3g  %8.3g  %8.3g  %8.3g  %8.3g\n",
              col,
              min(x), mean(x), median(x), max(x), sd(x)))
}

}

── 4. 因子 / 字符列 ───────────

cat_cols <- names(df)[sapply(df, function(x) is.factor(x) | is.character(x))]

if (length(cat_cols) > 0) { cat("\n── 分类列 ────────────────────\n") for (col in cat_cols) { x <- df[[col]] tbl <- sort(table(x, useNA = "no"), decreasing = TRUE) n_lv <- length(tbl) cat(sprintf("\n %s (%d unique values)\n", col, n_lv))

  show <- min(top_n_levels, n_lv)
  for (i in seq_len(show)) {
    lbl <- names(tbl)[i]
    cnt <- tbl[i]
    pct <- round(cnt / n_row * 100, 1)
    cat(sprintf("    %-25s  %5d  (%.1f%%)\n", lbl, cnt, pct))
  }
  if (n_lv > top_n_levels) {
    cat(sprintf("    ... and %d more levels\n", n_lv - top_n_levels))
  }
}

}

── 5. 重复行 ────────────────────────

cat("\n── 重复项 ─────────────────────────────\n") n_dup <- sum(duplicated(df)) if (n_dup == 0) { cat(" 没有重复行。\n") } else { cat(sprintf(" 发现 %d 个重复行（占数据的 %.1f%%）\n", n_dup, n_dup / n_row * 100)) }

cat("\n══════════════════════════════════════════\n") cat(" 报告结束\n") cat("══════════════════════════════════════════\n")

以不可见方式返回，供程序化使用

invisible(list( dims = c(rows = n_row, cols = n_col), na_counts = na_counts, n_dupes = n_dup )) } FILE:scripts/scaffold_analysis.R #!/usr/bin/env Rscript

scaffold_analysis.R — 生成一个入门分析脚本

用法（在终端中）：

Rscript scaffold_analysis.R myproject

Rscript scaffold_analysis.R myproject outcome_var group_var

用法（在 R 控制台中）：

source("scaffold_analysis.R")

scaffold_analysis("myproject", outcome = "score", group = "treatment")

输出：myproject_analysis.R（可直接编辑）

scaffold_analysis <- function(project_name, outcome = "outcome", group = "group", data_file = NULL) { if (is.null(data_file)) data_file <- paste0(project_name, ".csv") out_file <- paste0(project_name, "_analysis.R")

template <- sprintf( '# ============================================================

项目 : %s

创建 : %s

============================================================

── 0. 库 ─────────────────────────────────────────────

在此添加你需要的包

library(ggplot2)

library(haven) # 用于 .dta 文件

library(openxlsx) # 用于 Excel 输出

── 1. 加载数据 ─────────────────────────────────────────────

df <- read.csv("%s", stringsAsFactors = FALSE)

快速检查——务必先做这一步

cat("Dimensions:", dim(df), "\n") str(df) head(df)

── 2. 探索 / EDA ─────────────────────────────────────────

summary(df)

NA 检查

na_counts <- colSums(is.na(df)) na_counts[na_counts > 0]

关键变量分布

hist(df$%s, main = "Distribution of %s", xlab = "%s")

if ("%s" %%in%% names(df)) { table(df$%s) barplot(table(df$%s), main = "Counts by %s", col = "steelblue", las = 2) }

── 3. 清洗 / 转换 ──────────────────────────────────────

df <- df[complete.cases(df), ] # 删除含任何 NA 的行

df$%s <- as.factor(df$%s) # 转换为因子

── 4. 分析 ───────────────────────────────────────────────

按组描述性统计

tapply(df$%s, df$%s, mean, na.rm = TRUE) tapply(df$%s, df$%s, sd, na.rm = TRUE)

t 检验（两组）

t.test(%s ~ %s, data = df)

线性模型

fit <- lm(%s ~ %s, data = df) summary(fit) confint(fit)

方差分析（多组）

fit_aov <- aov(%s ~ %s, data = df)

summary(fit_aov)

TukeyHSD(fit_aov)

── 5. 可视化结果 ──────────────────────────────────────

par(mfrow = c(1, 2))

按组绘制箱线图

boxplot(%s ~ %s, data = df, main = "%s by %s", xlab = "%s", ylab = "%s", col = "lightyellow")

模型诊断

plot(fit, which = 1) # 残差 vs 拟合值

par(mfrow = c(1, 1))

── 6. 保存输出 ────────────────────────────────────────────

保存清洗后的数据

write.csv(df, "%s_clean.csv", row.names = FALSE)

将模型摘要保存为文本

sink("%s_results.txt")

summary(fit)

sink()

将图形保存到文件

png("%s_boxplot.png", width = 800, height = 600, res = 150)

boxplot(%s ~ %s, data = df, col = "lightyellow")

dev.off()

', project_name, format(Sys.Date(), "%%Y-%%m-%%d"), data_file, # 第 2 节 — 探索性数据分析 outcome, outcome, outcome, group, group, group, group, # 第 3 节 group, group, # 第 4 节 outcome, group, outcome, group, outcome, group, outcome, group, outcome, group, outcome, group, # 第 5 节 outcome, group, outcome, group, group, outcome, # 第 6 节 project_name, project_name, project_name, outcome, group )

writeLines(template, out_file) cat(sprintf("Created: %s\n", out_file)) invisible(out_file) }

── 从命令行运行 ─────────────────────────────────────

if (!interactive()) { args <- commandArgs(trailingOnly = TRUE)

if (length(args) == 0) { cat("Usage: Rscript scaffold_analysis.R [outcome_var] [group_var]\n") cat("Example: Rscript scaffold_analysis.R myproject score treatment\n") quit(status = 1) }

project <- args[1] outcome <- if (length(args) >= 2) args[2] else "outcome" group <- if (length(args) >= 3) args[3] else "group" scaffold_analysis(project, outcome = outcome, group = group) } FILE:README.md

base-r-skill

GitHub: https://github.com/iremaydas/base-r-skill

一个用于 base R 编程的 Claude Code 技能。

背景故事

我是一名政治学博士生，经常使用 R，但绝不会自称 R 圈内人。我需要一个用于 base R 的 Claude Code 技能——不要 tidyverse，不要 ggplot2，只要纯正的 R——可我到处都找不到。

于是我自己做了一个。在晚上 11 点。让 Claude 帮我为 Claude 构建一个技能。

如果你也是那种每次都要 Google how to drop NA rows in R 的人，那这个就是为你准备的。🫶

里面有什么

base-r/
├── SKILL.md                    # Main skill file
├── references/                 # Gotchas & non-obvious behaviors
│   ├── data-wrangling.md       # Subsetting traps, apply family, merge, factor quirks
│   ├── modeling.md             # Formula syntax, lm/glm/aov/nls, optim
│   ├── statistics.md           # Hypothesis tests, distributions, clustering
│   ├── visualization.md        # par, layout, devices, colors
│   ├── io-and-text.md          # read.table, grep, regex, format
│   ├── dates-and-system.md     # Date/POSIXct traps, options(), file ops
│   └── misc-utilities.md       # tryCatch, do.call, time series, utilities
├── scripts/
│   ├── check_data.R            # Quick data quality report for any data frame
│   └── scaffold_analysis.R     # Generates a starter analysis script
└── assets/
    └── analysis_template.R     # Copy-paste analysis template

参考文件是从官方 R 4.5.3 手册中精简而来的——19,518 行 → 945 行（减少 95%）。只保留了那些不那么显而易见的内容：坑点、出人意料的默认行为、棘手的交互。Claude 已经熟知的内容被删掉了。

如何使用

把这个技能添加到你的 Claude Code 配置中，指向这个仓库即可。之后当你在处理 R 任务时，Claude 会自动加载相关的参考文件。

最适合用于：

基础 R 数据操作（不用 tidyverse）
使用 lm、glm、aov 进行统计建模
使用 plot、par、barplot 进行基础绘图
理解你的 R 代码为什么会做出那种奇怪的行为

不适用于：tidyverse、ggplot2、Shiny 或 R 包开发。

`check_data.R` 脚本

这大概是这里最有用的独立工具。source 它，然后对任意数据框运行 check_data(df)，就能得到一份格式化报告，包含维度、NA 计数、数值摘要和分类明细。

source("scripts/check_data.R")
check_data(your_df)

在以下帮助下构建

Claude（显然）
官方 R 手册（全部 19,518 行）
轻微的挫败感和几杯咖啡

贡献

如果你发现遗漏的坑、错误的默认值，或者某些应该出现在参考资料里的内容——非常欢迎提交 PR。我也在学习。

由 @iremaydas 制作——博士生，偶尔用 R，全职谷歌搜索那些我大概早该知道的东西。

English original

Base R Programming Skill

A comprehensive reference for base R programming — covering data structures, control flow, functions, I/O, statistical computing, and plotting.

Quick Reference

Data Structures

# Vectors (atomic)
x <- c(1, 2, 3)              # numeric
y <- c("a", "b", "c")        # character
z <- c(TRUE, FALSE, TRUE)    # logical

# Factor
f <- factor(c("low", "med", "high"), levels = c("low", "med", "high"), ordered = TRUE)

# Matrix
m <- matrix(1:6, nrow = 2, ncol = 3)
m[1, ]       # first row
m[, 2]       # second column

# List
lst <- list(name = "ali", scores = c(90, 85), passed = TRUE)
lst$name      # access by name
lst[[2]]      # access by position

# Data frame
df <- data.frame(
  id = 1:3,
  name = c("a", "b", "c"),
  value = c(10.5, 20.3, 30.1),
  stringsAsFactors = FALSE
)
df[df$value > 15, ]    # filter rows
df$new_col <- df$value * 2  # add column

Subsetting

# Vectors
x[1:3]             # by position
x[c(TRUE, FALSE)]  # by logical
x[x > 5]           # by condition
x[-1]              # exclude first

# Data frames
df[1:5, ]                    # first 5 rows
df[, c("name", "value")]     # select columns
df[df$value > 10, "name"]    # filter + select
subset(df, value > 10, select = c(name, value))

# which() for index positions
idx <- which(df$value == max(df$value))

Control Flow

# if/else
if (x > 0) {
  "positive"
} else if (x == 0) {
  "zero"
} else {
  "negative"
}

# ifelse (vectorized)
ifelse(x > 0, "pos", "neg")

# for loop
for (i in seq_along(x)) {
  cat(i, x[i], "\n")
}

# while
while (condition) {
  # body
  if (stop_cond) break
}

# switch
switch(type,
  "a" = do_a(),
  "b" = do_b(),
  stop("Unknown type")
)

Functions

# Define
my_func <- function(x, y = 1, ...) {
  result <- x + y
  return(result)  # or just: result
}

# Anonymous functions
sapply(1:5, function(x) x^2)
# R 4.1+ shorthand:
sapply(1:5, \(x) x^2)

# Useful: do.call for calling with a list of args
do.call(paste, list("a", "b", sep = "-"))

Apply Family

# sapply — simplify result to vector/matrix
sapply(lst, length)

# lapply — always returns list
lapply(lst, function(x) x[1])

# vapply — like sapply but with type safety
vapply(lst, length, integer(1))

# apply — over matrix margins (1=rows, 2=cols)
apply(m, 2, sum)

# tapply — apply by groups
tapply(df$value, df$group, mean)

# mapply — multivariate
mapply(function(x, y) x + y, 1:3, 4:6)

# aggregate — like tapply for data frames
aggregate(value ~ group, data = df, FUN = mean)

String Operations

paste("a", "b", sep = "-")    # "a-b"
paste0("x", 1:3)              # "x1" "x2" "x3"
sprintf("%.2f%%", 3.14159)    # "3.14%"
nchar("hello")                # 5
substr("hello", 1, 3)         # "hel"
gsub("old", "new", text)      # replace all
grep("pattern", x)            # indices of matches
grepl("pattern", x)           # logical vector
strsplit("a,b,c", ",")        # list("a","b","c")
trimws("  hi  ")              # "hi"
tolower("ABC")                # "abc"

Data I/O

# CSV
df <- read.csv("data.csv", stringsAsFactors = FALSE)
write.csv(df, "output.csv", row.names = FALSE)

# Tab-delimited
df <- read.delim("data.tsv")

# General
df <- read.table("data.txt", header = TRUE, sep = "\t")

# RDS (single R object, preserves types)
saveRDS(obj, "data.rds")
obj <- readRDS("data.rds")

# RData (multiple objects)
save(df1, df2, file = "data.RData")
load("data.RData")

# Connections
con <- file("big.csv", "r")
chunk <- readLines(con, n = 100)
close(con)

Base Plotting

# Scatter
plot(x, y, main = "Title", xlab = "X", ylab = "Y",
     pch = 19, col = "steelblue", cex = 1.2)

# Line
plot(x, y, type = "l", lwd = 2, col = "red")
lines(x, y2, col = "blue", lty = 2)  # add line

# Bar
barplot(table(df$category), main = "Counts",
        col = "lightblue", las = 2)

# Histogram
hist(x, breaks = 30, col = "grey80",
     main = "Distribution", xlab = "Value")

# Box plot
boxplot(value ~ group, data = df,
        col = "lightyellow", main = "By Group")

# Multiple plots
par(mfrow = c(2, 2))  # 2x2 grid
# ... four plots ...
par(mfrow = c(1, 1))  # reset

# Save to file
png("plot.png", width = 800, height = 600)
plot(x, y)
dev.off()

# Add elements
legend("topright", legend = c("A", "B"),
       col = c("red", "blue"), lty = 1)
abline(h = 0, lty = 2, col = "grey")
text(x, y, labels = names, pos = 3, cex = 0.8)

Statistics

# Descriptive
mean(x); median(x); sd(x); var(x)
quantile(x, probs = c(0.25, 0.5, 0.75))
summary(df)
cor(x, y)
table(df$category)  # frequency table

# Linear model
fit <- lm(y ~ x1 + x2, data = df)
summary(fit)
coef(fit)
predict(fit, newdata = new_df)
confint(fit)

# t-test
t.test(x, y)                    # two-sample
t.test(x, mu = 0)               # one-sample
t.test(before, after, paired = TRUE)

# Chi-square
chisq.test(table(df$a, df$b))

# ANOVA
fit <- aov(value ~ group, data = df)
summary(fit)
TukeyHSD(fit)

# Correlation test
cor.test(x, y, method = "pearson")

Data Manipulation

# Merge (join)
merged <- merge(df1, df2, by = "id")                  # inner
merged <- merge(df1, df2, by = "id", all = TRUE)      # full outer
merged <- merge(df1, df2, by = "id", all.x = TRUE)    # left

# Reshape
wide <- reshape(long, direction = "wide",
                idvar = "id", timevar = "time", v.names = "value")
long <- reshape(wide, direction = "long",
                varying = list(c("v1", "v2")), v.names = "value")

# Sort
df[order(df$value), ]              # ascending
df[order(-df$value), ]             # descending
df[order(df$group, -df$value), ]   # multi-column

# Remove duplicates
df[!duplicated(df), ]
df[!duplicated(df$id), ]

# Stack / combine
rbind(df1, df2)    # stack rows (same columns)
cbind(df1, df2)    # bind columns (same rows)

# Transform columns
df$log_val <- log(df$value)
df$category <- cut(df$value, breaks = c(0, 10, 20, Inf),
                   labels = c("low", "med", "high"))

Environment & Debugging

ls()                  # list objects
rm(x)                 # remove object
rm(list = ls())       # clear all
str(obj)              # structure
class(obj)            # class
typeof(obj)           # internal type
is.na(x)              # check NA
complete.cases(df)    # rows without NA
traceback()           # after error
debug(my_func)        # step through
browser()             # breakpoint in code
system.time(expr)     # timing
Sys.time()            # current time

Reference Files

For deeper coverage, read the reference files in references/:

Function Gotchas & Quick Reference (condensed from R 4.5.3 Reference Manual)

Non-obvious behaviors, surprising defaults, and tricky interactions — only what Claude doesn't already know:

data-wrangling.md — Read when: subsetting returns wrong type, apply on data frame gives unexpected coercion, merge/split/cbind behaves oddly, factor levels persist after filtering, table/duplicated edge cases.
modeling.md — Read when: formula syntax is confusing (I(), * vs :, /), aov gives wrong SS type, glm silently fits OLS, nls won't converge, predict returns wrong scale, optim/optimize needs tuning.
statistics.md — Read when: hypothesis test gives surprising result, need to choose correct p.adjust method, clustering parameters seem wrong, distribution function naming is confusing (d/p/q/r prefixes).
visualization.md — Read when: par settings reset unexpectedly, layout/mfrow interaction is confusing, axis labels are clipped, colors don't look right, need specialty plots (contour, persp, mosaic, pairs).
io-and-text.md — Read when: read.table silently drops data or misparses columns, regex behaves differently than expected, sprintf formatting is tricky, write.table output has unwanted row names.
dates-and-system.md — Read when: Date/POSIXct conversion gives wrong day, time zones cause off-by-one, difftime units are unexpected, need to find/list/test files programmatically.
misc-utilities.md — Read when: do.call behaves differently than direct call, need Reduce/Filter/Map, tryCatch handler doesn't fire, all.equal returns string not logical, time series functions need setup.

Tips for Writing Good R Code

Use vapply() over sapply() in production code — it enforces return types
Prefer seq_along(x) over 1:length(x) — the latter breaks when x is empty
Use stringsAsFactors = FALSE in read.csv() / data.frame() (default changed in R 4.0)
Vectorize operations instead of writing loops when possible
Use stop(), warning(), message() for error handling — not print()
<<- assigns to parent environment — use sparingly and intentionally
with(df, expr) avoids repeating df$ everywhere
Sys.setenv() and .Renviron for environment variables FILE:references/misc-utilities.md

Miscellaneous Utilities — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

do.call

do.call(fun, args_list) — args must be a list, even for a single argument.
quote = TRUE prevents evaluation of arguments before the call — needed when passing expressions/symbols.
Behavior of substitute inside do.call differs from direct calls. Semantics are not fully defined for this case.
Useful pattern: do.call(rbind, list_of_dfs) to combine a list of data frames.

Reduce / Filter / Map / Find / Position

R's functional programming helpers from base — genuinely non-obvious.

Reduce(f, x) applies binary function f cumulatively: Reduce("+", 1:4) = ((1+2)+3)+4. Direction matters for non-commutative ops.
Reduce(f, x, accumulate = TRUE) returns all intermediate results — equivalent to Python's itertools.accumulate.
Reduce(f, x, right = TRUE) folds from the right: f(x1, f(x2, f(x3, x4))).
Reduce with init adds a starting value: Reduce(f, x, init = v) = f(f(f(v, x1), x2), x3).
Filter(f, x) keeps elements where f(elem) is TRUE. Unlike x[sapply(x, f)], handles NULL/empty correctly.
Map(f, ...) is a simple wrapper for mapply(f, ..., SIMPLIFY = FALSE) — always returns a list.
Find(f, x) returns the first element where f(elem) is TRUE. Find(f, x, right = TRUE) for last.
Position(f, x) returns the index of the first match (like Find but returns position, not value).

lengths

lengths(x) returns the length of each element of a list. Equivalent to sapply(x, length) but faster (implemented in C).
Works on any list-like object. Returns integer vector.

conditions (tryCatch / withCallingHandlers)

tryCatch unwinds the call stack — handler runs in the calling environment, not where the error occurred. Cannot resume execution.
withCallingHandlers does NOT unwind — handler runs where the condition was signaled. Can inspect/log then let the condition propagate.
tryCatch(expr, error = function(e) e) returns the error condition object.
tryCatch(expr, warning = function(w) {...}) catches the first warning and exits. Use withCallingHandlers + invokeRestart("muffleWarning") to suppress warnings but continue.
tryCatch finally clause always runs (like Java try/finally).
globalCallingHandlers() registers handlers that persist for the session (useful for logging).
Custom conditions: stop(errorCondition("msg", class = "myError")) then catch with tryCatch(..., myError = function(e) ...).

all.equal

Tests near equality with tolerance (default 1.5e-8, i.e., sqrt(.Machine$double.eps)).
Returns TRUE or a character string describing the difference — NOT FALSE. Use isTRUE(all.equal(x, y)) in conditionals.
tolerance argument controls numeric tolerance. scale for absolute vs relative comparison.
Checks attributes, names, dimensions — more thorough than ==.

combn

combn(n, m) or combn(x, m): generates all combinations of m items from x.
Returns a matrix with m rows; each column is one combination.
FUN argument applies a function to each combination: combn(5, 3, sum) returns sums of all 3-element subsets.
simplify = FALSE returns a list instead of a matrix.

modifyList

modifyList(x, val) replaces elements of list x with those in val by name.
Setting a value to NULL removes that element from the list.
Does add new names not in x — it uses x[names(val)] <- val internally, so any name in val gets added or replaced.

relist

Inverse of unlist: given a flat vector and a skeleton list, reconstructs the nested structure.
relist(flesh, skeleton) — flesh is the flat data, skeleton provides the shape.
Works with factors, matrices, and nested lists.

txtProgressBar

txtProgressBar(min, max, style = 3) — style 3 shows percentage + bar (most useful).
Update with setTxtProgressBar(pb, value). Close with close(pb).
Style 1: rotating |/-\, style 2: simple progress. Only style 3 shows percentage.

object.size

Returns an estimate of memory used by an object. Not always exact for shared references.
format(object.size(x), units = "MB") for human-readable output.
Does not count the size of environments or external pointers.

installed.packages / update.packages

installed.packages() can be slow (scans all packages). Use find.package() or requireNamespace() to check for a specific package.
update.packages(ask = FALSE) updates all packages without prompting.
lib.loc specifies which library to check/update.

vignette / demo

vignette() lists all vignettes; vignette("name", package = "pkg") opens a specific one.
demo() lists all demos; demo("topic") runs one interactively.
browseVignettes() opens vignette browser in HTML.

Time series: acf / arima / ts / stl / decompose

ts(data, start, frequency): frequency is observations per unit time (12 for monthly, 4 for quarterly).
acf default type = "correlation". Use type = "partial" for PACF. plot = FALSE to suppress auto-plotting.
arima(x, order = c(p,d,q)) for ARIMA models. seasonal = list(order = c(P,D,Q), period = S) for seasonal component.
arima handles NA values in the time series (via Kalman filter).
stl requires s.window (seasonal window) — must be specified, no default. s.window = "periodic" assumes fixed seasonality.
decompose: simpler than stl, uses moving averages. type = "additive" or "multiplicative".
stl result components: $time.series matrix with columns seasonal, trend, remainder. FILE:references/data-wrangling.md

Data Wrangling — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

Extract / Extract.data.frame

Indexing pitfalls in base R.

m[j = 2, i = 1] is m[2, 1] not m[1, 2] — argument names are ignored in [, positional matching only. Never name index args.
Factor indexing: x[f] uses integer codes of factor f, not its character labels. Use x[as.character(f)] for label-based indexing.
x[[]] with no index is always an error. x$name does partial matching by default; x[["name"]] does not (exact by default).
Assigning NULL via x[[i]] <- NULL or x$name <- NULL deletes that list element.
Data frame [ with single column: df[, 1] returns a vector (drop=TRUE default for columns), but df[1, ] returns a data frame (drop=FALSE for rows). Use drop = FALSE explicitly.
Matrix indexing a data frame (df[cbind(i,j)]) coerces to matrix first — avoid.

subset

Use interactively only; unsafe for programming.

subset argument uses non-standard evaluation — column names are resolved in the data frame, which can silently pick up wrong variables in programmatic use. Use [ with explicit logic in functions.
NAs in the logical condition are treated as FALSE (rows silently dropped).
Factors may retain unused levels after subsetting; call droplevels().

match / %in%

%in% never returns NA — this makes it safe for if() conditions unlike ==.
match() returns position of first match only; duplicates in table are ignored.
Factors, raw vectors, and lists are all converted to character before matching.
NaN matches NaN but not NA; NA matches NA only.

apply

On a data frame, apply coerces to matrix via as.matrix first — mixed types become character.
Return value orientation is transposed: if FUN returns length-n vector, result has dim c(n, dim(X)[MARGIN]). Row results become columns.
Factor results are coerced to character in the output array.
... args cannot share names with X, MARGIN, or FUN (partial matching risk).

lapply / sapply / vapply

sapply can return a vector, matrix, or list unpredictably — use vapply in non-interactive code with explicit FUN.VALUE template.
Calling primitives directly in lapply can cause dispatch issues; wrap in function(x) is.numeric(x) rather than bare is.numeric.
sapply with simplify = "array" can produce higher-rank arrays (not just matrices).

tapply

Returns an array (not a data frame). Class info on return values is discarded (e.g., Date objects become numeric).
... args to FUN are not divided into cells — they apply globally, so FUN should not expect additional args with same length as X.
default = NA fills empty cells; set default = 0 for sum-like operations. Before R 3.4.0 this was hard-coded to NA.
Use array2DF() to convert result to a data frame.

mapply

Argument name is SIMPLIFY (all caps) not simplify — inconsistent with sapply.
MoreArgs must be a list of args not vectorized over.
Recycles shorter args to common length; zero-length arg gives zero-length result.

merge

Default by is intersect(names(x), names(y)) — can silently merge on unintended columns if data frames share column names.
by = 0 or by = "row.names" merges on row names, adding a "Row.names" column.
by = NULL (or both by.x/by.y length 0) produces Cartesian product.
Result is sorted on by columns by default (sort = TRUE). For unsorted output use sort = FALSE.
Duplicate key matches produce all combinations (one row per match pair).

split

If f is a list of factors, interaction is used; levels containing "." can cause unexpected splits unless sep is changed.
drop = FALSE (default) retains empty factor levels as empty list elements.
Supports formula syntax: split(df, ~ Month).

cbind / rbind

cbind on data frames calls data.frame(...), not cbind.matrix. Mixing matrices and data frames can give unexpected results.
rbind on data frames matches columns by name, not position. Missing columns get NA.
cbind(NULL) returns NULL (not a matrix). For consistency, rbind(NULL) also returns NULL.

table

By default excludes NA (useNA = "no"). Use useNA = "ifany" or exclude = NULL to count NAs.
Setting exclude non-empty and non-default implies useNA = "ifany".
Result is always an array (even 1D), class "table". Convert to data frame with as.data.frame(tbl).
Two kinds of NA (factor-level NA vs actual NA) are treated differently depending on useNA/exclude.

duplicated / unique

duplicated marks the second and later occurrences as TRUE, not the first. Use fromLast = TRUE to reverse.
For data frames, operates on whole rows. For lists, compares recursively.
unique keeps the first occurrence of each value.

data.frame (gotchas)

stringsAsFactors = FALSE is the default since R 4.0.0 (was TRUE before).
Atomic vectors recycle to match longest column, but only if exact multiple. Protect with I() to prevent conversion.
Duplicate column names allowed only with check.names = FALSE, but many operations will de-dup them silently.
Matrix arguments are expanded to multiple columns unless protected by I().

factor (gotchas)

as.numeric(f) returns integer codes, not original values. Use as.numeric(levels(f))[f] or as.numeric(as.character(f)).
Only == and != work between factors; factors must have identical level sets. Ordered factors support <, >.
c() on factors unions level sets (since R 4.1.0), but earlier versions converted to integer.
Levels are sorted by default, but sort order is locale-dependent at creation time.

aggregate

Formula interface (aggregate(y ~ x, data, FUN)) drops NA groups by default.
The data frame method requires by as a list (not a vector).
Returns columns named after the grouping variables, with result column keeping the original name.
If FUN returns multiple values, result column is a matrix column inside the data frame.

complete.cases

Returns a logical vector: TRUE for rows with no NAs across all columns/arguments.
Works on multiple arguments (e.g., complete.cases(x, y) checks both).

order

Returns a permutation vector of indices, not the sorted values. Use x[order(x)] to sort.
Default is ascending; use -x for descending numeric, or decreasing = TRUE.
For character sorting, depends on locale. Use method = "radix" for locale-independent fast sorting.
sort.int() with method = "radix" is much faster for large integer/character vectors. FILE:references/dates-and-system.md

Dates and System — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

Dates (Date class)

Date objects are stored as integer days since 1970-01-01. Arithmetic works in days.
Sys.Date() returns current date as Date object.
seq.Date(from, to, by = "month") — "month" increments can produce varying-length intervals. Adding 1 month to Jan 31 gives Mar 3 (not Feb 28).
diff(dates) returns a difftime object in days.
format(date, "%Y") for year, "%m" for month, "%d" for day, "%A" for weekday name (locale-dependent).
Years before 1CE may not be handled correctly.
length(date_vector) <- n pads with NAs if extended.

DateTimeClasses (POSIXct / POSIXlt)

POSIXct: seconds since 1970-01-01 UTC (compact, a numeric vector).
POSIXlt: list with components $sec, $min, $hour, $mday, $mon (0-11!), $year (since 1900!), $wday (0-6, Sunday=0), $yday (0-365).
Converting between POSIXct and Date: as.Date(posixct_obj) uses tz = "UTC" by default — may give different date than intended if original was in another timezone.
Sys.time() returns POSIXct in current timezone.
strptime returns POSIXlt; as.POSIXct(strptime(...)) to get POSIXct.
difftime arithmetic: subtracting POSIXct objects gives difftime. Units auto-selected ("secs", "mins", "hours", "days", "weeks").

difftime

difftime(time1, time2, units = "auto") — auto-selects smallest sensible unit.
Explicit units: "secs", "mins", "hours", "days", "weeks". No "months" or "years" (variable length).
as.numeric(diff, units = "hours") to extract numeric value in specific units.
units(diff_obj) <- "hours" changes the unit in place.

system.time / proc.time

system.time(expr) returns user, system, and elapsed time.
gcFirst = TRUE (default): runs garbage collection before timing for more consistent results.
proc.time() returns cumulative time since R started — take differences for intervals.
elapsed (wall clock) can be less than user (multi-threaded BLAS) or more (I/O waits).

Sys.sleep

Sys.sleep(seconds) — allows fractional seconds. Actual sleep may be longer (OS scheduling).
The process yields to the OS during sleep (does not busy-wait).

options (key options)

Selected non-obvious options:

options(scipen = n): positive biases toward fixed notation, negative toward scientific. Default 0. Applies to print/format/cat but not sprintf.
options(digits = n): significant digits for printing (1-22, default 7). Suggestion only.
options(digits.secs = n): max decimal digits for seconds in time formatting (0-6, default 0).
options(warn = n): -1 = ignore warnings, 0 = collect (default), 1 = immediate, 2 = convert to errors.
options(error = recover): drop into debugger on error. options(error = NULL) resets to default.
options(OutDec = ","): change decimal separator in output (affects format, print, NOT sprintf).
options(stringsAsFactors = FALSE): global default for data.frame (moot since R 4.0.0 where it's already FALSE).
options(expressions = 5000): max nested evaluations. Increase for deep recursion.
options(max.print = 99999): controls truncation in print output.
options(na.action = "na.omit"): default NA handling in model functions.
options(contrasts = c("contr.treatment", "contr.poly")): default contrasts for unordered/ordered factors.

file.path / basename / dirname

file.path("a", "b", "c.txt") → "a/b/c.txt" (platform-appropriate separator).
basename("/a/b/c.txt") → "c.txt". dirname("/a/b/c.txt") → "/a/b".
file.path does NOT normalize paths (no .. resolution); use normalizePath() for that.

list.files

list.files(pattern = "*.csv") — pattern is a regex, not a glob! Use glob2rx("*.csv") or "\\.csv$".
full.names = FALSE (default) returns basenames only. Use full.names = TRUE for complete paths.
recursive = TRUE to search subdirectories.
all.files = TRUE to include hidden files (starting with .).

file.info

Returns data frame with size, isdir, mode, mtime, ctime, atime, uid, gid.
mtime: modification time (POSIXct). Useful for file.info(f)$mtime.
On some filesystems, ctime is status-change time, not creation time.

file_test

file_test("-f", path): TRUE if regular file exists.
file_test("-d", path): TRUE if directory exists.
file_test("-nt", f1, f2): TRUE if f1 is newer than f2.
More reliable than file.exists() for distinguishing files from directories. FILE:references/io-and-text.md

I/O and Text Processing — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

read.table (gotchas)

sep = "" (default) means any whitespace (spaces, tabs, newlines) — not a literal empty string.
comment.char = "#" by default — lines with # are truncated. Use comment.char = "" to disable (also faster).
header auto-detection: set to TRUE if first row has one fewer field than subsequent rows (the missing field is assumed to be row names).
colClasses = "NULL" skips that column entirely — very useful for speed.
read.csv defaults differ from read.table: header = TRUE, sep = ",", fill = TRUE, comment.char = "".
For large files: specifying colClasses and nrows dramatically reduces memory usage. read.table is slow for wide data frames (hundreds of columns); use scan or data.table::fread for matrices.
stringsAsFactors = FALSE since R 4.0.0 (was TRUE before).

write.table (gotchas)

row.names = TRUE by default — produces an unnamed first column that confuses re-reading. Use row.names = FALSE or col.names = NA for Excel-compatible CSV.
write.csv fixes sep = ",", dec = ".", and uses qmethod = "double" — cannot override these via ....
quote = TRUE (default) quotes character/factor columns. Numeric columns are never quoted.
Matrix-like columns in data frames expand to multiple columns silently.
Slow for data frames with many columns (hundreds+); each column processed separately by class.

read.fwf

Reads fixed-width format files. widths is a vector of field widths.
Negative widths skip that many characters (useful for ignoring fields).
buffersize controls how many lines are read at a time; increase for large files.
Uses read.table internally after splitting fields.

count.fields

Counts fields per line in a file — useful for diagnosing read errors.
sep and quote arguments match those of read.table.

grep / grepl / sub / gsub (gotchas)

Three regex modes: POSIX extended (default), perl = TRUE, fixed = TRUE. They behave differently for edge cases.
Name arguments explicitly — unnamed args after x/pattern are matched positionally to ignore.case, perl, etc. Common source of silent bugs.
sub replaces first match only; gsub replaces all matches.
Backreferences: "\\1" in replacement (double backslash in R strings). With perl = TRUE: "\\U\\1" for uppercase conversion.
grep(value = TRUE) returns matching elements; grep(value = FALSE) (default) returns indices.
grepl returns logical vector — preferred for filtering.
regexpr returns first match position + length (as attributes); gregexpr returns all matches as a list.
regexec returns match + capture group positions; gregexec does this for all matches.
Character classes like [:alpha:] must be inside [[:alpha:]] (double brackets) in POSIX mode.

strsplit

Returns a list (one element per input string), even for a single string.
split = "" or split = character(0) splits into individual characters.
Match at beginning of string: first element of result is "". Match at end: no trailing "".
fixed = TRUE is faster and avoids regex interpretation.
Common mistake: unnamed arguments silently match fixed, perl, etc.

substr / substring

substr(x, start, stop): extracts/replaces substring. 1-indexed, inclusive on both ends.
substring(x, first, last): same but last defaults to 1000000L (effectively "to end"). Vectorized over first/last.
Assignment form: substr(x, 1, 3) <- "abc" replaces in place (must be same length replacement).

trimws

which = "both" (default), "left", or "right".
whitespace = "[ \\t\\r\\n]" — customizable regex for what counts as whitespace.

nchar

type = "bytes" counts bytes; type = "chars" (default) counts characters; type = "width" counts display width.
nchar(NA) returns NA (not 2). nchar(factor) works on the level labels.
keepNA = TRUE (default since R 3.3.0); set to FALSE to count "NA" as 2 characters.

format / formatC

format(x, digits, nsmall): nsmall forces minimum decimal places. big.mark = "," adds thousands separator.
formatC(x, format = "f", digits = 2): C-style formatting. format = "e" for scientific, "g" for general.
format returns character vector; always right-justified by default (justify = "right").

type.convert

Converts character vectors to appropriate types (logical, integer, double, complex, character).
as.is = TRUE (recommended): keeps characters as character, not factor.
Applied column-wise on data frames. tryLogical = TRUE (R 4.3+) converts "TRUE"/"FALSE" columns.

Rscript

commandArgs(trailingOnly = TRUE) gets script arguments (excluding R/Rscript flags).
#! line on Unix: /usr/bin/env Rscript or full path.
--vanilla or --no-init-file to skip .Rprofile loading.
Exit code: quit(status = 1) for error exit.

capture.output

Captures output from cat, print, or any expression that writes to stdout.
file = NULL (default) returns character vector. file = "out.txt" writes directly to file.
type = "message" captures stderr instead.

URLencode / URLdecode

URLencode(url, reserved = FALSE) by default does NOT encode reserved chars (/, ?, &, etc.).
Set reserved = TRUE to encode a URL component (query parameter value).

glob2rx

Converts shell glob patterns to regex: glob2rx("*.csv") → "^.*\\.csv$".
Useful with list.files(pattern = glob2rx("data_*.RDS")). FILE:references/modeling.md

Modeling — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

formula

Symbolic model specification gotchas.

I() is required to use arithmetic operators literally: y ~ x + I(x^2). Without I(), ^ means interaction crossing.
* = main effects + interaction: a*b expands to a + b + a:b.
(a+b+c)^2 = all main effects + all 2-way interactions (not squaring).
- removes terms: (a+b+c)^2 - a:b drops only the a:b interaction.
/ means nesting: a/b = a + b %in% a = a + a:b.
. in formula means "all other columns in data" (in terms.formula context) or "previous contents" (in update.formula).
Formula objects carry an environment used for variable lookup; as.formula("y ~ x") uses parent.frame().

terms / model.matrix

model.matrix creates the design matrix including dummy coding. Default contrasts: contr.treatment for unordered factors, contr.poly for ordered.
terms object attributes: order (interaction order per term), intercept, factors matrix.
Column names from model.matrix can be surprising: e.g., factorLevelName concatenation.

glm

Default family = gaussian(link = "identity") — glm() with no family silently fits OLS (same as lm, but slower and with deviance-based output).
Common families: binomial(link = "logit"), poisson(link = "log"), Gamma(link = "inverse"), inverse.gaussian().
binomial accepts response as: 0/1 vector, logical, factor (second level = success), or 2-column matrix cbind(success, failure).
weights in glm means prior weights (not frequency weights) — for frequency weights, use the cbind trick or offset.
predict.glm(type = "response") for predicted probabilities; default type = "link" returns log-odds (for logistic) or log-rate (for Poisson).
anova(glm_obj, test = "Chisq") for deviance-based tests; "F" is invalid for non-Gaussian families.
Quasi-families (quasibinomial, quasipoisson) allow overdispersion — no AIC is computed.
Convergence: control = glm.control(maxit = 100) if default 25 iterations isn't enough.

aov

aov is a wrapper around lm that stores extra info for balanced ANOVA. For unbalanced designs, Type I SS (sequential) are computed — order of terms matters.
For Type III SS, use car::Anova() or set contrasts to contr.sum/contr.helmert.
Error strata for repeated measures: aov(y ~ A*B + Error(Subject/B)).
summary.aov gives ANOVA table; summary.lm(aov_obj) gives regression-style summary.

nls

Requires good starting values in start = list(...) or convergence fails.
Self-starting models (SSlogis, SSasymp, etc.) auto-compute starting values.
Algorithm "port" allows bounds on parameters (lower/upper).
If data fits too exactly (no residual noise), convergence check fails — use control = list(scaleOffset = 1) or jitter data.
weights argument for weighted NLS; na.action for missing value handling.

step / add1

step does stepwise model selection by AIC (default). Use k = log(n) for BIC.
Direction: direction = "both" (default), "forward", or "backward".
add1/drop1 evaluate single-term additions/deletions; step calls these iteratively.
scope argument defines the upper/lower model bounds for search.
step modifies the model object in place — can be slow for large models with many candidate terms.

predict.lm / predict.glm

predict.lm with interval = "confidence" gives CI for mean response; interval = "prediction" gives PI for new observation (wider).
newdata must have columns matching the original formula variables — factors must have the same levels.
predict.glm with type = "response" gives predictions on the response scale (e.g., probabilities for logistic); type = "link" (default) gives on the link scale.
se.fit = TRUE returns standard errors; for predict.glm these are on the link scale regardless of type.
predict.lm with type = "terms" returns the contribution of each term.

loess

span controls smoothness (default 0.75). Span < 1 uses that proportion of points; span > 1 uses all points with adjusted distance.
Maximum 4 predictors. Memory usage is roughly quadratic in n (1000 points ~ 10MB).
degree = 0 (local constant) is allowed but poorly tested — use with caution.
Not identical to S's loess; conditioning is not implemented.
normalize = TRUE (default) standardizes predictors to common scale; set FALSE for spatial coords.

lowess vs loess

lowess is the older function; returns list(x, y) — cannot predict at new points.
loess is the newer formula interface with predict method.
lowess parameter is f (span, default 2/3); loess parameter is span (default 0.75).
lowess iter default is 3 (robustifying iterations); loess default family = "gaussian" (no robustness).

smooth.spline

Default smoothing parameter selected by GCV (generalized cross-validation).
cv = TRUE uses ordinary leave-one-out CV instead — do not use with duplicate x values.
spar and lambda control smoothness; df can specify equivalent degrees of freedom.
Returns object with predict, print, plot methods. The fit component has knots and coefficients.

optim

Minimizes by default. To maximize: set control = list(fnscale = -1).
Default method is Nelder-Mead (no gradients, robust but slow). Poor for 1D — use "Brent" or optimize().
"L-BFGS-B" is the only method supporting box constraints (lower/upper). Bounds auto-select this method with a warning.
"SANN" (simulated annealing): convergence code is always 0 — it never "fails". maxit = total function evals (default 10000), no other stopping criterion.
parscale: scale parameters so unit change in each produces comparable objective change. Critical for mixed-scale problems.
hessian = TRUE: returns numerical Hessian of the unconstrained problem even if box constraints are active.
fn can return NA/Inf (except "L-BFGS-B" which requires finite values always). Initial value must be finite.

optimize / uniroot

optimize: 1D minimization on a bounded interval. Returns minimum and objective.
uniroot: finds a root of f in [lower, upper]. Requires f(lower) and f(upper) to have opposite signs.
uniroot with extendInt = "yes" can auto-extend the interval to find sign change — but can find spurious roots for functions that don't actually cross zero.
nlm: Newton-type minimizer. Gradient/Hessian as attributes of the return value from fn (unusual interface).

TukeyHSD

Requires a fitted aov object (not lm).
Default conf.level = 0.95. Returns adjusted p-values and confidence intervals for all pairwise comparisons.
Only meaningful for balanced or near-balanced designs; can be liberal for very unbalanced data.

anova (for lm)

anova(model): sequential (Type I) SS — order of terms matters.
anova(model1, model2): F-test comparing nested models.
For Type II or III SS use car::Anova(). FILE:references/statistics.md

Statistics — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

chisq.test

correct = TRUE (default) applies Yates continuity correction for 2x2 tables only.
simulate.p.value = TRUE: Monte Carlo with B = 2000 replicates (min p ~ 0.0005). Simulation assumes fixed marginals (Fisher-style sampling, not the chi-sq assumption).
For goodness-of-fit: pass a vector, not a matrix. p must sum to 1 (or set rescale.p = TRUE).
Return object includes $expected, $residuals (Pearson), and $stdres (standardized).

wilcox.test

exact = TRUE by default for small samples with no ties. With ties, normal approximation used.
correct = TRUE applies continuity correction to normal approximation.
conf.int = TRUE computes Hodges-Lehmann estimator and confidence interval (not just the p-value).
Paired test: paired = TRUE uses signed-rank test (Wilcoxon), not rank-sum (Mann-Whitney).

fisher.test

For tables larger than 2x2, uses simulation (simulate.p.value = TRUE) or network algorithm.
workspace controls memory for the network algorithm; increase if you get errors on large tables.
or argument tests a specific odds ratio (default 1) — only for 2x2 tables.

ks.test

Two-sample test or one-sample against a reference distribution.
Does not handle ties well — warns and uses asymptotic approximation.
For composite hypotheses (parameters estimated from data), p-values are conservative (too large). Use dgof or ks.test with exact = NULL for discrete distributions.

p.adjust

Methods: "holm" (default), "BH" (Benjamini-Hochberg FDR), "bonferroni", "BY", "hochberg", "hommel", "fdr" (alias for BH), "none".
n argument: total number of hypotheses (can be larger than length(p) if some p-values are excluded).
Handles NAs: adjusted p-values are NA where input is NA.

pairwise.t.test / pairwise.wilcox.test

p.adjust.method defaults to "holm". Change to "BH" for FDR control.
pool.sd = TRUE (default for t-test): uses pooled SD across all groups (assumes equal variances).
Returns a matrix of p-values, not test statistics.

shapiro.test

Sample size must be between 3 and 5000.
Tests normality; low p-value = evidence against normality.

kmeans

nstart > 1 recommended (e.g., nstart = 25): runs algorithm from multiple random starts, returns best.
Default iter.max = 10 — may be too low for convergence. Increase for large/complex data.
Default algorithm is "Hartigan-Wong" (generally best). Very close points may cause non-convergence (warning with ifault = 4).
Cluster numbering is arbitrary; ordering may differ across platforms.
Always returns k clusters when k is specified (except Lloyd-Forgy may return fewer).

hclust

method = "ward.D2" implements Ward's criterion correctly (using squared distances). The older "ward.D" did not square distances (retained for back-compatibility).
Input must be a dist object. Use as.dist() to convert a symmetric matrix.
hang = -1 in plot() aligns all labels at the bottom.

dist

method = "euclidean" (default). Other options: "manhattan", "maximum", "canberra", "binary", "minkowski".
Returns a dist object (lower triangle only). Use as.matrix() to get full matrix.
"canberra": terms with zero numerator and denominator are omitted from the sum (not treated as 0/0).
Inf values: Euclidean distance involving Inf is Inf. Multiple Infs in same obs give NaN for some methods.

prcomp vs princomp

prcomp uses SVD (numerically superior); princomp uses eigen on covariance (less stable, N-1 vs N scaling).
scale. = TRUE in prcomp standardizes variables; important when variables have very different scales.
princomp standard deviations differ from prcomp by factor sqrt((n-1)/n).
Both return $rotation (loadings) and $x (scores); sign of components may differ between runs.

density

Default bandwidth: bw = "nrd0" (Silverman's rule of thumb). For multimodal data, consider "SJ" or "bcv".
adjust: multiplicative factor on bandwidth. adjust = 0.5 halves the bandwidth (less smooth).
Default kernel: "gaussian". Range of density extends beyond data range (controlled by cut, default 3 bandwidths).
n = 512: number of evaluation points. Increase for smoother plotting.
from/to: explicitly bound the evaluation range.

quantile

Nine type options (1-9). Default type = 7 (R default, linear interpolation). Type 1 = inverse of empirical CDF (SAS default). Types 4-9 are continuous; 1-3 are discontinuous.
na.rm = FALSE by default — returns NA if any NAs present.
names = TRUE by default, adding "0%", "25%", etc. as names.

Distributions (gotchas across all)

All distribution functions follow the d/p/q/r pattern. Common non-obvious points:

n argument in r*() functions: if length(n) > 1, uses length(n) as the count, not n itself. So rnorm(c(1,2,3)) generates 3 values, not 1+2+3.
log = TRUE / log.p = TRUE: compute on log scale for numerical stability in tails.
lower.tail = FALSE gives survival function P(X > x) directly (more accurate than 1 - pnorm() in tails).
Gamma: parameterized by shape and rate (= 1/scale). Default rate = 1. Specifying both rate and scale is an error.
Beta: shape1 (alpha), shape2 (beta) — no mean/sd parameterization.
Poisson dpois: x can be non-integer (returns 0 with a warning for non-integer values if log = FALSE).
Weibull: shape and scale (no rate). R's parameterization: f(x) = (shape/scale)(x/scale)^(shape-1) exp(-(x/scale)^shape).
Lognormal: meanlog and sdlog are mean/sd of the log, not of the distribution itself.

cor.test

Default method: "pearson". Also "kendall" and "spearman".
Returns $estimate, $p.value, $conf.int (CI only for Pearson).
Formula interface: cor.test(~ x + y, data = df) — note the ~ with no LHS.

ecdf

Returns a function (step function). Call it on new values: Fn <- ecdf(x); Fn(3.5).
plot(ecdf(x)) gives the empirical CDF plot.
The returned function is right-continuous with left limits (cadlag).

weighted.mean

Handles NA in weights: observation is dropped if weight is NA.
Weights do not need to sum to 1; they are normalized internally. FILE:references/visualization.md

Visualization — Quick Reference

Non-obvious behaviors, gotchas, and tricky defaults for R functions. Only what Claude doesn't already know.

par (gotchas)

par() settings are per-device. Opening a new device resets everything.
Setting mfrow/mfcol resets cex to 1 and mex to 1. With 2x2 layout, base cex is multiplied by 0.83; with 3+ rows/columns, by 0.66.
mai (inches), mar (lines), pin, plt, pty all interact. Restoring all saved parameters after device resize can produce inconsistent results — last-alphabetically wins.
bg set via par() also sets new = FALSE. Setting fg via par() also sets col.
xpd = NA clips to device region (allows drawing in outer margins); xpd = TRUE clips to figure region; xpd = FALSE (default) clips to plot region.
mgp = c(3, 1, 0): controls title line (mgp[1]), label line (mgp[2]), axis line (mgp[3]). All in mex units.
las: 0 = parallel to axis, 1 = horizontal, 2 = perpendicular, 3 = vertical. Does not respond to srt.
tck = 1 draws grid lines across the plot. tcl = -0.5 (default) gives outward ticks.
usr with log scale: contains log10 of the coordinate limits, not the raw values.
Read-only parameters: cin, cra, csi, cxy, din, page.

layout

layout(mat) where mat is a matrix of integers specifying figure arrangement.
widths/heights accept lcm() for absolute sizes mixed with relative sizes.
More flexible than mfrow/mfcol but cannot be queried once set (unlike par("mfrow")).
layout.show(n) visualizes the layout for debugging.

axis / mtext

axis(side, at, labels): side 1=bottom, 2=left, 3=top, 4=right.
Default gap between axis labels controlled by par("mgp"). Labels can overlap if not managed.
mtext: line argument positions text in margin lines (0 = adjacent to plot, positive = outward). adj controls horizontal position (0-1).
mtext with outer = TRUE writes in the outer margin (set by par(oma = ...)).

curve

First argument can be an expression in x or a function: curve(sin, 0, 2*pi) or curve(x^2 + 1, 0, 10).
add = TRUE to overlay on existing plot. Default n = 101 evaluation points.
xname = "x" by default; change if your expression uses a different variable name.

pairs

panel function receives (x, y, ...) for each pair. lower.panel, upper.panel, diag.panel for different regions.
gap controls spacing between panels (default 1).
Formula interface: pairs(~ var1 + var2 + var3, data = df).

coplot

Conditioning plots: coplot(y ~ x | a) or coplot(y ~ x | a * b) for two conditioning variables.
panel function can be customized; rows/columns control layout.
Default panel draws points; use panel = panel.smooth for loess overlay.

matplot / matlines / matpoints

Plots columns of one matrix against columns of another. Recycles col, lty, pch across columns.
type = "l" by default (unlike plot which defaults to "p").
Useful for plotting multiple time series or fitted curves simultaneously.

contour / filled.contour / image

contour(x, y, z): z must be a matrix with dim = c(length(x), length(y)).
filled.contour has a non-standard layout — it creates its own plot region for the color key. Cannot use par(mfrow) with it. Adding elements requires the plot.axes argument.
image: plots z-values as colored rectangles. Default color scheme may be misleading; set col explicitly.
For image, x and y specify cell boundaries or midpoints depending on context.

persp

persp(x, y, z, theta, phi): theta = azimuthal angle, phi = colatitude.
Returns a transformation matrix (invisible) for projecting 3D to 2D — use trans3d() to add points/lines to the perspective plot.
shade and col control surface shading. border = NA removes grid lines.

segments / arrows / rect / polygon

All take vectorized coordinates; recycle as needed.
arrows: code = 1 (head at start), code = 2 (head at end, default), code = 3 (both).
polygon: last point auto-connects to first. Fill with col; border controls outline.
rect(xleft, ybottom, xright, ytop) — note argument order is not the same as other systems.

dev / dev.off / dev.copy

dev.new() opens a new device. dev.off() closes current device (and flushes output for file devices like pdf).
dev.off() on the last open device reverts to null device.
dev.copy(pdf, file = "plot.pdf") followed by dev.off() to save current plot.
dev.list() returns all open devices; dev.cur() the active one.

pdf

Must call dev.off() to finalize the file. Without it, file may be empty/corrupt.
onefile = TRUE (default): multiple pages in one PDF. onefile = FALSE: one file per page (uses %d in filename for numbering).
useDingbats = FALSE recommended to avoid issues with certain PDF viewers and pch symbols.
Default size: 7x7 inches. family controls font family.

png / bitmap devices

res controls DPI (default 72). For publication: res = 300 with appropriate width/height in pixels or inches (with units = "in").
type = "cairo" (on systems with cairo) gives better antialiasing than default.
bg = "transparent" for transparent background (PNG supports alpha).

colors / rgb / hcl / col2rgb

colors() returns all 657 named colors. col2rgb("color") returns RGB matrix.
rgb(r, g, b, alpha, maxColorValue = 255) — note maxColorValue default is 1, not 255.
hcl(h, c, l): perceptually uniform color space. Preferred for color scales.
adjustcolor(col, alpha.f = 0.5): easy way to add transparency.

colorRamp / colorRampPalette

colorRamp returns a function mapping [0,1] to RGB matrix.
colorRampPalette returns a function taking n and returning n interpolated colors.
space = "Lab" gives more perceptually uniform interpolation than "rgb".

palette / recordPlot

palette() returns current palette (default 8 colors). palette("Set1") sets a built-in palette.
Integer colors in plots index into the palette (with wrapping). Index 0 = background color.
recordPlot() / replayPlot(): save and restore a complete plot — device-dependent and fragile across sessions. FILE:assets/analysis_template.R

============================================================

Analysis Template — Base R

Copy this file, rename it, and fill in your details.

============================================================

Author :

Date :

Data :

Purpose :

============================================================

── 0. Setup ─────────────────────────────────────────────────

Clear environment (optional — comment out if loading into existing session)

rm(list = ls())

Set working directory if needed

setwd("/path/to/your/project")

Reproducibility

set.seed(42)

Libraries — uncomment what you need

library(haven) # read .dta / .sav / .sas

library(readxl) # read Excel files

library(openxlsx) # write Excel files

library(foreign) # older Stata / SPSS formats

library(survey) # survey-weighted analysis

library(lmtest) # Breusch-Pagan, Durbin-Watson etc.

library(sandwich) # robust standard errors

library(car) # Type II/III ANOVA, VIF

── 1. Load Data ─────────────────────────────────────────────

df <- read.csv("your_data.csv", stringsAsFactors = FALSE)

df <- readRDS("your_data.rds")

df <- haven::read_dta("your_data.dta")

First look — always run these

dim(df) str(df) head(df, 10) summary(df)

── 2. Data Quality Check ────────────────────────────────────

Missing values

na_report <- data.frame( column = names(df), n_miss = colSums(is.na(df)), pct_miss = round(colMeans(is.na(df)) * 100, 1), row.names = NULL ) print(na_report[na_report$n_miss > 0, ])

Duplicates

n_dup <- sum(duplicated(df)) cat(sprintf("Duplicate rows: %d\n", n_dup))

Unique values for categorical columns

cat_cols <- names(df)[sapply(df, function(x) is.character(x) | is.factor(x))] for (col in cat_cols) { cat(sprintf("\n%s (%d unique):\n", col, length(unique(df[[col]])))) print(table(df[[col]], useNA = "ifany")) }

── 3. Clean & Transform ─────────────────────────────────────

Rename columns (example)

names(df)[names(df) == "old_name"] <- "new_name"

Convert types

df$group <- as.factor(df$group)

df$date <- as.Date(df$date, format = "%Y-%m-%d")

Recode values (example)

df$gender <- ifelse(df$gender == 1, "Male", "Female")

Create new variables (example)

df$log_income <- log(df$income + 1)

df$age_group <- cut(df$age,

breaks = c(0, 25, 45, 65, Inf),

labels = c("18-25", "26-45", "46-65", "65+"))

Filter rows (example)

df <- df[df$year >= 2010, ]

df <- df[complete.cases(df[, c("outcome", "predictor")]), ]

Drop unused factor levels

df <- droplevels(df)

── 4. Descriptive Statistics ────────────────────────────────

Numeric summary

Cross-tabulation

table(df$group, df$category, useNA = "ifany")

prop.table(table(df$group, df$category), margin = 1) # row proportions

── 5. Visualization (EDA) ───────────────────────────────────

par(mfrow = c(2, 2))

Histogram of main outcome

hist(df$outcome_var, main = "Distribution of Outcome", xlab = "Outcome", col = "steelblue", border = "white", breaks = 30)

Boxplot by group

boxplot(outcome_var ~ group_var, data = df, main = "Outcome by Group", col = "lightyellow", las = 2)

Scatter plot

Correlation matrix (numeric columns only)

cor_mat <- cor(df[num_cols], use = "complete.obs") image(cor_mat, main = "Correlation Matrix", col = hcl.colors(20, "RdBu", rev = TRUE))

par(mfrow = c(1, 1))

── 6. Analysis ───────────────────────────────────────────────

·· 6a. Comparison of means ··

t.test(outcome_var ~ group_var, data = df)

·· 6b. Linear regression ··

fit <- lm(outcome_var ~ predictor1 + predictor2 + group_var, data = df) summary(fit) confint(fit)

Check VIF for multicollinearity (requires car)

car::vif(fit)

Robust standard errors (requires lmtest + sandwich)

lmtest::coeftest(fit, vcov = sandwich::vcovHC(fit, type = "HC3"))

·· 6c. ANOVA ··

fit_aov <- aov(outcome_var ~ group_var, data = df)

summary(fit_aov)

TukeyHSD(fit_aov)

·· 6d. Logistic regression (binary outcome) ··

fit_logit <- glm(binary_outcome ~ x1 + x2,

data = df,

family = binomial(link = "logit"))

summary(fit_logit)

exp(coef(fit_logit)) # odds ratios

exp(confint(fit_logit)) # OR confidence intervals

── 7. Model Diagnostics ─────────────────────────────────────

par(mfrow = c(2, 2)) plot(fit) par(mfrow = c(1, 1))

Residual normality

shapiro.test(residuals(fit))

Homoscedasticity (requires lmtest)

lmtest::bptest(fit)

── 8. Save Output ────────────────────────────────────────────

Cleaned data

write.csv(df, "data_clean.csv", row.names = FALSE)

saveRDS(df, "data_clean.rds")

Model results to text file

sink("results.txt")

cat("=== Linear Model ===\n")

print(summary(fit))

cat("\n=== Confidence Intervals ===\n")

print(confint(fit))

sink()

Plots to file

png("figure1_distributions.png", width = 1200, height = 900, res = 150)

par(mfrow = c(2, 2))

# ... your plots ...

par(mfrow = c(1, 1))

dev.off()

============================================================

END OF TEMPLATE

============================================================

FILE:scripts/check_data.R

check_data.R — Quick data quality report for any R data frame

Usage: source("check_data.R") then call check_data(df)

Or: source("check_data.R"); check_data(read.csv("yourfile.csv"))

check_data <- function(df, top_n_levels = 8) {

if (!is.data.frame(df)) stop("Input must be a data frame.")

n_row <- nrow(df) n_col <- ncol(df)

── 1. Column overview ──────────────────────

cat("── COLUMN OVERVIEW ────────────────────────\n")

for (col in names(df)) { x <- df[[col]] cls <- class(x)[1] n_na <- sum(is.na(x)) pct <- round(n_na / n_row * 100, 1) n_uniq <- length(unique(x[!is.na(x)]))

na_flag <- if (n_na == 0) "" else sprintf("  *** %d NAs (%.1f%%)", n_na, pct)
cat(sprintf("  %-20s  %-12s  %d unique%s\n",
            col, cls, n_uniq, na_flag))

}

── 2. NA summary ────────────────────────────

cat("\n── NA SUMMARY ─────────────────────────────\n")

na_counts <- sapply(df, function(x) sum(is.na(x))) cols_with_na <- na_counts[na_counts > 0]

if (length(cols_with_na) == 0) { cat(" No missing values. \n") } else { cat(sprintf(" Columns with NAs: %d of %d\n\n", length(cols_with_na), n_col)) for (col in names(cols_with_na)) { bar_len <- round(cols_with_na[col] / n_row * 20) bar <- paste0(rep("█", bar_len), collapse = "") pct_na <- round(cols_with_na[col] / n_row * 100, 1) cat(sprintf(" %-20s [%-20s] %d (%.1f%%)\n", col, bar, cols_with_na[col], pct_na)) } }

── 3. Numeric columns ───────────────────────

num_cols <- names(df)[sapply(df, is.numeric)]

if (length(num_cols) > 0) { cat("\n── NUMERIC COLUMNS ────────────────────────\n") cat(sprintf(" %-20s %8s %8s %8s %8s %8s\n", "Column", "Min", "Mean", "Median", "Max", "SD")) cat(sprintf(" %-20s %8s %8s %8s %8s %8s\n", "──────", "───", "────", "──────", "───", "──"))

for (col in num_cols) {
  x  <- df[[col]][!is.na(df[[col]])]
  if (length(x) == 0) next
  cat(sprintf("  %-20s  %8.3g  %8.3g  %8.3g  %8.3g  %8.3g\n",
              col,
              min(x), mean(x), median(x), max(x), sd(x)))
}

}

── 4. Factor / character columns ───────────

cat_cols <- names(df)[sapply(df, function(x) is.factor(x) | is.character(x))]

if (length(cat_cols) > 0) { cat("\n── CATEGORICAL COLUMNS ────────────────────\n")

for (col in cat_cols) {
  x    <- df[[col]]
  tbl  <- sort(table(x, useNA = "no"), decreasing = TRUE)
  n_lv <- length(tbl)
  cat(sprintf("\n  %s  (%d unique values)\n", col, n_lv))
  
  show <- min(top_n_levels, n_lv)
  for (i in seq_len(show)) {
    lbl <- names(tbl)[i]
    cnt <- tbl[i]
    pct <- round(cnt / n_row * 100, 1)
    cat(sprintf("    %-25s  %5d  (%.1f%%)\n", lbl, cnt, pct))
  }
  if (n_lv > top_n_levels) {
    cat(sprintf("    ... and %d more levels\n", n_lv - top_n_levels))
  }
}

}

── 5. Duplicate rows ────────────────────────

cat("\n── DUPLICATES ─────────────────────────────\n") n_dup <- sum(duplicated(df)) if (n_dup == 0) { cat(" No duplicate rows.\n") } else { cat(sprintf(" %d duplicate row(s) found (%.1f%% of data)\n", n_dup, n_dup / n_row * 100)) }

cat("\n══════════════════════════════════════════\n") cat(" END OF REPORT\n") cat("══════════════════════════════════════════\n")

Return invisibly for programmatic use

invisible(list( dims = c(rows = n_row, cols = n_col), na_counts = na_counts, n_dupes = n_dup )) } FILE:scripts/scaffold_analysis.R #!/usr/bin/env Rscript

scaffold_analysis.R — Generates a starter analysis script

Usage (from terminal):

Rscript scaffold_analysis.R myproject

Rscript scaffold_analysis.R myproject outcome_var group_var

Usage (from R console):

source("scaffold_analysis.R")

scaffold_analysis("myproject", outcome = "score", group = "treatment")

Output: myproject_analysis.R (ready to edit)

scaffold_analysis <- function(project_name, outcome = "outcome", group = "group", data_file = NULL) {

if (is.null(data_file)) data_file <- paste0(project_name, ".csv") out_file <- paste0(project_name, "_analysis.R")

template <- sprintf( '# ============================================================

Project : %s

Created : %s

============================================================

── 0. Libraries ─────────────────────────────────────────────

Add packages you need here

library(ggplot2)

library(haven) # for .dta files

library(openxlsx) # for Excel output

── 1. Load Data ─────────────────────────────────────────────

df <- read.csv("%s", stringsAsFactors = FALSE)

Quick check — always do this first

cat("Dimensions:", dim(df), "\n") str(df) head(df)

── 2. Explore / EDA ─────────────────────────────────────────

summary(df)

NA check

na_counts <- colSums(is.na(df)) na_counts[na_counts > 0]

Key variable distributions

hist(df$%s, main = "Distribution of %s", xlab = "%s")

if ("%s" %%in%% names(df)) { table(df$%s) barplot(table(df$%s), main = "Counts by %s", col = "steelblue", las = 2) }

── 3. Clean / Transform ──────────────────────────────────────

df <- df[complete.cases(df), ] # drop rows with any NA

df$%s <- as.factor(df$%s) # convert to factor

── 4. Analysis ───────────────────────────────────────────────

Descriptive stats by group

tapply(df$%s, df$%s, mean, na.rm = TRUE) tapply(df$%s, df$%s, sd, na.rm = TRUE)

t-test (two groups)

t.test(%s ~ %s, data = df)

Linear model

fit <- lm(%s ~ %s, data = df) summary(fit) confint(fit)

ANOVA (multiple groups)

fit_aov <- aov(%s ~ %s, data = df)

summary(fit_aov)

TukeyHSD(fit_aov)

── 5. Visualize Results ──────────────────────────────────────

par(mfrow = c(1, 2))

Boxplot by group

boxplot(%s ~ %s, data = df, main = "%s by %s", xlab = "%s", ylab = "%s", col = "lightyellow")

Model diagnostics

plot(fit, which = 1) # residuals vs fitted

par(mfrow = c(1, 1))

── 6. Save Output ────────────────────────────────────────────

Save cleaned data

write.csv(df, "%s_clean.csv", row.names = FALSE)

Save model summary to text

sink("%s_results.txt")

summary(fit)

sink()

Save plot to file

png("%s_boxplot.png", width = 800, height = 600, res = 150)

boxplot(%s ~ %s, data = df, col = "lightyellow")

dev.off()

', project_name, format(Sys.Date(), "%%Y-%%m-%%d"), data_file, # Section 2 — EDA outcome, outcome, outcome, group, group, group, group, # Section 3 group, group, # Section 4 outcome, group, outcome, group, outcome, group, outcome, group, outcome, group, outcome, group, # Section 5 outcome, group, outcome, group, group, outcome, # Section 6 project_name, project_name, project_name, outcome, group )

writeLines(template, out_file) cat(sprintf("Created: %s\n", out_file)) invisible(out_file) }

── Run from command line ─────────────────────────────────────

if (!interactive()) { args <- commandArgs(trailingOnly = TRUE)

if (length(args) == 0) { cat("Usage: Rscript scaffold_analysis.R [outcome_var] [group_var]\n") cat("Example: Rscript scaffold_analysis.R myproject score treatment\n") quit(status = 1) }

project <- args[1] outcome <- if (length(args) >= 2) args[2] else "outcome" group <- if (length(args) >= 3) args[3] else "group"

scaffold_analysis(project, outcome = outcome, group = group) } FILE:README.md

base-r-skill

GitHub: https://github.com/iremaydas/base-r-skill

A Claude Code skill for base R programming.

The Story

I'm a political science PhD candidate who uses R regularly but would never call myself an R person. I needed a Claude Code skill for base R — something without tidyverse, without ggplot2, just plain R — and I couldn't find one anywhere.

So I made one myself. At 11pm. Asking Claude to help me build a skill for Claude.

If you're also someone who Googles how to drop NA rows in R every single time, this one's for you. 🫶

What's Inside

base-r/
├── SKILL.md                    # Main skill file
├── references/                 # Gotchas & non-obvious behaviors
│   ├── data-wrangling.md       # Subsetting traps, apply family, merge, factor quirks
│   ├── modeling.md             # Formula syntax, lm/glm/aov/nls, optim
│   ├── statistics.md           # Hypothesis tests, distributions, clustering
│   ├── visualization.md        # par, layout, devices, colors
│   ├── io-and-text.md          # read.table, grep, regex, format
│   ├── dates-and-system.md     # Date/POSIXct traps, options(), file ops
│   └── misc-utilities.md       # tryCatch, do.call, time series, utilities
├── scripts/
│   ├── check_data.R            # Quick data quality report for any data frame
│   └── scaffold_analysis.R     # Generates a starter analysis script
└── assets/
    └── analysis_template.R     # Copy-paste analysis template

The reference files were condensed from the official R 4.5.3 manual — 19,518 lines → 945 lines (95% reduction). Only the non-obvious stuff survived: gotchas, surprising defaults, tricky interactions. The things Claude already knows well got cut.

How to Use

Add this skill to your Claude Code setup by pointing to this repo. Then Claude will automatically load the relevant reference files when you're working on R tasks.

Works best for:

Base R data manipulation (no tidyverse)
Statistical modeling with lm, glm, aov
Base graphics with plot, par, barplot
Understanding why your R code is doing that weird thing

Not for: tidyverse, ggplot2, Shiny, or R package development.

The `check_data.R` Script

Probably the most useful standalone thing here. Source it and run check_data(df) on any data frame to get a formatted report of dimensions, NA counts, numeric summaries, and categorical breakdowns.

source("scripts/check_data.R")
check_data(your_df)

Built With Help From

Claude (obviously)
The official R manuals (all 19,518 lines of them)
Mild frustration and several cups of coffee

Contributing

If you spot a missing gotcha, a wrong default, or something that should be in the references — PRs are very welcome. I'm learning too.

Made by @iremaydas — PhD candidate, occasional R user, full-time Googler of things I should probably know by now.