尧图建网站 尧图建网站 YAOTU WEB BUILD 免费咨询
ARTICLE DETAIL

资讯详情

深耕网站建设与建站编程的一线实战洞察。

R语言描述性统计:summary()的局限与扩展工具链

R语言描述性统计:summary()的局限与扩展工具链 简介这份PDF文档面向R语言初学者与数据分析入门者系统讲解描述性统计分析的完整方法体系帮助读者快速掌握从基础统计量到分组分析与可视化的核心技能。内容覆盖summary()、sapply()等内置函数偏度与峰度的自定义函数实现Hmisc与psych包的describe()函数以及aggregate()、by()、doBy包的summaryBy()分组统计方案并借助ggplot2绘制直方图与箱线图直观呈现数据分布特征。资源包共1个PDF文件约98KB篇幅精炼、代码示例完整便于随查随用。目前已有553人学习下载适合希望夯实数据探索基础、为后续数据挖掘与机器学习打底的学习者参考。1. 描述性统计为什么 summary() 给出的答案经常不够用很多人第一次在 R 里做描述性统计都是敲一句summary(df)就完事。最小值、最大值、四分位数、均值六个数一屏打完看起来数据已经被“看透”了。但真正踩过坑的人会发现这六个数字经常骗人均值被极端值拉走、四分位数掩盖双峰分布、缺失值悄悄吃掉样本量而summary()对此一言不发。描述性统计分析的本质不是“算几个指标”而是用尽量少的统计量把数据分布的中心趋势、离散程度和形状刻画清楚为后续建模和异常识别提供可信的起点。R 语言在这件事上提供了从内置函数、自定义函数到 Hmisc、psych 等扩展包的多层工具链理解它们各自的计算口径和边界比记住函数名重要得多。搞懂这套机制你才知道什么时候该信summary()什么时候必须自己动手。2. summary() 与 sapply()内置函数的计算口径和适用边界2.1 summary() 到底对每一列做了什么summary()是一个泛型函数它的行为完全取决于输入对象的类。对数据框调用时它会逐列分派数值型列走summary.default输出最小值、下四分位数、中位数、均值、上四分位数、最大值因子型列走频数统计逻辑型列统计 TRUE/FALSE 个数。这个分派机制意味着同一张表里不同类型的列会得到完全不同形态的输出这既是便利也是陷阱。data(mtcars) myvars - c(mpg, hp, wt) summary(mtcars[myvars]) # 输出形如 # mpg hp wt # Min. :10.40 Min. : 52.0 Min. :1.513 # 1st Qu.:15.43 1st Qu.: 96.5 1st Qu.:2.581 # Median :19.20 Median :123.0 Median :3.325 # Mean :20.09 Mean :146.7 Mean :3.217 # 3rd Qu.:22.80 3rd Qu.:180.0 3rd Qu.:3.610 # Max. :33.90 Max. :335.0 Max. :5.424关键点在于summary()默认不处理缺失值参数它内部对NA的处理是“统计非缺失部分”但不会告诉你缺失了多少个。也就是说如果一列有 30% 的值是NA你看到的均值只基于剩下 70% 的样本而输出里没有任何提示。这是最常见的翻车点之一。另外summary()不输出标准差、偏度、峰度对分布形状的刻画几乎为零——它只给你位置信息不给形状信息。2.2 sapply() 的向量化思维一次算一族统计量sapply()的作用是把一个函数应用到列表或向量的每个元素上并尝试把结果简化为向量或矩阵。在描述性统计场景里它的价值在于“对多列同时算同一个统计量”避免手写循环。# 对三列分别计算均值和标准差 sapply(mtcars[myvars], mean) sapply(mtcars[myvars], sd) # 一次算多个统计量传入自定义的匿名函数 sapply(mtcars[myvars], function(x) c( n length(x), mean mean(x), sd sd(x), median median(x) ))第二段代码返回的是一个矩阵行是统计量名列是变量名。这个结构非常适合直接喂给后续的报表或绘图。但要注意sapply()的“简化”行为是自动的——如果每次返回的长度不一致它会退化成列表这时候下游代码如果假设是矩阵就会报错。常见的实现思路是需要稳定结构时改用vapply()显式声明返回模板比如vapply(mtcars[myvars], mean, numeric(1))一旦某列算出非标量结果就直接报错而不是悄悄变形。还有一个容易被忽略的细节sapply()对数据框是按列迭代的如果某一列是字符型mean()会直接报错并中断整个调用。所以实际使用时要么先做类型筛选要么在函数内部加类型判断。这跟summary()的“逐列分派、各算各的”是两种哲学前者要求你保证输入同质后者替你做了分派但牺牲了可控性。3. 自定义统计函数偏度、峰度与缺失值处理的源码级拆解3.1 为什么内置函数不给你偏度和峰度偏度和峰度描述的是分布形状偏度衡量不对称性峰度衡量尾部厚度。R 的基础包stats里其实有mean、sd、var、quantile但没有直接暴露偏度和峰度的计算函数。原因在于这两个统计量的定义存在多种口径——样本偏度有有偏估计和无偏估计之分峰度有超额峰度减 3和非超额之分。基础包选择不替用户做这个决定把口径选择留给上层包或用户自己。原文给出的自定义函数是这样的mystats - function(x, na.omit FALSE) { if (na.omit) { x - x[!is.na(x)] } m - mean(x) n - length(x) s - sd(x) skew - sum((x - m)^3 / s^3) / n kurt - sum((x - m)^4 / s^4) / n - 3 return(c(n n, mean m, stdev s, skew skew, kurtosis kurt)) }逐行看na.omit参数控制是否先剔除缺失值这是一个显式的开关默认 FALSE 意味着有NA时结果就是NA——这比悄悄忽略更安全。skew用的是三阶中心矩除以标准差三次方再除以 n这是有偏的样本偏度估计kurt是四阶中心矩除以标准差四次方再减 3即超额峰度。注意这里的s是样本标准差分母 n-1而矩计算里除以的是 n这种混用在统计上并不严格一致但在描述性场景下足够用。如果要做推断应该换成基于无偏估计的公式或者直接用e1071包的skewness()和kurtosis()。3.2 自定义函数在 sapply 下的返回结构把mystats套进sapply返回的是一个 5 行、列数等于变量数的矩阵sapply(mtcars[myvars], mystats) # mpg hp wt # n 32.000000 32.00000 32.0000000 # mean 20.090625 146.68750 3.2172500 # stdev 6.026948 68.56287 0.9784574 # skew 0.610655 0.72602 0.4231465 # kurtosis -0.372766 -0.13556 -0.0227108这个矩阵的行名来自return(c(n..., mean...))里的命名列名来自数据框列名。机制上sapply检测到每次调用返回长度相同的命名向量就把它们按列拼成矩阵。这里有个血泪经验如果函数内部因为某列全是NA导致sd()返回NA整个矩阵对应位置就是NA不会报错但下游如果做排序或筛选就会静默出错。所以自定义函数里最好加一句if (anyNA(x) !na.omit) return(rep(NA_real_, 5))之类的防御或者干脆在调用前统一做缺失值审计。4. Hmisc 与 psych 的 describe()同名函数背后的不同设计取舍4.1 两个 describe() 不是一回事Hmisc::describe()和psych::describe()名字一样输出形态和计算口径差别很大。Hmisc的版本偏向“人读”它输出每个变量的样本量、缺失数、唯一值数、均值、分位数还会给出最频繁出现的几个值适合做数据质量初筛。psych的版本偏向“机读”返回一个数据框列包括 n、mean、sd、median、min、max、range、skew、kurtosis、se直接可以拿去做表或继续计算。library(Hmisc) describe(mtcars[myvars]) library(psych) describe(mtcars[myvars])机制上的差异来自包的设计目标Hmisc是通用数据描述工具箱describe()要兼顾数值、因子、字符多种类型psych是心理测量学背景describe()默认输出偏度和峰度因为量表分析里形状指标是刚需。这也解释了为什么psych::describe()的峰度默认是超额峰度而Hmisc的输出里峰度的呈现方式不同——口径不统一是跨包协作时最容易踩的坑。4.2 加载顺序导致的函数遮蔽两个包都导出describe谁后加载谁生效。如果先library(psych)再library(Hmisc)直接调describe()走的是 Hmisc 版本反过来则走 psych 版本。这种遮蔽不会报错只会让输出形态突然变化排查起来很费时间。稳妥的写法是始终带命名空间调用psych::describe()或Hmisc::describe()。这不是洁癖而是可复现性的基本要求——脚本换一台机器、换一个加载顺序结果就不一样属于典型的玄学问题根子就在命名空间冲突。另外psych::describe()有一个fast参数设为 TRUE 时跳过一些耗时计算适合大数据集快速扫一眼Hmisc::describe()则有digits参数控制输出精度。这些参数在原文里没展开但实际使用时直接影响输出可信度值得翻一下?describe的帮助页确认默认值。5. 分组描述性统计aggregate、by 与 summaryBy 的机制差异5.1 aggregate() 的公式接口与返回值形态aggregate()支持公式写法mpg ~ cyl表示按cyl分组算mpg。它的返回值是一个数据框分组列在前统计结果在后。机制上它内部先对分组变量做split再对每个子集应用函数最后rbind回来。aggregate(mpg ~ cyl, data mtcars, mean) # cyl mpg # 1 4 26.66364 # 2 6 19.74286 # 3 8 15.10000局限在于一次只能算一个统计量或者用cbind合并多个公式且返回的是数据框多统计量时列名不够直观。如果要对多个变量、多个统计量同时分组公式会变得很长。5.2 by() 的分组-应用模型by()是tapply()在数据框上的推广它把数据框按分组变量切成子数据框列表对每个子集调用函数返回一个列表。dstats - function(x) sapply(x, mystats) by(mtcars[myvars], mtcars$cyl, dstats)返回的是一个列表每个元素对应一个分组元素内部是mystats算出的矩阵。这个结构信息量大但不适合直接输出成表——需要再do.call(rbind, ...)或配合plyr/purrr整理。by()的优势是灵活函数可以是任意复杂逻辑返回任意结构代价是结果需要二次加工。5.3 summaryBy() 的公式化多统计量doBy::summaryBy()把公式接口和多统计量结合得最自然library(doBy) summaryBy(mpg hp wt ~ cyl, data mtcars, FUN mystats)它返回一个数据框分组列加上每个变量每个统计量的列列名形如mpg.mean、hp.skew。机制上它内部对FUN的返回向量做展开和命名拼接所以要求FUN返回命名向量——这正是mystats里c(n..., mean...)命名的意义所在。如果返回的是无名向量列名会变成mpg.1、mpg.2这种可读性骤降。三种方式的选择逻辑只要一个统计量、要数据框用aggregate()要复杂自定义逻辑、能接受列表用by()要多变量多统计量、要整齐数据框用summaryBy()。没有哪个绝对更好差别在返回结构和后续消费方式。6. 分组统计与可视化的避坑清单6.1 分组变量类型导致的静默错误现象aggregate(mpg ~ cyl, ...)结果里分组数比预期少。原因cyl如果是数值型但实际是分类含义aggregate按唯一值分组没问题但如果分组变量是连续数值会得到大量只有一两个样本的组均值毫无意义。解决分组前用factor()显式转换并检查table()的组样本量。6.2 缺失值在分组统计中的连锁反应现象某个分组的均值是NA其他分组正常。原因该组内目标变量存在NA而mean()默认na.rm FALSE。解决在FUN里显式传na.rm TRUE或者用mystats的na.omit开关并在结果里保留n以便核对有效样本量。6.3 ggplot2 直方图 binwidth 选错导致分布误判现象直方图看起来是单峰但箱线图显示有离群点。原因binwidth太大把双峰抹平了。解决多试几个binwidth或者改用geom_density()叠加箱线图的outlier点也要单独看。原文示例里binwidth 1对mpg是合理的但换一个量纲不同的变量就必须重调。6.4 箱线图分组因子顺序错乱现象箱线图 x 轴分组顺序不是 4、6、8 而是按字母或出现顺序。原因factor(cyl)默认按字典序或唯一值顺序。解决factor(cyl, levels c(4, 6, 8))显式指定否则图的可读性和结论方向都可能被误导。6.5 包加载顺序引发的 describe 遮蔽现象同一段代码昨天输出是数据框今天输出是文本块。原因Hmisc和psych的describe互相遮蔽取决于library()顺序。解决始终用psych::describe()或Hmisc::describe()带命名空间调用脚本里不依赖加载顺序。7. 从统计量到判断描述性分析的取舍与验证习惯描述性统计最容易给人的错觉是“算完了就等于看懂了”。summary()给你六个数字psych::describe()给你十几个summaryBy()给你一张宽表但数字多不等于信息多。真正的取舍在于你到底要回答什么问题。如果只是数据质量初筛Hmisc::describe()的缺失数、唯一值数、频次Top值比偏度峰度有用得多如果要做分布形状判断偏度和峰度必须配合直方图或密度图交叉验证单看数字容易被样本量小导致的波动带偏。一个可操作的验证习惯是任何描述性统计结果至少用两种口径交叉核对。均值用summary()和sapply(..., mean)各算一遍确认缺失值处理一致分组统计用aggregate()和summaryBy()各跑一次确认分组数和样本量对得上偏度用自定义函数和psych::describe()对比确认口径差异在可接受范围内。这种“双算”看起来笨但能拦住大部分静默错误——描述性统计的坑几乎都不报错只是悄悄给你一个错的数。需求首选工具注意点快速扫一眼summary()不报缺失数无形状指标多列同统计量sapply()/vapply()返回结构可能退化形状指标psych::describe()峰度为超额峰度数据质量初筛Hmisc::describe()输出为人读文本分组多统计量doBy::summaryBy()FUN 必须返回命名向量分布可视化ggplot2直方图箱线图binwidth 和因子顺序要显式设我自己的习惯是拿到任何新数据集先跑一遍summary()看量纲和缺失再跑psych::describe()看形状然后按分组变量画箱线图三件事做完才允许自己下“数据正常”的结论。有一次就是靠箱线图发现某个分组的样本量只有 2而summary()的均值看起来毫无异常——从那以后我每次做分组统计都强制先看table()的组样本量再谈均值。希望帮到你。本文还有配套的精品资源点击获取
返回列表