R Stats Skill

R's statistical modeling ecosystem centered on base stats with key extension packages: sandwich for robust/clustered variance estimators, lmtest for coefficient testing with custom covariance matrices, car for regression diagnostics and hypothesis tests, broom for tidy model output, and MASS for negative binomial regression. Covers OLS/WLS/GLS via lm(), GLMs via glm() (logit, probit, Poisson, negative binomial, Gamma), robust and clustered standard errors via sandwich+lmtest, diagnostics (VIF, heteroskedasticity tests, normality tests), tidy output via broom, and classical hypothesis tests (t, chi-squared, Wilcoxon, Fisher). Use when execution language is R and the model does not require absorbed fixed effects or panel structure. Python equivalent: statsmodels. For fixed effects regressions use fixest; for panel RE/between/Fama-MacBeth use plm; for complex survey designs use survey-r.

What is R Stats?

R is the original statistical computing environment. The stats package ships with every R installation and provides the core modeling functions that the entire R ecosystem builds upon:

• Formula interface: y ~ x1 + x2 — the original formula syntax that Python libraries (patsy, formulaic) later adopted
• Linear models: lm() for OLS, WLS, GLS; summary(), confint(), predict(), anova() for post-estimation
• Generalized linear models: glm() with family/link specifications (binomial/logit, Poisson/log, Gamma/inverse, etc.)
• Classical tests: t.test(), chisq.test(), wilcox.test(), fisher.test(), prop.test(), ks.test(), cor.test()
• Time series: ts(), arima(), acf(), pacf(), stl()

The extension packages covered by this skill complete the modeling toolkit:

Package	Version	Role
sandwich	3.1-1	Robust (HC0-HC4) and clustered variance-covariance estimators
lmtest	0.9-40	coeftest() with custom vcov; bptest(), dwtest(), bgtest()
car	3.1-5	vif(), linearHypothesis(), ncvTest(), Anova (type II/III)
broom	1.0.12	tidy(), glance(), augment() for tidy model output
MASS	7.3-65	glm.nb() for negative binomial; stepAIC()
modelsummary	2.6.0	Publication-quality regression tables (multiple models)
marginaleffects	0.32.0	Average marginal effects, contrasts, predictions

Version Notes

• R 4.5.3: Base R stats included. The |> native pipe (R 4.1+) is the preferred pipe style in DAAF scripts (not magrittr %>%).
• sandwich 3.1-1: vcovHC types HC0-HC4, vcovCL for clustering, vcovHAC for time series, vcovBS for bootstrap.
• lmtest 0.9-40: coeftest(), waldtest(), bptest(), dwtest(), bgtest(), resettest(), grangertest().
• car 3.1-5: vif(), linearHypothesis(), Anova() type II/III, ncvTest(), influencePlot().
• broom 1.0.12: Tidiers for 100+ model classes including lm, glm, nls, t.test, and many extension package models.
• MASS 7.3-65: Ships with R. glm.nb() for negative binomial, polr() for ordered logit/probit, rlm() for robust regression.
• modelsummary 2.6.0: Supports lm, glm, fixest, plm, and many other model classes. Output formats: markdown, LaTeX, HTML, Word, PNG.
• marginaleffects 0.32.0: avg_slopes(), avg_comparisons(), avg_predictions(), hypotheses() for delta-method inference.

How to Use This Skill

Reference File Structure

File	Purpose	When to Read
quickstart.md	lm() basics, summary(), confint(), predict(), formula syntax	Starting with R stats or need formula reference
glm.md	glm() families, MASS::glm.nb(), link functions, deviance, overdispersion	Non-linear models, binary/count outcomes
robust-se.md	sandwich vcovHC/vcovCL/vcovHAC, lmtest coeftest()	Robust or clustered standard errors
diagnostics.md	car::vif(), linearHypothesis(), residual analysis, heteroskedasticity tests	Checking model assumptions
reporting.md	broom tidy/glance/augment, modelsummary, marginaleffects	Tidy output, publication tables, marginal effects
tests.md	t.test(), chisq.test(), wilcox.test(), fisher.test(), correlation tests	Classical hypothesis tests
time-series.md	ts objects, arima(), acf/pacf, Durbin-Watson, ADF test	Time series analysis
gotchas.md	Formula pitfalls, factor contrasts, na.action, predict() type arg	Debugging common issues

Reading Order

1. New to R stats? Start with quickstart.md then robust-se.md
2. Need GLM or logit/probit? Read quickstart.md then glm.md
3. Need robust or clustered SEs? Read robust-se.md
4. Checking model assumptions? Read diagnostics.md
5. Building publication tables? Read reporting.md
6. Coming from Python? Read quickstart.md (R is the original formula syntax)

The reference-file routing in this skill applies to advisory and brainstorming turns as much as implementation. Recommending a model, reviewing an analysis plan, or answering a question that touches a routed topic calls for reading the routed reference file just as much as writing code does — the reference files carry curated caveats and environment-specific constraints (e.g., which packages are pre-installed, sandwich function signatures) that this overview and general knowledge lack.

Related Skills

• statsmodels: Python equivalent for OLS/GLM/time series. If execution language is Python, use statsmodels instead of r-stats. API patterns differ (statsmodels requires .fit() call; R returns fitted object from lm() directly)
• fixest: Use instead of base R stats when model needs absorbed fixed effects, IV with FE, or DiD. fixest is faster for FE models; base R stats is broader for GLMs and classical tests
• plm: Use for panel data models (FE within, RE, between, first difference, Fama-MacBeth). plm extends the lm/formula interface for panel structure
• survey-r: Use for survey-weighted regression with complex designs. Base R lm() with weights is NOT equivalent to survey-weighted regression — it does not account for stratification, clustering, or finite population corrections
• gt: Publication-quality tables with gt and modelsummary. Load for formatted regression output beyond what broom provides -- modelsummary (covered in this skill's reporting.md) uses gt as its default backend
• data-scientist: Provides methodology guidance (when to use which model, assumption checking protocol). Load alongside r-stats for the "why"; r-stats provides the "how"
• r-python-translation: Cross-language reference mapping R stats functions to Python equivalents. Load when annotating R code for Python-background users

Quick Decision Trees

"I need to fit a regression model"

What kind of regression?
+-- Linear (continuous outcome)
|   +-- Basic OLS --> ./references/quickstart.md
|   +-- Weighted least squares --> ./references/quickstart.md (WLS section)
|   |   (WLS != survey-weighted regression -- for complex surveys, use survey-r)
|   +-- Robust to outliers (M-estimator) --> MASS::rlm() in ./references/quickstart.md
|   +-- Quantile regression --> quantreg::rq() (external package)
|   +-- Need fixed effects?
|       +-- Use fixest instead (faster FE absorption)
+-- Binary outcome (0/1)
|   +-- Logit --> ./references/glm.md
|   +-- Probit --> ./references/glm.md
+-- Count outcome (0, 1, 2, ...)
|   +-- Poisson --> ./references/glm.md
|   +-- Negative binomial --> ./references/glm.md (MASS::glm.nb)
|   +-- Zero-inflated --> pscl package (external)
+-- Ordinal (ordered categories)
|   +-- MASS::polr() --> ./references/glm.md
+-- GLM (custom family/link)
|   +-- glm() framework --> ./references/glm.md
+-- Need robust/clustered SEs?
    +-- sandwich + lmtest --> ./references/robust-se.md

"I need to check model assumptions"

What assumption to check?
+-- Heteroskedasticity --> ./references/diagnostics.md
|   +-- Breusch-Pagan (lmtest::bptest)
|   +-- NCV test (car::ncvTest)
+-- Normality of residuals --> ./references/diagnostics.md
|   +-- Shapiro-Wilk (shapiro.test)
|   +-- Jarque-Bera (tseries::jarque.bera.test)
+-- Multicollinearity --> ./references/diagnostics.md
|   +-- VIF (car::vif)
+-- Specification / functional form --> ./references/diagnostics.md
|   +-- RESET test (lmtest::resettest)
+-- Serial correlation --> ./references/diagnostics.md
|   +-- Durbin-Watson (lmtest::dwtest)
|   +-- Breusch-Godfrey (lmtest::bgtest)
+-- Joint hypothesis test --> ./references/diagnostics.md
|   +-- car::linearHypothesis
+-- Influential observations --> ./references/diagnostics.md
    +-- Cook's distance, leverage, DFBETAS

"I need classical hypothesis tests"

What test?
+-- Compare two means --> ./references/tests.md
|   +-- t.test() (parametric)
|   +-- wilcox.test() (non-parametric)
+-- Compare proportions --> ./references/tests.md
|   +-- prop.test()
+-- Test independence --> ./references/tests.md
|   +-- chisq.test() (large sample)
|   +-- fisher.test() (exact, small sample)
+-- Test correlation --> ./references/tests.md
|   +-- cor.test() (Pearson, Spearman, Kendall)
+-- Goodness of fit --> ./references/tests.md
|   +-- chisq.test() with expected frequencies
|   +-- ks.test() (Kolmogorov-Smirnov)

"Something isn't working"

Common issues?
+-- Formula I() for arithmetic --> ./references/gotchas.md
+-- Factor contrasts (treatment vs sum) --> ./references/gotchas.md
+-- predict() type argument --> ./references/gotchas.md
+-- na.action behavior --> ./references/gotchas.md
+-- summary() vs anova() --> ./references/gotchas.md
+-- Convergence warnings (GLM) --> ./references/gotchas.md
+-- r-stats vs fixest boundary --> ./references/gotchas.md

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), R analyses are executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern:

1. Write model code to scripts/stage8_analysis/{step}_{model-name}.R
2. Execute via Bash with automatic output capture wrapper script
3. Validation results get automatically embedded in scripts as comments
4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol covering complete code file writing, output capture, and file versioning rules.

See:

• agent_reference/SCRIPT_EXECUTION_REFERENCE.md -- Script execution protocol and format with validation

The examples below show R syntax. In research workflows, wrap them in scripts following the file-first pattern.

---

Quick Reference

Essential Library Loading

# --- Config ---
library(sandwich)         # Robust variance estimators
library(lmtest)           # Coefficient testing with custom vcov
library(car)              # Diagnostics: VIF, linearHypothesis
library(broom)            # Tidy model output
library(MASS)             # glm.nb(), polr(), rlm()
library(modelsummary)     # Publication tables
library(marginaleffects)  # Marginal effects

Core Operations

Operation	Code
OLS	lm(y ~ x1 + x2, data = df)
WLS	lm(y ~ x1 + x2, data = df, weights = w)
Logit	glm(y ~ x1 + x2, data = df, family = binomial)
Probit	glm(y ~ x1 + x2, data = df, family = binomial(link = "probit"))
Poisson	glm(count ~ x1 + x2, data = df, family = poisson)
Neg. binomial	MASS::glm.nb(count ~ x1 + x2, data = df)
Gamma	glm(y ~ x1 + x2, data = df, family = Gamma(link = "log"))
Robust SE (HC1)	coeftest(fit, vcov = vcovHC(fit, type = "HC1"))
Clustered SE	coeftest(fit, vcov = vcovCL(fit, cluster = df$group))
Summary	summary(fit)
Tidy output	broom::tidy(fit, conf.int = TRUE)
Model stats	broom::glance(fit)
Predict	predict(fit, newdata = new_df)
Confidence intervals	confint(fit)
VIF	car::vif(fit)
Joint hypothesis	car::linearHypothesis(fit, c("x1 = 0", "x2 = 0"))
Marginal effects	marginaleffects::avg_slopes(fit)
Pub table	modelsummary::modelsummary(list(fit1, fit2))

Formula Syntax

# Additive terms
y ~ x1 + x2 + x3

# Interaction (with main effects)
y ~ x1 * x2           # equivalent to x1 + x2 + x1:x2

# Interaction only (no main effects)
y ~ x1 : x2

# Factor variable (explicit)
y ~ factor(region)

# Suppress intercept
y ~ x1 + x2 - 1       # or: y ~ 0 + x1 + x2

# Polynomial: I() protects arithmetic operators
y ~ x1 + I(x1^2)

# Log transformation
y ~ log(income) + age

# All pairwise interactions
y ~ (x1 + x2 + x3)^2

Topic Index

Topic	Reference File
lm() basics	./references/quickstart.md
Formula syntax	./references/quickstart.md
summary() output	./references/quickstart.md
confint()	./references/quickstart.md
predict()	./references/quickstart.md
WLS / weighted regression	./references/quickstart.md
Model matrix / design matrix	./references/quickstart.md
anova() vs summary()	./references/quickstart.md
Comparison to statsmodels	./references/quickstart.md
glm() framework	./references/glm.md
Logit / probit	./references/glm.md
Poisson regression	./references/glm.md
Negative binomial (MASS)	./references/glm.md
Ordered logit/probit (MASS)	./references/glm.md
GLM families and links	./references/glm.md
Deviance and overdispersion	./references/glm.md
Odds ratios and IRR	./references/glm.md
Marginal effects	./references/reporting.md
Robust SEs (HC0-HC4)	./references/robust-se.md
Clustered SEs	./references/robust-se.md
HAC / Newey-West SEs	./references/robust-se.md
Bootstrap SEs	./references/robust-se.md
coeftest()	./references/robust-se.md
waldtest()	./references/robust-se.md
VIF / multicollinearity	./references/diagnostics.md
Breusch-Pagan test	./references/diagnostics.md
NCV test	./references/diagnostics.md
Shapiro-Wilk test	./references/diagnostics.md
RESET test	./references/diagnostics.md
Durbin-Watson test	./references/diagnostics.md
Breusch-Godfrey test	./references/diagnostics.md
linearHypothesis()	./references/diagnostics.md
Cook's distance	./references/diagnostics.md
Influence plots	./references/diagnostics.md
Residual analysis	./references/diagnostics.md
broom::tidy()	./references/reporting.md
broom::glance()	./references/reporting.md
broom::augment()	./references/reporting.md
modelsummary tables	./references/reporting.md
marginaleffects	./references/reporting.md
t.test()	./references/tests.md
chisq.test()	./references/tests.md
wilcox.test()	./references/tests.md
fisher.test()	./references/tests.md
prop.test()	./references/tests.md
ks.test()	./references/tests.md
cor.test()	./references/tests.md
ts objects	./references/time-series.md
arima()	./references/time-series.md
acf / pacf	./references/time-series.md
Durbin-Watson (time series)	./references/time-series.md
ADF test	./references/time-series.md
Formula pitfalls (I())	./references/gotchas.md
Factor contrasts	./references/gotchas.md
na.action behavior	./references/gotchas.md
predict() type argument	./references/gotchas.md
Convergence warnings	./references/gotchas.md
r-stats vs fixest boundary	./references/gotchas.md

Citation

When base R stats functions are used as a primary analytical tool, include in the report's Software & Tools references:

R Core Team (2026). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

For extension packages, cite when they contribute substantially to the analysis:

Zeileis, A. (2004). "Econometric Computing with HC and HAC Covariance Matrix Estimators." Journal of Statistical Software, 11(10), 1-17. (sandwich package)

Zeileis, A. & Hothorn, T. (2002). "Diagnostic Checking in Regression Relationships." R News, 2(3), 7-10. (lmtest package)

Fox, J. & Weisberg, S. (2019). An R Companion to Applied Regression, Third Edition. Sage, Thousand Oaks CA. (car package)

Robinson, D., Hayes, A., & Couch, S. (2023). broom: Convert Statistical Objects into Tidy Tibbles. R package. (broom package)

Arel-Bundock, V. (2022). "modelsummary: Data and Model Summaries in R." Journal of Statistical Software, 103(1), 1-23. (modelsummary package)

Arel-Bundock, V., Greifer, N., & Heiss, A. (2024). "How to Interpret Statistical Models Using marginaleffects for R and Python." Journal of Statistical Software, 111(9), 1-32. (marginaleffects package)

Cite when: The package's functions are central to the analytical approach (e.g., sandwich for robust inference, car for diagnostics). Do not cite when: Only used incidentally for output formatting.