摘要：发表于《Geochimica et Cosmochimica Acta》，题为“新元古代沉积地层中异常重黄铁矿的微量元素证据”，作者们通过机器学习算法，对大量黄铁矿微量元素数据库进行分析，以区分不同来源的黄铁矿，特别是在新元古代沉积地层中同沉积期和后沉积期热液
大家好，我是python_ml。
发表于《Geochimica et Cosmochimica Acta》，题为“新元古代沉积地层中异常重黄铁矿的微量元素证据”，作者们通过机器学习算法，对大量黄铁矿微量元素数据库进行分析，以区分不同来源的黄铁矿，特别是在新元古代沉积地层中同沉积期和后沉积期热液过程对黄铁矿的影响。研究结合详细的岩相分析，对新元古代两个沉积序列中的黄铁矿起源进行了有效性测试，强调了黄铁矿微量元素在解析其起源方面的重要性。
沉积黄铁矿长期以来一直被用作地球历史上海洋环境的档案。然而，为了捕捉可靠的古环境信号，我们需要首先评估沉积地层中的黄铁矿，因为它可能被后来的成岩和/或热液过程改变和掩盖。在此，我们在一个大的激光剥蚀电感耦合等离子体质谱（LA-ICP-MS）黄铁矿微量元素数据库上训练了两种监督机器学习算法，以区分不同来源的黄铁矿。分析验证了基于12种微量元素（Co、Ni、Cu、Zn、As、Mo、Ag、Sb、Te、Au、Tl和Pb）的共行为构建的两个模型可以准确预测黄铁矿的来源。
总之，这是首次应用监督机器学习来理解新元古代黄铁矿起源的研究。测试的两种模型——随机森林和XGBoost，在利用原位微量元素数据区分沉积和热液黄铁矿方面表现出色。我们的分析表明，新元古代黄铁矿的起源可能是多样的：一些黄铁矿可能与当代海洋环境有关，而另一些可能经历了热液改造。未来的研究需要进一步探讨异常硫同位素信号与硫生物地球化学循环之间的关系。对于 exceptional fossil preservation（异常化石保存），我们的分析揭示了比原先想象更为复杂的成岩历史。
1 数据读取和准备将训练数据集（training.csv）划分为训练集和样本内局部性测试数据集。从每个标签（沉积型（0）、同沉积期热液型（1）以及后沉积期热液型（2））的黄铁矿分析数据中随机选取 240 条，构成一个平衡的训练数据集。剩余的分析数据则组成了样本内局部性测试数据集。另一个独立的数据集（OSL_Test.csv）被用于样本外局部性测试。
这一部分主要是读取和准备数据:
导入需要的库
读取训练数据、域外测试数据和Xiao & Cui数据
将不同类型的数据(如沉积型、热液角砾岩型等)按照一定比例分为训练集和域内测试集
提取训练数据、域内外测试数据以及Xiao & Cui数据的特征和标签
对所有数据进行标准化处理
import pandas as pdimport matplotlib.pyplot as pltimport matplotlibimport numpy as npfrom sklearn.model_selection import train_test_splitfrom sklearn.preprocessing import StandardScalerfrom sklearn.ensemble import RandomForestClassifierfrom sklearn.model_selection import GridSearchCVfrom sklearn.metrics import classification_reportfrom sklearn.metrics import confusion_matrixfrom sklearn.metrics import ConfusionMatrixDisplayimport xgboost as xgb#确保matplotlib导出可编辑的文本matplotlib.rcParams['pdf.fonttype'] = 42#读取数据为dataframefp1 = 'Training.csv' # 训练数据文件路径fp2 = 'OSL_Test.csv' # 域外测试数据文件路径 fp3 = 'Xiao_Cui.csv' # Xiao和Cui数据文件路径Data = pd.read_csv(fp1) # 读取训练数据OSL_Test = pd.read_csv(fp2) # 读取域外测试数据XC = pd.read_csv(fp3) # 读取Xiao和Cui数据#将沉积数据分为训练集和域内测试集Sedimentary_data = Data[Data['Deposit style']=='Sedimentary'].copy # 获取沉积数据Sedimentary_deposit = np.unique(Sedimentary_data['Deposit']) # 获取沉积矿床名称ISL_Test_data = pd.DataFrame # 创建空的dataframe存储域内测试数据Train_data = pd.DataFrame # 创建空的dataframe存储训练数据#对每个矿床,随机选择15个数据作为训练数据for i in Sedimentary_deposit: Train, Test = train_test_split(Data[Data['Deposit']==i], train_size=15, random_state=0) Train_data = pd.concat([Train_data, Train]) ISL_Test_data = pd.concat([ISL_Test_data, Test])#将热液角砾岩数据分为训练集和域内测试集 Train, Test = train_test_split(Data[Data['Deposit']=='hydrothermal breccia Menninnie Dam'], train_size=30, random_state=0)Train, Test = train_test_split(Data[Data['Deposit']=='hydrothermal breccia Telephone Dam'], train_size=30, random_state=0)#将IOCG数据分为训练集和域内测试集Train, Test = train_test_split(Data[Data['Deposit']=='IOCG Manxman'], train_size=30, random_state=0)Train, Test = train_test_split(Data[Data['Deposit']=='IOCG Punt Hill'], train_size=30, random_state=0) #将造山金数据分为训练集和域内测试集Au_data = Data[Data['Deposit style']=='Orogenic Au'].copy # 获取造山金数据Au_deposit = np.unique(Au_data['Deposit']) # 获取造山金矿床名称#对每个矿床,随机选择4个数据作为训练数据 for i in Au_deposit: Train, Test = train_test_split(Data[Data['Deposit']==i], train_size=4, random_state=0)#将斑岩数据分为训练集和域内测试集Train, Test = train_test_split(Data[Data['Deposit']=='porphyry Chalkidiki'], train_size=30, random_state=0)Train, Test = train_test_split(Data[Data['Deposit']=='porphyry Cadia'], train_size=30, random_state=0)#将SEDEX数据分为训练集和域内测试集SEDEX_data = Data[Data['Deposit style']=='SEDEX'].copy # 获取SEDEX数据SEDEX_deposit = np.unique(SEDEX_data['Deposit']) # 获取SEDEX矿床名称#对每个矿床,随机选择24个数据作为训练数据for i in SEDEX_deposit: Train, Test = train_test_split(Data[Data['Deposit']==i], train_size=24, random_state=0)#将VHMS数据分为训练集和域内测试集 Train, Test = train_test_split(Data[Data['Deposit']=='VHMS Kutlular'], train_size=15, random_state=0)Train, Test = train_test_split(Data[Data['Deposit']=='VHMS Yaman-Kasy'], train_size=46, random_state=0)ISL_Test_data = pd.concat([ISL_Test_data, Test], ignore_index=True)Train_data = pd.concat([Train_data, Data[Data['Deposit']=='VHMS Golden Grove']])Train_data = pd.concat([Train_data, Data[Data['Deposit']=='VHMS Jaguar Mine']])Train_data = pd.concat([Train_data, Data[Data['Deposit']=='VHMS Kyzilkaya']])Train_data = pd.concat([Train_data, Data[Data['Deposit']=='VHMS Lahanos']]) Train_data = pd.concat([Train_data, Data[Data['Deposit']=='VHMS Scuddles']], ignore_index=True)#获取数据和标签X_train = Train_data.iloc[:, 3:].copy # 训练数据特征y_train = Train_data.iloc[:, 2].copy # 训练数据标签X_ISL_test = ISL_Test_data.iloc[:, 3:].copy # 域内测试数据特征 y_ISL_test = ISL_Test_data.iloc[:, 2].copy # 域内测试数据标签X_OSL_test = OSL_Test.iloc[:, 3:].copy # 域外测试数据特征y_OSL_test = OSL_Test.iloc[:, 2].copy # 域外测试数据标签XC_data = XC.iloc[:, 2:].copy # Xiao和Cui数据特征feature_names = list(X_train.columns) # 存储特征名称#标准化数据 scaler = StandardScaler scaler.fit(X_train) # 用X_train计算均值和标准差作为scaler#标准化所有数据,然后为特征添加列名X_train = pd.DataFrame(scaler.transform(X_train), columns = feature_names) X_ISL_test = pd.DataFrame(scaler.transform(X_ISL_test), columns = feature_names)X_OSL_test = pd.DataFrame(scaler.transform(X_OSL_test), columns = feature_names)XC_data = pd.DataFrame(scaler.transform(XC_data), columns = feature_names)2 随机森林模型的训练和评估#寻找最佳参数params = {'n_estimators': list(range(100,1100,100)), 'min_samples_leaf': list(range(1,10,1)), 'criterion': ['gini','entropy']}gs_rf = GridSearchCV(RandomForestClassifier(random_state = 0), param_grid = params, scoring = 'f1_macro')gs_rf.fit(X_train, y_train) # 在训练数据上进行网格搜索gs_rf.best_params_ # 输出最佳参数#训练随机森林clf = RandomForestClassifier(criterion = 'gini', n_estimators = 300, min_samples_leaf = 2, oob_score = True, random_state = 0)clf.fit(X_train, y_train) # 用最佳参数在训练数据上训练随机森林#用分类器预测域内测试数据y_ISL_pred = clf.predict(X_ISL_test) #域内测试数据的分类报告print(classification_report(y_ISL_test, y_ISL_pred))#域内测试数据的混淆矩阵cm_ISL = confusion_matrix(y_ISL_test, y_ISL_pred)cm_display = ConfusionMatrixDisplay(cm_ISL).plot(cmap='YlGnBu') plt.savefig("ConfusionMatrix_RD_ISL_Test.pdf", dpi=150) # 保存混淆矩阵图#计算不同特征的重要性(基于不纯度的均值减少)和标准差importances = clf.feature_importances_#计算标准差std = np.std([tree.feature_importances_ for tree in clf.estimators_], axis=0)clf_importances = pd.Series(importances, index = clf.feature_names_in_)fig, ax = plt.subplots clf_importances.plot.bar(yerr=std, ax=ax)ax.set_title("Feature importances using MDI")ax.set_ylabel("Mean decrease in impurity")fig.tight_layoutplt.savefig("Feature_Importance_RD.pdf", dpi=150) # 保存特征重要性图#用分类器预测域外测试数据y_OSL_pred = clf.predict(X_OSL_test)#域外测试数据的分类报告 print(classification_report(y_OSL_test, y_OSL_pred))#域外测试数据的混淆矩阵cm_OSL = confusion_matrix(y_OSL_test, y_OSL_pred)cm_OSL_display = ConfusionMatrixDisplay(cm_OSL).plot(cmap='YlGnBu')plt.savefig("ConfusionMatrix_RD_OSL_Test.pdf", dpi=150) # 保存混淆矩阵图#用分类器预测Xiao和Cui数据y_XC_RD_pred = clf.predict(XC_data)#复制Xiao和Cui数据XC_RD = XC.copy#在XC_RD中添加标签作为新列XC_RD['Label'] = y_XC_RD_pred#将预测概率转为三列的dataframe XC_RD_prob = pd.DataFrame(clf.predict_proba(XC_data), columns = ['0', '1', '2'])#合并XC_RD和概率dataframeXC_RD = pd.concat([XC_RD, XC_RD_prob], axis=1)XC_RD.to_csv("XC_RD.csv") # 保存预测结果到csv文件主要是随机森林模型的训练和评估:
用网格搜
索寻找随机森林的最佳参数
用最佳参数在训练数据上训练随机森林模型
用训练好的模型预测域内测试数据
,输出分类报告和混淆矩阵
计算特征重要性并绘图
用模型预测域外测试数据,输出分类报告和混淆矩阵
用模型预测Xiao & Cui数据,输出预测标签和概率
将Xiao & Cui数据的预测结果保存到csv文件
3 XGBoost模型的训练和评估
#将X_train转换为DMatrix对象dtrain = xgb.DMatrix(X_train, label=y_train, feature_names=feature_names)#设置初始参数并检查最优树的数量params = {'eta': 0.1, 'max_depth': 5,'min_child_weight': 1,'subsample': 0.8,'colsample_bytree': 0.8,'objective': 'multi:softprob','num_class': 3}xgb.cv(params, dtrain, num_boost_round=1000, early_stopping_rounds=100, metrics='auc', seed=0)#调整max_depth和min_child_weightparam_test1 = { 'max_depth': range(3,10,2), 'min_child_weight': range(1,6,2)}gsearcg1 = GridSearchCV(xgb.XGBClassifier(learning_rate =0.1, n_estimators=492, max_depth=5, min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8, objective= 'multi:softprob', num_class=3, seed=0), param_grid = param_test1, scoring='f1_macro')gsearcg1.fit(X_train, y_train)gsearcg1.best_params_, gsearcg1.best_score_#微调max_depth和min_child_weightparam_test2 = { 'max_depth': [1,2,3,4], 'min_child_weight': [0.5,1,1.5]}gsearcg2 = GridSearchCV(xgb.XGBClassifier(learning_rate =0.1, n_estimators=492, max_depth=5, objective= 'multi:softprob', num_class=3, seed=0), param_grid = param_test2, scoring='f1_macro')gsearcg2.fit(X_train, y_train)gsearcg2.best_params_, gsearcg2.best_score_#调整gammaparam_test3 = { 'gamma': np.arange(0,0.5,0.1)}gsearcg3 = GridSearchCV(xgb.XGBClassifier(learning_rate =0.1, n_estimators=492, max_depth=2, min_child_weight=0.5, gamma=0, subsample=0.8, colsample_bytree=0.8, objective= 'multi:softprob', num_class=3, seed=0), param_grid = param_test3, scoring='f1_macro')gsearcg3.fit(X_train, y_train)gsearcg3.best_params_, gsearcg3.best_score_#重新校准boosting轮数'max_depth': 2,'min_child_weight': 0.5,'gamma': 0.3,'subsample': 0.8,'num_class': 3}#调整subsample和colsample_bytreeparam_test4 = { 'subsample': np.arange(0.5,1,0.05), 'colsample_bytree': np.arange(0.5,1,0.05)}gsearcg4 = GridSearchCV(xgb.XGBClassifier(learning_rate =0.1, n_estimators=261, max_depth=2, min_child_weight=0.5, gamma=0.3, subsample=0.8, colsample_bytree=0.8, objective= 'multi:softprob', num_class=3, seed=0), param_grid = param_test4, scoring='f1_macro')gsearcg4.fit(X_train, y_train)gsearcg4.best_params_, gsearcg4.best_score_#选择learning_rate=0.1, n_estimators=300clf = xgb.XGBClassifier(learning_rate =0.1, n_estimators=300, max_depth=2, min_child_weight=0.5, gamma=0.3, subsample=0.65, colsample_bytree=0.5, objective= 'multi:softprob', num_class=3, seed=0, importance_type = 'gain')clf.fit(X_train, y_train)y_ISL_pred = clf.predict(X_ISL_test)from sklearn.metrics import classification_reportprint(classification_report(y_ISL_test, y_ISL_pred))#计算不同特征的重要性(基于增益)importances = clf.feature_importances_clf_importances = pd.Series(importances, index = clf.feature_names_in_)fig, ax = plt.subplotsclf_importances.plot.barax.set_title("Feature importances using gain") ax.set_ylabel("Gain")fig.tight_layoutplt.savefig("Feature_Importance_XG.pdf", dpi=150) # 保存特征重要性图#域内测试数据的混淆矩阵cm_ISL = confusion_matrix(y_ISL_test, y_ISL_pred)cm_display = ConfusionMatrixDisplay(cm_ISL).plot(cmap='YlGnBu')plt.savefig("ConfusionMatrix_XG_ISL_Test.pdf", dpi=150) # 保存混淆矩阵图#用分类器预测域外测试数据y_OSL_pred = clf.predict(X_OSL_test)#域外测试数据的混淆矩阵cm_OSL = confusion_matrix(y_OSL_test, y_OSL_pred)cm_OSL_display = ConfusionMatrixDisplay(cm_OSL).plot(cmap='YlGnBu')plt.savefig("ConfusionMatrix_XG_OSL_Test.pdf", dpi=150) # 保存混淆矩阵图print(classification_report(y_OSL_test, y_OSL_pred))#用分类器预测Xiao和Cui数据y_XC_XG_pred = clf.predict(XC_data)#复制Xiao和Cui数据XC_XG = XC.copy#在XC_XG中添加标签作为新列XC_XG['Label'] = y_XC_XG_pred#将预测概率转为三列的dataframeXC_XG_prob = pd.DataFrame(clf.predict_proba(XC_data), columns = ['0', '1', '2'])#合并XC_RD和概率dataframeXC_XG = pd.concat([XC_XG, XC_XG_prob], axis=1)XC_XG.to_csv("XC_XG.csv") # 保存预测结果到csv文件XGBoost模型的训练和评估:
将训练数据转换为DMatrix对象,设置初始参数
用网格搜索调整max_depth、min_child_weight、gamma、subsample、colsample_bytree等参数
用最佳参数在训练数据上训练XGBoost模型
用训练好的模型预测域内测试数据,输出分类报告和混淆矩阵
计算特征重要性(基于增益)并绘图
用模型预测域外测试数据,输出分类报告和混淆矩阵
用模型预测Xiao & Cui数据,输出预测标签和概率
将Xiao & Cui数据的预测结果保存到csv文件
4 XGBoost模型的训练和评估
1. 导入SHAP库并初始化：python
importshapshap.initjs
这行代码导入了SHAP库并初始化了JavaScript，JavaScript的初始化是为了能够在后续的可视化中显示交互式图形。2. 创建SHAP解释器对象：
python
explainer = shap.TreeExplainer(clf)这行代码使用训练好的clf（XGBoost分类器）创建了一个TreeExplainer对象。TreeExplainer专门为树模型（如XGBoost、LightGBM等）设计，能够高效计算每个特征的SHAP值。3. 计算训练数据的SHAP值：
python
shap_values = explainer.shap_values(X_train)这行代码计算了训练数据的SHAP值。shap_values是一个列表，包含每个类别的SHAP值。对于每个类别，shap_values[i]会给出一个矩阵，每一行对应一个样本，每一列对应一个特征。4. 绘制特征重要性条形图：
python
shap.summary_plot(shap_values, X_train, plot_type='bar')python
shap.dependence_plot('Mo95', shap_values[0], X_train, feature_names=feature_names)'Mo95', shap_values[1
'Cu65', shap_values[0
第一行显示了类别0下Mo95特征的依赖关系。
第二行显示了类别1下Mo95特征的依赖关系。
第三行显示了类别0下Cu65特征的依赖关系。
6. 单个样本的SHAP值分析：python
one_sample_shap_values = shap_values[0][0]one_sample_base_value = explainer.expected_value[0]
one_sample_data = X_train.iloc[0, :]
explanation = shap.Explanation(
values=one_sample_shap_values,
base_values=one_sample_base_value,
data=one_sample_data,
feature_names=X_train.columns
)
shap.plots.waterfall(explanation, max_display=10)
这段代码选择了第一个类别（shap_values[0]）中的第0个样本，计算该样本的SHAP值、基准值和原始特征值，并使用waterfall图进行可视化。one_sample_shap_values
是第一个样本的SHAP值。
one_sample_base_value是该样本的基准值，即模型的输出值（未考虑任何特征的情况下的预测值）。
one_sample_data是该样本的特征数据。
python
base_value_for_class_0 = explainer.expected_value[0]one_sample_shap_values = shap_values[0][0]
one_sample_features = X_train.iloc[0, :]
shap.plots.force(
base_value_for_class_0,
one_sample_shap_values,
features=one_sample_features,
)
红色表示特征推动预测值增加（正贡献）。
蓝色表示特征推动预测值减少（负贡献）。
8. 绘制多个样本的Force图：python
shap.plots.force(explainer.expected_value[0],
shap_values[0],
features=X_train,
matplotlib=False
)
python
shap.decision_plot(explainer.expected_value[0],
shap_values[0][::10],
X_train.columns
)
python
shap.summary_plot(shap_values[1], X_train, plot_type="dot")0"dot")
2"dot")
表示使用散点图。每个点对应一个样本。
这段代码主要使用SHAP库来解释XGBoost模型的预测过程。通过SHAP，我们能够深入理解模型如何基于输入特征做出分类决策，特别是在涉及复杂的非线性模型（如XGBoost）时。通过SHAP值和可视化图表，我们不仅能够评估哪些特征最重要，还能洞察每个特征对每个预测样本的贡献，甚至可以针对单个样本进行详细分析。这有助于模型的可解释性，提升信任度并为进一步的模型改进提供洞察。Githup项目地址:
来源：九焰山灰太狼
