R语言在分析化学教学中的应用

谭巧国

厦门大学环境与生态学院，福建 厦门 361102

分析化学由于其定量属性，涉及大量的计算和统计分析，并常需借助可视化手段展示复杂概念、过程和方法。因此，在分析化学的教学中，兼具计算和绘图功能的 Excel软件成为广泛采用的教学辅助工具。分析化学的经典教材也普遍结合Excel进行教学[1,2]，而且有专著专门讲解Excel在分析化学中的应用[3,4]。此外，也有文献零星地介绍其他计算机软件在分析化学中应用，所介绍软件包括Matlab[5-7]，Origin[8]等。

R语言是近年来兴起的一门计算机语言，在生物学、生态环境科学、经济学、社会科学等各领域都得到了广泛的应用，成为当前最受欢迎的编程语言之一[9,10]。R有比Excel更卓越的数据处理和可视化能力，理论上非常适用于分析化学的教学。然而，与其在科学研究中的热度形成鲜明对比，目前尚未见论文和书籍介绍将R用于包括分析化学在内的课程教学。本文列举数例，说明R可以在分析化学教学的多个方面得到应用。

1 R语言简介

R语言是一款开源免费的编程语言，主要用于统计分析和数据可视化。R最初是由两位统计学家罗斯·伊哈卡和罗伯特·杰特曼基于S语言开发而成。由于两位开发者名字的首字母都是R，R因此得名。目前R由“R开发核心团队”负责维护和开发[9]。

R的功能可以通过程序包得到扩展。例如本文中绘图用到了ggplot2程序包，数据整理用到了reshape2程序包，而这两个程序包与其他多个常用的程序包一起又包括在tidyverse程序包中。初次使用时，需运行以下代码安装：

install.packages("tidyverse")

除了以上提及的统计分析和绘图功能，R用于分析化学教学还具有以下几方面的优势：

R以代码形式实现计算和作图，操作过程完整地体现在代码中。代码具有自释性，理论上无需另外的解释性文字。代码可以纯文本的形式储存和传播，操作过程能够通过代码精确重现。与此相比，Excel计算过程的展示则必须用到图片(甚至视频)结合繁琐的步骤说明[11]，不仅占用大量储存空间，而且也不利于教学中的传播、分享和交流。

R有广阔的应用场景。R作为以数据分析和可视化见长的编程语言，是科学研究的有力工具。R借助扩展程序包，不断发展出新的功能满足科研需求。在分析化学的教学中使用R，为学习者接触掌握这项科研利器提供良好契机，可提高学习者的知识获得感和实际收益。

R软件免费，轻便。在写作本文时，最新版本R 4.0.0的安装文件仅84 MB。因此便于在不同的场合(如教室、师生的个人电脑等)随时安装使用R。与此相比，Matlab等收费软件不仅价格昂贵，而且需要5-8 GB的安装空间。

2 R应用于分析化学教学举例

2.1 R应用于分析化学计算

在学习化学平衡和滴定分析时，常涉及复杂的计算。例如，一元弱酸溶液 pH值的计算式是一元三次方程，多元弱酸和其他溶液还会涉及更高次的方程。高次方程不一定存在公式解，即使存在其公式也极为复杂。因此，直接解高次方程缺乏可操作性。目前教材中常采用忽略次要组分，忽略方程中的次要项等做法，以损失一定的准确度为代价，将高次方程简化成二次或一次方程进行求解[12]。简化中涉及多个判别式，判别式因溶液体系而异。这部分内容因此记忆负担较大，常成为学习分析化学的障碍。

R语言可看成一个功能全面的计算器，能执行各类科学计算，能用于解决分析化学中的高次方程求解问题。下面举一例说明求解过程。例题取自柳青和王海水[13]的论文，关于磷酸溶液pH值的计算。柳青和王海水提出了简化多元弱酸溶液 pH值计算的新方法，但思路仍然是将高次方程简化为二次或一次方程。其方法改善了通用性，但过程仍然较为复杂，仍有较大的记忆负担。

例 1 已知磷酸(H3A)的 Ka1= 7.6 × 10−3，Ka2= 6.3 × 10−8，Ka3= 4.4 × 10−13，计算 0.20 mol·L−1磷酸溶液的pH。

以上例题可通过执行以下 R代码进行计算。需说明的是#号后是代码注释，用于说明代码的用途，注释不会被作为代码执行。本文中的代码在关键处添加了注释，便于读者理解。

以上代码中，起核心作用的是解方程函数uniroot()。在使用uniroot()时，所输入的f是H+浓度精确计算式，因溶液体系而异；参数interval是粗略预估的方程根的上下限；参数tol是所允许的误差。运行以上代码后，得到结果：1.45134，因此可知溶液pH值为1.45，与参考文献中的计算值一致。

运用该方法，除了计算准确，还可以将所有酸碱溶液的 pH值计算都转变为同一种方法。计算过程可分为三步：(1) 根据质子得失平衡写出质子条件式；(2) 将式中酸(碱)各型体的浓度都替换成酸(碱)的分析浓度与分布分数之积，得到关于[H+]([OH−])的一元高次方程，即以上所用的f；(3) 解高次方程。

此方法只用到质子条件式和分布分数这两个容易掌握的基本方法，避免了复杂的判别式和简化过程，便于理解和记忆。从而使学习者从复杂的公式和繁琐的计算中解脱出来，将更多精力用于学习分析化学知识本身。

2.2 R应用于分析化学绘图

R语言另一项特色功能是科学绘图。图在分析化学的教学中发挥着不可或缺的作用，可以将复杂的概念、溶液体系、反应过程等可视化，降低理解难度。例如，分布分数图之于化学平衡，滴定曲线图之于滴定原理，都是极为有效的教学工具。学习者若能自己绘制这些图，则更能提升学习效果，增加学习乐趣。

以下以多元酸的分布分数图和络合滴定图为例，演示如何用将R应用于分析化学绘图，并帮助理解其背后的化学原理。需要说明的是，为了节省篇幅，以下作图的R代码仅保留主干部分，图的细节修饰所需代码未给出。读者运行以下代码，可看到图的主体信息，但格式细节与文中所示的图略有区别。完整作图代码获取方式见文末说明。

2.2.1 分布分数图的绘制

三元弱酸磷酸经三级解离，得到四种型体。根据磷酸的三个Ka值，可得到四种型体分布分数的表达式，并据此绘制分布分数图。以下是绘图的R代码：

运行以上代码，可得到磷酸的分布分数图(图1)。该图可用于判断不同pH值条件下，磷酸的主要型体。例如，天然海水的pH约为8.2，因此可判断海水中的磷酸盐主要以型体存在。该图亦可反向用于估算磷酸或磷酸盐溶液的pH值。例如，可估算KH2PO4和K2HPO4溶液的pH值分别在5和10左右。例1中的问题亦可使用该图帮助简化。假设0.20 mol∙L−1的H3PO4溶于水后第一级电离发生了50% (此比例据经验粗略估算，无需准确)，则溶液pH值为1左右。据图1判断，溶液中几乎不存在，因此可忽略H3PO4的第二、三级电离，可直接将磷酸当成一元弱酸来处理，大为简化了其pH值的计算。

图1 用R绘制磷酸的分布分数图

除了一元、多元酸碱的分布分数图，其他分布分数图，如Cu(NH3)n逐级络合物的分布分数图，EDTA各型体的分布分数图等，亦可按照以上方法绘制。学习者若改变作图参数(如Ka值)，便可即时看到图的变化。这样的过程可以帮助理解稳定常数、解离常数等如何决定了曲线的形状，进而理解这些常数的内涵。

2.2.2 滴定曲线的绘制

以下以EDTA络合滴定金属离子M为例，演示如何用R绘制络合滴定曲线。络合滴定曲线是滴定过程中金属离子浓度(pM')与滴定分数(a)之间关系的曲线，由络合滴定曲线方程描述：

在此例中，假设用0.1000 mol∙L−1的EDTA滴定同浓度的金属离子M，EDTA络合M的条件稳定常数K'MY为1.0 × 1014L∙mol−1。绘制其滴定曲线图的R代码如下：

运行以上代码，得到络合滴定曲线(图2)。曲线中包含了络合滴定的多项信息，这些信息在曲线绘制过程中均得到求解。例如，化学计量点的pM'即滴定分数a等于1时的pM'，为7.65；滴定突跃范围是滴定分数a为(100.0 ± 0.1)%对应的pM'范围，为4.30-11.00。通过调整K'MY和cM的取值，观察曲线形状的变化，有助于理解络合稳定常数以及金属离子浓度对于滴定突跃的影响。

图2 用R绘制的络合滴定曲线

此例中所用的方程可以定量描述络合滴定的各个阶段。对于其他滴定，也能建立类似的方程用于绘制滴定曲线，英文教材中将这类方程称为“master equation”[14]。此外，还可将滴定过程分为滴定前、化学计量点前、化学计量点、化学计量点后四个阶段，将问题简化，分别计算 pM'和滴定分数的关系，并进而绘制曲线。滴定曲线绘制充分运用了化学平衡与滴定的知识，同时使用了R的解方程和绘图两项主要功能，适用于教学中的知识讲解，但也因此步骤较多。若在日常教学科研中需频繁绘制滴定曲线，还可所采用R扩展程序包“titrationCurves”一步完成曲线绘制[15]。

2.3 R应用于分析化学中的统计分析

定量分析中的误差不可避免，计算和评估误差则需要用到统计分析方法。统计分析的理论和方法较为抽象，对其理解有赖于扎实的数学和统计学基础，因此常成为教和学的难点。分析化学中统计分析的教学应有别于数学、统计学类课程，目标应当是准确理解，正确使用，可以淡化数学推导。R最初是由统计学家开发，其统计学功能齐全而严谨，且容易使用。将R运用于分析化学中统计方法的教学，不仅可以帮助学习者掌握数据处理的正确方法，还有助于更好地理解相关的理论和方法。

以下以t检验为例，演示如何将R用于分析化学中的统计分析，并结合可视化的方法帮助学习者理解统计学概念。例题取自武汉大学《分析化学(上册)》[12]。

例2 采用一种新方法测定基准明矾中铝的质量分数，9次测定结果为10.74%，10.77%，10.77%，10.77%，10.81%，10.82%，10.73%，10.86%，10.81%。已知明矾中铝含量的标准值(以理论值代替)为10.77%。采用新方法后，是否引起系统误差(置信度95%)？

教材上介绍的是常规方法：由于测定次数较少，判断随机误差符合t分布。计算t值得到1.43，查t值表得到临界值t(0.05,8)= 2.31。t值小于临界值，从而判断无显著系统误差。

运用R语言，无需查表，进行单样本t检验直接得到结果。具体操作如下：

输出结果如下(为节省篇幅，部分信息省略)：

结果显示p值为0.2631，大于显著性水平0.05，因而可判断新方法的测定值与理论值无显著差异，亦即新方法未导致显著的系统误差。结果中还给出了测定值的95%置信区间，即(10.755，10.819)，此区间内包含理论值10.77，也表明新方法无显著系统误差。

在此基础上，还可将以上结果可视化，显示9次测定值、置信区间、理论值在坐标轴上的位置(图3)，帮助学习者理解t检验及置信区间的含义。以下是绘制图3a的代码：

从图上可看到各次测定值的位置决定了置信区间的位置和宽度；置信区间若将理论值包含在内，则表明测定值的均值与理论值无显著差异。为了进一步理解，另外生成一组有系统误差的数据作为对比。将例题中的测定值均减去0.05，并针对这组新数据做t检验：

结果显示p值为0.04266，小于0.05；95%置信区间为(10.705, 10.768)，不包含理论值10.77。表明存在显著的系统误差(图3b)。若将显著性水平调为99%，再进行t检验：

图3 用R绘图对比测定值及其均值的95%置信区间与理论值

得到结果 p值仍为 0.04266，虽小于 0.05，但大于新采用的显著性水平 0.01；99%置信区间为(10.690, 10.783)，比95%置信区间更宽，包含理论值10.77，表明不存在显著系统误差。因此可知，关于是否存在系统误差的结论与置信度的选择有密切相关。

除以上的t检验外，分析化学教学中涉及正态分布、t分布、置信区间的估计、F检验、可疑值的取舍、回归分析等统计分析内容均可借助R进行教学。

2.4 R数值模拟用于解释复杂概念

R作为编程语言，有生成符合特定分布的随机数，进行随机抽样等功能，因而可以进行数值模拟，帮助解释分析化学中一些复杂的概念和方法，以下举一例说明。

例3 测量值m1= 0.00 ± 0.10，测量值m2= 2.00 ± 0.10，则m2与m1之差m的值和标准差是多少？

按照误差传递公式，可以算得m的标准差是(0.102+ 0.102)1/2，即0.14，因此m = 2.00 ± 0.14。误差传递计算并不难，但是误差传递公式的推导过程较为复杂，学习和理解的难度较高。因此，学习者虽然能够按照公式算出误差，但未必能理解误差传递的真正含义。借助R的数值模拟和可视化，可以帮助学习者克服这一困难。以下针对例题3，进行数值模拟，得到图4。

图4 用R做数值实验帮助理解误差传递的原理

从图4中可知，差值(m2− m1)的离散程度大于m1和m2，但是小于r2，而与r1接近。从R代码中可知，r2的标准差为m1和m2标准差之和(即0.20)；r1的标准差为0.1414，与误差公式的计算结果吻合。通过数值模让学习者“眼见为实”，验证了误差传递公式。类似的模拟还可用于其他概念的解释，如置信区间、标准偏差的计算公式等[16]。

3 结语

以上举例展示了R在分析化学教学中的应用。R兼具计算、绘图、统计、模拟等功能，因此，分析化学中凡是涉及这些功能的内容，不限于以上所举的例子，均可用R辅助教学。R通过代码实现各项功能，操作过程均体现在文本化的代码中，便于信息的传播交流。因此R不仅适用于课堂演示，还便于学习者课后自学。相比于Excel，R在初学时有较为陡峭的学习曲线。考虑到本科生通常在低年级学习过C语言等更为复杂的编程语言，学习R已经有了较好的基础。目前在高校和学术界有人数众多气氛活跃的R学习群体，互联网上有丰富且优质的R学习资料，这也为自学R营造了良好的环境。将R用于分析化学的教与学，也让学习者有机会在应用实践中学习R，反过来促进R的学习。从分析化学课上学到R的技能，也将有助于其他课程的学习和日后的科学研究。

本文中计算、绘图等所涉及R代码的完整版本均可在此链接(http://tinyurl.com/y9hh36zq)下载或通过联系作者获取。