第2章 数据处理与简单计算

2.1 加载和处理数据

在 R 中启动一个新会话时，无论是否使用 GCDkit，首要任务都是导入 data。在纯 R 环境中，常见的火成岩地球化学表格 data 可通过 read.table 命令最方便地加载。另一种方法是使用 data 命令访问内置的 dataset 之一。在 GCDkit 中，通常使用 loadData 函数读取 analyses，或通过剪贴板从任意 Windows 应用程序（例如电子表格）复制 data。

数据加载后，最常见的任务包括显示、子集操作（基本的数据库功能）以及简单的算术运算。接下来，我们将针对向量和数据框这两种最常见的数据类型练习这些技能。

练习2.1：子集化数值向量，向量运算

GCDkit 包含多个内置 data，通常仅用于确保系统正常运行。其中之一是原子量，存储在名为 mw 的向量中。我们将使用它来展示一些简单的向量运算。

a) 显示整个向量。Rb的原子量是多少？
b) 整个向量的平均值是多少？
c) 哪些原子的原子量高于170？
d) 显示 data 集中最轻的六个和最重的六个元素的名称。

a) GCDkit‐> mw
银 铝 砷 金 硼 钡 107.86820 26.98154 74.92160 196.96650 10.81100 137.33000
铍 铋 溴 碳 钙 镉 9.01218 208.98040 79.90400 12.01100 40.07800 112.41000 …

GCDkit‐> mw[“Rb”]
[1] 85.4678

b) GCDkit‐> mean(mw)
[1] 107.9206

c) GCDkit‐> names(mw[mw > 170])
[1] “金” “铋” “铪” “汞” “铱” “镥” “镎” “锇” “铅” “铂” “钚”
[12] “铼” “钽” “钍” “铊” “铀” “钨” “镱”

d) GCDkit‐> names(sort(mw)[1:6])
[1] “H” “Li” “Be” “B” “C” “N”
1.00797 6.94100 9.01218 10.81100 12.01100 14.00670

GCDkit‐> names(rev(sort(mw))[1:6])
[1] “Pu” “U” “Np” “Th” “Bi” “Pb”
244.0640 238.0290 237.0482 232.0381 208.9804 207.2000

练习 2.2：加载文件，矩阵/数据框操作

文件 sazava.data 包含来自捷克共和国波希米亚地块中波希米亚深成杂岩（CBPC）萨扎瓦岩套（~354 Ma）的主量和微量元素分析数据（雅努舍克等，2000，2004）。

a) 将制表符分隔的数据文件中存储的分析结果读取到一个数据框WR中。
b) 查看可用变量（=列名）的名称。
c) 样品Sa‐1的MgO含量是多少？
d) 显示样品Po‐1和Po‐4的所有可用数值数据。
e) 计算“Na2O”列的总和。
f) 显示SiO2含量最低和最高的三个样品的名称。
g) 计算所有变量的平均值。
h) 显示一个包含三列的表格：SiO2, MgO和 Na2O/K2O。

a) > sazava <- read.table(“sazava.data”, sep=”\t”)
GCDkit‐> loadData(“sazava.data”) # GCDkit 中的替代方法
GCDkit‐> sazava <- cbind(labels, whole.rock.chemistry)

b) > colnames(sazava)
[1] “侵入体” “产地” “岩石学” “露头”
[5] “符号” “颜色” “SiO2” “TiO2”
[9] “Al2O3” “FeO” “Fe2O3” “MnO” …

c) > sazava[“Sa‐1”, “MgO”]
[1] 3.21

d) > sazava[c(“Po‐1”,”Po‐4”),]
SiO2 TiO2 Al2O3 FeO Fe2O3 MnO MgO CaO Na2O K2O
Po-1 62.95 0.28 20.02 1.65 0.67 0.05 0.55 6.61 3.91 1.99
Po-4 71.09 0.30 15.09 2.12 0.38 0.06 0.52 3.75 3.68 1.87 …

e) > sum(sazava[,”Na2O”])
[1] 39.13

f) > silica.content <- sazava[,”SiO2”]

names(silica.content) <- rownames(sazava)
names(sort(silica.content))[1:5]
[1] “Gbs-1” “Gbs-20” “Sa-4” “Gbs-2” “SaD-1”
names(rev(sort(silica.content)))[1:5]
[1] “Po-4” “Po-3” “Po-5” “Po-1” “Sa-1”

g) > apply(sazava, 2, mean, na.rm=TRUE)
SiO2 TiO2 Al2O3 FeO
57.95285714 0.63928571 16.94285714 4.73071429
Fe2O3 MnO MgO CaO
1.74642857 0.13785714 3.57000000 8.16000000 …

h) > var.x <- cbind(sazava[,”SiO2”], sazava[,”MgO”], sazava[,”Na2O”]/sazava[,”K2O”])

colnames(var.x) <- c(“SiO2”,”MgO”,”Na2O/K2O”)
rownames(var.x) <- rownames(sazava)
var.x
SiO2 MgO Na2O/K2O
Sa-1 59.98 3.21 1.008000
Sa-2 55.17 3.67 1.976471 …

2.2 将全岩化学成分与矿物化学计量关联

加载后，一项常见任务是将整体地球化学分析结果转化为多个与岩石矿物学相关的指数。这通常紧接着

1 注意，此情况下将创建两个变量。WR 将仅包含数值，所有文本信息将被转移到数据框 labels 中。详见附录B。

2.2.1 基性指数

在R语言中，通过简单的地球化学指数示例可以最好地演示计算过程。其中使用最广泛且功能强大的是将Fe重新计算为全亚铁或高铁氧化物，以及两种类型的镁数：

$$
FeOt = FeO + 0.89981 \times Fe_2O_3 \quad [\text{wt.%}]
$$
(2.1)

$$
mg# = \frac{MgO}{FeO + MgO} \times 100 \quad [\text{mol.%}]
$$
(2.2)

$$
Mg# = \frac{MgO}{FeOt + MgO} \times 100 \quad [\text{mol.%}]
$$
(2.3)

许多基于主量元素的图解是通过将氧化物重量百分比转换为阳离子数来构建的；这使得能够轻松地与矿物化学式进行比较。例如，常用的铝饱和指数 A/CNK（尚德 1943）（有时也简称为ASI）模仿了长石的化学计量关系：

$$
A/CNK = \frac{Al_2O_3}{CaO + Na_2O + K_2O} \quad [\text{mol.%}]
$$
(2.4)

一个类似的指数，用于区分过碱性岩石（含有过量碱金属），其表达式为：

$$
A/NK = \frac{Al_2O_3}{Na_2O + K_2O} \quad [\text{mol.%}]
$$
(2.5)

如果A/CNK > 1，则铝含量超过形成长石所需的量。这类岩石被称为过铝质，而A/CNK < 1且A/NK > 1的岩石为偏铝质，A/CNK ~ 1的为亚铝质（Shand 1943）。A/CNK值与实测矿物学有直接关系—出现钙和/或碱金属富集的矿物相，如角闪石和辉石，表明铝不足，其寄主岩石为偏铝质。黑云母呈弱过铝质，因此其出现表明岩石具弱过铝质特征（Miller 1985）。强过铝质花岗岩类（按Miller 1985定义）含有更多过铝质矿物相，如白云母，甚至铝硅酸盐（蓝晶石、硅线石或红柱石）、堇青石、石榴子石、电气石、黄玉或刚玉（克拉克 1981）。然而，

铁数（弗罗斯特等人 2001；弗罗斯特和弗罗斯特 2008）的定义类似。

不同作者对铝饱和指数（ASI）的定义存在差异——有时钙（Ca）会针对磷灰石进行校正，甚至其定义具有误导性，完全不能反映长石的化学计量比（弗罗斯特等人 2001）。

练习2.3：计算简单指数

在萨扎瓦数据集上，我们可以演示如何定义一个计算地球化学指数的函数。通过这种方式，系统可以快速且高效地得到扩展。

a) 根据以下分子量，设计一个函数来计算mg数量。
FeO MgO Al2O3 CaO Na2O K2O
71.85 40.31 101.96 56.08 61.98 94.20

b) 编写一个返回尚德指数（A/CNK 和 A/NK）的函数。
c) 为萨扎瓦数据集计算所有这些值。
d) 将主要元素氧化物重新换算到100%无挥发物基准。

sazava <- read.table(“sazava.data”, sep=”\t”)
MW <- c(71.85, 40.31, 101.96, 56.08, 61.98, 94.20)
oxides <- c(“FeO”,”MgO”,”Al2O3”,”CaO”,”Na2O”,”K2O”)
names(MW) <- oxides
mol <- t(sazava[, oxides]) / MW[oxides]

a) > mgno <- function() {

mg <- 100 * mol[“MgO”,] / (mol[“FeO”,] + mol[“MgO”,])
return(mg)
}

b) > ank <- function() {

ANK <- mol[“Al2O3”,] / (mol[“Na2O”,] + mol[“K2O”,])
return(ANK)
}
acnk <- function() {
ACNK <- mol[“Al2O3”,] / (mol[“Na2O”,] + mol[“K2O”,] + mol[“CaO”,])
return(ACNK)
}

c) > x <- cbind(mgno(), acnk(), ank())

colnames(x) <- c(“mg.no”,”A/CNK”,”A/NK”)
x
mg.no A/CNK A/NK
Sa-1 51.16987 0.8355806 2.396569
Sa-2 55.42955 0.7619109 2.307463
Sa-3 51.92059 0.8079150 2.562820 …

d) > major <- c(“SiO2”,”TiO2”,”Al2O3”,”Fe2O3”,”FeO”,”MnO”, “MgO”,”CaO”,”Na2O”,”K2O”,”P2O5”)

sums <- apply(sazava[, major], 1, sum)
anh <- sazava[, major] / sums * 100
anh
SiO2 TiO2 Al2O3 Fe2O3 FeO
Sa-1 60.30565 0.6334205 16.50915 1.3573296 5.489644
Sa-2 56.25000 0.7238989 17.33279 2.7120718 5.362969
Sa-3 56.03133 0.7628153 17.89056 2.1663954 5.909276 …

在GCDkit中计算分子量使用的是函数molecularWeight。该函数还会返回每个化学式的阳离子数和氧原子数，因此需要对结果进行子集提取，例如如下方式：

GCDkit‐> molecularWeight(“Al2O3”)[1]
[1] 101.9613

GCDkit中的简单地球化学指数
加载新数据后，会自动计算若干有用的岩石学指数，包括FeOt、mg#和Mg#，以及A/CNK以及可供绘图等用途的A/NK值。此外，一个矩阵WRanh包含了换算为100%无水基的常量元素氧化物。

因此，练习 a–c 有一个简单的 GCDkit 解决方案：
GCDkit‐> loadData(“sazava.data”)
GCDkit‐> WR[, c(“mg#”,”A/CNK”,”A/NK”)]
GCDkit‐> WRanh

2.2.2 阳离子参数

尼格利（1948）强调了简单阳离子值在火成岩成因解释中的重要性。尼格利的一些阳离子值，si, al, f m, c, alk, k, mg, ti, p, c/f m, 和 qz 仍被使用。该概念后来由法国学者进一步发展为基于毫当量的多阳离子参数：

$$
D = 1000 \times \frac{C}{n \times MW}
$$
(2.6)

其中MW为分子量，n为氧化物公式中的原子数（例如2O，以及Al2O3中的2）。请注意，即使原始分析（氧化物重量百分比）的总和为100，因此没有理由要求总量达到某个特定值。对于某些应用（例如锆石饱和度，第13.1.1节），需要将此总量归一化为1。

出于分类目的，德·拉·罗什等人（1980）采用了两个参数的投影：R1= 4Si – 11(Na+ K) – 2(Fe+ Ti) 和 R2= 6Ca+ 2Mg+ Al，从而涵盖了所有主要元素的氧化物。此外，R1–R2图解具有成因的和大地构造的意义（巴彻勒和鲍登，1985）。

在德邦和勒福尔（1983, 1988）复杂的分类体系中，或许最有用的参数是A= Al –(K+ Na+ 2Ca)（反映过铝性），B= Fe+ Mg+ Ti（镁铁性），P= K –(Na+ Ca)（长石中钾长石的比例）以及Q= Si/3 –(K+ Na+ 2Ca/3)（石英含量）。

毫当量及相关分类方案

将新 data 加载到 GCDkit 后，所有 analyses 都会重新计算为毫当量，并存储在 data 矩阵 milli 中。此外还提供了一些函数计算基于毫当量的指数，并生成一些相关图件（De La Roche 等人，1980； Debon 和 Le Fort，1983，1988；巴彻勒和鲍登，1985；比利亚塞卡等人，1998）。请参阅 LaRocheCalc、LaRoche、DebonCalc、Debon、Batchelor 和 Villaseca 的帮助以了解更多信息。

2.2.3 火成岩的规范性计算与分类

在GCDkit中计算标准矿物

GCDkit的大部分标准值重新计算方案均沿用自其前身NORMAN（亚努舍克，2001）。可用的包括 用于CIPW标准矿物计算的模块，包括含黑云母和角闪石的修正方法（哈奇森，1974，1975）、 Catanorm（哈奇森，1974及其中引用的文献）以及改进型花岗岩Mesonorm（米尔克和温克勒，1979）。感兴趣的读者可在GCDkit提示符下输入CIPW（不带括号！），查看该函数的代码。有关详细信息，请参阅CIPW, CIPWhb, Catanorm和Mesonorm的帮助文档。

尽管标准矿物（norms）在火成岩岩石学发展早期就已提出，但至今并未过时。例如，克罗斯等人，1902设计的CIPW标准矿物仍是火山岩TAS分类中的重要组成部分，用于区分某些岩石类型（利巴斯等人，1986；勒梅特，2002）。该计算过程涉及一系列层级规则，这些规则往往存在歧义，并在不同程度上导致不同的结果（哈奇森，1974，1975；Verma等，2002，2003）。CIPW标准矿物的另一个缺点是未包含含水矿物，因此其得出的矿物相通常与所研究火成岩（尤其是酸性岩）的实际矿物学存在差异。

处理GCDkit计算的结果

最新计算的值始终存储在变量 results 中。结果（即同名变量）可以复制到剪贴板，可追加到数据中（追加到数据矩阵WR），或通过右键点击R-Console窗口后出现的菜单保存为多种格式。

2.3 统计

在初步解释新获得的地球化学数据集时，检查所选元素或氧化物的描述性统计是很有帮助的。R语言包含了大量的统计工具，这些工具要么是内置的，要么以附加模块（包）的形式提供。然而，在此阶段，像mean、median、sd（标准差）和summary（统计概览）这样的简单函数已足够使用。此外，一些简单的图形工具也很有启发性，例如箱线图（箱须图；函数boxplot）和直方图（hist）。散点矩阵（pairs）可用于发现潜在的显著相关性。有关语法和更多细节，请参见附录B；对于包含多个数据组（使用因子）的更复杂数据集的特定问题，将在第2.4节中专门讨论。

GCDkit中的统计学
更复杂的工具超出了当前文本的范围，感兴趣的读者可参考R/S文档或专门出版物（例如，钱伯斯和黑斯蒂1992；韦纳布尔斯和里普利1999；梅因登阿尔德和布朗2003；Reimann 等 2008；van den Bogaard 和 Tolosana‐Delgado 2013）。

练习2.4：简单统计学

a) 计算文件 sazava.data中所有列（变量）的均值。
b) 绘制锶的箱线图，并找出表征其分布的所有主要统计参数（范围、中位数、观测数量和未测定情况等）。
c) 绘制以下氧化物所有可能组合的二元图（散点图矩阵）：SiO2, MgO, CaO, Na2O, K2O, 和 P2O5。

sazava <- read.table(“sazava.data”, sep=”\t”)
sazava <- sazava[, - (1:6)] # 或 sazava[,7:ncol(x)] 以仅获取数值数据

a) > result <- apply(sazava, 2, mean, na.rm=TRUE)

round(result, 2)
SiO2 TiO2 Al2O3 FeO
57.95 0.64 16.94 4.73
Fe2O3 MnO MgO CaO
1.75 0.14 3.57 8.16 …

b) > boxplot(sazava[,”Sr”], xlab=”Sr”, ylab=”ppm”)

summary(sazava[,”Sr”])
最小值 第一四分位数 中位数 均值 第三四分位数 最大值 NA的原子数
278.0 392.5 430.0 443.0 537.5 599.0 2.0

c) > oxides <- c(“SiO2”,”MgO”,”CaO”,”Na2O”,”K2O”,”P2O5”)

pairs(sazava[, oxides])

第2章 数据处理与简单计算（续）

2.4 分类与分组—使用因子

假设所研究的深成岩复合体由若干个火成岩系组成。那么，每个样品的全岩分析均可附带指明其所属的岩系。通过定义一个包含该分类信息的因子，便可以分别计算各个岩系的平均A/CNK值。我们应首先演示因子的定义方法，然后将其用于日益复杂的统计和分类任务。

2.4.1 复杂数据集的统计检验

对包含多个侵入体分析等复杂地球化学data集进行统计检查是繁琐的。幸运的是，R语言中的因子与函数tapply相结合，提供了一种非常灵活且优雅的解决方案。

练习 2.5：使用因子处理复杂数据集 I

a) 对于萨扎瓦数据集，根据“Intrusion”列中给出的说明定义一个因子 intrusion，将该岩石序列分为三组：basic（石英闪长岩至角闪辉长岩的多个小型侵入体）、Sazava（萨扎瓦侵入体主体：主要为角闪‐黑云花岗闪长岩至石英闪长岩）和Pozary（波扎里奥长花岗岩）。
b) 显示该因子的所有可能取值（水平）。
c) 利用因子intrusion，计算萨扎瓦数据集中各岩石组的 SiO2 平均含量。
d) 类似地，计算钡的平均浓度。

sazava <- read.table(“sazava.data”, sep=”\t”)

a) > intrusion <- factor(sazava[,”Intrusion”])

intrusion
[1] Sazava Sazava Sazava Sazava Sazava basic basic basic
[9] basic basic Pozary Pozary Pozary Pozary
Levels: basic Pozary Sazava

b) > levels(intrusion)
[1] “basic” “Pozary” “sazava”

c) > tapply(sazava[,”SiO2”], intrusion, mean)
basic Pozary sazava
51.778 68.440 55.738

d) > tapply(sazava[,”Ba”], intrusion, mean)
basic Pozary Sazava
NA 1291.25 NA

在上一个命令中，三个组中有两个组给出了NA，因为存在一些缺失值：

tapply(sazava[,”Ba”], intrusion, is.na)
$basic
[1] FALSE FALSE TRUE FALSE FALSE
$Pozary
[1] FALSE FALSE FALSE FALSE
$Sazava
[1] FALSE FALSE TRUE FALSE FALSE

如果要忽略缺失值并计算其余分析的均值，我们可以将参数 na.rm=TRUE 传递给函数 mean ：

tapply(sazava[,”Ba”], intrusion, mean, na.rm=TRUE)
basic Pozary sazava
676.25 1291.25 682.25

R语言提供了额外的、甚至可以说更强大的工具。例如， aggregate 会根据分组（由因子或因子列表定义）对数值矩阵或数据框x的每个变量（即列）应用指定的函数。类似地，函数 by 可根据因子（或因子列表）将数据框拆分为多个较小的数据框。

练习2.6：使用因子处理复杂数据集 II

a) 利用函数 summary ，计算萨扎瓦岩套各岩石组（因子 intrusion）中SiO2分布的基本统计参数。
b) 各侵入体中选定的微量元素（Ba、Rb、Sr和 Zr）的均值是多少？
c) 使用函数 by ，打印各岩石组中主要元素氧化物的基本统计摘要。

sazava <- read.table(“sazava.data”, sep=”\t”)
intrusion <- factor(sazava$”Intrusion”)

a) > tapply(sazava[,”SiO2”], intrusion, summary)
$basic
最小值 第一四分位数 中位数 均值 第三四分位数 最大值
48.84 49.63 51.72 51.78 52.90 55.80
$Pozary
最小值 第一四分位数 中位数 均值 第三四分位数 最大值
62.95 66.96 69.69 68.44 71.17 71.42
$Sazava
最小值 第一四分位数 中位数 均值 第三四分位数 最大值
50.72 55.09 55.17 55.74 57.73 59.98

b) > trace <- c(“Rb”,”Sr”,”Ba”,”Zr”)

aggregate(sazava[, trace], list(Rock = intrusion), mean, na.rm = TRUE)
Rock Rb Sr Ba Zr
1 basic 34.5 346.25 676.25 65.75
2 Pozary 25.0 595.00 1291.25 157.60
3 sazava 52.2 522.00 682.25 61.00

c) > by(sazava[, 7:17], list(Rock = intrusion), summary)
Rock: basic
SiO2 TiO2 Al2O3 …
最小值:48.84 最小值:0.340 最小值:13.34
第一四分位数:49.63 第一四分位数:0.670 第一四分位数:14.17
…
Rock: Pozary
SiO2 TiO2 Al2O3 …
最小值:62.95 最小值:0.280 最小值:15.09
第一四分位数:66.96 第一四分位数:0.300 第一四分位数:15.38
…
Rock: sazava
SiO2 TiO2 Al2O3 …
最小值:50.72 最小值:0.590 最小值:16.48
第一四分位数:55.09 第一四分位数:0.630 第一四分位数:16.94
…

GCDkit中的组

GCDkit 解决方案

GCDkit‐> loadData(“sazava.data”)
GCDkit‐> groupsByLabel(“Intrusion”)
Sa-1 Sa-2 Sa-3 Sa-4 Sa-7 SaD-1 “Sazava” “Sazava” “Sazava” “Sazava” “Sazava” “basic” …
已分配组：
basic Pozary sazava
5 4 5

复杂数据集中的统计学——GCDkit方式

两个函数为具有多个组的复杂数据集（如火成岩系）提供基本的统计参数（并可选择绘制直方图）： summarySingleByGroup （针对单个变量）和 summaryByGroup （多个元素/氧化物）。如果仅需要某些变量的范围，请使用函数 summaryRangesByGroup 。

继续前面的例子⋯⋯（输出省略）

GCDkit‐> summarySingleByGroup(“SiO2”)
GCDkit‐> trace <- c(“Rb”,”Sr”,”Ba”,”Zr”)
GCDkit‐> summaryByGroup(trace)
GCDkit‐> summaryRangesByGroup(trace)
GCDkit‐> summaryByGroup(major.elements)

2.4.2 数值向量转换为因子

函数 cut 将数值向量 x 切分为给定数量的区间，并根据其所属的秩对各个元素进行编码。因此，此函数可用于简单的分类目的。

练习 2.7：使用因子进行分类

根据萨扎瓦样品集中SiO2含量（wt.%）将样品分为四组：U (< 45)、B (45–52)、I (52–63) 和 A (> 63)。

sazava <- read.table(“sazava.data”, sep=”\t”)
silica <- cut(sazava[,”SiO2”], breaks = c(0, 45, 52, 63, 100),
labels = c(“U”, “B”, “I”, “A”))
silica
[1] I I I B I I B B B I I A A A
Levels: U B I A

注意，那些在数据中完全未出现的水平（此处为超基性岩石，U）不会被删除。

如果我们想了解各个样品的分类情况，可以将因子 silica 转换为字符向量：

acidity <- as.vector(silica)
names(acidity) <- rownames(sazava)
acidity
Sa-1 Sa-2 Sa-3 Sa-4 Sa-7 SaD-1 Gbs-1 Gbs-20
“I” “I” “I” “B” “I” “I” “B” “B”
Gbs-2 Gbs-3 Po-1 Po-3 Po-4 Po-5
“B” “I” “I” “A” “A” “A”

根据单个数值变量进行分组

类似的任务，即根据数值变量的值将样品分类到若干组中，可通过GCDkit函数 cutMy 实现。

GCDkit 解决方案

GCDkit‐> loadData(“sazava.data”)
GCDkit‐> cutMy(“SiO2”, c(0, 45, 52, 63, 100), c(“U”,”B”,”I”,”A”))
SiO2 interval
Sa-1 59.98 I
Sa-2 55.17 I
Sa-3 55.09 I
Sa-4 50.72 B
Sa-7 57.73 I
…

2.4.3 频数（列联）表

因子的一个巧妙应用是能够创建频数表。

练习2.8：频数表

承接上一个练习，我们演示如何制作频数表。

a) 使用因子Intrusion，统计萨扎瓦数据集中每个岩石组的分析次数。
b) 类似地，统计超基性、基性、中性和酸性岩石（来自上一练习的因子silica）的分析次数。
c) 建立一个频数表，显示二氧化硅含量与岩石类型之间的依赖关系。

intrusion <- factor(sazava[,”Intrusion”])

a) > table(intrusion)
intrusion
basic Pozary sazava
5 4 5

b) > table(silica)
silica
U B I A
0 4 7 3

c) > table(intrusion, silica)
silica
intrusion U B I A
basic 0 3 2 0
Pozary 0 0 1 3
sazava 0 1 4 0