引言
我们要对植物内元素含量进行测定，因此就直接上了ICP。实验中植物使用了微波消解法与干灰化法两种方法（具体试验步骤见同一文集的推文）。得到的机读数据如下。
library(xlsx)
library(readxl)
setwd("C:/Users/Administrator/Desktop/已办/土壤植物环境下/土壤植物环境分析实验-ICP/微波消解法")
shuju_ICP<-read_xls("20210512-微波消解.xls")
names(shuju_ICP)[1:2]<-c("处理","样品号")
shuju_ICP处理 样品号 B 249.677 Ca 317.933 Cu 327.393 Fe 238.204 K 766.490 Mg 285.213 Mn 257.610 Na 589.592 P 213.617 Zn 206.200
1 ck1 Sample011 0.005 mg/L 8.215 mg/L 0.064 mg/L 0.182 mg/L 0.714 mg/L 0.996 mg/L 0.005 mg/L 1.025 mg/L -0.145 mg/L 0.139 mg/L
2 ck2 Sample012 0.009 mg/L 5.202 mg/L 0.061 mg/L 0.168 mg/L 0.442 mg/L 0.479 mg/L 0.004 mg/L 1.050 mg/L -0.158 mg/L 0.040 mg/L
3 1 Sample013 0.072 mg/L 64.31 mg/L 0.151 mg/L 13.05 mg/L 156.4 mg/L 23.16 mg/L 0.620 mg/L 4.375 mg/L 32.23 mg/L 0.507 mg/L
4 2 Sample014 0.085 mg/L 70.85 mg/L 0.151 mg/L 5.475 mg/L 117.9 mg/L 18.36 mg/L 0.422 mg/L 2.577 mg/L 28.99 mg/L 0.538 mg/L
5 3 Sample015 0.072 mg/L 66.55 mg/L 0.054 mg/L 4.136 mg/L 119.3 mg/L 18.37 mg/L 0.419 mg/L 2.319 mg/L 31.56 mg/L 0.390 mg/L
6 4 Sample016 0.070 mg/L 61.66 mg/L 0.123 mg/L 13.34 mg/L 153.8 mg/L 23.78 mg/L 0.692 mg/L 3.008 mg/L 33.31 mg/L 0.479 mg/L
7 5 Sample017 0.073 mg/L 63.93 mg/L 0.126 mg/L 13.77 mg/L 156.5 mg/L 24.11 mg/L 0.634 mg/L 2.993 mg/L 32.91 mg/L 0.477 mg/L
8 6 Sample018 0.085 mg/L 72.04 mg/L 0.148 mg/L 5.476 mg/L 121.8 mg/L 18.36 mg/L 0.419 mg/L 4.261 mg/L 28.50 mg/L 0.621 mg/L
9 7 Sample019 0.072 mg/L 61.74 mg/L 0.140 mg/L 13.67 mg/L 156.7 mg/L 24.06 mg/L 0.591 mg/L 3.870 mg/L 31.92 mg/L 0.448 mg/L
10 8 Sample020 0.073 mg/L 61.38 mg/L 0.144 mg/L 12.64 mg/L 146.8 mg/L 22.76 mg/L 0.666 mg/L 3.046 mg/L 33.00 mg/L 0.537 mg/L
11 9 Sample021 0.090 mg/L 69.35 mg/L 0.101 mg/L 5.636 mg/L 124.6 mg/L 19.29 mg/L 0.433 mg/L 2.150 mg/L 29.62 mg/L 0.398 mg/L
12 10 Sample022 0.067 mg/L 69.30 mg/L 0.176 mg/L 13.61 mg/L 153.6 mg/L 24.18 mg/L 0.608 mg/L 3.862 mg/L 32.34 mg/L 0.503 mg/L
13 11 Sample023 0.068 mg/L 59.15 mg/L 0.110 mg/L 12.85 mg/L 152.8 mg/L 23.42 mg/L 0.604 mg/L 3.103 mg/L 32.49 mg/L 0.429 mg/L
14 12 Sample024 0.067 mg/L 61.64 mg/L 0.137 mg/L 13.21 mg/L 155.4 mg/L 23.95 mg/L 0.595 mg/L 3.238 mg/L 31.65 mg/L 0.462 mg/L
15 17 Sample025 0.086 mg/L 66.45 mg/L 0.112 mg/L 5.685 mg/L 122.8 mg/L 18.11 mg/L 0.429 mg/L 2.005 mg/L 28.42 mg/L 0.406 mg/L
16 18 Sample026 0.071 mg/L 75.88 mg/L 0.215 mg/L 13.46 mg/L 156.8 mg/L 24.84 mg/L 0.611 mg/L 6.337 mg/L 31.99 mg/L 0.513 mg/L
17 19 Sample027 0.069 mg/L 61.68 mg/L 0.123 mg/L 13.14 mg/L 158.0 mg/L 23.98 mg/L 0.629 mg/L 3.025 mg/L 33.54 mg/L 0.462 mg/L
18 20 Sample028 0.089 mg/L 78.73 mg/L 0.164 mg/L 5.250 mg/L 122.7 mg/L 18.97 mg/L 0.424 mg/L 3.248 mg/L 28.65 mg/L 0.511 mg/L
19 21 Sample029 0.071 mg/L 59.89 mg/L 0.117 mg/L 13.22 mg/L 151.1 mg/L 23.20 mg/L 0.635 mg/L 2.743 mg/L 32.82 mg/L 0.449 mg/L
20 22 Sample030 0.070 mg/L 66.57 mg/L 0.164 mg/L 12.67 mg/L 158.5 mg/L 24.35 mg/L 0.605 mg/L 3.423 mg/L 33.06 mg/L 0.455 mg/L
21 23 Sample031 0.075 mg/L 70.88 mg/L 0.195 mg/L 4.250 mg/L 134.5 mg/L 21.61 mg/L 0.432 mg/L 4.883 mg/L 34.66 mg/L 0.447 mg/L
22 24 Sample032 0.071 mg/L 70.40 mg/L 0.209 mg/L 13.25 mg/L 161.7 mg/L 25.22 mg/L 0.599 mg/L 7.407 mg/L 32.50 mg/L 0.470 mg/L
23 25 Sample033 0.074 mg/L 63.17 mg/L 0.121 mg/L 13.32 mg/L 161.8 mg/L 24.87 mg/L 0.590 mg/L 3.567 mg/L 32.09 mg/L 0.452 mg/L
24 26 Sample034 0.101 mg/L 69.16 mg/L 0.112 mg/L 5.264 mg/L 126.7 mg/L 20.21 mg/L 0.417 mg/L 2.975 mg/L 29.20 mg/L 0.411 mg/L数据提取
观察干灰化法与微波分解法的ICP原始数据，发现数据采用类似的格式，即为“数值+空格+ mL/L”的数据格式。根据此项数据格式，以空格划分数据，提取数值部分。
shuju2_ICP<-data.frame(shuju_ICP[1])
for(a in c(3:12)){
c<-c()
for(i in shuju_ICP[[a]]){
string<-strsplit(i,split = " ")
string<-unlist(string)[1]
string<-as.numeric(string)
c<-c(c,string)
}
shuju2_ICP<-data.frame(shuju2_ICP,c)
}
names(shuju2_ICP)[-1]<-names(shuju_ICP)[3:12]
print(shuju2_ICP)数据输出如下
> head(shuju2_ICP)
处理 B 249.677 Ca 317.933 Cu 327.393 Fe 238.204 K 766.490 Mg 285.213 Mn 257.610 Na 589.592 P 213.617 Zn 206.200
1 ck1 0.005 8.215 0.064 0.182 0.714 0.996 0.005 1.025 0.00 0.139
2 ck2 0.009 5.202 0.061 0.168 0.442 0.479 0.004 1.050 0.00 0.040
3 1 0.072 64.310 0.151 13.050 156.400 23.160 0.620 4.375 32.23 0.507
4 2 0.085 70.850 0.151 5.475 117.900 18.360 0.422 2.577 28.99 0.538
5 3 0.072 66.550 0.054 4.136 119.300 18.370 0.419 2.319 31.56 0.390
6 4 0.070 61.660 0.123 13.340 153.800 23.780 0.692 3.008 33.31 0.479观察数据的两个CK部分，发现存在负值的数据值，将负值的数据值人为调整为0。
for(a in c(1:2)){
for(b in c(2:11)){
if(shuju2_ICP[a,b]<0){
shuju2_ICP[a,b]<-0
}}}
shuju_CK<-shuju2_ICP[1:2,]
print(shuju_CK)输出CK组数据
> shuju_CK
处理 B 249.677 Ca 317.933 Cu 327.393 Fe 238.204 K 766.490 Mg 285.213 Mn 257.610 Na 589.592 P 213.617 Zn 206.200
1 ck1 0.005 8.215 0.064 0.182 0.714 0.996 0.005 1.025 0 0.139
2 ck2 0.009 5.202 0.061 0.168 0.442 0.479 0.004 1.050 0 0.040依照土植环实验的经验，将所有数据按如下方式进行分组：1组（1、4、7、10、19、22、25号）、2组（2、5、8、11、20、23、26号）、3组（3、6、9、12、21、24号）。其中以下的数据分析只使用微波消解法数据，原因在于缺少各组样品质量的数据，而干灰化法的称量处于0.3000-0.3500g的范围；微波消解的称量范围为0.2990-0.3010g，相比而言用微波消解法的溶液浓度数据当做植物元素浓度数据具有更小的误差并具有分析意义，因此，我们假设微波消解法各组样品的称样量为0.3000g，同一组内具有多个重复。
group_1<-seq(from=1,to=26,by=3)
group_2<-seq(from=2,to=26,by=3)
group_3<-seq(from=3,to=26,by=3)
shuju2_ICP<-within(shuju2_ICP[-1:-2,],{
组别<-NA
组别[处理%in% group_1] <- 1
组别[处理%in% group_2] <- 2
组别[处理%in% group_3 ]<- 3
})
shuju2_ICP<-data.frame(shuju2_ICP[12],shuju2_ICP[1],shuju2_ICP[2:11])
print(shuju2_ICP) 提取的样品数据根据1-12号样品扣除CK1,17-26号样品扣除CK2得到扣除空白的数据列表
shuju2_ICP$处理<-as.numeric(shuju2_ICP$处理)
for(b in c(3:12)){
for(a in c(1:22)){
if(shuju2_ICP[a,2]<=14){
shuju2_ICP[a,b]<-shuju2_ICP[a,b]-shuju_CK[1,b-1]
}else{
shuju2_ICP[a,b]<-shuju2_ICP[a,b]-shuju_CK[2,b-1]
}}}
print(shuju2_ICP)进行到这里，数据提取与整理完毕，如下显示最后提取的数据
> shuju2_ICP
组别 处理 B.249.677 Ca.317.933 Cu.327.393 Fe.238.204 K.766.490 Mg.285.213 Mn.257.610 Na.589.592 P.213.617 Zn.206.200
3 1 1 0.067 56.095 0.087 12.868 155.686 22.164 0.615 3.350 32.23 0.368
4 2 2 0.080 62.635 0.087 5.293 117.186 17.364 0.417 1.552 28.99 0.399
5 3 3 0.067 58.335 -0.010 3.954 118.586 17.374 0.414 1.294 31.56 0.251
6 1 4 0.065 53.445 0.059 13.158 153.086 22.784 0.687 1.983 33.31 0.340
7 2 5 0.068 55.715 0.062 13.588 155.786 23.114 0.629 1.968 32.91 0.338
8 3 6 0.080 63.825 0.084 5.294 121.086 17.364 0.414 3.236 28.50 0.482
9 1 7 0.067 53.525 0.076 13.488 155.986 23.064 0.586 2.845 31.92 0.309
10 2 8 0.068 53.165 0.080 12.458 146.086 21.764 0.661 2.021 33.00 0.398
11 3 9 0.085 61.135 0.037 5.454 123.886 18.294 0.428 1.125 29.62 0.259
12 1 10 0.062 61.085 0.112 13.428 152.886 23.184 0.603 2.837 32.34 0.364
13 2 11 0.063 50.935 0.046 12.668 152.086 22.424 0.599 2.078 32.49 0.290
14 3 12 0.062 53.425 0.073 13.028 154.686 22.954 0.590 2.213 31.65 0.323
15 2 17 0.077 61.248 0.051 5.517 122.358 17.631 0.425 0.955 28.42 0.366
16 3 18 0.062 70.678 0.154 13.292 156.358 24.361 0.607 5.287 31.99 0.473
17 1 19 0.060 56.478 0.062 12.972 157.558 23.501 0.625 1.975 33.54 0.422
18 2 20 0.080 73.528 0.103 5.082 122.258 18.491 0.420 2.198 28.65 0.471
19 3 21 0.062 54.688 0.056 13.052 150.658 22.721 0.631 1.693 32.82 0.409
20 1 22 0.061 61.368 0.103 12.502 158.058 23.871 0.601 2.373 33.06 0.415
21 2 23 0.066 65.678 0.134 4.082 134.058 21.131 0.428 3.833 34.66 0.407
22 3 24 0.062 65.198 0.148 13.082 161.258 24.741 0.595 6.357 32.50 0.430
23 1 25 0.065 57.968 0.060 13.152 161.358 24.391 0.586 2.517 32.09 0.412
24 2 26 0.092 63.958 0.051 5.096 126.258 19.731 0.413 1.925 29.20 0.371异常值的剔除与各组数据集中程度的衡量
数据异常值的剔除
这里我们使用箱型图来检测异常值，通过对微波消解法各组数据进行箱型图的绘制，以小于Q1-1.5IQR或大于Q3+1.5IQR作为异常值的划定标准，检测到四个异常值并在箱型图上进行标注。对于检测出的异常值，我们较为保守的进行了保留。
par(mfrow=c(2,5))
for(i in c(3:12)){
box<-boxplot(shuju2_ICP[[i]]~shuju2_ICP[[1]],data = shuju2_ICP,
main=names(shuju2_ICP)[i],
xlab = "组别",
ylab = "溶液浓度")
out.vals=box$out
for (a in out.vals) {
text(a,adj = -0.2,labels = a)
}
}
在Mn和Na组分别检测到了1个和2个异常值。我们发现报告异常值大多来自于Na的测定，我们猜测是由于Na的污染（如汗液）导致Na组异常值偏多。
通过箱型图可以得到在大部分的元素测定数据中，1组数据具有较其他两组数据更好的集中程度
正态性检验
library(car)
par(mfrow=c(2,2))
shuju2_ICP$组别<-as.character(shuju2_ICP$组别)
qqPlot(lm(shuju2_ICP$B.249.677~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[3],labels=F)
qqPlot(lm(shuju2_ICP$Ca.317.933~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[4],labels=F)
qqPlot(lm(shuju2_ICP$Cu.327.393~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[5],labels=F)
qqPlot(lm(shuju2_ICP$Fe.238.204~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[6],labels=F)
qqPlot(lm(shuju2_ICP$K.766.490~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[7],labels=F)
qqPlot(lm(shuju2_ICP$Mg.285.213~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[8],labels=F)
qqPlot(lm(shuju2_ICP$Mn.257.610~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[9],labels=F)
qqPlot(lm(shuju2_ICP$Na.589.592~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[10],labels=F)
qqPlot(lm(shuju2_ICP$P.213.617~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[11],labels=F)
qqPlot(lm(shuju2_ICP$Zn.206.200~shuju2_ICP$组别),
simulate = T,main=names(shuju2_ICP)[12],labels=F
方差齐次分析我们这里进行了偷懒，因为单因素方差分析对于方差是否齐次并不敏感，即使不通过方差齐次性检验，还是建议通过aov方法进行单因素方差分析
单因素方差分析
attach(shuju2_ICP)
data_result1<-aov(B.249.677~组别)
summary(data_result1)
data_result2<-aov(Ca.317.933~组别)
summary(data_result2)
data_result3<-aov(Cu.327.393~组别)
summary(data_result3)
data_result4<-aov(Fe.238.204~组别)
summary(data_result4)
data_result5<-aov(K.766.490~组别)
summary(data_result5)
data_result6<-aov(Mg.285.213~组别)
summary(data_result6)
data_result7<-aov(Mn.257.610~组别)
summary(data_result7)
data_result8<-aov(Na.589.592~组别)
summary(data_result8)
data_result9<-aov(P.213.617~组别)
summary(data_result9)
data_result10<-aov(Zn.206.200~组别)
summary(data_result10)| 各元素数据方差分析结果（只显示最终结果） | B | Ca | Cu | Fe | K | Mg | Mn | Na | P | Zn | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Pf(>F) | 0.0736 | 0.391 | 0.988 | 0.0348 | 0.0217 | 0.0562 | 0.0572 | 0.37 | 0.202 | 0.986 | |
| * | * |
备注：标注*为差异显著（95%）；标注**为差异极显著（99%）
多重比较
使用TukeyHSD方法
par(mfrow=c(1,2))
tuk_Fe<-aov(shuju2_ICP$Fe.238.204~shuju2_ICP$组别)
tuk_Fe<-TukeyHSD(tuk_Fe)
plot(tuk_Fe)
tuk_k<-aov(shuju2_ICP$K.766.490~shuju2_ICP$组别)
tuk_k<-TukeyHSD(tuk_k)
plot(tuk_k)