如何在R中绘制地球化学分析的对数累积概率图？

我理解你遇到的头疼问题：qualityTools::ppPlot 在处理部分元素数据时，居然因为内部拟合对数正态分布时意外生成负值报错，而base R的qqnorm又完全不是你要的累积概率图样式。下面给你两个靠谱的解决方案，尤其是手动构建的方法，完全可控，特别适配地球化学数据的特点。

方法一：手动绘制累积概率图（强烈推荐） #

地球化学的累积概率图说白了就是排序后的元素浓度（对数刻度）对应经验累积概率，手动计算完全能避开包的内部拟合坑，步骤也很清晰：

代码示例（用你提供的Mg数据） #

# 你的Mg数据
mg <- c(51.400, 149.000, 276.000, 135.000, 179.000, 81.000, 116.000, 8.150, 7.770, 7.870, 8.840, 15.600, 13.400, 57.400, 7.440, 14.800, 40.800, 15.100, 21.400, 5.550, 3.390, 18.800, 20.100, 19.600, 11.600, 11.700, 12.200, 12.500, 11.700, 12.100, 13.000, 12.300, 13.300, 13.200, 12.600, 29.700, 25.400, 21.000, 11.100, 11.500, 11.000, 32.600, 17.500, 16.500, 18.100, 27.200, 21.200, 26.400, 18.800, 19.900, 32.000, 28.600, 29.400, 30.700, 2.370, 2.070, 1.850, 1.970, 24.900, 19.100, 17.400, 23.100, 50.100, 48.800, 18.000, 15.800, 27.100, 43.500, 4.820, 13.400, 14.600, 24.100, 22.700, 22.500, 43.500, 41.300, 43.700, 41.100, 40.800, 63.700, 7.700, 8.360, 60.000, 58.400, 63.100, 65.100, 219.000, 25.800, 4.940, 3.670, 13.800, 5.190, 14.700, 15.000, 13.100, 12.300, 10.700, 10.700, 11.100, 10.100, 10.600, 63.200, 19.800, 22.200, 17.600, 11.500, 10.600, 9.380, 3.190, 9.180, 10.800, 189.000, 190.000, 152.000, 119.000, 194.000, 56.100)

# 1. 先把数据从小到大排序
mg_sorted <- sort(mg)
# 2. 计算经验累积概率（地球化学圈常用 (排名-0.5)/总样本数，比简单的1:n/(n+1)更贴合实际）
n <- length(mg_sorted)
cum_prob <- (seq_along(mg_sorted) - 0.5)/n

# 3. 画图！直接把x轴设为对数刻度
plot(x = mg_sorted, y = cum_prob, 
     log = "x",  # 开启x轴对数刻度，完美符合你的需求
     pch = 16, col = "steelblue",  # 用实心点更清晰
     xlab = "Mg 浓度（对数刻度）",
     ylab = "累积概率",
     main = "Mg 元素累积概率图")

# 要是需要加对数正态拟合参考线，也可以自己加（可选）
library(MASS)
# 先对排序后的浓度取对数，拟合正态分布
ln_fit <- fitdistr(log(mg_sorted), "normal")
# 生成拟合用的x序列
x_seq <- seq(min(mg_sorted), max(mg_sorted), length.out = 100)
# 计算拟合的累积概率
y_fit <- pnorm(log(x_seq), mean = ln_fit$estimate[1], sd = ln_fit$estimate[2])
# 画拟合线
lines(x_seq, y_fit, col = "red", lwd = 2)
# 加个图例
legend("topleft", legend = c("实测数据", "对数正态拟合"), 
       col = c("steelblue", "red"), pch = c(16, NA), lty = c(NA, 1))

为啥推荐这个方法？ #

• 完全避开了ppPlot内部的拟合bug，所有数据都是你提供的正值，绝对不会出现负值报错。
• 累积概率的计算用的是地球化学领域常用的公式，结果更专业。
• 对数刻度直接通过log="x"设置，一步到位，不需要额外调整。

方法二：绕开ppPlot的对数正态拟合坑 #

如果你非要用qualityTools包，其实可以换个思路：对数正态分布的PP图，等价于原始数据取对数后的正态PP图，这样就能绕开包直接拟合对数正态时的问题：

library(qualityTools)
# 先对Mg数据取对数
mg_log <- log(mg)
# 绘制正态分布的PP图
ppPlot(mg_log, distribution = "normal")
# 手动把x轴标签改成原始浓度的对数刻度，不然别人以为是对数后的值
title(xlab = "Mg 浓度（对数刻度）")

注意事项 #

• 这个方法的本质是把对数正态的问题转换成正态问题，避开了包的内部报错，但x轴显示的是对数后的值，必须手动修改标签说明，不然容易误解。

为啥ppPlot会报错？ #

你猜的没错！qualityTools::ppPlot在拟合对数正态分布时，内部的参数估计可能生成了一个左侧尾部延伸到负值区域的分布（比如均值和标准差的组合导致的），而你的数据都是正值，就触发了报错。手动方法完全不需要包去拟合分布来生成曲线（除非你自己主动加），所以根本不会有这个问题。

内容的提问来源于stack exchange，提问作者Nicola Gambaro