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

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

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

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

例1  已知磷酸(H₃A)的K_a1 = 7.6 × 10⁻³，K_a2 = 6.3 × 10⁻⁸，K_a3= 4.4 × 10⁻¹³，计算0.20 mol∙L⁻¹磷酸溶液的pH。

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

#输入已知条件

Ka1 = 7.6E-3

Ka2 = 6.3E-8

Ka3 = 4.4E-13

c = 0.20

#定义三元酸H+浓度的精确计算式

f = function(H) {

    #将分母命名为D

    D = H^3 + H^2*Ka1 + H*Ka1*Ka2 + Ka1*Ka2*Ka3

    H - 1E-14/H - (c*H^2*Ka1 + c*2*H*Ka1*Ka2 + c*3*Ka1*Ka2*Ka3)/D

}

#运用uniroot函数解方程

H = uniroot(f, interval = c(0, 0.6), tol = 10^-20)$root

#根据H+浓度计算

pH = -log10(H)

pH

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

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

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

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

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

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

library(ggplot2) #加载用于绘图的程序包

Ka1 = 7.6E-3

Ka2 = 6.3E-8

Ka3 = 4.4E-13

#选取0~14范围内的100个pH值

pH = seq(from = 0, to = 14, length.out = 100)

H = 10^(-pH)

#根据H+浓度计算分布分数的分母

D = H^3 + H^2*Ka1 + H*Ka1*Ka2 + Ka1*Ka2*Ka3

#计算磷酸各型体的分布分数

delta_0 = Ka1*Ka2*Ka3/D

delta_1 = H*Ka1*Ka2/D

delta_2 = H^2*Ka1/D

delta_3 = H^3/D

#将pH值和各型体分布分数组合成数据表，用于绘图

df = data.frame(pH, delta_0, delta_1, delta_2, delta_3)

#绘制分布分数图

ggplot(df, aes(x = pH))+ #x轴为pH

    geom_line(aes(y = delta_0), linetype = 1)+ #PO4对应的曲线

    geom_line(aes(y = delta_1), linetype = 2)+ #HPO4对应的曲线

    geom_line(aes(y = delta_2), linetype = 3)+ #H2PO4对应的曲线

    geom_line(aes(y = delta_3), linetype = 4) #H3PO4对应的曲线

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

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

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

(1)$ {K}_{\mathrm{M}\mathrm{Y}}^{\text{'}}·{\left[\mathrm{M}\mathrm{\text{'}}\right]}^{2}+[{K}_{\mathrm{M}\mathrm{Y}}^{\text{'}}·{c}_{\mathrm{M}}·\left(a-1\right)+1]·\left[\mathrm{M}\mathrm{\text{'}}\right]-{c}_{\mathrm{M}}=0 $

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

#设定金属离子和EDTA溶液的分析浓度

cM = 0.1

#设定条件稳定常数值

KMY = 1E14

#设定滴定分数，在0~2之间取200个点

a_seq = seq(0, 2, length.out = 200)

#设定名为pM，长度为200的向量，用于储存计算所得的pM值

pM = numeric(200)

#针对每一个滴定分数值a，解络合滴定方程计算对应的pM值，此处运用循环语句

for (i in 1:200) {

    #定义络合滴定曲线方程

    f = function(x) {

    #将已知条件导入方程

    a < < - a_seq[i]

    KMY < < - KMY

    cM < < - cM

    KMY*x^2+(KMY*cM/(1+a)*(a-1)+1)*x-cM/(1+a)

    }

    #将方程的根存入名为M的变量

    M = uniroot(f, interval = c(1E-20, 0.1), tol=1E-20)$root

    #计算pM值并保存

    pM[i] = -log10(M)

}

#将滴定分数值与对应的pM值合并成数据表，命名为df

df = data.frame(a = a_seq, pM = pM)

#调用作图所需软件包ggplot2

library(ggplot2)

#作图

ggplot(df, aes(x = a, y = pM))+

    geom_line()+

    labs(x = "滴定分数", y = "pM'")

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

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

定量分析中的误差不可避免，计算和评估误差则需要用到统计分析方法。统计分析的理论和方法较为抽象，对其理解有赖于扎实的数学和统计学基础，因此常成为教和学的难点。分析化学中统计分析的教学应有别于数学、统计学类课程，目标应当是准确理解，正确使用，可以淡化数学推导。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检验直接得到结果。具体操作如下：

#将9次测定值存入变量x

x = c(10.74, 10.77, 10.77, 10.77, 10.81, 10.82, 10.73, 10.86, 10.81)

#用t.test函数进行t检验，参数mu是理论值，参数conf.level是置信度

t.test(x, mu = 10.77, conf.level = 0.95)

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

##            One Sample t-test

## t = 1.2039, df = 8, p-value = 0.2631

## alternative hypothesis: true mean is not equal to 10.77

## 95 percent confidence interval:

## 10.75474 10.81859

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

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

library(ggplot2)

ggplot(data.frame(x), aes(x, y=0))+

    geom_rect(xmin=10.75474, xmax=10.81859, ymin=-0.01, ymax=0.01)+ #置信区间

    geom_point(alpha = 1/3)+ #各测定值对应的点

    geom_text(label=x, nudge_y=0.01)+ #各测定值

    geom_point(x=10.77, y=0, shape= 1, size =4)+ #理论值

    ylim(-0.1, 0.1)

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

x_2 = x - 0.05 #原数据各减去0.05

t.test(x_2, mu = 10.77, conf.level = 0.95)

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

t.test(x_2, mu = 10.77, conf.level = 0.99)

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

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

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

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

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

library(ggplot2)

library(reshape2)

#随机生成4组包含100个符合正态分布的数，参数mean为均值，sd为标准偏差

m1 = rnorm(100, mean = 0, sd = 0.10) #均值为0，标准差为0.10

m2 = rnorm(100, mean = 2, sd = 0.10) #均值为2，标准差为0.10

r1 = rnorm(100, mean = 2, sd = 0.20) #均值为2，标准差为0.20

r2 = rnorm(100, mean = 2, sd = 0.1414) #均值为2，标准差为0.1414

#计算m2与m1之差

m = m2 - m1

#运用data.frame函数将上述5组数合组合成数据表d

d = data.frame(m1, m2, r1, r2, m)

#运用gather函数排列的数据表d的转化成纵向排列的数据表d_1

d1 < - melt(d)

#运用ggplot函数作图

ggplot(d1, aes(x = variable, y = value))+

    geom_violin()+

    geom_jitter(width = 0.1, shape = 1)

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

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

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