kaggle入门三部曲的尾声，终于迎来了一个难搞的数据集。之前的Titanic用到的特征工程和数据清洗技巧放在这里已经完全不够用了，于是这个比赛的主要目的在于学习对数据的处理。

首先，我们读入数据并观察一下数据的形式：

train = pd.read_csv(r'C:\Users\Sidney Nash\Desktop\data\House\train.csv')
test = pd.read_csv(r'C:\Users\Sidney Nash\Desktop\data\House\test.csv')

train.head()

5 行 81 列的数据表中除了我们关注的 target 值 SalePrice 外还有 80 个 feature ，看起来是好事，毕竟我们的资料非常齐全，可供提取的信息更多，但麻烦随之而来，这样大规模的数据必然是杂乱和冗余的，我们需要对数据进行清洗和分析。

首先，Id 这一列不属于特征的范畴，我们可以将其删除。

train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)

接下来我们对数据进行预处理。

1、异常值

fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

可以看到，右下角的两个数据明显是和整个趋势不符的，而且偏离很大，因此我们将这两个数据点移除。

train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index, 
           inplace=True)

这里需要注意的是表示与运算要用 & 而不是 and，用 and 会导致报错。

此时我们再次画图，验证一下异常点是否已经删除。

在训练数据中可能还有其他的异常值。然而，如果在测试数据中也有异常值，那么删除所有这些异常值可能会严重影响我们的模型。这就是为什么我们不把它们全部删除，而只是设法使我们的模型对异常值更加稳健。

2、错误值

数据集中的错误往往难以发现，即使是大的离谱或小的离谱的值，也可能只是特殊情况下的异常值罢了，除非有非常明显的错误（一般是由于人工录入的失误），否则我们并不对原始数据进行修改。

巧的是，在这个问题的数据集中真的出现了一个错误值，而且是可以看出如何修改的错误，如下：

test.describe()

可以看到，GarageYrBlt的最大值是2207，而数据集只收集至2010年，因此这个数字显然是错误的，不难猜出，这个数字应该是2007，因此，我们对这个数值进行修改。

test[test['GarageYrBlt'] == 2207]

test.loc[1132,'GarageYrBlt'] = 2007

修改后可以再次列出GarageYrBlt的统计表，发现数值范围正常了。

3、目标值处理

因为我们希望用线性回归的方法来拟合数据，因此我们希望目标值的分布接近正态分布，这基于各个变量是正态分布的假设，多个正态分布加起来还是正态分布。

sns.distplot(train['SalePrice'] , fit=norm)

可以看到我们的 Target 的分布相比正态分布整体是右偏的。画出qq-plot如下：

fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)

为了让 Target 更接近正态分布，我们可以对 Target 进行 log-transformation。

train["SalePrice"] = np.log1p(train["SalePrice"])
sns.distplot(train['SalePrice'] , fit=norm)

看起来更接近正态分布了，我们再看一下qq-plot：

res = stats.probplot(train['SalePrice'], plot=plt)

可以看到蓝色的样本曲线更接近理论的正态分布红线了。

4、特征工程

4.1、缺失值

缺失值是实际的数据集中经常出现的问题，而造成数据缺失的原因也各不相同，有的是因为各种原因无法采集到数据，有的缺失值则代表某个默认值，比如数据文档中写到的关于泳池质量特征的缺失值说明：

在考察缺失值时，比较方便的处理方式是将训练集和测试集合并考虑，当然这里要注意的是，当我们想删除某个有缺失值的 item 的时候，只能删除训练集中的 item，而不能删除测试集中的 item，缺失值的替换可以在两者中同时进行。

ntrain = train.shape[0]
ntest = test.shape[0]
y_train = train.SalePrice.values

all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data[:20]

上表列出了缺失值较多的特征，接下来我们要看一下如何对这些缺失值进行处理。

文档显示，缺失值最多的前 5 个特征的 NA 都表示没有相应的设施，因此可以直接用 None 来填充。

all_data["PoolQC"] = all_data["PoolQC"].fillna("None")
all_data["MiscFeature"] = all_data["MiscFeature"].fillna("None")
all_data["Alley"] = all_data["Alley"].fillna("None")
all_data["Fence"] = all_data["Fence"].fillna("None")
all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")

LotFrontage（房屋到街道的距离）的缺失值可以用此房屋所在区的其它房屋的此指标的中位数来填充，这是合理的。

all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
                          lambda x: x.fillna(x.median()))

文档显示，一系列的 GarageX 特征中的缺失值均 NA 表示当前房屋没有车库。

观察上表可以发现，GarageType 特征的缺失值比其它 GarageX 的缺失值要少。这说明有部分房屋只能确认是有车库的，但车库的信息一概不知……不知道为什么会出现这样的情况。

总之我们可以把一系列的非数值型GarageX特征的缺失值用 None 填充。

for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    all_data[col] = all_data[col].fillna('None')

把一系列的数值型 GarageX 特征的缺失值用 0 填充。

for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)

同理，一系列的 BsmtX 特征的缺失值是没有地下室，因此也用 None 和 0 来填充。

for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    all_data[col] = all_data[col].fillna('None')
    
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    all_data[col] = all_data[col].fillna(0)

MasVnrX 特征的缺失值是没有 Masonry veneer，因此也用 None 和 0 来填充。

all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)

接下来是 MSZoning，这个特征文档给出的解释大概就是地段的意思……这个如何填充呢？我们先来看下地段分布情况。

all_data["MSZoning"].value_counts()

可以看到，地段为 RL 的数量最多，因此我们就把缺失值用 RL 来填充。类似地，我们把 Utilities 特征用其众数进行填充。

all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])
all_data['Utilities'] = all_data['Utilities'].fillna(all_data['Utilities'].mode()[0])

除了上述缺失值较多的特征，还有以下缺失值较少的特征需要处理：

由文档得知，Functional 这个特征默认值是 Typical ，因此我们用 Typ 进行填充。

all_data["Functional"] = all_data["Functional"].fillna("Typ")

接下来几个特征缺失值非常少，我们均用该特征对应的众数进行填充。

all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])
all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])
all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])
all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])

4.2、特征的进一步处理

首先，我们把一些数值型的代表类别的特征转换成字符串形式，以方便我们之后用 dummy categorical features 来代替这些特征：

#MSSubClass=The building class
all_data['MSSubClass'] = all_data['MSSubClass'].astype(str)
#Changing OverallCond into a categorical variable
all_data['OverallCond'] = all_data['OverallCond'].astype(str)
#Year and month sold are transformed into categorical features.
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)

接下来我们把这些代表类别的特征（包括上面我们转成字符串形式的特征）进行类别编码，具体如下：

from sklearn.preprocessing import LabelEncoder
cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond', 'ExterQual',
        'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1', 'BsmtFinType2',
        'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope', 'LotShape',
        'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold')
# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl = LabelEncoder() 
    lbl.fit(list(all_data[c].values)) 
    all_data[c] = lbl.transform(list(all_data[c].values))

这样看或许我们不太清楚 LabelEncoder 做了什么事，简单来说，它就是把代表类别的各个值按顺序标成 0，1……numClass - 1，比如：

le = preprocessing.LabelEncoder()
le.fit([6, 2, 2, 1])
le.transform([1, 1, 2, 6]) 

out:
array([0, 0, 1, 2], dtype=int64)

处理好类别特征，接下来我们要关注那些数值型的特征了，我们希望这部分特征都是服从正态分布的，但实际情况并非这么理想，我们可以通过计算各数值特征的偏度来寻找偏离正态分布较大的特征并对其加以修正。

numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

# Check the skew of all numerical features
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\nSkew in numerical features: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(10)

可以观察到，偏度最大的几个特征的偏度非常大，因此我们要对它们进行正态化处理。

skewness = skewness[abs(skewness) > 0.75]

from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
    #all_data[feat] += 1
    all_data[feat] = boxcox1p(all_data[feat], lam)

这里使用的 box-cox transform 是一种正态化处理的变换，具体见这里。

4.3、去除冗余特征

其实我们往往很难判断一个特征是否是冗余的，除非多个特征之间有明显的相关关系或某个特征的值大体都相同。

这个问题中我们重点关注后者，因为前者可以通过某些手段（比如有特征选择功能的L1正则化）来解决。

上面我们已经对数值特征进行了转换，现在我们关注类别特征。

首先筛选出类别特征并考察类别数：

objects = []
for i in all_data.columns:
    if all_data[i].dtype == object:
        objects.append(i)
        
sums_features = all_data[objects].apply(lambda x: len(np.unique(x)))
sums_features.sort_values(ascending=False)

接下来我们看一下各类别众数对应的占比多少：

all_data_balance = pd.Series((all_data.values == all_data.mode().values).sum(axis=0) / len(all_data)) * 100
all_data_balance = all_data_balance.sort_values(ascending=False)
mode_data = pd.DataFrame({'Mode Ratio' :all_data_balance})
mode_data.index = all_data.columns[mode_data.index]
for c in mode_data.index:
    if c not in objects:
        mode_data.drop(c,inplace=True)
mode_data

可以看到 Utilities 特征的类别数很少且分布非常集中，因此对最终结果影响不大，我们可以将其去除。

all_data = all_data.drop(['Utilities'], axis=1)

5、模型搭建

这部分内容我参考了很多大神的kernel，发现表现最好的基本上都用了一种叫做stacking的技术，stacking和bagging与boosting一样，都是模型的聚合方法，具体实现方法一会再说。

首先，我们直接用 Lasso 回归来解决问题，因为这是一个特征非常多的问题，而 Lasso 可以有效解决特征的选择问题。

具体实现也很简单：

from sklearn.linear_model import Lasso
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler

lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))
Lasso_model = lasso.fit(train.values,y_train)
result = np.expm1(Lasso_model.predict(test.values))

sub = pd.DataFrame()
sub['Id'] = test_ID
sub['SalePrice'] = result
sub.to_csv(r'C:\Users\Sidney Nash\Desktop\data\House\submission.csv',index=False)

提交后的分数是0.11958，排名top 20%。作为一个base model，这样的表现算是不错了，接下来我尝试了将多个模型的预测结果进行简单平均，结果有不小的提升。

具体实现如下：

KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))
ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))
GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,
                                   max_depth=4, max_features='sqrt',
                                   min_samples_leaf=15, min_samples_split=10, 
                                   loss='huber', random_state =5)

LassoMd = lasso.fit(train.values,y_train)
ENetMd = ENet.fit(train.values,y_train)
KRRMd = KRR.fit(train.values,y_train)
GBoostMd = GBoost.fit(train.values,y_train)

finalMd = (np.expm1(LassoMd.predict(test.values)) + np.expm1(ENetMd.predict(test.values)) + np.expm1(KRRMd.predict(test.values)) + np.expm1(GBoostMd.predict(test.values)) ) / 4

sub = pd.DataFrame()
sub['Id'] = test_ID
sub['SalePrice'] = result
sub.to_csv(r'C:\Users\Sidney Nash\Desktop\data\House\submission.csv',index=False)

排名提升到了top 13%，可见对多个模型的结果进行聚合是有效的做法，那么如果我们不用简单平均，而是使用更高级的聚合方法，模型的表现应该会得到进一步提升。下面我们看一下 Stacking 方法。

要理解 Stacking 方法的思路，我们先要明白它的算法流程（图自Here）。

考虑 A、B、C 三个数据集，其中 A 和 B 数据集我们知道目标变量 y，C 则不知道。

stacking可以分为以下四步：

（1）在数据集 A 中使用多种机器学习算法（回归器或分类器）进行训练得到多个模型。

（2）使用（1）中得到的模型对数据集 B 和 C 的 y 进行预测，并创建仅包含这些预测的新数据集 B1 和 C1 。如果我们运行 10 个模型那么 B1 和 C1 各有 10 列。

（3）使用 B1 作为特征，B 数据集的 y 作为目标变量训练一个新的机器学习算法（通常被称为元学习算法）得到元模型。

（4）使用元模型在 C1 上进行预测。

在实际使用 Stacking 算法的时候，我们的 C 就是测试集，而 A 和 B 是由训练集划分得到的，A 用于多种算法的模型训练，B 则作为留出集对多种模型进行融合。

那么这个融合的过程究竟是如何发生的呢？直觉来看，不同算法在数据中的不同样本上的表现是有差异的，举个极端的例子来说，A 算法训练得到的模型在前半部分的样本上正确率100%而在后半部分样本上正确率只有50%，而 B 算法则刚好相反，此时我们希望用 A 算法去对前半部分样本做预测，用 B 算法去对后半部分样本做预测，这就是 A 和 B 的融合。

当然，在实际操作中，我们不可能具体地去考虑各个样本上哪个模型更好，我们只需要把各个样本所做的预测再作为新模型（元模型）的输入，让这个元模型来帮助我们分辨各个模型在各个样本分类时所占的权重，这样一来就达到了各个模型取长补短的效果。

在本问题中，我们使用多个自身表现已经很好的回归模型来进行Stacking，以达到超越任何单个模型表现的新模型。

首先规定误差函数：

#Validation function
n_folds = 5

def rmsle_cv(model):
    kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)
    rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))
    return(rmse)

然后建立多个用于回归的模型，超参选取可以使用网格搜索的方法。

lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))
ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=1))
KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,
                                   max_depth=4, max_features='sqrt',
                                   min_samples_leaf=15, min_samples_split=10, 
                                   loss='huber', random_state =5)
model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, 
                             learning_rate=0.05, max_depth=3, 
                             min_child_weight=1.7817, n_estimators=2200,
                             reg_alpha=0.4640, reg_lambda=0.8571,
                             subsample=0.5213, silent=1,
                             random_state =7, nthread = -1)
model_lgb = lgb.LGBMRegressor(objective='regression',num_leaves=5,
                              learning_rate=0.05, n_estimators=720,
                              max_bin = 55, bagging_fraction = 0.8,
                              bagging_freq = 5, feature_fraction = 0.2319,
                              feature_fraction_seed=9, bagging_seed=9,
                              min_data_in_leaf =6, min_sum_hessian_in_leaf = 11)

输出各模型的误差大小：

score = rmsle_cv(lasso)
print("\nLasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(ENet)
print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(KRR)
print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(GBoost)
print("Gradient Boosting score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(model_xgb)
print("Xgboost score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
score = rmsle_cv(model_lgb)
print("LGBM score: {:.4f} ({:.4f})\n" .format(score.mean(), score.std()))

out:
Lasso score: 0.1123 (0.0074)
ElasticNet score: 0.1123 (0.0075)
Kernel Ridge score: 0.1145 (0.0073)
Gradient Boosting score: 0.1167 (0.0070)
Xgboost score: 0.1153 (0.0076)
LGBM score: 0.1154 (0.0071)

可以看到，表现较好的模型是Lasso和ENet。接下来我们采用Stacking的方式来将各个模型进行融合：

class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds
   
    # We again fit the data on clones of the original models
    def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)
        
        # Train cloned base models then create out-of-fold predictions
        # that are needed to train the cloned meta-model
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X[train_index], y[train_index])
                y_pred = instance.predict(X[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred
                
        # Now train the cloned  meta-model using the out-of-fold predictions as new feature
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self
   
    #Do the predictions of all base models on the test data and use the averaged predictions as 
    #meta-features for the final prediction which is done by the meta-model
    def predict(self, X):
        meta_features = np.column_stack([
            np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
            for base_models in self.base_models_ ])
        return self.meta_model_.predict(meta_features)

stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),
                                                 meta_model = lasso)

score = rmsle_cv(stacked_averaged_models)
print("Stacking Averaged models score: {:.4f} ({:.4f})".format(score.mean(), score.std()))

out:
Stacking Averaged models score: 0.1081 (0.0070)

可以看到，Stacking后的模型表现优于任何单个模型。

最后，我们将Stacking模型和xgb以及lgb线性组合，进一步提升模型的表现，至于各模型系数，可以用网格搜索来得出。

stacked_averaged_models.fit(train.values, y_train)
stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))

model_xgb.fit(train, y_train)
xgb_pred = np.expm1(model_xgb.predict(test))

model_lgb.fit(train, y_train)
lgb_pred = np.expm1(model_lgb.predict(test.values))

ensemble = stacked_pred*0.70 + xgb_pred*0.15 + lgb_pred*0.15

最终得到的模型精度相比融合前有了显著提高。rank 也终于到了top 10%。

以上思路及代码多数参考自serigne大神的Notebook，受益匪浅，在此表示感谢。