【Kaggle】练习赛《洪水数据集的回归预测》(上)

前言

关于 kaggle 月赛也不多说明，前面两篇《肥胖风险的多类别预测》和 《鲍鱼年龄预测》 已做详细说明。分别是一个分类模型和一个回归模型。本期是2024年5月份的题目《Regression with a Flood Prediction Dataset
》即《洪水数据集的回归预测》，本以为回归模型，与前一篇差不多，没有什么新意，也想写有没有写这篇文章的必要。可随着参与这竞赛后，发现与我之前参与的题目不一样，颠覆了我的认知，使我重认识机器学习的魅力。
因此，也用二篇博客来完成这次题目的分享，按照惯例，上篇讲解数据探索（ EDA） 方面，下篇讲解建模优化方面。

题目简介

洪水探测是指识别、监测和提醒当局或个人特定地区存在或可能发生洪水的过程。它涉及使用各种技术和方法来检测、预测和减轻洪水的影响。
先上一张关于洪水的图片。

数据集介绍

本次比赛的数据集（训练和测试）由在洪水预测因子数据集上训练的深度学习模型生成。特征分布与原始分布接近，但并不完全相同。作为这场比赛的一部分，可以随意使用原始数据集，既可以探索差异，也可以看看在训练中加入原始数据集是否能提高模型性能。
注意：此数据集特别适合于可视化、集群和通用EDA。展示你的技能！

以上是官方的说明，着重说明了这个数据集是深度学习生成的，还有原始数据集。原始数据集 链接

以下是数据集中各列的描述（包括功能名称的含义）：

MonsoonIntensity（季风强度）：这一特征可能衡量该地区季风降雨的强度和频率，较高的值表示降雨强度更大，可能更频繁，这可能会导致更高的洪水风险。
TopographyDrainage（地形排水）： 反映地形的自然排水能力。更好的排水能力（可能由更高的分数表示）可能表明，由于人口稠密或关键地区的水流更好，洪水风险更低。
RiverManagement（河流管理）： 评估河流在防洪方面的管理情况，包括河岸、大坝和其他基础设施的维护。更高的分数可能意味着更好的管理，有可能降低洪水风险。
Deforestation（森林砍伐）：衡量影响土壤稳定性和吸水性的森林砍伐率或程度。更高的森林砍伐分数可能表明森林覆盖的损失更大，更容易受到洪水的影响。
Urbanization（城市化）：表示城市发展水平，由于不透水表面的增加，通常会降低土地吸收降雨的自然能力。城市化程度的提高可能与洪水风险的增加有关。
ClimateChange（气候变化）：　评估气候变化的影响，如降雨量增加或海平面上升，这可能会加剧洪水。得分越高可能表示更容易受到这些变化的影响。
DamsQuality（水坝质量）：　研究大坝在防洪中的状况和有效性。大坝质量差或维护不善可能导致更高的洪水风险。
Siltation（淤积）：　测量水体中淤泥堆积的程度，这会降低其管理水流的能力，增加洪水风险。
AgriculturalPractices （农业实践）：　评估农业活动对洪水风险的影响，考虑灌溉方法和土地使用等可能影响径流和土壤侵蚀的因素。
Encroachments（侵扰）：　对人类入侵洪水易发地区的程度进行评分，这可能会加剧洪水的严重性。
IneffectiveDisasterPreparedness（无效的灾难准备）：　反映了备灾计划及其实施的不足。得分越高可能表明准备工作越差，潜在的破坏和洪水恢复时间越长。
DrainageSystems（排水系统）：　评估城市和农村地区排水系统处理强降雨和水流的效率和容量。
CoastalVulnerability（海岸脆弱性）：　评估沿海地区因风暴潮、海平面上升和气旋活动等因素而发生洪水的风险。
Landslides（滑坡）：　表示山体滑坡的风险和历史，当这些天然大坝决堤时，山体滑坡会堵塞河流，并在下游造成突发洪水。
Watersheds（流域）：　评估流域的健康和管理，流域在管理水资源和减轻洪水风险方面发挥着关键作用。
DeterioratingInfrastructure（不断恶化的基础设施）：　查看与洪水管理相关的基础设施的总体状况，如下水道、桥梁和道路。恶化会阻碍有效的洪水应对。
PopulationScore（人口得分）：　测量洪水易发地区的人口密度或增长，这可能会影响洪水对人类社区的影响。
WetlandLoss（湿地流失）：　量化湿地的减少，湿地通过吸收洪水起到天然缓冲作用。
InadequatePlanning（计划不足）：　评估区域和城市规划在多大程度上整合了洪水风险管理，包括分区和土地利用政策。
PoliticalFactors（政治因素）：　考虑政治决策、治理和政策实施如何影响洪水管理实践。
FloodProbability（洪水概率）：　该结果变量基于上述因素预测洪水的可能性，可能表示为0到1之间的概率。

加载库

# 加载库
import pandas as pd
import numpy as np
from numpy import mean
from numpy import std
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
from datetime import datetime
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import StratifiedKFold
from xgboost import XGBClassifier
from sklearn.metrics import f1_score, classification_report, confusion_matrix, roc_auc_score, roc_curve, accuracy_score
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

设置样式

sns.set()
sns.set_palette('mako')
SNS_CMAP = 'mako'

colors = sns.palettes.color_palette(SNS_CMAP)

加载数据

# 加载所有数据
original_data = pd.read_csv('/kaggle/input/flood-prediction-factors/flood.csv')
train_data = pd.read_csv('/kaggle/input/playground-series-s4e5/train.csv')
test_data = pd.read_csv('/kaggle/input/playground-series-s4e5/test.csv')

这里说明一下，这个original_data 是原始数据集。参考链接为：https://www.kaggle.com/datasets/brijlaldhankour/flood-prediction-factors

数据探索（EDA）

original_missing = original_data.isnull().sum()
train_missing = train_data.isnull().sum()
test_missing = test_data.isnull().sum()

original_missing, train_missing, test_missing

以上给出没有缺失值。

# The summary statistics for original_data
original_data.describe().style.background_gradient(cmap=SNS_CMAP)

给出原始数据统计情况

# The summary statistics for train_data
train_data_drop_id = train_data.drop(columns='id')
train_data_drop_id.describe().style.background_gradient(cmap=SNS_CMAP)

给出训练数据统计情况

# The summary statistics for test_data
test_data_drop_id = test_data.drop(columns='id')
test_data_drop_id.describe().style.background_gradient(cmap=SNS_CMAP)

给出测试数据统计情况

数据集大小情况

original_data.shape,train_data_drop_id.shape,test_data_drop_id.shape

((50000, 21), (1117957, 21), (745305, 20))

数据分布

import matplotlib.pyplot as plt
import seaborn as sns

# Set the style and size for the plots
sns.set(style='whitegrid')
sns.set_palette('mako')
plt.figure(figsize=(16, 24))

# Create a list of datasets and titles
datasets = [original_data, train_data_drop_id, test_data_drop_id]
titles = ['Flood Original Data', 'Train Data', 'Test Data']

# Loop through the datasets and create a boxplot for each
for i, dataset in enumerate(datasets):
    plt.subplot(3, 1, i + 1)
    sns.boxplot(data=dataset.drop(['id'], axis=1, errors='ignore'))
    plt.title(titles[i])
    plt.xticks(rotation=90)

plt.tight_layout()
plt.show()

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from math import ceil, log2

cont_cols = [f for f in train_data_drop_id.columns if train_data_drop_id[f].dtype in [float, int] and train_data_drop_id[f].nunique() > 2 and f not in ["FloodProbability"]]

# Calculate the number of rows needed for the subplots
num_rows = (len(cont_cols) + 2) // 3

# Create subplots for each continuous column
fig, axs = plt.subplots(num_rows, 3, figsize=(18, num_rows * 5), constrained_layout=True)

# Loop through each continuous column and plot the histograms
for i, col in enumerate(cont_cols):
    # Determine the range of values to plot
    max_val = max(train_data_drop_id[col].max(), test_data_drop_id[col].max(), original_data[col].max())
    min_val = min(train_data_drop_id[col].min(), test_data_drop_id[col].min(), original_data[col].min())
    range_val = max_val - min_val
    
    # Determine the bin size and number of bins
    bin_size = range_val / 20
    num_bins_train = round(range_val / bin_size)
    num_bins_test = round(range_val / bin_size)
    num_bins_original = round(range_val / bin_size)
    
    # Calculate the subplot position
    row = i // 3
    col_pos = i % 3
    
    # Plot the histograms
    sns.histplot(train_data_drop_id[col], ax=axs[row][col_pos], color='darkturquoise', kde=True, label='Train', bins=num_bins_train)
    sns.histplot(test_data_drop_id[col], ax=axs[row][col_pos], color='salmon', kde=True, label='Test', bins=num_bins_test)
    sns.histplot(original_data[col], ax=axs[row][col_pos], color='orange', kde=True, label='Original', bins=num_bins_original)
    axs[row][col_pos].set_title(col)
    axs[row][col_pos].set_xlabel('Value')
    axs[row][col_pos].set_ylabel('Frequency')
    axs[row][col_pos].legend()

# Remove any empty subplots
if len(cont_cols) % 3 != 0:
    for col_pos in range(len(cont_cols) % 3, 3):
        axs[-1][col_pos].remove()

plt.tight_layout()
plt.show()

# Calculate the correlation matrix for the original data
correlation_matrix = original_data.corr()

# Plot the correlation heatmap
plt.figure(figsize=(14, 12))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f", cmap='flare', cbar=True)
plt.title('Correlation Matrix for Original Data')
plt.show()

# Calculate the correlation matrix for train_data
correlation_matrix = train_data_drop_id.corr()

# Plot the correlation heatmap
plt.figure(figsize=(14, 12))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f", cmap='crest', cbar=True)
plt.title('Correlation Matrix for Train Data')
plt.show()

import seaborn as sns
import matplotlib.pyplot as plt

# Assuming 'train_data_drop_id' is your dataset and it has been preprocessed appropriately
features = train_data_drop_id.columns.drop('FloodProbability')  # Exclude target variable

for feature in features:
    sns.set_theme(style="white", palette=None)
    g = sns.jointplot(x=feature, y="FloodProbability", data=train_data_drop_id, kind="scatter", color="darkturquoise")
    g.plot_joint(plt.scatter, c="black", s=30, linewidth=1, marker="+")
    g.ax_joint.collections[0].set_alpha(0)  # Make the kde plot transparent
    g.set_axis_labels(f"${feature}$", "$FloodProbability$")
    g.fig.suptitle(f'Distribution of Flood Probability vs {feature}')
    g.fig.subplots_adjust(top=0.95)  # Adjust the title to not overlap with plots
    
    plt.show()

数据特征小节

小伙伴们，如果你看到这里，那就恭喜你，你快要Get到写这篇文章的初衷。有没有发现什么？之所以我把这些（这么多）图都展示出来，主要我看不出各特征的区别，每个特征非常相似，特征与特征之间相关性都是接近零，这是我之前没有碰到的过的事，到现在你可能没有意识到问题的严重性，所以呢，我先卖个关子，请大伙耐心地往下看。

目标值特征

import matplotlib.pyplot as plt
import seaborn as sns

def plot_distribution_with_stats(data, column, title, color):
    plt.figure(figsize=(10, 5))
    sns.histplot(data[column], kde=True, color=color, bins=30)
    mean_value = data[column].mean()
    median_value = data[column].median()
    plt.axvline(mean_value, color='r', linestyle='--', label=f'Mean: {mean_value:.2f}')
    plt.axvline(median_value, color='g', linestyle='-', label=f'Median: {median_value:.2f}')
    plt.title(title)
    plt.xlabel('Flood Probability')
    plt.ylabel('Frequency')
    plt.legend()
    plt.show()

# For flood_data
plot_distribution_with_stats(original_data, 'FloodProbability', 'Original Data Flood Probability Distribution', 'blue')

# For train_data
plot_distribution_with_stats(train_data, 'FloodProbability', 'Train Data Flood Probability Distribution', 'red')

Shapiro-Wilk 检验

from scipy.stats import shapiro

# Perform Shapiro-Wilk test
flood_shapiro_test = shapiro(original_data['FloodProbability'])
train_shapiro_test = shapiro(train_data['FloodProbability'])

print("Flood Data Shapiro-Wilk Test:")
print(f"Statistic: {flood_shapiro_test.statistic}, p-value: {flood_shapiro_test.pvalue}")

print("\nTrain Data Shapiro-Wilk Test:")
print(f"Statistic: {train_shapiro_test.statistic}, p-value: {train_shapiro_test.pvalue}")

Flood Data Shapiro-Wilk Test:
Statistic: 0.9986612796783447, p-value: 2.831748137029014e-19

Train Data Shapiro-Wilk Test:
Statistic: 1.0006940364837646, p-value: 1.0

Shapiro-Wilk检验是一种统计方法，用于检验数据集是否符合正态分布。在这种检验中，统计量（S值）接近1表示数据与正态分布非常接近，而接近0则表示数据与正态分布差异很大。p-value值表示观察到这种统计量或更极端统计量在数据符合正态分布的假设下的概率。
根据上述结果，两组数据的统计量都接近1，这表明这两组数据都与正态分布非常接近。但是，original_data 数据的p-value接近0，这表明在正态分布的假设下，观察到这种统计量或更极端统计量的概率非常低，因此我们拒绝 original_data 数据符合正态分布的零假设。相反，train_data 数据的p-value接近1，这表明在正态分布的假设下，观察到这种统计量或更极端统计量的概率非常高，因此我们不能拒绝第二组数据符合正态分布的零假设。
综上所述，尽管两组数据都与正态分布非常接近，但 original_data 数据与正态分布的差异在统计上是显著的，而 train_data 数据与正态分布的差异在统计上不显著。

不同数据分布

分别用无界约翰逊分布，正态分布，对数分布进行可视化

import scipy.stats as st
y = original_data['FloodProbability']
plt.figure(1); plt.title('Johnson SU')
sns.distplot(y, kde=False, fit=st.johnsonsu)
plt.figure(2); plt.title('Normal')
sns.distplot(y, kde=False, fit=st.norm)
plt.figure(3); plt.title('Log Normal')
sns.distplot(y, kde=False, fit=st.lognorm)

y = train_data['FloodProbability']
plt.figure(1); plt.title('Johnson SU')
sns.distplot(y, kde=False, fit=st.johnsonsu)
plt.figure(2); plt.title('Normal')
sns.distplot(y, kde=False, fit=st.norm)
plt.figure(3); plt.title('Log Normal')
sns.distplot(y, kde=False, fit=st.lognorm)

这里更进一步验证目标值的正态分布情况。

PCA作图

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# Prepare data for PCA
features = train_data_drop_id.columns[1:-1]  # Exclude 'id' and 'FloodProbability'
X = train_data_drop_id[features]
y = train_data_drop_id['FloodProbability'] > 0.5  # Create binary labels based on the threshold

# Standardize the data
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Apply PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# Plotting the PCA result
plt.figure(figsize=(10, 6))
plt.scatter(X_pca[y, 0], X_pca[y, 1], label='FloodProbability >= 0.5', alpha=0.5, color='red')
plt.scatter(X_pca[~y, 0], X_pca[~y, 1], label='FloodProbability < 0.5', alpha=0.5, color='blue')
plt.title('PCA of Dataset (2 Components)')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.legend()
plt.grid(True)
plt.show()

从上图，可以看出，数据非常集中，分布非常理想。

写到这里，我们去数据集有了一个总体的了解，这个数据集应该是非常理想。接下来就建立模型了。

建立模型

建模的过程全留给下篇吧，这里先写个Baseline的结果。
数据集以7：3的方式划分训练集和验证集，在不做任何优化条件下。【小伙伴想重现结果的话，可以用以下参数train_test_split(X,y,test_size=0.3,random_state=100)】

模型	验证集$R^2$
LinearRegression	0.845460
Lasso	-5.51093
Ridge	0.845460
ElasticNet	-5.51093
SVR	0.69667
RandomForestRegressor	0.65394
XGBRegressor	0.80942
LGBMRegressor	0.767198
CatBoostRegressor	0.8466935

以上表上名称直接用代码的原名，方便小伙伴直接使用，参考代码如下，只要替换LinearRegression成上表的名称即可。
特别提醒，在执行 RandomForestRegressor 需要花费较长时间待。

from sklearn.metrics import mean_absolute_error
from sklearn.metrics import make_scorer, r2_score

line = LinearRegression()
line.fit(X_train,y_train)
preds = line.predict(X_valid)
mse=mean_absolute_error(y_valid, preds)
r2=r2_score(y_valid, preds)
mse,r2

	feature	importance
8	AgriculturalPractices	0.050861
11	DrainageSystems	0.050823
4	Urbanization	0.050602
19	PoliticalFactors	0.050539
12	CoastalVulnerability	0.050538
0	MonsoonIntensity	0.050519
3	Deforestation	0.050474
13	Landslides	0.050451
9	Encroachments	0.050377
17	WetlandLoss	0.050323
5	ClimateChange	0.050207
18	InadequatePlanning	0.050175
7	Siltation	0.050089
14	Watersheds	0.050029
16	PopulationScore	0.049935
10	IneffectiveDisasterPreparedness	0.049464
15	DeterioratingInfrastructure	0.049378
6	DamsQuality	0.048836
1	TopographyDrainage	0.048542
2	RiverManagement	0.047839

总结

1. 从上表可以看出线
2. 性结果比较理想
3. 树性模型 除了 CatBoostRegressor 之外均不太理想。
4. 随机森林重要性排序结果可以看出各个特征的重要接近，无法进行筛选。
   总之，如果直接用上述的数据建模的话，$R^2$ 的上限不太可能会突破0.85，必须寻求新的思路。
   这个问题，就留给下一篇文章《洪水数据集的回归预测》(下)。