简介面向数据分析初学者与希望夯实 R 语言基础的中级读者这份 PDF 指南围绕描述性统计分析展开帮助读者理解数据的中心趋势、离散程度与分布形状并衔接基础统计量、分组比较与图形呈现等环节。包内共 1 个 PDF 文件大小约 98KB篇幅紧凑便于在电脑或移动设备上随时查阅与检索。内容覆盖 summary() 与 sapply() 等基础函数、自定义偏度峰度计算以及 Hmisc、psych 中 describe() 的扩展用法同时结合 aggregate()、by()、doBy 的 summaryBy() 说明分组统计思路并用 ggplot2 的直方图、箱线图展示分布与异常值。已有 554 人学习适合作为课程作业、数据探索和项目复盘的案头速查材料。从基础统计量到分组汇总再到图形表达形成连贯路径可据此整理自己的统计分析函数与图表模板减少在多个包与函数之间反复查找的成本。1. 一份陌生数据框到手为什么我建议先看 summary() 而不是直接建模接手一份别人整理好的数据最容易犯的错是直接lm()或者t.test()跑起来等结果不对劲再回头翻数据。我的习惯是先在 R 语言里把描述性统计这一步做扎实——用summary()扫一遍列的分布骨架用sapply()把均值方差批量拉出来再用手写函数把偏度和峰度补齐。这套流程的好处很直接量纲错位、缺失扎堆、常量列、离群点几乎都会在描述性统计里露头。这篇笔记拆的是 R 里那组高频函数和包summary()、sapply()、自定义mystats()、Hmisc与psych的describe()、分组用的aggregate()/by()/summaryBy()最后用ggplot2把统计量画回图上去对账。适合刚拿到数据、要在半天内交出一份可信体检报告的人。2. 三行代码立基线summary()、sapply() 与手写 mystats() 分别管什么2.1 先看清列类型和缺失str() 与 colSums(is.na())很多人跳过这一步直接算均值结果某个本来该是数值的列被读成了字符mean()报错或者返回NA再花两小时排查。数据类型是描述性统计的前置条件不是可选项。str()看结构colSums(is.na())看缺失密度两行就够。# 用内置数据集做样本换成你自己的 read.csv() 结果即可 data(mtcars) myvars - c(mpg, hp, wt) # 1) 看类型数值列才能进均值方差因子/字符列要单独处理 str(mtcars[myvars]) # 2) 看缺失密度整列缺失的变量后面算什么都不可信 colSums(is.na(mtcars[myvars]))str()的输出里num表示数值型chr是字符型Factor带水平数。这一步真正的价值在于提前发现两种脏数据一是本该是数字的列带了单位或逗号读成chr二是缺失比例高到没有任何统计意义。colSums(is.na(df))比is.na(df)好用因为前者返回每列一个整数一眼就能排序。参数上没有什么可调的唯一需要注意的是对超大表先做head()抽样看结构别在几百万行上反复str()。2.2 summary() 给的是分布骨架不是全部summary()是描述性统计的起点输出的是数值列的六个数最小值、下四分位、中位数、均值、上四分位、最大值因子列给频数。它便宜、快、不需要装包所以适合做第一轮扫描。# 只对关心的列做摘要避免几十列一起刷屏 summary(mtcars[myvars])在mtcars上跑出来mpg的均值 20.09、中位数 19.2均值被右侧几个高油耗车型往上拽了一点这种均值和中位数的偏离就是最廉价的偏态信号。hp的极差从 52 到 335跨度接近 6 倍说明这个变量在后续建模里大概率需要变换或者标准化。summary()有几个容易忽略的行为边界。第一它对缺失值默认忽略但会在列尾打一个NAs计数看到这个数字不为零就得决定是删还是补。第二从 R 4.0 起读入数据默认stringsAsFactors FALSE字符列在summary()里只给Length、Class、Mode三行不再给频数表如果你的报告需要类别分布得先factor()再用table()或者用janitor::tabyl()。第三summary()不给你标准差、偏度、峰度这三样要自己补。2.3 sapply() 批量跑单统计量FUN 后面能跟参数summary()不出的量用sapply()一行一个地批量算。它的本质是对数据框每一列调用一次函数再把结果简化成向量或矩阵。# 一列一个统计量返回命名向量名字就是列名 sapply(mtcars[myvars], mean) sapply(mtcars[myvars], sd) sapply(mtcars[myvars], var) # 分位数会返回矩阵行是分位数列是变量 sapply(mtcars[myvars], quantile, probs c(0.25, 0.50, 0.75)) # FUN 后面跟的参数会被原样传给函数本体 sapply(mtcars[myvars], mean, na.rm TRUE) # 一次算多个量用匿名函数返回命名向量 sapply(mtcars[myvars], function(v) c(n length(v), mean mean(v), sd sd(v)))参数上最关键的是...sapply(x, mean, na.rm TRUE)里的na.rm不是sapply的参数而是通过...透传给mean()的。这一点不理解就会写出sapply(x, mean, na.rm TRUE, trim 0.1)却不知道哪些参数归谁。第二个注意点是返回结构quantile()返回长度 3 的向量时sapply会把它拼成 3 行 3 列的矩阵而不是列表后续取数要用[, mpg]这种矩阵下标。第三个是sapply的简化行为有不确定性——如果某次调用返回长度不一致它会静默退化成 list你的下游代码就炸了。需要稳定返回结构时换成vapply(x, mean, numeric(1))把返回类型和长度写死出错当场报不会拖到后面。2.4 手写 mystats()一次拿 n、均值、标准差、偏度、峰度基础函数凑齐偏度和峰度要写一长串封装成函数更实际。核心公式是中心矩的标准写法。# 一次返回 n、均值、标准差、偏度、峰度超额峰度 # 参数 na.rm是否先剔除缺失这里刻意不叫 na.omit避免遮蔽 base::na.omit mystats - function(x, na.rm TRUE) { if (na.rm) x - x[!is.na(x)] n - length(x) m - mean(x) s - sd(x) # 样本标准差分母 n-1 skew - sum((x - m)^3 / s^3) / n # 三阶矩 / s^3 kurt - sum((x - m)^4 / s^4) / n - 3 # 四阶矩 / s^4再减 3 c(n n, mean m, stdev s, skew skew, kurtosis kurt) } # 每列一列结果转置后更好读 round(t(sapply(mtcars[myvars], mystats)), 3)逻辑说明先剔缺失再算矩否则sum()里任何一个NA会把整列结果污染成NA。n放在第一位是有意的——分组分析时样本量小于 30 的组偏度峰度的解释力很弱得先看n再决定要不要写进报告。参数说明na.rm默认TRUE是因为描述性统计的场景下一个缺失值不该让整列统计量报废如果要做严格复现把它设成FALSE让缺失直接暴露。函数名和参数名都可以改但别用na.omit当参数名——它会遮蔽stats::na.omit()函数体里一旦还要调用na.omit()就会报attempt to apply non-function这个错误信息不会提示你名字冲突只能靠经验认出来。round(t(...), 3)里的t()是必须的sapply返回的是 5 行统计量× 3 列变量的矩阵转置之后每行一个变量才是能直接贴进报告的形态。3. 偏度与峰度的口径之争手写公式、Hmisc 与 psych 三套结果怎么对齐3.1 原文公式属于哪一套口径偏度峰度没有唯一写法不同软件默认口径不同这是跨工具对不上数的头号原因。手写公式sum((x-m)^3/s^3)/n用样本标准差s做尺度、用n做分母对应的是三阶矩除以s^3峰度同理四阶矩除以s^4再减 3减 3 是为了让正态分布的峰度落在 0 附近这叫超额峰度。按e1071包对偏度的三分类这套写法对应type 3也称 b1 口径和 MINITAB、BMDP 一致type 1用的是总体二阶矩做尺度type 2带样本量校正是 SPSS、SAS 常报的 G1。三者在 n 很小时差异明显n 上百之后差异会收窄但不会消失。# 用同一列数据把三种口径摊开看确认自己的手写函数落在哪一格 v - mtcars$mpg m - mean(v); s - sd(v); n - length(v) c(hand_skew sum((v - m)^3 / s^3) / n, e1071_type1 e1071::skewness(v, type 1), e1071_type2 e1071::skewness(v, type 2), e1071_type3 e1071::skewness(v, type 3), hand_kurt sum((v - m)^4 / s^4) / n - 3, e1071_kurt3 e1071::kurtosis(v, type 3))逻辑说明e1071的skewness()/kurtosis()都带type参数1/2/3把三个都跑出来和手写值并排放哪一列相等就说明手写实现是哪个口径。参数说明type默认是 3所以不写type的时候手写值和包值应该一致一旦有人把默认值当成唯一答案换到 SPSS 口径的报告里就会差出小数点后两位评审时很难解释。这个自校验建议做成脚本固化下来每次换数据集跑一遍比记公式可靠。3.2 Hmisc::describe 与 psych::describe 的输出差异两个包都提供describe()但产物形态完全不同选型要看下游是给人看还是给程序读。对比项Hmisc::describepsych::describe返回形态控制台文本报告为主data.frame可直接取列缺失信息单独列missing、distinctn列缺失需自行推算分位数5%、10%、25%、50%、75%、90%、95%中位数为主配合trim截尾均值极值打印最低 5 个和最高 5 个观测值min、max、range偏度峰度有有受type参数控制适合场景人工翻阅、写报告正文批量汇总、拼表、导出 CSV# 人工阅读Hmisc 会把分位数、极值、缺失一起打出来 Hmisc::describe(mtcars[myvars]) # 程序消费psych 返回 data.frame能直接 index 取列 d - psych::describe(mtcars[myvars]) str(d) # 先看你这台机器上的字段名版本间会增删 d[, c(n, mean, sd, skew, kurtosis)]逻辑说明Hmisc::describe()的输出信息密度高但不好机读适合贴进分析报告的第一节psych::describe()返回的是规整数据框适合往下游传比如拼成多组对比表。参数说明psych::describe()的trim控制截尾均值比例默认 0.1输出里的trimmed列就是它type控制偏度峰度口径默认 3与手写公式同口径的前提是你没动过它。至于kurtosis这一列是否已经减过 3各实现的习惯不完全统一别靠猜——拿一列已知正态性的模拟数据跑一次数值落在 0 附近就说明减过了落在 3 附近就是原始峰度。这一步花 30 秒能省掉后面被追问时的心虚。3.3 装包、加载与命名空间describe 到底调的是谁两个包都有describe()这是最典型的命名空间覆盖场景。install.packages(c(Hmisc, psych)) library(Hmisc) library(psych) # 后加载同名的 describe 被 psych 覆盖 describe(mtcars[myvars]) # 实际调用的是 psych::describe Hmisc::describe(mtcars[myvars]) # 显式指定永远是 Hmisc 的 # 查当前会话里有哪些冲突 conflicts(detail TRUE)逻辑说明R 的搜索路径是后进先出后library()的包排在前面同名函数后者胜出。conflicts()会把所有被遮蔽的函数名列出来describe这种高频名字一定要查一次。参数说明conflicts(detail TRUE)会额外列出冲突双方的包名比默认输出更直观。更要紧的是加载成本。Hmisc的依赖较重library(Hmisc)之后搜索路径会明显变长你的脚本里如果还有别的包同名冲突排查难度会翻倍。我的做法是全程用::显式调用只在交互式探索时才library()脚本里一律写全名——多打几个字符换来的是半年后回看还能知道当时调的是谁的函数。4. 分组描述性统计aggregate、by、summaryBy 三条路的取舍4.1 aggregate 的公式接口分组算均值是最常见的需求aggregate()的公式接口写起来最短左边是要算的变量右边是分组键。# 单变量分组均值 aggregate(mpg ~ cyl, data mtcars, FUN mean) # 多变量左侧用 cbind 包起来 aggregate(cbind(mpg, hp, wt) ~ cyl, data mtcars, FUN mean) # 多个统计量分开算再按分组键 merge别让 FUN 一次返回多个值 agg_mean - aggregate(cbind(mpg, hp, wt) ~ cyl, data mtcars, FUN mean) agg_sd - aggregate(cbind(mpg, hp, wt) ~ cyl, data mtcars, FUN sd) merge(agg_mean, agg_sd, by cyl, suffixes c(_mean, _sd))逻辑说明aggregate()的FUN必须返回单个值。如果硬塞一个返回长度大于 1 的函数比如quantile或range它不会报错退出而是给出警告并做循环填充结果列数对不上、数值错位这是最难发现的一类错误。所以多统计量的正确姿势是分次调用再merge()。参数说明na.action默认继承options(na.action)一般是na.omit意味着分组键带NA的行会被整行丢掉drop TRUE控制是否丢弃未使用的因子组合因子有空水平时这两个参数直接决定行数对不对。merge()的by要和分组列同名命名不一致时用by.x/by.y指定别指望它自动猜。4.2 by() 返回 list整形才是成本by()的定位不是出表而是把数据框按分组切开、逐块套函数返回一个 list。它比aggregate()灵活因为你可以对每组做任意复杂的事代价是结果得自己整形。# 每组返回一个 5 行统计量× 3 列变量的矩阵 dstats - function(df) sapply(df, mystats) res - by(mtcars[myvars], mtcars$cyl, dstats) str(res, max.level 1) # 先看清结构list of matrix # 整形成长表每组转置成每行一个变量再补上分组键和变量名 out - do.call(rbind, lapply(names(res), function(g) { m - as.data.frame(t(res[[g]])) m$var - rownames(m) m$cyl - as.numeric(g) m })) out逻辑说明by()输出的 list 每个元素继承输入的结构矩阵被转置后行名是变量名这一步要在整形时显式取出来存成var列否则变量信息就丢了。do.call(rbind, ...)比循环rbind快一个量级分组多的时候差别明显。参数说明第二个参数INDICES会被转成因子所以数值型分组键跑完之后names(res)是字符型的4、6、8转回数值要用as.numeric()这一步漏了后面所有按数值比较的代码都会失败。by()没有na.rm之类的参数缺失要在dstats内部处理这也是把手写函数设计成默认剔缺失的原因之一。4.3 summaryBy 一步出宽表doBy包的summaryBy()把分组 多统计量 拼宽表合成一步输出直接是 data.frame少一次整形。install.packages(doBy) library(doBy) # FUN 返回命名向量时列名由变量名和统计量名拼成 s - summaryBy(mpg hp wt ~ cyl, data mtcars, FUN mystats) head(s) names(s) # 先看列名长什么样再决定要不要重命名逻辑说明summaryBy()的列名拼接规则和FUN的返回值命名强相关——返回命名向量时列名可预期返回未命名向量时按位置编号可读性差很多。所以mystats()里那几个c(n ..., mean ...)的名字不是装饰是给这一步用的。参数说明FUN可以传自定义函数也可以传函数列表keep.names控制分组列是否保留原列名order.groups控制输出是否按分组键排序。不同版本的列名细节可能微调稳妥做法是跑完先names()看一眼再写下游代码别凭记忆写死列名。4.4 分组键的 NA、因子空水平与 drop 行为三个函数对分组键的处理并不一致这是分组分析最容易翻车的地方。第一条分组键里的NA。aggregate()走的是na.action默认整行丢弃by()走的是split行为取决于参数summaryBy()又是另一套。所以在分组之前先看清楚到底有多少行分组键是空值table(mtcars$cyl, useNA always) # 把 NA 显式算作一类useNA always是必须写的默认的useNA no会把NA藏起来让你误以为分组数是对的。第二条因子空水平。分组变量是因子、又经历过筛选之后未使用的水平仍然留在levels里某些函数会给你输出一行全NA的空组。解决办法是筛选后立刻droplevels()把空水平清掉。第三条分组列的数值/字符类型。分组键是字符时输出顺序按字母排10会排在2前面要按数值排就先转因子并显式指定levels。三个函数在这点上表现一致但后果一样难查——表格行序错了看图的人很容易理解反。5. 避坑与排查五条踩过的坑从现象倒推到修法5.1 数据类型与缺失值相关的坑坑一字符列的频数表凭空消失。现象是summary(df)里本该有的类别计数没了只剩Length、Class、Mode三行。原因是从 R 4.0 起读入数据默认stringsAsFactors FALSE字符向量交给summary()只给这三个属性频数统计是因子才有的待遇。解决是显式转换或者换工具summary(factor(df$col))、table(df$col, useNA always)或者用janitor::tabyl()一步出频数和占比。这一步别偷懒报告里缺了类别分布评审第一个问题就答不上。坑二sapply 算出来整列 NA。现象是sapply(df, mean)里某一列明明是数值结果却是NA。原因是这列中间夹了一个缺失值mean()默认na.rm FALSE一个NA污染整列。解决是sapply(df, mean, na.rm TRUE)把参数透传进去或者更保守一点先colSums(is.na(df))确认缺失密度缺失超过三成的列直接排除在统计之外别硬算。顺带一句sapply的简化行为不稳定某列返回长度不一致时会静默退化成 list下游取数就错位了用vapply(df, mean, numeric(1), na.rm TRUE)把返回类型写死更安全。5.2 分组聚合与手写函数相关的坑坑三mystats 在某些组返回 NaN。现象是分组跑完某一组的skew和kurtosis是NaN其余组正常。原因有两个一是该组样本量极小n 小于 3sd()本身就不稳二是该组某个变量是常量sd 0公式里出现0/0。解决是在函数里加护栏if (length(unique(x)) 2) return(c(n n, mean m, stdev NA, skew NA, kurtosis NA))让不成立的统计量显式返回NA而不是NaN下游过滤时能区分算不出来和算出来是空。另外报告里凡是偏度峰度都要和n一起呈现n 小的组别单独标注。坑四aggregate 的 FUN 返回多个值结果列错位。现象是aggregate(cbind(a, b) ~ g, data df, FUN quantile)只给警告不报错出来的列数对不上数值还挺像那么回事。原因是aggregate要求FUN返回长度 1返回向量时会做循环填充。解决是拆成多次调用再merge()或者换summaryBy()/dplyr::summarise()。判断方法是看列数cbind里两个变量、每组该出两列多出来或对不上就是踩了这条。血泪经验是这类错误不会中断脚本只会安静地把错误数字写进报告。5.3 包加载与命名空间相关的坑坑五describe 调的不是你以为的那个。现象是同一个describe()调用昨天输出是文本报告今天变成 data.frame。原因是脚本里先library(Hmisc)后library(psych)后者的同名函数把前者覆盖了。解决是全程用Hmisc::describe()和psych::describe()显式限定并在脚本开头跑一次conflicts(detail TRUE)把冲突列出来。更隐蔽的一层是Hmisc依赖较重加载后会带进一串包如果项目里还有别的统计包冲突面会扩大。我的处理原则是脚本里永远写全名只有交互式探索时才library()这样半年后回看还能确认当时调的是谁。顺手提一句手写偏度峰度和包结果对不上九成不是算错是口径不同type 1/2/3 的差别先按第 3 章的对照表校验一遍再改代码。6. 把统计量画回去对账ggplot2 直方图、箱线图与一个收尾习惯数值表和图必须能对上否则报告里两个部分会互相打脸。ggplot2在这里的价值不只是好看是能让你用图形反向验证summary()的数字。library(ggplot2) # 直方图 均值参考线binwidth 决定形状换一个值形状就变 p1 - ggplot(mtcars, aes(x mpg)) geom_histogram(binwidth 2, fill steelblue, color white) geom_vline(xintercept mean(mtcars$mpg), linetype dashed, color red) geom_vline(xintercept median(mtcars$mpg), linetype dotted) # 箱线图分组变量必须 factor 化否则被当成连续变量映射 p2 - ggplot(mtcars, aes(x factor(cyl), y mpg)) geom_boxplot(width 0.5, outlier.colour red) stat_summary(fun mean, geom point, shape 23, size 3, fill yellow) # 把离群点抽出来和箱线图对照 boxplot.stats(mtcars$mpg)$out逻辑说明binwidth 2是手动指定的组距不指定的话不同版本会有不同的默认分箱同一份数据两次出图可能长得不一样报告里就别谈分布形状了。geom_vline把均值和单位数画在同一张图上两条线离得越远偏态越明显和第 3 章算出来的偏度符号方向应该一致——正偏度对应均值在右、长尾在右。stat_summary(fun mean, ...)是把分组均值直接叠到箱线图上箱体位置和均值菱形点是否贴合一眼能看出组内是否有严重偏斜。参数说明outlier.colour只控制离群点颜色不影响判定规则判定规则固定为 1.5 倍四分位距和boxplot.stats()的$out一致这两个对上才算图数一致。factor(cyl)这一步不能省数值型x会被当成连续变量处理箱线图的横轴间距会变成 4、6、8 的真实距离视觉上误导性很强。对账的具体动作是三项summary(mtcars$mpg)里的 1st Qu 和 3rd Qu 应当等于箱体上下沿箱线图上的红点应当和boxplot.stats()$out一一对应直方图的分箱形状换binwidth后如果剧烈变化说明样本量不足以支撑形状结论报告里要写明。三项都过一遍再往下写建模部分。再补一个进阶手法把第 4 章的分组统计量和图绑在一起用stat_summary(fun.data mean_sdl, fun.args list(mult 1), geom errorbar, width 0.2)在箱线图上加均值 ± 1 标准差误差棒mean_sdl来自Hmisc用之前记得library(Hmisc)不然会报找不到函数。误差棒和箱体对不上通常就是第 4 章分组键类型或者drop行为出了问题图会第一时间提醒你。从那以后我每次交描述性统计表都强制走一遍三方对账summary()出骨架、手写mystats()出偏度峰度、psych::describe()出交叉验证三者数值差超过小数点后三位就停下来查口径绝不带着疑问往下写。图出完再回头核一遍分位数和离群点这套流程跑顺了后面建模阶段被数据问题打断的次数会明显少。希望帮到你。本文还有配套的精品资源点击获取