简介面向R语言初学者与中级数据分析学习者的描述性统计分析实操指南聚焦如何用基础函数与扩展包快速完成数据特征探查与结果呈现。内容覆盖 summary() 的基础汇总、sapply() 批量计算均值与标准差、自定义函数求偏度与峰度以及 Hmisc、psych 两个包中 describe() 输出的详尽统计量分组场景给出 aggregate()、by() 与 doBy 包 summaryBy() 三种实现路径可视化部分则借 ggplot2 绘制直方图与箱线图用于观察分布形态与异常值形成从数值汇总到图形呈现的完整链路。整份资料为1个PDF文件约98KB篇幅紧凑适合随查随用。已有554人学习下载可作为课程作业、数据探索报告或正式分析前的例行检查参考帮助读者建立规范的描述性统计流程少走函数挑选与分组聚合的弯路。1. 从一份 CSV 到一张统计底表R 描述性统计先解决哪三个问题拿到一份几十列的表第一反应不该是建模而是把每列的分布摸一遍。R 语言里的描述性统计分析正好干这件事中心趋势、离散程度、分布形状三组指标配上分组分析和 ggplot2 的一眼可视半小时内就能判断哪些列是废列、哪些列在分组之后明显漂移。这套流程适合已经会用read.csv()读数据、但每次出统计表还得手敲mean()、sd()的从业者。真正麻烦的从来不是“算不算得出来”而是算完之后怎么拼成一张能交付的表。summary()吐的是一段文本sapply()给的是矩阵describe()返回的是对象aggregate()一遇到多函数就变成矩阵列——每种返回结构背后都对应一次翻车。下面从最基础的函数一路推到分组和可视化重点放在返回结构、参数边界和排错上。2. summary 与 sapply把散落的统计量拼成一张可交付的表2.1 summary() 的输出不是数据框先认清它的类型summary()是最顺手的入口但它对data.frame的返回值并不是你以为的那种表格。对数值列给最小值、四分位数、均值、最大值对因子列给频数——看起来像表实际是一个table类的矩阵对象。直接用rbind()往报告里塞或者丢进write.csv()十有八九会发现列名全乱。data(mtcars) myvars - c(mpg, hp, wt) s - summary(mtcars[myvars]) class(s) # table dim(s) # 6 3行是统计量列是变量 # 转成规整的数据框列是变量行是统计量 m - as.data.frame.matrix(as.matrix(s)) m$stat - rownames(m) m逻辑上分两步走先as.matrix()把 table 还原成纯矩阵再as.data.frame.matrix()保持列结构转成数据框最后把行名单独存成一列stat这样后面和别的统计量合并才不会丢索引。参数上要注意summary()自带的两个开关digits控制小数位默认跟着getOption(digits)走通常是 3maxsum控制因子水平最多展示多少个默认 100水平一多就只给个计数摘要。这两个参数只影响打印结果不影响返回对象里的内容所以想改精度得在转完之后自己round()。summary()的定位是快速体检不是出报表。它给不出标准差、偏度、峰度也不告诉你缺失值有多少。我一般拿它做第一眼筛查哪一列的最小值和最大值量级明显不对哪一列的均值和中位数差得很远先记下来再去细查。2.2 sapply 加自定义函数一次拿到 n、mean、sd、偏度、峰度summary()缺的那几项用sapply()逐列套函数就能补齐。内置的mean()、sd()、var()、median()、quantile()都是标量返回套一层sapply()直接得到命名向量。但偏度和峰度没有基础函数必须自己写。mystats - function(x, na.omit TRUE, digits 3) { if (na.omit) x - x[!is.na(x)] n - length(x) # 样本量不足两个标准差无意义直接返回 NA 占位 if (n 2) return(c(n n, mean NA, sd NA, skew NA, kurt NA)) m - mean(x) s - sd(x) # 常量列 sd 为 0除以 s^3 会得到 Inf/NaN必须先拦 if (s .Machine$double.eps) return(c(n n, mean m, sd 0, skew NA, kurt NA)) skew - sum((x - m)^3) / (n * s^3) # 总体偏度正态为 0 kurt - sum((x - m)^4) / (n * s^4) - 3 # 超额峰度正态为 0 round(c(n n, mean m, sd s, skew skew, kurt kurt), digits) } res - sapply(mtcars[myvars], mystats) res # 行是统计量列是变量 t(res) # 转置成“一行一个变量”更适合出表这段代码有几个地方值得说清。na.omit默认给了TRUE因为只要列里混进一个NAsum((x - m)^3)整项就变NA后面的偏度峰度全废sd()的NA传播同理。s .Machine$double.eps这个判断是拿浮点数当零来比比写s 0稳wt这类精度高的列不会因为尾数误差漏判。偏度用的三阶矩除以s^3是总体口径如果需要样本无偏估计那就得换e1071::skewness()或psych::describe()两者给的数会差一点报告里必须注明用的是哪套公式否则复核的人一算对不上就麻烦了。sapply()的返回结构取决于函数返回值返回长度为 5 的命名向量结果就是 5 行 N 列的矩阵行名就是n、mean、sd、skew、kurt。所以t()之后才是常见的“变量在行、指标在列”。这一步转置不做后面write.csv()出来的表是横的业务方看不懂。2.3 类型和缺失值的预处理顺序比选哪个函数更影响结果函数选错顶多重算一次类型判错会静默出错。最常见的是数值列被读成字符带千分位逗号、带百分号、带单位后缀、夹杂N/A或-作为缺失标记。这类列sapply(..., mean)会直接报argument is not numeric or logical或者更糟——被自动转成因子后mean()在旧版 R 里对它返回NA加一个 warning很容易被刷屏的日志盖过去。probe - function(df) { data.frame( var names(df), cls sapply(df, function(x) class(x)[1]), n_na sapply(df, function(x) sum(is.na(x))), n_uniq sapply(df, function(x) length(unique(x))), row.names NULL, stringsAsFactors FALSE ) } probe(mtcars)cls那一列是判断要不要动手的第一依据出现character或factor而业务上应该是数值的一律先转换再统计。n_na用来交叉验证sum(is.na(x))会把NaN也算进去但空字符串不算——如果一列里n_na是 0 而n_uniq里混着一个空串那就是读取时没设na.strings。n_uniq等于行数的列基本是主键或流水号这种列算均值毫无意义出表前直接剔掉。读入阶段的写法比事后修补省事得多raw - read.csv(data.csv, na.strings c(, NA, N/A, -, null), stringsAsFactors FALSE, check.names FALSE) raw[] - lapply(raw, function(x) if (is.character(x)) type.convert(x, as.is TRUE) else x)na.strings一次列全常见的缺失标记比读进来再gsub干净stringsAsFactors FALSE保证字符列不会被悄悄转成因子type.convert()负责把“看起来是数字的字符”转回数值转不动的保持原样这时候再用probe()看一遍就能定位剩余问题列。check.names FALSE是为了保住中文列名和带空格的列名不被改写成mpg.x这种形式代价是后续用列名索引时得加反引号。3. Hmisc 与 psych 的两个 describe()字段、参数与转数据框的差别3.1 Hmisc::describe() 的输出是给人的不是给程序的Hmisc包的describe()信息最全样本量、缺失数、唯一值数、均值、五个分位点、标准差、偏度、峰度还额外给出最小和最大的几个实际观测值排查异常值的时候特别省事。但它的返回对象是describe类print()出来是分节的可读文本想拿它做后续处理必须自己抽。library(Hmisc) d - describe(mtcars[myvars]) class(d) # describe length(d) # 3每个变量一个元素 str(d[[1]]) # 先看清楚字段名再取值关键就在str(d[[1]])这一步。describe对象本质是 list每个变量对应一个子 list里面装着n、missing、distinct、mean、sd、各分位点、skewness、kurtosis等字段。字段名在不同版本之间有过调整所以别凭记忆硬写x$skew先str()确认一遍再写抽取代码。res - do.call(rbind, lapply(names(d), function(v) { x - d[[v]] data.frame( var v, n x$n, missing x$missing, mean x$mean, sd x$sd, skew x$skewness, kurt x$kurtosis, row.names NULL ) }))do.call(rbind, lapply(...))这套组合是处理“每个元素结构相同但返回 list”的标准解法比for循环加rbind快得多也不会出现行名冲突。参数方面describe()支持desc、omit、trans等几个开关omit TRUE会跳过缺失值统计适合数据已经清洗干净的场景trans可以传一个变换函数比如log让分位点在变换后的尺度上展示。注意Hmisc加载时会连带挂上一批依赖包和dplyr同时加载容易撞函数名summarize、src之类项目里我一般只在需要出详细体检表的那段代码里library()用完不卸载也没关系但要在脚本开头写清楚加载顺序。3.2 psych::describe() 直接给矩阵参数更值得调psych的describe()走的是另一条路返回现成的data.frame一行一个变量列是n、mean、sd、median、trimmed、mad、min、max、range、skew、kurtosis、se。省掉抽取这一步代价是默认不带分位点需要显式打开。library(psych) desc - describe(mtcars[myvars], trim 0.1, # 计算 10% 截尾均值输出 trimmed 列 skew TRUE, # 保留偏度 ranges TRUE, # 输出 min/max/range IQR TRUE, # 输出四分位距 quant c(.05, .25, .5, .75, .95)) write.csv(desc, desc_psych.csv, row.names TRUE)参数里最该说清的是trim。它算的是截尾均值先把两端各 10% 的观测砍掉再求平均。均值和中位数打架、又怀疑尾部有异常值时trimmed和mean差得越大说明尾部拉动力越强这个信号比单看偏度直观。IQR TRUE输出的四分位距是箱线图那套异常值判定的基础后面第 6 章直接用得上。quant是自定义分位点向量业务上关心 P95 响应时间就把它写进去别每次都事后quantile()补。fast参数值得单独提一句数据量上到千万行级别时describe()默认还会老老实实算偏度峰度耗时会明显上来fast TRUE会跳过这两项只留基础指标速度差出好几倍。日常几万行以内的表没必要动它。还有一个容易忽略的点psych的describe()对非数值列不是报错而是直接跳过。所以一列明明在数据里输出表里却找不到八成是类型问题回到 2.3 的probe()查一遍就行。3.3 因子列、逻辑列、日期列不会自动进表两个describe()都只对数值型负责。因子列的频数、逻辑列的 TRUE/FALSE 占比、日期列的时间跨度都得另起一段统计否则报告里就只剩半张表。freq_tab - function(x, name deparse(substitute(x))) { x - as.factor(x) t - table(x, useNA ifany) # ifany有缺失才显示 NA 行 data.frame( var name, level names(t), n as.integer(t), pct round(100 * as.numeric(t) / sum(t), 1), row.names NULL ) } freq_tab(factor(mtcars$cyl, levels c(4, 6, 8), labels c(四缸, 六缸, 八缸)))useNA ifany比默认的no有用得多缺失本身经常就是业务信号——某个渠道的订单来源字段大面积空着往往意味着埋点漏了。pct分母用的是sum(t)也就是含缺失的总数这样各水平占比加起来是 100%如果想让占比只在有效值里算把分母换成sum(t[!is.na(names(t))])即可两种口径在报告里必须写明白不然业务方对不上账。日期列我一般单独抽三个数最早、最晚、跨度天数再加一句“缺失率”。跨度直接决定了这份数据能不能支撑同比分析比任何分布指标都重要。4. 分组统计三条路aggregate、by、summaryBy 的返回结构对比4.1 aggregate() 的多函数写法会把列变成矩阵aggregate()的公式接口最像 SQL 的GROUP BY一行就能出结果aggregate(mpg ~ cyl, data mtcars, FUN mean)但它一次只能套一个函数。想同时要样本量、均值、标准差就得把 FUN 换成返回命名向量的匿名函数这时候返回结构就开始变形了。agg - aggregate(mtcars[myvars], by list(cyl mtcars$cyl), FUN function(x) c(n length(x), mean mean(x), sd sd(x))) str(agg) # mpg、hp、wt 三列的类型是 matrix不是 numericstr()会明确告诉你mpg这一列是个 3 列的矩阵。这就是最常见的翻车点表在控制台打印出来长得挺正常一write.csv()或者丢进ggplot()就报错因为下游拿到的是矩阵列而不是扁平列。必须再摊一层flat - do.call(data.frame, agg) names(flat) # cyl, mpg.n, mpg.mean, mpg.sd, hp.n, hp.mean, ...do.call(data.frame, ...)会把每个矩阵列拆成独立的向量列名按“变量名.统计量名”拼接。这套命名规则省事也坑人列名里带点后续用$取值时 R 会把它当成 S3 方法分派的一部分得写成flat$mpg.mean或者flat[[mpg.mean]]。列一多我宁愿重命名一遍。aggregate()还有一个必须记住的默认行为na.action na.omit是整行删除。也就是说只要你按cyl分组、而hp列里有缺失这一整行会被丢掉结果是mpg的均值也跟着少算了几条观测。解决办法是在 FUN 内部显式na.rm TRUE别指望外层参数aggregate(mtcars[myvars], by list(cyl mtcars$cyl), FUN function(x) c(n sum(!is.na(x)), mean mean(x, na.rm TRUE)))n用sum(!is.na(x))而不是length(x)这样表里的样本量是有效样本量和各组均值的实际计算基数对得上。这一点在出交付表时特别重要样本量对不上是业务方最先质疑的地方。4.2 by() 返回的是列表转表要自己动手by()的定位是“按组切开对每一块跑任意函数”灵活性最高代价是返回一个by类的 list不能直接出表。res - by(mtcars[myvars], mtcars$cyl, function(d) { m - t(sapply(d, mystats)) # 行变量列n/mean/sd/skew/kurt data.frame(var rownames(m), m, row.names NULL) }) out - do.call(rbind, lapply(names(res), function(g) { cbind(grp g, res[[g]], row.names NULL) })) rownames(out) - NULL结构上分两层内层对每个分组的数据框套sapply加t()得到“变量在行、指标在列”的小表外层用lapply遍历names(res)取出每个分组名cbind上去最后rbind合并。row.names NULL连着写两次是为了压掉 R 自动生成的1、1.1这种带后缀的行名否则合并之后行名会非常难看。需要留意的两个参数行为by()的INDICES会被强制转成因子分组顺序按因子水平排不是按原始出现顺序想让“高价值客户”排在“低价值客户”前面就得显式传factor(..., levels ...)另外分组键为NA的行会被直接丢弃而且不给任何提示。如果缺失分组本身有意义比如“未知渠道”得先把NA替换成一个明确的水平再跑。4.3 doBy::summaryBy() 的多函数写法与结果列名doBy包的summaryBy()走的是公式接口但比aggregate()更照顾“多个统计量并列”这个需求。library(doBy) sb - summaryBy(mpg hp wt ~ cyl, data mtcars, FUN function(x) c(n sum(!is.na(x)), mean mean(x, na.rm TRUE), sd sd(x, na.rm TRUE)), keep.names FALSE)返回的是扁平数据框列名形如mpg.n、mpg.mean、mpg.sd跟aggregate()摊平后的结果一致但不用再套do.call。keep.names FALSE是默认值保持“变量名.统计量名”的列名设成TRUE时列名只用统计量名多变量场景下会重名一般不推荐。分组键为NA的行同样会被丢掉处理方式和by()一样。三条路的选型我一般这么定只算一两个指标、要贴进现有脚本的用aggregate()分组内部要做多步计算、或者每组结果结构不一样用by()要一次性出多变量多指标的宽表用summaryBy()最省事。如果项目里已经引入了dplyrgroup_by()加summarise()的组合在列名控制和后续管道上更顺手这也是现在更常见的做法但换工具意味着所有老脚本都得动中小项目里我倾向于沿用已经在跑的那套别为了写法好看重写一遍分组逻辑。5. 避坑与排查数值悄悄变 NA 的五个现场5.1 现象summary() 里整列显示 Length/Mode均值返回 NA这是最典型的一类。控制台summary()输出里本该是数值的列显示成Length:32 Class:character Mode:character接着sapply(..., mean)直接报argument is not numeric or logical。原因几乎总是读取阶段列里混了N/A、-、空串、带千分位的1,234或者带单位的120kg。只要有一个非数字字符串read.csv()的自动类型推断就会把整列判成字符。解决上先别急着as.numeric()——直接转会得到一串NA加一条 warning原来的值全丢这是最亏的操作。正确顺序是先看unique(x[!grepl(^[0-9.]$, x)])捞出不规范的值确认它们代表什么是缺失还是单位再决定用na.strings重读、还是gsub()去掉单位后type.convert()。提示as.numeric(1,234)返回NA而不是 1234也不是报错。这种静默失败在几千行的表里几乎不可能靠肉眼发现。5.2 现象mean 有值但 sd 是 NA两个指标对不上同一列mean(x)给出了一个正常数字sd(x)却是NA。原因是mean()里写了na.rm TRUE而sd()里忘了写。只要列里有一个缺失sd()就返回NA然后由它算出来的偏度峰度全跟着变NA。解决方式不是逐处补参数而是把na.rm统一收口到自定义函数的一个开关上就是 2.2 里mystats()的na.omit参数全项目只走这一个入口。另外要区分NA和NaNNaN常常是0/0算出来的比如对全零列求标准差再相除is.na()会把两者都判为真但排查时得分开看is.nan()才是判NaN的那个。5.3 现象各分组的均值加权回去和总体均值对不上把aggregate()出来的各组均值按样本量加权平均结果跟直接对全表求mean()差了一点。差距不大所以最容易被放过。原因就是aggregate()默认na.action na.omit整行删除。某组里有一条观测的hp是缺失这条观测的mpg也被一起扔了该组的有效样本量和length(x)记的样本量不一致。解决FUN 内一律写na.rm TRUE并且样本量统计用sum(!is.na(x))而不是length(x)。出表之后做一次交叉校验——各组n之和应该等于该列的非缺失总数不相等就说明还有别的地方在偷偷删数据。5.4 现象偏度算出 NaN 或 ±Inf自定义函数跑完某一列的skew是NaN或者Inf而mean、sd都正常。sd正常说明这列不是完全常量但可能极度集中——绝大部分观测是同一个值只有零星几个不同。这种情况下s极小(x - m)^3 / s^3会直接溢出。另一个来源是全零列s 0除掉之后就是NaN。解决在函数里加if (s .Machine$double.eps) return(...)把常量列和近常量列拦住返回NA并在表里明确标出来。这类列的偏度本身就没有解释意义硬算出来的数字只会误导人。如果确实需要换e1071::skewness()这类带类型判断的实现它内部会给出NA而不是溢出值。5.5 现象describe() 结果表里少了几列或者报错找不到字段Hmisc::describe()跑完抽取代码报x$skew是NULL或者psych::describe()的输出里明明数据中有的列不见了。前者是版本差异describe对象的字段名在不同版本间调整过凭记忆写死的字段名一旦对不上就静默返回NULLdata.frame()里那一列就变空。解决是先str(d[[1]])打出真实结构再取值别照抄别人的代码。后者是类型问题psych::describe()对非数值列静默跳过不报错。回到probe()看cls那一列字符型和因子列先转换再统计或者用 3.3 的频数表单独出。6. ggplot2 收尾binwidth 定不准直方图会骗人数值表出完还得画两张图确认分布形状。ggplot2的geom_histogram()有个默认行为不指定binwidth时会给一条stat_bin() using bins 30的提示30 个分箱对 32 行的mtcars来说细得离谱画出来全是锯齿看起来像多峰实际只是样本太少。library(ggplot2) # Freedman-Diaconis 规则2 * IQR / n^(1/3) fd_bin - function(x) { x - x[!is.na(x)] 2 * IQR(x) / length(x)^(1/3) } bw - fd_bin(mtcars$mpg) # 约 4.6 ggplot(mtcars, aes(x mpg)) geom_histogram(binwidth bw, boundary 0, fill steelblue, colour white) geom_vline(xintercept mean(mtcars$mpg), linetype dashed) geom_vline(xintercept median(mtcars$mpg), linetype dotted) labs(x mpg, y count)binwidth用 FD 规则算出来的是数据驱动的值比默认的 30 箱靠谱得多当样本量少于 30 时更适合用 3 到 5 个分箱或者改用密度图叠加。boundary 0让分箱边界从 0 开始对齐否则同样的数据换个binwidth会因为起点偏移而出现“柱子左移一格”的假象。均值虚线和中位数虚线一叠偏态方向立刻就能看出来虚线离得越远尾部拉动越明显。箱线图看分组差异异常值判定直接复用前面的四分位距ggplot(mtcars, aes(x factor(cyl), y mpg)) geom_boxplot(width 0.5, outlier.shape 1, outlier.colour grey40) stat_summary(fun mean, geom point, shape 23, size 2.5, fill white) labs(x cyl, y mpg)outlier.shape 1用的是空心圆比默认的实心点更容易和stat_summary()叠上去的均值菱形区分开。箱线图默认按1.5 * IQR判定异常值也就是说超出上下四分位 1.5 倍四分位距的点会被单独画出来——这个阈值和psych::describe()里IQR TRUE输出的那列是同一套口径图和表能对上复核时省很多口舌。分组变量记得套factor()cyl是数值型不转因子的话 x 轴会变成连续刻度箱体会被挤成一团。血泪经验是从那以后我每次出图前都强制走一遍“binwidth 显式给定 分组变量转因子”这两步宁可多写一行fd_bin()也不接受一张看着挺顺眼、实际在骗人的直方图。图形是给别人看的第一印象一旦分布形状被分箱方式扭曲后面所有基于这张图的判断都跟着歪。希望帮到你。本文还有配套的精品资源点击获取