R语言描述性统计分析：从summary到分组与可视化核对

简介：一份面向R语言初学者与数据分析入门者的描述性统计实战指南，聚焦如何用R快速完成数据集中趋势、离散程度与分布形状的探索，覆盖从内置函数到分组分析与可视化的完整链路。压缩包内含1个PDF文档，大小约98KB，篇幅精炼，便于随查随用；已有554人学习，适合作为课程作业、数据分析项目前置探索及统计函数速查材料。文档以mtcars等数据集为例，系统梳理summary()与sapply()的用法，并给出自定义偏度、峰度函数的实现，同时介绍Hmisc、psych包的describe()输出。分组统计部分讲解aggregate()、by()和doBy::summaryBy()的实现思路，可视化部分结合ggplot2演示直方图与箱线图的绘制，帮助读者把数值结果与图形分布对应起来。整体内容兼顾基础操作与进阶包应用，可作为R语言描述性统计的案头参考。

1. 描述性统计分析在 R 语言里：四条命令跑通再往下谈

接手一份陌生的数据表，最先要回答的问题不是该用哪个模型，而是每一列长什么样：中心在哪、离散多大、有没有离群点、缺失值占不占坑。R 语言在描述性统计分析上的工具链是这个环节里最顺手的一档——summary()一行给出四分位数，sapply()把任意函数批量套到每一列，psych 与 Hmisc 的describe()补齐偏度和峰度，ggplot2 再把分布画出来和数字对照。数据科学项目里常说的 EDA，第一层砖就是这套描述性统计。刚接手数据清洗、要给团队交一份字段体检报告的从业者，或者想把零散几行 R 代码沉淀成可复用模板的人，往下读会有收获。下面从基础函数开始，一路走到分组统计和图形核对。

2. 基础函数层：summary() 与 sapply() 的输出结构怎么读

这一章解决的是「跑得出结果但看不懂结果」的问题。很多人第一次用summary()和sapply()就卡在输出形态上——summary()打出来的是一张看着像表格、实际上是特殊类对象的玩意，sapply()有时候返回矩阵有时候返回列表，接write.csv()就报错。先把返回结构的类型、维度、行名列名摸清，后面所有扩展包的分析结果才能稳定拼成可交付的表。

2.1 summary() 返回的不是 data.frame，而是 table

先跑个最小例子，看清它到底返回什么：

# 加载内置数据集，只用三列做演示 data(mtcars) myvars <- c("mpg", "hp", "wt") s <- summary(mtcars[myvars]) # 关键：看清类型和维度 class(s) # "table" dim(s) # 6 3 rownames(s) # Min, 1st Qu, Median, Mean, 3rd Qu, Max colnames(s) # mpg, hp, wt

summary.data.frame()会对每一列分别调用summary.default()，如果所有列都能凑成同样长度的结果，它就拼成一个 6 行 N 列的table对象；只要有一列返回结构不同，它就退化成 list 打印。这解释了为什么有的数据框summary()出来整整齐齐，加一列字符型之后就变得七零八落。

行序固定是六个数：最小值、下四分位、中位数、均值、上四分位、最大值。想单独取某一行做后续计算，直接s["Mean", ]或as.numeric(s["Mean", ])就行，不必再逐列mean()。列名的顺序跟mtcars[myvars]传入的列顺序一致，所以选择变量的向量最好显式写出来，别依赖列位置。

提示：summary()对因子列返回的是各水平频数，对字符列在 R 4.0 之后只会给Length / Class / Mode三个值，看不到频数。要频数得先factor()化再summary()，或者干脆table()。

参数层面summary()没什么可调的，它的价值就在「一眼看全」。但它有个常被忽略的行为：默认不做任何 NA 过滤，mean那一行在含缺失值的列上会直接给出NA。所以看到均值是NA，第一反应应该是回去查缺失值数量，而不是怀疑函数坏了。

2.2 sapply() 批量套函数：名字对齐是隐形的坑

sapply()的定位是把一个函数套到列表或向量的每个元素上，并尽量把结果简化成向量或矩阵。描述性统计里最常见的写法是：

# 单一统计量：返回命名向量，名字来自列名 sapply(mtcars[myvars], mean) # 多个统计量：函数返回命名向量，sapply 拼成矩阵 sapply(mtcars[myvars], function(x) { c(mean = mean(x), sd = sd(x), median = median(x)) })

第二种写法返回的是一个 3 行 3 列的矩阵，行是统计量名，列是变量名。这个方向容易记反，写报告时把行列搞混是新手常见翻车点。记住口诀：sapply沿列方向遍历变量，所以每个函数的返回值在行方向堆叠。

真正的坑在函数返回值的长度上。如果某个变量全是 NA，mean()返回NaN，长度还是 1，不影响；但如果函数内部有条件分支导致返回长度不一致，sapply会静默退化成 list，后续as.data.frame()直接报错。稳妥一点的替代是vapply()，强制指定返回类型和长度：

# vapply 要求显式声明返回值形状，长度不对立刻报错 vapply(mtcars[myvars], function(x) c(mean(x), sd(x)), FUN.VALUE = numeric(2))

FUN.VALUE是一个模板，告诉 R「我这个函数一定返回 2 个 double」。一旦某个变量触发了不同分支，R 会在那一列报values must be length 2，定位速度比sapply的沉默退化快得多。我一般写生产脚本用vapply，写一次性探索用sapply，这个习惯是从一次线上报告里某列莫名变成字符串之后养成的。

另一个细节是na.rm。sapply(mtcars[myvars], mean)默认na.rm = FALSE，只要有一个 NA 就整列返回 NA。批量计算时应该显式传参：

sapply(mtcars[myvars], mean, na.rm = TRUE) sapply(mtcars[myvars], sd, na.rm = TRUE)

2.3 把宽表结果整理成可交付的数据框

矩阵形态适合人眼看，不适合写进 CSV 交出去。一个合格的文件清单应该是一行一个变量、一列一个统计量。下面这段把前面几次调用合成一张表：

# 单变量统计量集合，返回命名向量，长度固定为 8 stats_one <- function(x) { x <- x[!is.na(x)] # 先剔除缺失，避免污染 c(n = length(x), mean = mean(x), sd = sd(x), min = min(x), q25 = unname(quantile(x, 0.25)), # unname 去掉分位数的名字 median = median(x), q75 = unname(quantile(x, 0.75)), max = max(x)) } # 对每列应用，得到 8 行 3 列的矩阵 m <- sapply(mtcars[myvars], stats_one) # 转置成 3 行 8 列，再挂上变量名列 out <- as.data.frame(t(m)) out$variable <- rownames(out) rownames(out) <- NULL out[, c("variable", setdiff(names(out), "variable"))]

三个细节值得说清：

• quantile()默认返回带25%这种名字的向量，直接塞进c()会让整行名字变成q25.25%这类拼接结果，所以用unname()剥掉。
• t(m)转置后行名是原来的列名，正好可以拿来当变量名列，但as.data.frame()之后这些行名会挂在rownames上，得手动提升成列。
• out$variable放到第一列纯粹为了阅读顺序，setdiff那一步保证不会重复。

跑完之后得到的是一张可以直接write.csv(out, row.names = FALSE)的表，后面接 psych 或 Hmisc 的输出也沿用同样的「转置 + 提行名」套路即可。这套整理逻辑一旦固定下来，整个描述性统计环节的产出格式就统一了。

3. 自定义统计函数：偏度、峰度公式实现与 na 处理的两个边界

内置的mean、sd、median覆盖了中心趋势和离散程度，但分布形状的指标——偏度和峰度——不在基础包里。要么装扩展包，要么自己写。自己写的好处是完全可控、不引入依赖、能塞进任何环境；坏处是公式细节一旦写歪，结果就会系统性偏离，而且不报错、不报警，只是悄悄错着。这一章把公式实现、参数命名、缺失值边界三件事拆开讲。

3.1 偏度与峰度公式：分母用 n 还是 n-1

先摆出原始版本：

# 基础版：矩法偏度和超额峰度 mystats <- function(x, na.omit = FALSE) { if (na.omit) x <- x[!is.na(x)] m <- mean(x) n <- length(x) s <- sd(x) # 样本标准差，分母 n-1 skew <- sum((x - m)^3 / s^3) / n # 三阶矩 kurt <- sum((x - m)^4 / s^4) / n - 3 # 四阶矩减 3 c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt) }

这段代码能跑，但里面混用了两套分母：s用的是n-1无偏估计，而三阶、四阶矩求和之后除以的是n。这在统计上叫「矩法估计」，严格来说和教科书里样本偏度g1、样本峰度g2不是同一个量。R 的e1071::skewness()默认用的是带校正的G1 = n/((n-1)(n-2)) * sum((x-m)^3) / s^3，样本量小的时候两者能差出百分之十几。

这意味着：如果你拿mystats()算的偏度去和 SPSS、SAS 或者psych::describe()的输出对照，数字对不上属于正常，不是谁算错了，是口径不同。要做交叉验证，得先把口径对齐——要么都用手写公式，要么都用同一个包的默认值。

# 交叉验证：手写值 vs e1071 默认值 x <- mtcars$mpg mystats(x)["skew"] e1071::skewness(x) # 带校正，样本量小时差异明显 e1071::skewness(x, type = 2) # type = 2 是未校正版本，更接近手写

e1071::skewness的type参数有 1、2、3 三档，type = 3是默认的G1，type = 2就是除以n的矩法。写模板函数时把口径注释清楚，比追求「哪个是标准答案」更实际。

峰度那一行同样要注意，-3是把正态分布的峰度归零，得到的是「超额峰度」。如果下游报告要的是原始峰度，去掉-3，别让读者自己加减。

3.2 参数叫 na.omit 会盖住同名函数

上面那段mystats()的参数名叫na.omit，一旦函数内部需要调用stats::na.omit(x)，就会因为参数遮蔽而报could not find function，或者更可怕——把参数本身当成函数调，报一个看不懂的错。这是抄网上代码时最容易带进来的地雷。

# 改掉参数名，并加长度下限检查 mystats_safe <- function(x, remove_na = TRUE) { if (remove_na) x <- x[!is.na(x)] if (length(x) < 3) { # 样本太短，四阶矩没有意义，返回结构保持对齐 return(c(n = length(x), mean = NA_real_, stdev = NA_real_, skew = NA_real_, kurtosis = NA_real_)) } m <- mean(x); n <- length(x); s <- sd(x) c(n = n, mean = m, stdev = s, skew = sum((x - m)^3 / s^3) / n, kurtosis = sum((x - m)^4 / s^4) / n - 3) }

两点改进都很关键：

• 参数改名为remove_na，避免和na.omit()撞车。
• 加length(x) < 3的短路返回。sd()在长度 1 时返回NA，s^3就是NA，整个偏度变成NA，不报错，直接污染下游。提前返回一个长度固定、名字固定、值全NA的向量，sapply拼矩阵时才不会因为长度不一致退化成 list。

下面这句能说明为什么长度对齐这么重要：

# 假设一列全是 NA，直接算会破坏 sapply 的形状 bad <- rep(NA_real_, 10) sapply(list(a = mtcars$mpg, b = bad), mystats) # 行列可能错位 sapply(list(a = mtcars$mpg, b = bad), mystats_safe) # 稳定返回 5 行 2 列

3.3 封装成可复用的字段体检模板

零散函数凑在一起是没法复用的，把它包成一个接受数据框、返回数据框的函数，团队里谁拿到都能跑：

# 输入数据框和变量名向量，输出一行一变量的描述统计表 desc_report <- function(df, vars = names(df)) { num_vars <- vars[vapply(df[vars], is.numeric, logical(1))] if (length(num_vars) == 0) stop("没有可分析的数值列") m <- vapply(df[num_vars], mystats_safe, FUN.VALUE = numeric(5)) out <- as.data.frame(t(m)) out$variable <- rownames(out) rownames(out) <- NULL out[, c("variable", "n", "mean", "stdev", "skew", "kurtosis")] } desc_report(mtcars, c("mpg", "hp", "wt", "cyl"))

里面有两个工程化的取舍：用vapply替掉sapply，让任何一列返回值不符合模板时立刻报错，而不是静默变 list；用is.numeric先过滤列类型，避免字符型或因子列进入公式计算时抛出难以定位的报错。

跑出来cyl是 0 偏度、峰度约 -1.8 左右（这一列只有 4、6、8 三个值），看起来像个连续变量但其实是分组变量。这也是描述性统计阶段必须人工过一眼的原因——函数不会替你判断某一列该不该参与数值分析，只会老老实实把数字算出来。

4. 扩展包 describe()：Hmisc 与 psych 的选型对比和加载冲突

手写函数可控但需要维护，扩展包的describe()一次性给出均值、标准差、截尾均值、中位数、绝对中位差、极值、偏度、峰度。两个主流的包是 Hmisc 和 psych，都提供同名函数，但返回类型和适用场景差别很大，混着用很容易在打印环节翻车。

4.1 Hmisc::describe() 输出的是给人看的报告

install.packages("Hmisc") library(Hmisc) describe(mtcars[c("mpg", "hp", "wt")])

打印出来的格式是这样一组一组：

• 顶部：n、missing、distinct、Info（有效比例）、Mean、Gmd（基尼均差）、0.05 / 0.10 / 0.25 / 0.50 / 0.75 / 0.90 / 0.95分位、lowest / highest五个值。
• 下方：连续变量的直方图字符画，一眼看形状。

对初学者很友好，但它是describe类对象，print()时有专门方法，as.data.frame()不一定能直接拿。想抽数值得用Hmisc::describe(...)的结果里逐项取，或者干脆sapply手写。所以 Hmisc 更适合人眼粗看，不适合塞进自动化流水线。

另一个实际问题是加载成本。library(Hmisc)会 attach 一堆依赖包，而且会覆盖一些同名函数——summarize是重灾区，如果你后面写 dplyr 的summarise()切换环境时容易搞混。我通常只在交互式探索时加载它，生产脚本里写清楚Hmisc::describe()带命名空间调用。

4.2 psych::describe() 返回数据框，更适合流水线

install.packages("psych") library(psych) d <- psych::describe(mtcars[c("mpg", "hp", "wt")]) class(d) # "data.frame" names(d) # vars n mean sd median trimmed mad min max range skew kurtosis se

psych 的describe()返回标准数据框，一行一个变量，列名从上到下覆盖了 n、均值、标准差、中位数、截尾均值（默认 10%）、绝对中位差、极值、极差、偏度、峰度、标准误。这直接可以write.csv()或rbind()，跟前面手写的模板输出结构几乎能拼在一起，这是它最大的优势。

可调参数主要有三个：

用法示例：

psych::describe(mtcars[c("mpg", "hp", "wt")], trim = 0.2, quant = c(.25, .5, .75), IQR = TRUE)

trim = 0.2之后trimmed那列会变，跟mean对比一下就能看出尾部对均值的影响有多大——如果两者差得远，说明这一列有长尾，后面建模时可能要做变换。这是 psych 比 Hmisc 更合适的地方：它给你的是可比较的数字列，不是一堆格式化文本。

4.3 两个包同时加载时的输出冲突

Hmisc::describe()和psych::describe()同名，谁后library()谁生效。这意味着你在一份脚本里先library(psych)再library(Hmisc)，后面写describe(df)用的就是 Hmisc 版本，返回类型从 data.frame 变成describe对象，$mean之类取不到，直接翻车。

三条处理方式，按稳妥程度排序：

1. 全程加命名空间前缀：Hmisc::describe()、psych::describe()。一行多敲几个字符，换来永远可复现，这是我现在的默认写法。
2. 只library()一个，另一个全部前缀调用。脚本头部显式library(psych)，需要 Hmisc 的地方写全名。
3. 用conflicts()自检：脚本开头跑一次conflicts()，把自己用的包全部列进去，看到describe出现两次就知道要处理。

# 自检：看一下当前环境里有哪些同名遮蔽 library(psych) library(Hmisc) conflicts() # 输出里 describe 会出现在 "Hmisc" 与 "psych" 两栏 # 稳妥写法：全程前缀 h <- Hmisc::describe(mtcars$mpg) p <- psych::describe(mtcars$mpg)

还有一个不太容易注意的点：psych 的describe()默认不接受字符列，遇到非数值会直接报错。传进去之前应该先过滤：

num_df <- mtcars[, vapply(mtcars, is.numeric, logical(1))] psych::describe(num_df)

这一句其实是第 2 章vapply套路的复用——用is.numeric做列筛选，把类型问题挡在函数入口外。以后接 csv、接数据库导出的表，这行前置过滤可以一直挂着。

5. 分组描述性统计的常见问题排查：aggregate、by、summaryBy 五处翻车

按某一列分组、对每个分组的数值列做统计，这是业务场景里出现频率最高的需求。R 里至少有三条主流路线：aggregate()、by()、doBy::summaryBy()。三条路都能出结果，但输出结构、NA 处理策略、多列多返回值的展开方式都不一样，选错一个后面就要花半天整理格式。这一章把三条路的用法和五个最常踩的坑摆在一起。

5.1 aggregate() 的公式接口与多返回值展开规则

# 单列单统计量：返回干净的 3 行 2 列 aggregate(mpg ~ cyl, data = mtcars, mean) # 多列单统计量：cbind 把多个 y 捆起来 aggregate(cbind(mpg, hp, wt) ~ cyl, data = mtcars, mean) # 多列多统计量：FUN 返回命名向量 aggregate(cbind(mpg, hp) ~ cyl, data = mtcars, FUN = function(x) c(m = mean(x), s = sd(x)))

前两种结果规整，第三种列名会被拼成mpg.m、mpg.s这种形式，字段一多就难以辨认。这不是 bug，是aggregate()把 FUN 的返回值展平之后的自然结果。如果统计量固定，直接用FUN = length、FUN = mean分多次调用再merge()，比一次返回多值更好整理。

aggregate()还有一个隐蔽行为：na.action默认是na.omit，丢的是整行，不是整列。分组变量里出现 NA 时，那一整组数据会被静默丢弃，输出少一组，你根本无从察觉。数据干净时看不出问题，一旦上游有缺失，结果就是错的。

# 让 NA 单独成组，而不是被静默丢弃 aggregate(mpg ~ cyl, data = mtcars, mean, na.action = na.pass)

5.2 by() 返回 list，后续拼接要显式处理形状

by()的语义是「按因子把数据框切成子集，对每个子集跑同一段函数」，返回值是一个by类的 list，打印时自带分组标题，人眼舒服，但接rbind()或write.csv()得手动展开。

dstats <- function(x) sapply(x, mystats_safe) r <- by(mtcars[c("mpg", "hp", "wt")], mtcars$cyl, dstats) class(r) # "by" length(r) # 3 r[[1]] # 第一组：5 行 3 列（统计量 × 变量）

r[[1]]是一个 5 行 3 列的矩阵，行是n/mean/stdev/skew/kurtosis，列是变量名。要拼回一张宽表，得把每一组的矩阵转置、加分组列，再rbind()：

long <- do.call(rbind, lapply(names(r), function(g) { m <- r[[g]] data.frame(cyl = g, variable = colnames(m), as.data.frame(t(m)), row.names = NULL, check.names = FALSE) })) head(long)

这里check.names = FALSE是为了让列名保持原样，不至于被 R 自动转成kurtosis.1这类合法变量名。row.names = NULL则是防止不同分组的行号重复。

5.3 doBy::summaryBy() 与自定义函数的组合

install.packages("doBy") library(doBy) summaryBy(mpg + hp + wt ~ cyl, data = mtcars, FUN = mystats_safe)

结果是一张 data.frame，一行一个分组，列名形如mpg.n、mpg.mean、mpg.skew。这是三条路里输出最规整的，直接可以写盘。两点需要注意：

• FUN 如果返回命名向量，summaryBy()会用「列名.函数名」的规则拼列名；FUN 如果返回未命名向量，则会按mpg.1、mpg.2顺序编号，可读性差。
• FUN 也可以传列表：FUN = list(mean = mean, sd = sd, median = median)，此时输出列名是mpg.mean、mpg.sd这一套，比手写命名向量更省心。但列表形式要求每个函数都返回标量，不适合偏度这种需要共用中间变量的计算。

所以我的选择习惯是：单统计量用 aggregate，多统计量且需要干净输出用 summaryBy，需要自定义组合逻辑用 by。三者可以混用，但输出结构必须先统一，否则后面merge()就是灾难。

5.4 五个最常见翻车记录

现象一：summary(df)里字符列只看到 Length / Class / Mode，没有频数。原因：R 4.0 之后stringsAsFactors默认FALSE，字符列不是因子，summary()按字符处理。 解决：summary(factor(df$col))或直接table(df$col)；数据读入时如果想保留因子语义，用read.csv(..., stringsAsFactors = TRUE)。

现象二：sapply(df, mean)报参数不是数值也不是逻辑值。原因：数据框里混了字符列或因子列，mean()在非数值上直接抛错。 解决：df[, vapply(df, is.numeric, logical(1))]先筛列，再做批量统计。别用位置索引，列顺序一变就全乱。

现象三：aggregate()输出的行数比分组数少。原因：分组变量里有 NA，默认na.action = na.omit把整组丢了。 解决：加上na.action = na.pass，或在建因子时用addNA(cyl, ifany = TRUE)把 NA 显式当作一个水平。

现象四：手写偏度和psych::describe()算的对不上，差了一个量级。原因：手写公式里标准化用的是sd()（分母n-1），除以的是n；psych 默认使用带校正的样本偏度，样本少时差异明显。 解决：固定口径，用e1071::skewness(x, type = 2)验证手写值，或者全程只用一端的结果，绝不混用。

现象五：by()之后do.call(rbind, ...)报numbers of columns of arguments do not match。原因：某一组只有 1 行或 2 行，sapply在退化成向量时行数不足，函数返回的形状和其他组不一致。 解决：mystats_safe里做长度下限判断，统一返回 5 行固定结构；实在没法保证时改成vapply加FUN.VALUE，形状不对立刻报错，而不是等到rbind()才炸。

注意：这五条里有三条是「不报错、结果错」，比报错更需要警惕。分组统计交付之前，先拿table(df$group)核对一遍分组计数，再拿nrow(result)和length(unique(df$group))对一下，能提前拦掉大部分问题。

6. ggplot2 可视化核对：直方图、箱线图与统计量的交叉验证

数值统计表只能告诉你「均值是多少、四分位在哪」，看不出分布是双峰还是长尾。可视化这一步的真正目的不是出图好看，是拿眼睛核对刚才那张表——直方图的分箱峰位应该和均值、中位数吻合，箱线图的箱体和须长应该和四分位、极值对得上。不一致的地方就是需要回去复查的地方。

6.1 直方图分箱与 summary() 四分位的对照

library(ggplot2) # binwidth = 1 让每箱宽度 1，boundary = 0 让分箱从整数开始 ggplot(mtcars, aes(x = mpg)) + geom_histogram(binwidth = 1, boundary = 0, fill = "steelblue", color = "black") + geom_vline(xintercept = quantile(mtcars$mpg, c(.25, .5, .75)), linetype = "dashed") + labs(x = "mpg", y = "count")

三条虚线分别是 Q1、中位数、Q3。看图时盯三件事：

• 中线是否落在分布的主体上。如果中位数线偏在一侧，说明尾巴长，后面建模要考虑稳健统计量。
• Q1 和 Q3 之间的箱体宽度和直方图里的「主体峰」宽度是否接近。差距太大，说明分箱粒度选得不对，binwidth需要调整。
• 有没有明显脱离主体的孤立柱子。这是离群点的第一层信号，和boxplot.stats()$out的结果要能对上。

boundary = 0常被省掉，但它决定了分箱的起点：不指定时 ggplot2 从数据最小值开始分，同一个变量换一份数据后分箱位置就会变，图形没法跨数据集对比。指定整数边界能让多次出图保持口径一致。

6.2 箱线图与 aggregate() 结果的互证

ggplot(mtcars, aes(x = factor(cyl), y = mpg)) + geom_boxplot(outlier.colour = "red", outlier.size = 2) + stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "darkgreen")

红色点是离群值，深绿色星号是分组均值。三件事可以立刻验证：

• 箱体上下沿是否和aggregate(mpg ~ cyl, mtcars, quantile)算出的 Q1、Q3 一致。R 的箱线图用的是 hinges 算法，小样本上跟quantile(type = 7)会差一点，差一条线的宽度内属于正常，差得多要回头查分组。
• 上下须的端点是否落在Q1 - 1.5*IQR和Q3 + 1.5*IQR之内。落在外面又被画成点的，就是被判定为离群值。
• 绿色星号（均值）和中位数线的水平位置是否错开。错开越大，偏度越明显，可以回第 3 章手写的skew值上对一下方向。

# 与图形一一对应的数值 aggregate(mpg ~ cyl, data = mtcars, FUN = function(x) c(q1 = unname(quantile(x, .25)), med = median(x), q3 = unname(quantile(x, .75)), mean = mean(x))) boxplot.stats(mtcars$mpg[mtcars$cyl == 4])$out # 离群点原始值

6.3 一页核对表：数字和图形怎么互相印证

把整条流水线的检查点固定成一张表，每次交付之前照着过一遍，比凭记忆靠谱：

这张表的用法不是逐一打勾，而是「有任何一格对不上就停」。描述性统计环节的错误会一路传到建模和结论，越早拦越省事。

我吃过一次亏：一份分组报告里 Q3 算出来 78.5，箱线图箱体上沿明显更高，当时以为是 hinge 差异没管，交付后才发现是分组变量里有个前导空格，被当成了独立一组，数量少但均值高，硬生生把上沿抬上去。从那以后我每次跑完分组统计，都强制走一遍table()对分组计数、aggregate()对分位、ggplot()看图的三步核对，任何一步对不上就先停下查数据而不是改代码。这套习惯救回来的返工时间，远比每次多花的十分钟多。希望帮到你。

本文还有配套的精品资源，点击获取