python实战（二）——房屋价格回归建模

一、任务背景

本章将使用一个经典的Kaggle数据集——House Prices - Advanced Regression Techniques进行回归建模的讲解。这是一个房价数据集，与我们熟知的波士顿房价数据集类似，但是特征数量要更多，数据也要更为复杂一些。下面，我们将使用这个房价数据进行机器学习中的回归建模，依次分解机器学习建模的各个步骤并给出详细的注解。

二、机器学习建模

1、数据获取

为了展示一个完整的机器学习建模流程，我们把数据集下载到了本地并使用python进行读取：

df = pd.read_csv('./data/housePrices.csv')
print('数据量：', len(df))
print(df.head())

结果如下，数据总量是1460条，可以看到最后一列SalePrice就是我们要预测的标签了，此外除Id列之外还有80个特征列。

2、探索性数据分析

通过常规的探索性分析，大概看一下这个数据集长什么样子，比如有多少行、有多少列、每一列中是否有空值、哪些列是数值类型的、哪些列是文本类型的，以及最重要的——我们要预测的标签是哪一列。

（1）查看空值

由于特征列较多，所以全都打出来比较难阅读，这里通过筛选的方式，仅打印有空值的特征列：

for col in list(df.columns):sum_na = df[col].isna().sum()if sum_na > 0:print("{:<20} {}".format(col, sum_na))

结果如下，有空值的列数为19列，其中一些列含有较多的空值，比如最后三列。

（2）查看特征取值类型

我们可以通过筛选数据字段类型的方式分别统计出有多少是数值型特征，有多少是文本特征：

# 选择数值型的列
numeric_df = df.select_dtypes(include=['int64', 'float64'])
# 选择文本型的列
text_df = df.select_dtypes(include=['object'])  # 注意：这里可能包括了混合类型的列
# 打印结果
print("数值型列：")
print('数量：{}， 列名：{}'.format(len(numeric_df.columns), numeric_df.columns.tolist()))
print("文本型列：")
print('数量：{}， 列名：{}'.format(len(text_df.columns), text_df.columns.tolist()))

结果如下：

或者，还有个一步到位的方式：

print(df.info())

结果如下，可以看到，使用.info()方法既可以打印出数据量，也能显示每一列非空值的数量以及字段类型等信息。

（3）查看标签分布

通过可视化的方式把标签分布可视化出来

plt.figure(figsize=(9, 8))
sns.histplot(df['SalePrice'], color='b', bins=100)
plt.show()

可以看到标签取值范围集中在100K到250K的区间内。

3、数据预处理

在进行了初步的数据探索之后，我们可以开始数据的预处理了。首先，部分特征列存在大量的空值，虽然可以取当前列有效数据的平均值或者中值进行填充，但是由于空值列比例太大，笔者在这里直接去掉这些列。对于数据为文本类型的列，我们进行一个简单的转换，确保每一列的数据都是数值类型，以方面后续的相关性计算以及建模。由于文本列也可能存在空值，这里统一编码为“missing”。

# 去掉空值数量超过600的列，不同数据集相应空值数量临界值的界定视具体情况调整
df = df.drop(['Alley', 'MasVnrType', 'FireplaceQu', 'PoolQC', 'Fence', 'MiscFeature'], axis=1)
# 对于文本数据列首先填充空值，再进行数值化转换
label_encoder = LabelEncoder()
for col in text_df.columns:if col not in df:continuedf[col].fillna('missing', inplace=True)# 转换成数值类型的标号df[col] = label_encoder.fit_transform(df[col])

其次，对于包含少量空值的数值类型列，我们使用中值进行填充。

# 计算每一列的中值
median_values = df.median()
# 使用中值填充空值
df = df.fillna(value=median_values)

4、特征工程

由于总共有将近80列的特征，实际上我们可能不需要这么多特征，所以我们进行一个简单的相关性分析（这里使用Pearson相关系数，衡量线性相关度，也可以使用斯皮尔曼系数等来衡量非线性相关度），只保留跟标签列较相关的特征（相关性大于0.5）。在本文中，笔者只进行特征工程中的“特征选择方法”的展示，至于特征构造等内容，视建模结果而定，若建模结果不佳，则可能需要人工构造新特征了。

correlation_matrix = df.corr(method='pearson')
# 计算所有特征列与标签列之间的相关系数
correlations = df.drop('SalePrice', axis=1).corrwith(df['SalePrice'])
# 筛选出相关度大于0.5的特征列，实际取多少根据数据规律确定，例如这个数据集超过0.8的有0列，自然无法指定高于0.8的取值
highly_correlated_features = correlations[abs(correlations) > 0.5].index.tolist()
print("特征列与标签列相关度大于0.5的特征：", highly_correlated_features)

5、训练集/测试集划分

在进行建模之前划分好训练集的测试集，训练集用于模型训练，测试集测试模型的效果。

# 选定训练特征和标签
X = df[highly_correlated_features].values.tolist()
y = df['SalePrice'].tolist()
# 数据划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=2024)
print('训练集数据量：', len(X_train))
print('测试集数据量：', len(X_test))

这里讲一下几个比较重要的参数：

test_size：指定训练集和测试集划分过程中的比例，0-1开区间之间的小数。
random_state：指定随机数，以保证结果可复现。

6、模型训练

这里我们使用一个GBDT回归器作为当前任务的模型。我们使用R2作为评价模型的指标，R2显示了一个回归模型的可解释方差占总方差的比例，R2得分越接近1代表模型效果越好。此外，MSE和MAE是两个常见的指标，不赘述了。

gbr = GradientBoostingRegressor(loss="squared_error", n_estimators=100, criterion="friedman_mse", random_state=2024)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
# 评估模型表现
print('r2:', r2_score(y_test, y_pred))
print('mae:', mean_absolute_error(y_test, y_pred))
print('mse:', mean_squared_error(y_test, y_pred))
print('rmse:', mean_squared_error(y_test, y_pred)**0.5)

这里同样有几个需要注意的参数：

n_estimators：指定要构建的决策树的数量，默认值是100。更多的树可以提高模型的复杂度和拟合能力，但也会增加过拟合的风险及计算量。
max_depth：指定每个决策树的最大深度，树的深度太浅会欠拟合，树的深度太深则会过拟合，需要凭经验调整，默认值是3。
loss：指定优化目标函数，可选squared_error（默认）、absolute_error、huber、quantile。squared _ error是回归的平方误差；absolute_error是回归的绝对误差；huber是两者的结合；quantile允许使用分位数回归(使用alpha指定分位数)。
criterion：指定衡量分裂质量的准则，可以选friedman_mse或mse。默认是friedman_mse，在某些情况下它可以得到更好的近似值。
random_state：随机数种子，用于控制树的抽样和特征选择的随机性，默认是None，即不指定。

指标结果如下，由于我们并没有对标签进行归一化，所以MAE和MSE数值会非常大，但这是正常的，毕竟我们的标签房价基本都是10万元以上的：

可视化模型预测结果和真实值的差异，可以见到预测效果还不错（这里我们的数据并不是时序性的，但是为了观察预测值和真实值的拟合情况，可以把数据看作是时序的并使用曲线进行可视化，这时候的X轴是标签对应的序号，Y轴是房价）：

plt.plot(y_test, label='real')
plt.plot(y_pred, label='pred')
plt.legend()
plt.show()

三、完整代码

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_errordf = pd.read_csv('./data/housePrices.csv')
print('数据量：', len(df))
print(df.head())# 查看包含空值的列
for col in list(df.columns):sum_na = df[col].isna().sum()if sum_na > 0:print("{:<20} {}".format(col, sum_na))# 选择数值型的列
numeric_df = df.select_dtypes(include=['int64', 'float64'])
# 选择文本型的列
text_df = df.select_dtypes(include=['object'])  # 注意：这里可能包括了混合类型的列
# 打印结果
print("数值型列：")
print('数量：{}， 列名：{}'.format(len(numeric_df.columns), numeric_df.columns.tolist()))
print("文本型列：")
print('数量：{}， 列名：{}'.format(len(text_df.columns), text_df.columns.tolist()))
print(df.info())# 可视化标签分布
print(df['SalePrice'].describe())
plt.figure(figsize=(9, 8))
sns.histplot(df['SalePrice'], color='b', bins=100)
plt.show()# 去掉空值数量超过600的列，不同数据集相应空值数量临界值的界定视具体情况调整
df = df.drop(['Alley', 'MasVnrType', 'FireplaceQu', 'PoolQC', 'Fence', 'MiscFeature'], axis=1)
# 对于文本数据列首先填充空值，再进行数值化转换
label_encoder = LabelEncoder()
for col in text_df.columns:if col not in df:continuedf[col].fillna('missing', inplace=True)# 转换成数值类型的标号df[col] = label_encoder.fit_transform(df[col])# 计算每一列的中值
median_values = df.median()
# 使用中值填充空值
df = df.fillna(value=median_values)correlation_matrix = df.corr(method='pearson')
# 计算所有特征列与标签列之间的相关系数
correlations = df.drop('SalePrice', axis=1).corrwith(df['SalePrice'])
# 筛选出相关度大于0.5的特征列，实际取多少根据数据规律确定，例如这个数据集超过0.8的有0列，自然无法指定高于0.8的取值
highly_correlated_features = correlations[abs(correlations) > 0.5].index.tolist()
print("特征列与标签列相关度大于0.5的特征：", highly_correlated_features)# 选定训练特征和标签
X = df[highly_correlated_features].values.tolist()
y = df['SalePrice'].tolist()
# 数据划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=2024)
print('训练集数据量：', len(X_train))
print('测试集数据量：', len(X_test))gbr = GradientBoostingRegressor(loss="squared_error", n_estimators=100, criterion="friedman_mse", random_state=2024)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
# 评估模型表现
print('r2:', r2_score(y_test, y_pred))
print('mae:', mean_absolute_error(y_test, y_pred))
print('mse:', mean_squared_error(y_test, y_pred))# 预测效果可视化
plt.plot(y_test, label='real')
plt.plot(y_pred, label='pred')
plt.legend()
plt.show()