R's built-in stats package covers the full range of classical statistics: linear/generalized models, hypothesis tests, probability distributions, time series, and clustering.
Load: library(stats) (base R auto-loaded)
| Function | Signature | Description |
|---|---|---|
| lm() | lm(y ~ x1 + x2, data=df) → lm object | Fit linear regression model |
| summary(lm) | summary(model) → summary.lm | Coefficients, R², p-values, F-statistic |
| predict() | predict(model, newdata=df) → vector | Predict from fitted model |
| glm() | glm(y~x, family=binomial, data=df) | Generalized linear model |
| t.test() | t.test(x, y, paired=FALSE) → htest | One or two-sample t-test |
| wilcox.test() | wilcox.test(x, y) → htest | Non-parametric Wilcoxon test |
| chisq.test() | chisq.test(table) → htest | Chi-squared test of independence |
| aov() | aov(y ~ group, data=df) → aov | Analysis of variance (ANOVA) |
| cor() | cor(x, y, method='pearson') → numeric | Correlation coefficient |
| cor.test() | cor.test(x, y) → htest | Correlation test with p-value |
| shapiro.test() | shapiro.test(x) → htest | Normality test (Shapiro-Wilk) |
| ks.test() | ks.test(x, 'pnorm') → htest | Kolmogorov-Smirnov test |
| dnorm() | dnorm(x, mean=0, sd=1) | Normal distribution density |
| pnorm() | pnorm(q, mean=0, sd=1) | Normal CDF — P(X <= q) |
| qnorm() | qnorm(p, mean=0, sd=1) | Normal inverse CDF — quantile |
| rnorm() | rnorm(n, mean=0, sd=1) | Random normal samples |
| pt/qt/dt | pt(q,df) / qt(p,df) / dt(x,df) | Student's t distribution |
| pchisq/qchisq | pchisq(q,df) / qchisq(p,df) | Chi-squared distribution |
| pbinom | pbinom(q, size, prob) | Binomial distribution CDF |
| rbinom | rbinom(n, size, prob) | Random binomial samples |
| quantile() | quantile(x, probs=c(.25,.5,.75)) | Sample quantiles |
| prcomp() | prcomp(df, scale.=TRUE) | Principal component analysis (PCA) |
| kmeans() | kmeans(df, centers=3) | K-means clustering |
| hclust() | hclust(dist(df)) | Hierarchical clustering |
| dist() | dist(x, method='euclidean') | Distance matrix |
| pairwise.t.test() | pairwise.t.test(x, g, p.adjust.method='bonferroni') | Post-hoc pairwise t-tests |
| TukeyHSD() | TukeyHSD(aov_obj) | Tukey HSD post-hoc test after ANOVA |
Run with Rscript script.R.
set.seed(2026)
math <- rnorm(50, mean=82, sd=10)
science <- rnorm(50, mean=78, sd=12)
cat("=== Descriptive Statistics ===\n")
cat("Math: mean=", round(mean(math),2), " sd=", round(sd(math),2), "\n")
cat("Science: mean=", round(mean(science),2), " sd=", round(sd(science),2), "\n")
cat("\n=== Two-Sample T-Test ===\n")
t_res <- t.test(math, science)
cat("t=", round(t_res$statistic,4), " df=", round(t_res$parameter,2),
" p=", round(t_res$p.value,4), "\n")
cat("Significant?", ifelse(t_res$p.value < 0.05,"YES","NO"), "\n")
cat("\n=== Linear Regression ===\n")
n <- 100
x <- runif(n, 0, 10)
y <- 2.5*x + 5 + rnorm(n, 0, 2)
df <- data.frame(x=x, y=y)
mdl <- lm(y ~ x, data=df)
sm <- summary(mdl)
cat("Intercept:", round(coef(mdl)[1],4), "\n")
cat("Slope: ", round(coef(mdl)[2],4), "\n")
cat("R-squared:", round(sm$r.squared,4), "\n")
cat("P(slope): ", round(sm$coefficients[2,4],6), "\n")
preds <- predict(mdl, data.frame(x=c(2,5,8)), interval="confidence")
cat("\n=== Predictions ===\n")
print(round(cbind(x=c(2,5,8), preds), 3))
cat("\n=== Normality Test ===\n")
sw <- shapiro.test(math)
cat("Shapiro-Wilk: W=", round(sw$statistic,4), " p=", round(sw$p.value,4), "\n")
cat("Normal?", ifelse(sw$p.value>0.05,"YES","NO"), "\n")
# Rscript stats_demo.R