生物统计必备：用R语言分析螃蟹数据中的正态分布特征

生物统计实战：用R语言透视螃蟹数据中的正态分布奥秘

在生态学和生物统计学的研究中，我们常常需要面对一个基础却至关重要的问题：手头的数据是否符合正态分布？这不仅仅是一个形式上的检验，它直接关系到后续应该选择哪种统计方法——是参数检验还是非参数检验？结论的可靠性有多高？对于许多刚接触数据分析的研究生或领域内的实践者来说，理解如何操作并解读正态性检验，是迈出科学分析的第一步。

今天，我们就以一个经典且有趣的案例——螃蟹形态测量数据，来深入探讨这个问题。我们将超越简单的“是”或“否”的判断，聚焦于如何利用R语言强大的可视化与统计检验工具，特别是lattice包的分组绘图能力，来细致地审视数据在不同组别（例如公蟹与母蟹）下的分布特征。你会发现，掌握这些技巧，不仅能让你更自信地处理自己的数据，更能从数据中挖掘出更深层次的生物学故事。

1. 理解数据与正态分布的核心逻辑

在动手写代码之前，我们有必要先厘清几个核心概念。正态分布，这个在统计学中无处不在的“钟形曲线”，其重要性在于许多经典的统计推断方法（如t检验、方差分析、线性回归）都建立在数据服从或近似服从正态分布的假设之上。如果这个前提不成立，盲目使用这些方法可能会导致错误的结论。

那么，我们如何判断数据是否正态呢？主要有三种相辅相成的方法：

1. 可视化检验：通过图形直观感受数据的分布形态，如直方图、密度曲线图、Q-Q图。
2. 统计量检验：计算偏度和峰度等描述性统计量，量化分布与正态的偏离程度。
3. 统计假设检验：执行正式的统计检验（如Shapiro-Wilk检验），给出一个概率值（p-value）作为客观依据。

对于生物数据，尤其是像螃蟹的甲壳宽度（Carapace Width）这类连续型测量数据，其分布往往受到性别、年龄、种群等多种因素的影响。因此，分组检验显得尤为重要。我们绝不能将公蟹和母蟹的数据混在一起检验，因为性别差异本身就可能导致混合后的分布呈现双峰或其他非正态形态，从而掩盖了各组内部真实的分布情况。

提示：在生物统计学中，分组分析是避免“辛普森悖论”等统计陷阱的关键。始终问自己：我的数据是否需要按某个分类变量（如性别、地点、处理组）进行分层审视？

2. 数据准备与lattice包的分组绘图魔法

我们的实战对象是R语言MASS包中的crabs数据集。它包含了200只螃蟹（100只雄性，100只雌性）的多种形态测量值。首先，让我们准备好战场。

# 加载必要的程序包
library(MASS)      # 包含crabs数据集
library(lattice)   # 用于高级条件绘图
library(dplyr)     # 用于数据操作，可选但推荐

# 查看数据结构
data(crabs)
str(crabs)
head(crabs)

运行str(crabs)后，你会看到数据框中有sp（物种）、sex（性别）、index（编号）以及FL（前叶长度）、RW（后叶宽度）等测量指标。我们重点关注CW（甲壳宽度）和sex。

lattice包是R语言中用于创建条件（或称“栅格”）图表的强大系统。它的核心思想是：基于一个或多个条件变量，将数据分割成多个子集，并在同一个图形框架中并排绘制。这正好完美契合了我们按性别分组分析的需求。

其基本绘图函数，如histogram()，使用一种**公式接口（formula interface）**来指定绘图关系。公式的写法是理解其魔法的钥匙：

~ y | condition

• ~：表示“由……描述”的关系。
• y：放在波浪号右边，是我们要绘图的数值变量。
• condition：放在竖线|右边，是作为分组条件的因子变量。这个|符号可以读作“给定”或“基于”。

所以，~ CW | sex 的意思就是：绘制变量CW的分布图，并基于sex的不同水平（M和F）进行分组展示。

3. 构建分组直方图：第一眼洞察

现在，让我们绘制第一张分组直方图，对公蟹和母蟹的甲壳宽度分布有一个直观的印象。

# 基础分组直方图
histogram(~ CW | sex, data = crabs,
          main = "螃蟹甲壳宽度（CW）按性别分布直方图",
          xlab = "甲壳宽度 (mm)",
          layout = c(1, 2)) # 设置布局为1列2行

执行这行代码，你将得到并排的两个直方图。仅仅一眼，你就能回答几个关键问题：两个分布的集中趋势（中心位置）相同吗？离散程度（宽度）呢？形状大致对称吗？有没有明显的异常值？

通常，我们会给直方图叠加上核密度估计曲线和正态分布参考线，以便更精确地对比。

# 增强版分组直方图：添加密度曲线和正态参考线
histogram(~ CW | sex, data = crabs,
          main = "螃蟹甲壳宽度分布（含密度曲线）",
          xlab = "甲壳宽度 (mm)",
          type = c("density", "bar"), # 同时绘制条形和密度曲线
          panel = function(x, ...) {
              panel.histogram(x, ...) # 绘制直方图条形
              panel.densityplot(x, col = "red", lwd = 2, ...) # 添加红色密度曲线
              panel.mathdensity(dmath = dnorm,
                                args = list(mean = mean(x), sd = sd(x)),
                                col = "blue", lwd = 2, lty = 2) # 添加蓝色正态参考线
          },
          layout = c(1, 2))

在这段代码中，panel参数允许我们自定义每个子图的绘制内容。我们依次画出了直方图条形、根据数据计算出的核密度曲线（红色实线）、以及基于该组数据均值和标准差拟合的理论正态分布曲线（蓝色虚线）。

如何解读？

• 如果数据的密度曲线（红色）与正态参考线（蓝色）贴合紧密，则正态性较好。
• 观察尾部：如果红色曲线在两端比蓝色虚线更高，说明数据具有“厚尾”特征；如果更低，则是“薄尾”。
• 观察峰值：如果红色曲线的峰比蓝色虚线更尖，称为“尖峰态”；更平缓则称为“平峰态”。

通过这幅图，你可能已经对数据的正态性有了初步判断。但视觉判断具有一定主观性，我们需要更定量的工具。

4. 量化检验：描述性统计与Shapiro-Wilk检验

在深入假设检验之前，先计算一些描述性统计量能给我们更扎实的感觉。我们可以快速计算每组的关键统计量。

# 使用dplyr分组计算描述性统计量
library(dplyr)
stats_summary <- crabs %>%
  group_by(sex) %>%
  summarise(
    Count = n(),
    Mean_CW = mean(CW, na.rm = TRUE),
    SD_CW = sd(CW, na.rm = TRUE),
    Median_CW = median(CW, na.rm = TRUE),
    Skewness = moments::skewness(CW), # 需要安装并加载moments包
    Kurtosis = moments::kurtosis(CW) - 3 # 超额峰度
  )
print(stats_summary)

计算出的表格可能类似这样（数值为示例）：

sex	Count	Mean_CW	SD_CW	Median_CW	Skewness	Kurtosis
F	100	28.1	3.2	28.0	0.15	-0.20
M	100	32.5	4.1	32.3	0.08	0.05

• 偏度（Skewness）：衡量分布不对称性。接近0表示对称；大于0为右偏（长尾在右）；小于0为左偏。
• 峰度（Kurtosis）：这里计算的是超额峰度。接近0表示与正态分布峰度相当；大于0为尖峰厚尾；小于0为平峰薄尾。

从示例数据看，公蟹和母蟹的偏度都接近0，峰度也接近0，这提示数据可能服从正态分布，但还需要假设检验来确认。

接下来是重头戏：Shapiro-Wilk检验。这是在小样本（n < 50）情况下检验正态性的首选方法。

# 对公蟹和母蟹分别进行Shapiro-Wilk检验
# 方法1：手动分组应用
shapiro_male <- shapiro.test(crabs$CW[crabs$sex == "M"])
shapiro_female <- shapiro.test(crabs$CW[crabs$sex == "F"])

cat("雄性螃蟹CW的Shapiro-Wilk检验结果:\n")
print(shapiro_male)
cat("\n雌性螃蟹CW的Shapiro-Wilk检验结果:\n")
print(shapiro_female)

# 方法2：使用by函数或tapply进行分组计算（更简洁）
by(crabs$CW, crabs$sex, shapiro.test)

检验结果会输出W统计量和对应的p值。我们需要关注p值：

• 原假设（H0）：数据服从正态分布。
• p值 > 0.05：没有足够证据拒绝原假设，即不能认为数据偏离正态分布（注意，这不等于“证明”数据是正态的）。
• p值 ≤ 0.05：有足够证据拒绝原假设，认为数据显著偏离正态分布。

注意：Shapiro-Wilk检验对于样本量非常敏感。当样本量很大时（比如n>100），即使数据对正态分布的偏离非常微小，检验也可能给出显著的p值（p<0.05）。此时，应结合图形（如Q-Q图）和描述性统计量进行综合判断，不能单纯依赖p值。

5. Q-Q图：与正态分布的精确对标

当我们需要更精细地查看数据在哪个分位数上偏离了正态分布时，Q-Q图（分位数-分位数图）是最佳工具。它的原理是将数据的实际分位数与理论正态分布的分位数进行散点图对比。如果数据完全正态，这些点应该大致落在一条直线上。

我们可以用lattice的qqmath()函数方便地绘制分组Q-Q图。

# 使用lattice绘制分组Q-Q图
qqmath(~ CW | sex, data = crabs,
       main = "螃蟹甲壳宽度按性别的Q-Q图",
       xlab = "理论正态分位数",
       ylab = "样本分位数 (CW)",
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
           panel.qqmathline(x, ...) # 添加参考线
           panel.qqmath(x, ...)     # 绘制散点
       },
       layout = c(1, 2))

解读Q-Q图的要点：

1. 关注散点与直线的整体贴合度。大部分点落在直线附近是好的迹象。
2. 系统性地观察偏离模式：
   • 如果散点呈“S”形曲线，说明分布有偏态。
   • 如果两端的点向上偏离直线，表示分布有厚尾；向下偏离则表示薄尾。
   • 如果只有头部或尾部的一两个点偏离，可能是异常值。

结合我们之前的分组直方图和Shapiro-Wilk检验，Q-Q图提供了第三维度的证据。如果三者结论一致，那么我们对数据正态性的判断就会非常稳健。

6. 当数据非正态时：应对策略与稳健方法

经过上述分析，如果发现某个组的数据显著偏离正态分布（例如p值远小于0.05，且Q-Q图显示系统性偏离），我们该怎么办？直接放弃参数检验吗？别急，还有几条路可以走。

策略一：数据变换 对于轻度到中度的非正态性（特别是偏态），对数据进行数学变换可能使其更接近正态。常见的变换有：

• 对数变换：适用于右偏（正偏）数据，且数据值全为正数。log(x)或log1p(x)（处理含0的数据）。
• 平方根变换：适用于轻度右偏的计数数据。
• Box-Cox变换：一种更通用的幂变换，可以自动寻找最优的变换参数λ。

# 示例：对数据尝试对数变换后再次检验
crabs$log_CW <- log(crabs$CW)
# 然后对log_CW重复上述的分组直方图、Q-Q图和Shapiro检验
histogram(~ log_CW | sex, data = crabs, ...)
shapiro.test(crabs$log_CW[crabs$sex == "M"])

策略二：使用非参数检验 如果变换无效或数据性质不适合变换，则应选择不依赖于正态分布假设的非参数检验方法。

• 替代t检验：曼-惠特尼U检验（Wilcoxon秩和检验），用于比较两个独立组。
• 替代配对t检验：威尔科克森符号秩检验。
• 替代方差分析：克鲁斯卡尔-沃利斯H检验。

# 示例：使用非参数方法比较公蟹和母蟹的CW中位数是否存在差异
wilcox.test(CW ~ sex, data = crabs, conf.int = TRUE)

策略三：使用稳健统计方法 一些现代统计方法对偏离正态性不那么敏感，例如：

• 自助法（Bootstrap）：通过重抽样来估计参数的置信区间，不依赖正态假设。
• 稳健回归：使用M估计量等方法，降低异常值的影响。

策略四：重新审视样本量与检验效力 有时，轻微的非正态性在大样本下可能不会对t检验等方法的结论产生实质性影响（得益于中心极限定理）。但这是一个需要谨慎评估的领域，最好咨询统计专家。

处理完一组数据，从数据导入、清洗、探索性可视化，到分组统计描述、正式的正态性检验，再到根据结果选择后续分析路径，这构成了一个完整的生物统计数据分析微循环。掌握这个流程，并理解每个工具背后的“为什么”，远比记住几个函数命令更重要。下次当你拿到一批新的实验测量数据时，不妨就从一幅清晰的分组直方图开始你的探索之旅。