TOC (Total Organic Carbon, 总有机碳)：用于评估水体中有机物的总量，广泛应用于饮用水处理、工业废水处理以及自然水体质量监测。
COD (Chemical Oxygen Demand, 化学需氧量)：测量水中可被氧化的物质（主要是有机物）的数量，是衡量水污染程度的重要指标之一，常用于工业废水和城市污水处理效果的评价。
电导率：反映水中离子浓度的间接指标，可用于判断水的纯度或含盐量，适用于饮用水、农业灌溉用水、工业用水的质量控制。
pH：表示溶液酸碱性强弱的程度，对于维持生态系统平衡至关重要，广泛应用于水质检测及各类水处理过程中的监控。
盐度：指水中溶解盐类的总量，通常用来描述海水或其他高盐度水体的特性，在海洋学研究、水产养殖等领域有重要应用。
Na, K, Ba, Sr, Ca, Li, Fe, Mg, B, Al, Mn：这些都是常见的金属元素或非金属元素（硼），它们的存在形式和浓度可以影响水质的安全性和适用性，比如饮用水安全标准、工业废水排放限制等。
阳离子、阴离子：分别指的是带正电荷和负电荷的离子，通过分析特定的阴阳离子可以帮助了解污染物来源及其对环境的影响。
F-, Cl-, NO2-, Br-, NO3-, SO42-：具体的阴离子，这些离子的存在与否及其浓度可以提供关于水源污染状况的重要信息，例如硝酸盐含量过高可能指示农业面源污染。
26PAHs (多环芳烃)、17PESTs (农药)、24ZGWT (重金属)、8NPEOs (壬基酚聚氧乙烯醚)、8 PAEs (邻苯二甲酸酯)、16 QACs (季铵化合物)、BPA (双酚A)、有机物：这些都是特定类型的有机污染物或者化学物质，它们可能来源于工业排放、农业活动或是日常生活用品。检测这些物质有助于评估环境污染水平，并采取相应的环境保护措施。
水质毒性检测数据预处理和分析源码：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
from datetime import datetime
from scipy.interpolate import CubicSpline
from scipy import stats
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity, calculate_kmo
from matplotlib.patches import Ellipse
from statsmodels.tsa.stattools import acf
import matplotlib.transforms as transforms
# 设置支持中文的字体
plt.rcParams['font.sans-serif'] = ['WenQuanYi Micro Hei'] # 使用黑体
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
# 创建保存结果的文件夹
current_time = datetime.now().strftime('%Y-%m-%d_%H-%M-%S')
output_folder = f'分析结果_{current_time}'
os.makedirs(output_folder, exist_ok=True)
# 创建子文件夹
sub_folders = [
'异常值箱型图',
'折线对比图',
'自相关对比图',
'分布直方图对比',
'QQ图',
'PCA相关分析结果图' # 添加新的子文件夹
]
for folder in sub_folders:
os.makedirs(os.path.join(output_folder, folder), exist_ok=True)
outliers_boxplot_folder = os.path.join(output_folder, '异常值箱型图')
line_chart_folder = os.path.join(output_folder, '折线对比图')
autocorrelation_plot_folder = os.path.join(output_folder, '自相关对比图')
histogram_folder = os.path.join(output_folder, '分布直方图对比')
qq_plot_folder = os.path.join(output_folder, 'QQ图')
pca_analysis_folder = os.path.join(output_folder, 'PCA相关分析结果图') # 新的子文件夹路径
print(f"结果将保存到文件夹: {output_folder}")
print(f"异常值箱型图将保存到子文件夹: {outliers_boxplot_folder}")
print(f"折线对比图将保存到子文件夹: {line_chart_folder}")
print(f"自相关对比图将保存到子文件夹: {autocorrelation_plot_folder}")
print(f"分布直方图将保存到子文件夹: {histogram_folder}")
print(f"QQ图将保存到子文件夹: {qq_plot_folder}")
print(f"PCA相关分析结果图将保存到子文件夹: {pca_analysis_folder}")
# 从 Excel 文件中读取指定列的数据
input_file = "toxicity_prediction.xlsx"
# 检查文件是否存在
if not os.path.exists(input_file):
raise FileNotFoundError(f"文件 {input_file} 未找到，请检查文件路径和文件名。")
df = pd.read_excel(input_file)
print("\n原始数据预览：")
print(df.head())
# 数据预处理（先对原始数据的缺失值、异常值进行检验、插值，再插入空白日期进行插值，最后再进行一次异常值检验）
# 处理原始数据缺失值和异常值
print("\n数据预处理开始...")
# 处理原始数据缺失值
print("\n处理原始数据缺失值...")
# 检查日期列是否存在
if '日期' not in df.columns:
raise ValueError("输入数据中缺少 '日期' 列。请检查文件内容和列名。")
df['日期'] = pd.to_datetime(df['日期'], format='%y.%m.%d') # 将日期列转换为日期时间类型
df.set_index('日期', inplace=True) # 将日期列设置为索引
# 使用三次样条插值填充缺失值
def cubic_spline_interpolation(data):
"""
使用三次样条插值填充缺失值
:param data: 输入数据（DataFrame）
:return: 插值后的数据（DataFrame）
"""
for col in data.columns:
if data[col].isnull().sum() > 0: # 检查是否有缺失值
valid_idx = data[col].notnull() # 获取非缺失值的索引
if valid_idx.sum() >= 2: # 至少需要两个非缺失值才能进行插值
x = data.index[valid_idx].astype(np.int64) // 10**9 # 将时间戳转换为秒
y = data.loc[valid_idx, col]
cs = CubicSpline(x, y) # 创建三次样条插值对象
x_missing = data.index[~valid_idx].astype(np.int64) // 10**9 # 缺失值的时间戳
data.loc[~valid_idx, col] = cs(x_missing) # 填充缺失值
return data
# 执行三次样条插值
df_missing_filled = cubic_spline_interpolation(df.copy()) # 复制一份数据并填充缺失值
print("\n原始数据缺失值处理后数据预览：")
print(df_missing_filled.head())
# 处理原始数据异常值（基于箱线图的IQR方法）
print("\n原始数据处理异常值...")
def detect_outliers_iqr(data, threshold=2.0):
"""
基于IQR方法检测异常值
:param data: 输入数据（DataFrame或Series）
:param threshold: IQR倍数阈值（默认2.0）
:return: 异常值掩码（True表示异常值）
"""
if data.empty:
print("输入数据为空，无法检测异常值。")
return pd.DataFrame()
Q1 = data.quantile(0.25)
Q3 = data.quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - threshold * IQR
upper_bound = Q3 + threshold * IQR
return (data < lower_bound) | (data > upper_bound)
# 检测异常值
outliers_mask = detect_outliers_iqr(df.select_dtypes(include=[np.number]), threshold=2.0)
if outliers_mask.empty:
print("没有检测到异常值。")
else:
df_outliers_marked = df.copy() # 复制一份原始数据用于标记异常值
df_outliers_marked[outliers_mask] = np.nan # 将异常值标记为NaN
# 先对最尾的缺失值使用前向填充
for col in df_outliers_marked.columns:
if df_outliers_marked[col].isnull().sum() > 0: # 检查是否有缺失值
# 找到最尾的缺失值位置
last_valid_index = df_outliers_marked[col].last_valid_index()
if last_valid_index is not None:
# 使用 loc 索引避免链式索引
df_outliers_marked.loc[last_valid_index:, col] = df_outliers_marked.loc[last_valid_index:, col].ffill()
# 使用三次样条插值填充剩余的缺失值
df_outliers_filled = cubic_spline_interpolation(df_outliers_marked.copy()) # 复制一份数据并填充异常值
print("\n原始数据异常值处理后数据预览（最尾缺失值前向填充后，剩余缺失值三次样条插值）：")
print(df_outliers_filled.head())
# 绘制箱线图展示原始数据异常值分布（每个指标单独绘制）
print("\n绘制原始数据箱线图展示异常值分布（每个指标单独绘制）...")
plt.figure(figsize=(12, 8))
# 获取数值列
numeric_columns = df.select_dtypes(include=[np.number]).columns
# 检测原始数据的异常值（包括IQR和负值）
outliers_mask_original = (df.select_dtypes(include=[np.number]) < df.select_dtypes(include=[np.number]).quantile(
0.25) - 2 * (df.select_dtypes(include=[np.number]).quantile(0.75) - df.select_dtypes(include=[np.number]).quantile(
0.25))) | \
(df.select_dtypes(include=[np.number]) > df.select_dtypes(include=[np.number]).quantile(
0.75) + 2 * (df.select_dtypes(include=[np.number]).quantile(0.75) - df.select_dtypes(
include=[np.number]).quantile(0.25))) | \
(df.select_dtypes(include=[np.number]) < 0)
# 绘制每个指标的箱线图
for i, column in enumerate(numeric_columns, 1):
plt.subplot(5, 7, i) # 5行7列的子图布局
# 绘制箱线图（统一颜色为lightblue）
sns.boxplot(y=df[column], whis=2, color='lightblue')
# 标记异常值（统一颜色为红色）
outliers = df[column][outliers_mask_original[column]]
if not outliers.empty:
plt.scatter(x=[0] * len(outliers), y=outliers, color='red', marker='o', facecolors='none', edgecolors='red',linewidths=1.5, zorder=5)
plt.title(column, fontsize=12)
plt.suptitle('原始数据异常值分布箱线图（分指标显示）', fontsize=14)
plt.tight_layout()
# 保存图片到异常值箱型图的子文件夹
output_file_path = os.path.join(outliers_boxplot_folder, '0_原始数据异常值分布箱线图_分指标显示.png')
plt.savefig(output_file_path)
print(f"箱线图已保存到路径: {output_file_path}")
plt.close()
# 插入空白日期
print("\n插入空白日期...")
# 获取数据的开始和结束日期
start_date = df_outliers_filled.index.min()
end_date = df_outliers_filled.index.max()
# 获取所有需要的日期范围
date_range = pd.date_range(start=start_date, end=end_date, freq='D')
# 将原来的日期索引转换为 DataFrame 以便操作
df_temp = df_outliers_filled.reset_index()
df_temp.columns = ['日期'] + df_temp.columns[1:].tolist()
# 合并完整的日期范围和原始数据
df_full_date = pd.DataFrame(date_range, columns=['日期'])
df_merged = pd.merge(df_full_date, df_temp, on='日期', how='left')
# 将日期列设置为索引
df_merged.set_index('日期', inplace=True)
# 对空白日期进行三次样条插补
df_final = df_merged.copy() # 复制一份数据用于插补
# 对每一列进行三次样条插值
for col in df_final.columns:
if df_final[col].isnull().sum() > 0: # 只对有缺失值的列进行插补
# 使用三次样条插值（B-spline，order=3）
df_final[col] = df_final[col].interpolate(method='spline', order=3, limit_direction='both', s=0)
print("\n插入空白日期数据填充预览：")
print(df_final.head())
# 空白日期填充数据进行异常值检验
print("\n再次进行空白日期填充数据异常值检验...")
# 使用IQR方法检测异常值
outliers_mask_final = detect_outliers_iqr(df_final.select_dtypes(include=[np.number]), threshold=2.0)
# 检查填充值是否小于零，并将其标记为异常值
negative_values_mask = df_final.select_dtypes(include=[np.number]) < 0
# 将小于零的值也标记为异常值
outliers_mask_final = outliers_mask_final | negative_values_mask
# 检查是否存在异常值
if outliers_mask_final.empty:
print("无异常值，数据质量良好。")
else:
print("检测到异常值，进行标注和填充处理...")
# 复制一份数据用于标记异常值
df_final_outliers_marked = df_final.copy()
df_final_outliers_marked[outliers_mask_final] = np.nan # 将异常值标记为NaN
# 使用指数移动平均（EMA）填充缺失值
def ema_fill(data, span=3):
"""
使用指数移动平均（EMA）填充缺失值
:param data: 输入数据（DataFrame）
:param span: EMA的跨度（默认3）
:return: 填充后的数据（DataFrame）
"""
for col in data.columns:
if data[col].isnull().sum() > 0: # 检查是否有缺失值
# 计算指数移动平均值
ema_values = data[col].ewm(span=span, adjust=False).mean()
# 仅用EMA值填充缺失值的位置
data[col] = data[col].fillna(ema_values) # 直接赋值，避免链式赋值
return data
# 执行指数移动平均填充
df_final_outliers_filled = ema_fill(df_final_outliers_marked.copy()) # 复制一份数据并填充异常值
# 使用前向填充和后向填充处理剩余的缺失值（最头和最尾的缺失值）
df_final_outliers_filled = df_final_outliers_filled.ffill() # 前向填充
df_final_outliers_filled = df_final_outliers_filled.bfill() # 后向填充
print("\n异常值处理后数据预览（指数移动平均填充后，并补充前向和后向填充）：")
print(df_final_outliers_filled.head())
# 绘制箱线图展示空白日期填充数据异常值检验后的异常值分布（每个指标单独绘制）
print("\n绘制箱线图展示空白日期填充数据异常值检验后的异常值分布（每个指标单独绘制）...")
plt.figure(figsize=(12, 8))
# 获取数值列
numeric_columns = df_final.select_dtypes(include=[np.number]).columns
# 绘制每个指标的箱线图
for i, column in enumerate(numeric_columns, 1):
plt.subplot(5, 7, i) # 5行7列的子图布局
# 绘制箱线图（统一颜色为lightblue）
sns.boxplot(y=df_final[column], whis=2, color='lightblue')
# 标记异常值（统一颜色为红色）
outliers = df_final[column][outliers_mask_final[column]]
if not outliers.empty:
plt.scatter(x=[0] * len(outliers), y=outliers, color='red', marker='o', facecolors='none', edgecolors='red',linewidths=1.5, zorder=5)
plt.title(column, fontsize=12)
plt.suptitle('空白日期填充数据异常值检验后的异常值分布箱线图（分指标显示）', fontsize=14)
plt.tight_layout()
# 保存图片到异常值箱型图的子文件夹
output_file_path = os.path.join(outliers_boxplot_folder,'1_空白日期填充数据异常值检验后的异常值分布箱线图_分指标显示.png')
plt.savefig(output_file_path)
print(f"箱线图已保存到路径: {output_file_path}")
plt.close()
# 填充数据评估
# 统计指标对比
print("\n开始进行统计指标对比...")
def compare_filled_vs_benchmark(benchmark_data, filled_data):
"""
对比填充数据与基准数据的统计指标（针对所有列）
:param benchmark_data: 基准数据（异常值处理后的数据）
:param filled_data: 填充后的数据
"""
benchmark_data = benchmark_data.reindex(filled_data.index)
filled_data = filled_data.reindex(benchmark_data.index)
# 获取所有列名
columns = benchmark_data.columns
for column in columns:
# 均值
benchmark_mean = benchmark_data[column].mean()
filled_mean = filled_data[column].mean()
# 标准差
benchmark_std = benchmark_data[column].std()
filled_std = filled_data[column].std()
# 中位数
benchmark_median = benchmark_data[column].median()
filled_median = filled_data[column].median()
# 众数
try:
benchmark_mode = benchmark_data[column].mode().iloc[0]
except IndexError:
benchmark_mode = np.nan
try:
filled_mode = filled_data[column].mode().iloc[0]
except IndexError:
filled_mode = np.nan
# 最大值和最小值
benchmark_max = benchmark_data[column].max()
filled_max = filled_data[column].max()
benchmark_min = benchmark_data[column].min()
filled_min = filled_data[column].min()
# 偏度
try:
benchmark_skew = stats.skew(benchmark_data[column])
except:
benchmark_skew = np.nan
try:
filled_skew = stats.skew(filled_data[column])
except:
filled_skew = np.nan
# 峰度
try:
benchmark_kurtosis = stats.kurtosis(benchmark_data[column])
except:
benchmark_kurtosis = np.nan
try:
filled_kurtosis = stats.kurtosis(filled_data[column])
except:
filled_kurtosis = np.nan
# 计算相关系数
try:
correlation = np.corrcoef(benchmark_data[column], filled_data[column])[0, 1]
except ValueError:
correlation = np.nan # 如果列中存在缺失值或数据不匹配，返回NaN
print(f"列名: {column}")
print(f"基准均值: {benchmark_mean:.4f}, 填充均值: {filled_mean:.4f}")
print(f"基准标准差: {benchmark_std:.4f}, 填充标准差: {filled_std:.4f}")
print(f"基准中位数: {benchmark_median:.4f}, 填充中位数: {filled_median:.4f}")
print(f"基准众数: {benchmark_mode}, 填充众数: {filled_mode}")
print(f"基准最大值: {benchmark_max:.4f}, 填充最大值: {filled_max:.4f}")
print(f"基准最小值: {benchmark_min:.4f}, 填充最小值: {filled_min:.4f}")
print(f"基准偏度: {benchmark_skew:.4f}, 填充偏度: {filled_skew:.4f}")
print(f"基准峰度: {benchmark_kurtosis:.4f}, 填充峰度: {filled_kurtosis:.4f}")
print(f"相关系数: {correlation:.4f}")
print("-" * 80) # 分割线，便于区分不同列的结果
# 调用
compare_filled_vs_benchmark(df_outliers_filled, df_final_outliers_filled)
# 折线对比图
def plot_filled_vs_benchmark(benchmark_data, filled_data, column, output_path):
"""
绘制填充数据与基准数据的对比图
:param benchmark_data: 基准数据（原始数据异常值插值）
:param filled_data: 填充后的数据（空白日期数据填充后异常值插值）
:param column: 需要对比的列名
:param output_path: 图片保存路径
"""
plt.figure(figsize=(12, 6))
plt.plot(benchmark_data.index, benchmark_data[column], label='基准数据（原始数据异常值插值）', marker='o', linestyle='-', color='blue')
plt.plot(filled_data.index, filled_data[column], label='填充数据（空白日期数据填充后异常值插值）', marker='x', linestyle='--', color='red')
plt.title(f'{column} 填充数据与基准数据对比')
plt.xlabel('日期')
plt.ylabel(column)
plt.legend()
plt.grid(True)
plt.savefig(output_path)
plt.close()
# 获取所有数值型列作为指标
numeric_columns = df_final.select_dtypes(include=[np.number]).columns
# 遍历所有指标列，绘制对比图
print("绘制折线对比图...")
for column in numeric_columns:
output_path = os.path.join(line_chart_folder, f'{column}_填充与基准对比.png') # 保存到折线对比图子文件夹
plot_filled_vs_benchmark(df_outliers_filled, df_final_outliers_filled, column, output_path)
print(f"所有指标的填充数据与基准数据对比图已保存到子文件夹 '{line_chart_folder}' 中。")
# 时间序列特性对比
def compare_all_autocorrelation(benchmark_data, filled_data, output_folder):
"""
对比填充数据与基准数据的所有指标的自相关性
:param benchmark_data: 基准数据（异常值处理后的数据）
:param filled_data: 填充后的数据
:param output_folder: 图片保存文件夹路径
"""
# 确保输出文件夹存在
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# 检查两个数据框的列名是否一致
if set(benchmark_data.columns) != set(filled_data.columns):
raise ValueError("基准数据和填充数据的列名不一致，请检查数据！")
# 遍历数据框的所有列
print("绘制自相关性对比图...")
for column in benchmark_data.columns:
# 计算自相关系数
benchmark_acf = acf(benchmark_data[column], nlags=20, fft=True)
filled_acf = acf(filled_data[column], nlags=20, fft=True)
# 绘制自相关性对比图
plt.figure(figsize=(12, 6))
plt.plot(benchmark_acf, label='基准数据', marker='o', linestyle='-', color='blue')
plt.plot(filled_acf, label='填充数据', marker='x', linestyle='--', color='red')
plt.title(f'{column} 自相关性对比')
plt.xlabel('滞后阶数')
plt.ylabel('自相关系数')
plt.legend()
plt.grid(True)
# 保存图片
output_path = os.path.join(output_folder, f'{column}_自相关性对比.png')
plt.savefig(output_path)
plt.close()
print(f"所有指标的自相关性对比图已保存到 {output_folder}")
# 调用
compare_all_autocorrelation(df_outliers_filled, df_final_outliers_filled, autocorrelation_plot_folder)
# 绘制分布直方图
def plot_distribution_histograms(benchmark_data, filled_data, column, output_path):
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
sns.histplot(benchmark_data[column].dropna(), bins=30, kde=True, color='skyblue')
plt.title(f'{column} 基准数据分布直方图（原始数据异常值插值）')
plt.xlabel(column)
plt.ylabel('频数')
plt.subplot(2, 1, 2)
sns.histplot(filled_data[column].dropna(), bins=30, kde=True, color='salmon')
plt.title(f'{column} 填充数据分布直方图（空白日期数据填充后异常值插值）')
plt.xlabel(column)
plt.ylabel('频数')
plt.tight_layout()
plt.savefig(output_path)
plt.close()
print("绘制分布直方图对比...")
for column in numeric_columns:
output_path = os.path.join(histogram_folder, f'{column}_分布直方图对比.png')
plot_distribution_histograms(df_outliers_filled, df_final_outliers_filled, column, output_path)
print(f"所有指标的分布直方图对比图已保存到子文件夹 '{histogram_folder}' 中。")
# 绘制QQ图
def plot_qq_comparison(benchmark_data, filled_data, column, output_folder):
"""
绘制QQ图对比（上下子图布局）
:param benchmark_data: 基准数据（原始数据异常值插值）
:param filled_data: 填充后的数据（空白日期数据填充后异常值插值）
:param column: 需要对比的列名
:param output_folder: 图片保存路径
"""
plt.figure(figsize=(10, 8))
# 上子图：基准数据QQ图
plt.subplot(2, 1, 1)
stats.probplot(benchmark_data[column].dropna(), dist="norm", plot=plt)
plt.title(f'基准数据QQ图 - {column}', fontsize=10)
plt.grid(True)
# 下子图：填充数据QQ图
plt.subplot(2, 1, 2)
stats.probplot(filled_data[column].dropna(), dist="norm", plot=plt)
plt.title(f'填充数据QQ图 - {column}', fontsize=10)
plt.grid(True)
# 调整布局并保存
plt.tight_layout()
output_path = os.path.join(output_folder, f'{column}_QQ图对比.png')
plt.savefig(output_path, dpi=150)
plt.close()
# 获取所有数值型列作为指标
numeric_columns = df_final.select_dtypes(include=[np.number]).columns
# 确保数据时间索引对齐
benchmark_data = df_outliers_filled.reindex(df_final_outliers_filled.index)
filled_data = df_final_outliers_filled.copy()
# 遍历所有指标生成QQ图
print("绘制QQ图对比...")
for column in numeric_columns:
try:
# 过滤无效数据
valid_benchmark = benchmark_data[column].dropna()
valid_filled = filled_data[column].dropna()
if len(valid_benchmark) > 1 and len(valid_filled) > 1:
plot_qq_comparison(benchmark_data, filled_data, column, qq_plot_folder)
else:
print(f"跳过 {column}：有效数据点不足（基准数据:{len(valid_benchmark)}，填充数据:{len(valid_filled)}）")
except Exception as e:
print(f"生成 {column} 的QQ图时出错: {str(e)}")
print(f"\n所有QQ图对比已保存到 {qq_plot_folder}")
'''
数据扩充、预处理完成
以下进行PCA分析
'''
# 提取特征数据（仅提取指定的数值列）
# 定义需要提取的特征列名列表
feature_columns = ['TOC', 'COD', '电导率', 'pH', '盐度', '阳离子', '阴离子', '有机物']
# 检查这些列是否存在于df_final_outliers_filled中
missing_columns = [col for col in feature_columns if col not in df_final_outliers_filled.columns]
if missing_columns:
raise ValueError(f"处理后的数据中缺少以下列：{missing_columns}。请检查数据预处理是否正确。")
# 提取指定的特征列
features = df_final_outliers_filled[feature_columns].copy() # 复制一份数据以进行分析
print("\n特征数据预览：")
print(features.head())
# 标准化数据
print("\n标准化数据...")
scaler = StandardScaler()
numeric_columns = features.select_dtypes(include=[np.number])
if numeric_columns.empty:
raise ValueError("特征数据中没有数值列，无法进行标准化。")
features_scaled = scaler.fit_transform(numeric_columns)
df_features_scaled = pd.DataFrame(features_scaled, columns=numeric_columns.columns) # 将标准化后的数据转换为 DataFrame
print("标准化后数据预览：")
print(df_features_scaled.head())
# 相关性分析
print("\n进行相关性分析...")
correlation_matrix = features.corr() # 计算相关性矩阵
print("\n相关性矩阵：")
print(correlation_matrix)
# 绘制相关性矩阵热力图
print("\n绘制相关性矩阵热力图...")
plt.figure(figsize=(8, 6))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f', linewidths=0.5)
plt.title('相关性矩阵热力图', fontsize=14)
# 调整布局
plt.tight_layout()
# 保存并关闭
plt.savefig(os.path.join(pca_analysis_folder, '相关性矩阵热力图.png'), dpi=300, bbox_inches='tight')
plt.close()
# 检验
print("\n进行PCA分析...")
# Bartlett球形检验
chi_square_value, p_value = calculate_bartlett_sphericity(numeric_columns)
print("\nBartlett球形检验结果：")
print(f"卡方值: {chi_square_value}, p值: {p_value}")
if p_value < 0.05:
print("拒绝球形假设，数据适合进行PCA。")
else:
print("不能拒绝球形假设，数据可能不适合进行PCA。")
# KMO检验
kmo_all, kmo_model = calculate_kmo(numeric_columns)
print("\nKMO检验结果：")
print(f"KMO值: {kmo_model}")
if kmo_model > 0.6:
print("KMO值大于0.6，数据适合进行PCA。")
else:
print("KMO值小于0.6，数据可能不适合进行PCA。")
# 执行PCA
pca = PCA(n_components=None) # 提取所有特征值
pca_scores = pca.fit_transform(features_scaled)
# 将PCA结果转换为DataFrame
pca_df = pd.DataFrame(data=pca_scores, columns=[f'PC{i + 1}' for i in range(pca_scores.shape[1])])
pca_df['日期'] = df_final_outliers_filled.index # 将日期添加到PCA结果中
print("PCA得分预览：")
print(pca_df.head())
# 结果解释
print("解释方差比：", pca.explained_variance_ratio_)
# 获取PCA成分矩阵
components = pca.components_
components_df = pd.DataFrame(components, columns=numeric_columns.columns,
index=[f'PC{i + 1}' for i in range(components.shape[0])])
print("\nPCA成分矩阵：")
print(components_df)
# 获取特征值
explained_variance = pca.explained_variance_
print("\n特征值：", explained_variance)
# 分析每个主成分的主要贡献者（保留正负号）
for pc in components_df.index:
print(f"\n{pc}的主要贡献者：")
# 按绝对值降序排序，但保留原始符号
sorted_contributions = components_df.loc[pc].sort_values(key=lambda x: abs(x), ascending=False)
print(sorted_contributions)
# 将 PCA 结果与日期和毒性值合并
result_df = pd.merge(pca_df, df_final_outliers_filled[['LC50']], left_on='日期', right_index=True, how='left')
print("\n合并后的结果预览：")
print(result_df.head())
#PCA结果绘图
# 绘制PCA散点图
# 散点图1：颜色按毒性大小分
print("\n绘制PCA散点图(带变量箭头)...")
plt.figure(figsize=(8, 6)) # 创建画布
scatter = plt.scatter(
result_df['PC1'], # X轴：PC1
result_df['PC2'], # Y轴：PC2
c=result_df['LC50'], # 颜色：毒性值
cmap='viridis', # 颜色映射
s=100, # 点的大小
alpha=0.7,
edgecolor='k', # 点的边框颜色
linewidth=1 # 点的边框宽度
)
# 添加颜色条
cbar = plt.colorbar(scatter)
cbar.set_label('LC50', fontsize=12) # 颜色条标签
# 添加变量箭头（修改比例因子部分）
scale_factor = 3 # 缩放因子，数值越大箭头越长
for i, (x, y) in enumerate(zip(components[0], components[1])):
# 缩放箭头长度
scaled_x = x * scale_factor
scaled_y = y * scale_factor
plt.arrow(0, 0, scaled_x, scaled_y, color='r', alpha=0.5, head_width=0.05, head_length=0.1)
plt.text(scaled_x, scaled_y, features.columns[i], color='b', fontsize=12, ha='center', va='center')
# 添加标题和标签
plt.title('PCA散点图：毒性值分布', fontsize=14)
plt.xlabel('PC1 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[0])) # 添加解释方差比
plt.ylabel('PC2 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[1])) # 添加解释方差比
plt.grid(True)
# 保存图片
plt.savefig(os.path.join(pca_analysis_folder, 'PCA散点图_毒性变化.png'), dpi=300, bbox_inches='tight')
plt.close()
# 散点图2：毒性、含95%置信区间圆
print("\n绘制PCA散点图（颜色表示毒性高低，添加置信椭圆和变量箭头，仅2022-11-02至2022-12-04）...")
# 确保日期列是datetime格式
result_df['日期'] = pd.to_datetime(result_df['日期']) # 转换为日期时间格式，确保正确比较
# 根据日期范围过滤数据
date_filtered_df = result_df[
(result_df['日期'] >= pd.to_datetime('2022-10-21')) &
(result_df['日期'] <= pd.to_datetime('2022-11-01'))
]
# 检查数据是否为空
if date_filtered_df.empty:
print("注意：2022-11-02至2022-12-04期间无数据，无法绘制散点图。")
else:
# 根据毒性值分类
high_toxicity = date_filtered_df[date_filtered_df['LC50'] < 0.2]
low_toxicity = date_filtered_df[date_filtered_df['LC50'] >= 0.2]
# 创建画布
plt.figure(figsize=(8, 6))
# 绘制散点图
plt.scatter(high_toxicity['PC1'], high_toxicity['PC2'], c='darkred', label='毒性较高 (LC50 < 0.2)', s=100, edgecolor='k', linewidth=1)
plt.scatter(low_toxicity['PC1'], low_toxicity['PC2'], c='lightblue', label='毒性较低 (LC50 >= 0.2)', s=100, edgecolor='k', linewidth=1)
# 添加95%置信区间椭圆
def confidence_ellipse(x, y, ax, n_std=2.0, facecolor='none', **kwargs):
cov = np.cov(x, y)
pearson = cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])
ell_radius_x = np.sqrt(1 + pearson)
ell_radius_y = np.sqrt(1 - pearson)
ellipse = Ellipse((0, 0), width=ell_radius_x * 2, height=ell_radius_y * 2, facecolor=facecolor, **kwargs)
scale_x = np.sqrt(cov[0, 0]) * n_std
scale_y = np.sqrt(cov[1, 1]) * n_std
transf = transforms.Affine2D().rotate_deg(45).scale(scale_x, scale_y).translate(np.mean(x), np.mean(y))
ellipse.set_transform(transf + ax.transData)
return ax.add_patch(ellipse)
# 为高毒性和低毒性数据添加椭圆
confidence_ellipse(high_toxicity['PC1'], high_toxicity['PC2'], plt.gca(), edgecolor='darkred', linewidth=2, label='高毒性95%置信区间')
confidence_ellipse(low_toxicity['PC1'], low_toxicity['PC2'], plt.gca(), edgecolor='lightblue', linewidth=2, label='低毒性95%置信区间')
# 添加变量箭头
for i, (x, y) in enumerate(zip(components[0], components[1])):
plt.arrow(0, 0, x, y, color='r', alpha=0.5, head_width=0.05, head_length=0.1)
plt.text(x, y, features.columns[i], color='b', fontsize=12)
# 添加图例
plt.legend()
# 添加标题和标签
plt.title('PCA散点图：毒性值分布 22.11.02-22.12.04（带95%置信椭圆）', fontsize=14)
plt.xlabel('PC1 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[0]))
plt.ylabel('PC2 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[1]))
plt.grid(True)
# 调整布局并保存
plt.tight_layout()
plt.savefig(os.path.join(pca_analysis_folder, 'PCA散点图_毒性变化_带置信圆.png'), dpi=300, bbox_inches='tight')
plt.close()
# 散点图3：颜色代表时间变化
print("\n绘制 PCA 散点图...")
# 将日期转换为时间序列的索引值，用于颜色映射
pca_df['时间索引'] = pd.to_datetime(pca_df['日期']).map(datetime.toordinal)
# 绘制散点图
plt.figure(figsize=(8, 6))
scatter = plt.scatter(
pca_df['PC1'],
pca_df['PC2'],
c=pca_df['时间索引'],
cmap='Accent',
s=100,
alpha=0.7,
edgecolor='k',
linewidth=1
)
plt.title('PCA 散点图：时间分布', fontsize=12)
plt.xlabel('PC1 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[0])) # 添加解释方差比
plt.ylabel('PC2 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[1])) # 添加解释方差比
plt.grid(True)
# 创建颜色条
cbar = plt.colorbar(scatter)
cbar.set_label('日期（时间顺序）')
# 自定义时间点
custom_dates = ['22.10.21', '22.12.04']
# 将自定义时间点转换为时间戳
custom_date_ticks = [pd.to_datetime(date, format='%y.%m.%d').toordinal() for date in custom_dates]
# 设置颜色条的标签
cbar.set_ticks(custom_date_ticks)
cbar.set_ticklabels(custom_dates)
# 添加变量箭头（修改比例因子部分）
scale_factor = 3 # 缩放因子，数值越大箭头越长
for i, (x, y) in enumerate(zip(components[0], components[1])):
# 缩放箭头长度
scaled_x = x * scale_factor
scaled_y = y * scale_factor
plt.arrow(0, 0, scaled_x, scaled_y, color='r', alpha=0.5, head_width=0.05, head_length=0.1)
plt.text(scaled_x, scaled_y, features.columns[i], color='b', fontsize=12, ha='center', va='center')
# 保存图片
plt.savefig(os.path.join(pca_analysis_folder, 'PCA散点图_时间变化.png'), dpi=300, bbox_inches='tight')
plt.close()
# 散点图4：时间、含95%置信区间圆
print("\n绘制PCA散点图（颜色表示时间段，添加置信椭圆和变量箭头）...")
plt.figure(figsize=(8, 6)) # 创建画布
# 根据日期分类
period1 = result_df[(result_df['日期'] >= '2022-10-21') & (result_df['日期'] <= '2022-11-02')]
period2 = result_df[(result_df['日期'] >= '2022-11-03') & (result_df['日期'] <= '2022-12-04')]
# 绘制散点图
plt.scatter(period1['PC1'], period1['PC2'], c='blue', label='2022-10-21-2022-11-02', s=100, edgecolor='k', linewidth=1)
plt.scatter(period2['PC1'], period2['PC2'], c='orange', label='2022-11-03-2022-12-04', s=100, edgecolor='k', linewidth=1)
# 添加95%置信区间椭圆
confidence_ellipse(period1['PC1'], period1['PC2'], plt.gca(), edgecolor='blue', linewidth=2, label='2022-10-21-2022-11-02 95%置信区间')
confidence_ellipse(period2['PC1'], period2['PC2'], plt.gca(), edgecolor='orange', linewidth=2, label='2022-11-03-2022-12-04 95%置信区间')
# 添加变量箭头（修改比例因子部分）
scale_factor = 3 # 缩放因子，数值越大箭头越长
for i, (x, y) in enumerate(zip(components[0], components[1])):
# 缩放箭头长度
scaled_x = x * scale_factor
scaled_y = y * scale_factor
plt.arrow(0, 0, scaled_x, scaled_y, color='r', alpha=0.5, head_width=0.05, head_length=0.1)
plt.text(scaled_x, scaled_y, features.columns[i], color='b', fontsize=12, ha='center', va='center')
# 添加图例
plt.legend(loc='upper left') # 将图例放在左上方
# 添加标题和标签
plt.title('PCA散点图：时间分布（带95%置信椭圆）', fontsize=14)
plt.xlabel('PC1 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[0])) # 添加解释方差比
plt.ylabel('PC2 (解释方差比: {:.2%})'.format(pca.explained_variance_ratio_[1])) # 添加解释方差比
plt.grid(True)
# 调整布局并保存
plt.tight_layout()
plt.savefig(os.path.join(pca_analysis_folder, 'PCA散点图_时间变化_带置信圆.png'), dpi=300, bbox_inches='tight')
plt.close()
# 主成分的方差贡献率图和碎石图合并图
print("\n绘制主成分的方差贡献率图和碎石图...")
plt.figure(figsize=(12, 6))
# 绘制方差贡献率柱状图
plt.bar(range(1, len(pca.explained_variance_ratio_) + 1), pca.explained_variance_ratio_, alpha=0.6, label='方差贡献率')
# 绘制碎石图（累计方差贡献率）
cumulative_variance = np.cumsum(pca.explained_variance_ratio_)
plt.plot(range(1, len(cumulative_variance) + 1), cumulative_variance, marker='o', linestyle='--', color='red', label='累计方差贡献率')
plt.title('主成分的方差贡献率图和碎石图', fontsize=14)
plt.xlabel('主成分编号', fontsize=12)
plt.ylabel('方差贡献率', fontsize=12)
plt.grid(True)
plt.legend(fontsize=10)
plt.savefig(os.path.join(pca_analysis_folder, '主成分的方差贡献率图和碎石图.png'))
plt.close()
# 将结果保存为Excel文件
with pd.ExcelWriter(os.path.join(output_folder, '结果.xlsx')) as writer:
df_missing_filled.to_excel(writer, sheet_name='原始数据缺失值处理后数据')
df_outliers_marked.to_excel(writer, sheet_name='原始数据异常值标记')
df_outliers_filled.to_excel(writer, sheet_name='原始数据异常值插值')
df_merged.to_excel(writer, sheet_name='插入空白日期')
df_final.to_excel(writer, sheet_name='插入空白日期数据填充')
df_final_outliers_marked.to_excel(writer, sheet_name='空白日期数据填充后异常值标注')
df_final_outliers_filled.to_excel(writer, sheet_name='空白日期数据填充后异常值插值')
df_features_scaled.to_excel(writer, sheet_name='标准化后的数据')
correlation_matrix.to_excel(writer, sheet_name='相关性矩阵')
pca_df.to_excel(writer, sheet_name='PCA得分')
components_df.to_excel(writer, sheet_name='PCA成分矩阵')
result_df.to_excel(writer, sheet_name='PCA与毒性值数据')
print(f"结果已保存到文件夹 '{output_folder}' 中。")
机器学习建模预测LOC:
import pandas as pd
import os
from datetime import datetime
from sklearn.model_selection import KFold
from sklearn.preprocessing import MinMaxScaler, StandardScaler, PolynomialFeatures
from sklearn.svm import SVR
from xgboost import XGBRegressor, DMatrix
from lightgbm import LGBMRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.linear_model import LinearRegression
import numpy as np
# 创建保存结果的文件夹
current_time = datetime.now().strftime('%Y-%m-%d_%H-%M-%S')
output_folder = f'分析结果_{current_time}'
os.makedirs(output_folder, exist_ok=True)
# 读取Excel文件
# data = pd.read_excel("填充数据.xlsx")
data = pd.read_excel("toxicity_prediction.xlsx")
data = data[['TOC', '电导率', '盐度', 'LC50']]
# data = data[['TOC', 'COD', '电导率', 'pH', '盐度', 'Na', 'K', 'Ba', 'Sr', 'Ca', 'Li', 'Fe', 'Mg', 'B', 'Al', 'Mn', '阳离子', 'F-', 'Cl-', 'NO2-', 'Br-', 'NO3-', 'SO42-', '阴离子', '26PAHs', '17PESTs', '24ZGWT', '8NPEOs', '8 PAEs', '16 QACs', 'BPA', '有机物', 'LC50']]
print("预览数据:")
print(data.head())
# 提取特征和目标变量
X = data[['TOC', '电导率', '盐度']]
# X = data[['TOC', 'COD', '电导率', 'pH', '盐度', 'Na', 'K', 'Ba', 'Sr', 'Ca', 'Li', 'Fe', 'Mg', 'B', 'Al', 'Mn', '阳离子', 'F-', 'Cl-', 'NO2-', 'Br-', 'NO3-', 'SO42-', '阴离子', '26PAHs', '17PESTs', '24ZGWT', '8NPEOs', '8 PAEs', '16 QACs', 'BPA', '有机物']]
y = data['LC50']
# breakpoint()
# 特征工程：多项式扩展
poly = PolynomialFeatures(degree=3, include_bias=False)
X_poly = poly.fit_transform(X)
poly_feature_names = poly.get_feature_names_out(['TOC', '电导率', '盐度'])
# poly_feature_names = poly.get_feature_names_out(['TOC', 'COD', '电导率', 'pH', '盐度', 'Na', 'K', 'Ba', 'Sr', 'Ca', 'Li', 'Fe', 'Mg', 'B', 'Al', 'Mn', '阳离子', 'F-', 'Cl-', 'NO2-', 'Br-', 'NO3-', 'SO42-', '阴离子', '26PAHs', '17PESTs', '24ZGWT', '8NPEOs','8 PAEs', '16 QACs', 'BPA', '有机物'])
print(f"poly_feature_names is {poly_feature_names}")
# features = ['TOC', 'COD', '电导率', 'pH', '盐度', 'Na', 'K', 'Ba', 'Sr', 'Ca', 'Li', 'Fe', 'Mg', 'B', 'Al', 'Mn', '阳离子', 'F-', 'Cl-', 'NO2-', 'Br-', 'NO3-', 'SO42-', '阴离子', '26PAHs', '17PESTs', '24ZGWT', '8NPEOs','8 PAEs', '16 QACs', 'BPA', '有机物']
X_poly_df = pd.DataFrame(X_poly, columns=poly_feature_names)
# X_poly_df = pd.DataFrame(X, columns=['TOC', '电导率', '盐度'])
# X_poly_df = pd.DataFrame(X, columns=features)
# 数据标准化
# scaler = StandardScaler()
scaler = MinMaxScaler()
X_scaled = scaler.fit_transform(X_poly_df)
X_scaled_df = pd.DataFrame(X_scaled, columns=X_poly_df.columns)
print("\n多项式扩展后标准化数据预览:")
print(X_scaled_df.head())
# 交叉验证设置
kf = KFold(n_splits=5, shuffle=True, random_state=42)
# kf = KFold(n_splits=5, shuffle=False)
# 模型定义
models = {
"Support Vector Regression (SVR)": SVR(
kernel='rbf',
C=10.0,
gamma='auto',
epsilon=0.1
),
"XGBoost": XGBRegressor(
n_estimators=100,
max_depth=5,
learning_rate=0.1,
subsample=0.8,
colsample_bytree=0.6,
random_state=42,
),
"Random Forest": RandomForestRegressor(
n_estimators=100,
max_depth=5,
min_samples_split=3,
min_samples_leaf=2,
max_features='sqrt',
random_state=42
),
"Decision Tree": DecisionTreeRegressor(
max_depth=5,
min_samples_split=3,
min_samples_leaf=2,
random_state=42
),
"Kernel Ridge Regression": KernelRidge(
kernel='polynomial',
alpha=0.1,
gamma=0.1
),
"LightGBM": LGBMRegressor(
n_estimators=100,
max_depth=5,
learning_rate=0.1,
num_leaves=31,
min_child_samples=5,
reg_alpha=0.1,
reg_lambda=0.1,
random_state=42,
),
"LinearModel": LinearRegression(),
}
# models = {
# "Support Vector Regression (SVR)": SVR(
# kernel='rbf',
# C=2.0,
# gamma='auto',
# epsilon=0.1
# ),
# "XGBoost": XGBRegressor(
# n_estimators=1,
# max_depth=2,
# learning_rate=0.1,
# subsample=0.8,
# colsample_bytree=0.6,
# random_state=42,
# ),
# "Random Forest": RandomForestRegressor(
# n_estimators=1,
# max_depth=2,
# min_samples_split=3,
# min_samples_leaf=2,
# max_features='sqrt',
# random_state=42
# ),
# "Decision Tree": DecisionTreeRegressor(
# max_depth=1,
# min_samples_split=3,
# min_samples_leaf=2,
# random_state=42
# ),
# "Kernel Ridge Regression": KernelRidge(
# kernel='polynomial',
# alpha=0.1,
# gamma=0.1
# ),
# "LightGBM": LGBMRegressor(
# n_estimators=1,
# max_depth=3,
# learning_rate=0.1,
# num_leaves=11,
# min_child_samples=5,
# reg_alpha=0.1,
# reg_lambda=0.1,
# random_state=42,
# ),
# "LinearModel": LinearRegression(),
# }
# 模型训练与评估
results = {}
for name, model in models.items():
print(f"\n正在训练模型: {name}")
mse_scores = []
mae_scores = []
r2_scores = []
train_r2_scores = []
for train_idx, val_idx in kf.split(X_scaled_df):
# 划分数据集
X_train, X_val = X_scaled_df.iloc[train_idx], X_scaled_df.iloc[val_idx]
y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]
# 模型训练
if name == "XGBoost":
model.fit(X_train, y_train, verbose=True)
elif name == "LightGBM":
model.fit(X_train, y_train)
else:
model.fit(X_train, y_train)
# 预测与评估
y_pred_val = model.predict(X_val)
y_pred_train = model.predict(X_train)
# 计算指标
mse = mean_squared_error(y_val, y_pred_val)
mae = mean_absolute_error(y_val, y_pred_val)
val_r2 = r2_score(y_val, y_pred_val)
train_r2 = r2_score(y_train, y_pred_train)
mse_scores.append(mse)
mae_scores.append(mae)
r2_scores.append(val_r2)
train_r2_scores.append(train_r2)
# 计算平均指标
results[name] = {
"MSE": np.mean(mse_scores),
"MAE": np.mean(mae_scores),
"Val R²": np.mean(r2_scores),
"Train R²": np.mean(train_r2_scores),
"Overfit Risk": np.mean(train_r2_scores) - np.mean(r2_scores)
}
# 结果输出
results_df = pd.DataFrame(results).T
results_df_sorted = results_df.sort_values(by='Val R²', ascending=False)
print("\n模型性能对比:")
print(results_df_sorted)
# 保存结果到Excel
output_file = os.path.join(output_folder, '分析结果.xlsx')
with pd.ExcelWriter(output_file) as writer:
data.to_excel(writer, sheet_name='原始特征', index=False)
X_poly_df.to_excel(writer, sheet_name='多项式扩展特征', index=False)
X_scaled_df.to_excel(writer, sheet_name='标准化特征', index=False)
results_df_sorted.to_excel(writer, sheet_name='模型性能对比', index=True)
print(f"\n分析结果已保存到 {output_folder}")