Sequential Ensemble Model (Voting, Stacking generalization) of Precipitation Downscaling¶

# reading data
data= pd.read_csv("sample_data.csv", index_col=["time"], parse_dates=[0])
y = data["Precipitation"]
X = data.drop(["Precipitation"], axis=1)
print(X.shape, y.shape)
data.head()

# defining training and testing time 
train_period = pd.date_range(start="1958-01-01", end="2000-12-31", freq="MS")
test_period = pd.date_range(start="2001-01-01", end="2019-12-31", freq="MS")
full_time = pd.date_range(start="1958-01-01", end="2019-12-31", freq="MS")
print(train_period.shape, test_period.shape)

#transforming data after feature selection 
rfcv= rfcv.fit(X_train, y_train)
X_train_new = rfcv.transform(X_train)
X_test_new = rfcv.transform(X_test)
print(X_train_new.shape)
print("Optimal number of predictors:", rfcv.n_features_)
print("Selected predictors:", X_train.columns[rfcv.support_])
print("Feature selection score:", rfcv.cv_results_["mean_test_score"].mean())

# checking the performance of the individual models
estimators = [("lassoCV", lassoCV), ("ARD", ARD), ("Gamma", gamma), ("MLP", MLP), 
                                        ("RandomForest", randomforest), 
             ("Dense", keras_reg), ("Lars",lassoLarsCV)]


for reg in (lassoCV, ARD, gamma, MLP, randomforest, lassoLarsCV, bayesianRidge, poissonRegressor, 
           decisiontree, extratree):
    reg.fit(X_train_new, y_train)
    yhat = reg.predict(X_test_new)
    print(reg.__class__.__name__, mean_squared_error(y_test, yhat, squared=False))
    print("R² score: {:.2f}".format(r2_score(y_test, yhat)))

#combinging level 1 output in level 2  GradientBoostingRegressor,   GradientBoostingRegressor, StackingRegressor
# AdaBoostRegressor, ExtraTreesRegressor, RandomForestRegressor, VotingRegressor, XGBoost

#set display of models 
set_config(display="diagram")

# ensemble types (based on boasting ) voting or stacking generalization
# different types of final estimator 
# tree-based 
extra_tree = ExtraTreesRegressor(random_state=42, criterion="squared_error")
#gradient boosting 
gradient_boost = GradientBoostingRegressor(n_estimators=25, subsample=0.5, min_samples_leaf=25, random_state=42)
#HistGradientBoosting
histgradient_boost = HistGradientBoostingRegressor(random_state=42)
#XGBoosting
xgboost_reg = XGBRegressor(random_state=42)
#adaboost
adaboost_reg = AdaBoostRegressor(random_state=42, n_estimators=100)

#voting regressor
#voting_reg = VotingRegressor(estimators=estimators) # default estimator: RidgeCV

#Random Forest
randomforest_reg = RandomForestRegressor(random_state=42)

for estimator,name in zip([extra_tree, gradient_boost, histgradient_boost, xgboost_reg, adaboost_reg, 
                          randomforest_reg], ["ExtraTree", "GradientBoosting", "HistGradientBoosting", "XGBoosting",
                                              "AdaBoost", "RandomForest"]):
    stacked_regressor = StackingRegressor(estimators=estimators, final_estimator=estimator)
    stacked_regressor.fit(X_train_new, y_train)
    stacked_regressor
    yhat = stacked_regressor.predict(X_test_new)
    print(name, mean_squared_error(y_test, yhat, squared=False))
    print("R² score: {:.2f}".format(r2_score(y_test, yhat)))
    

stacked_regressor

#prediciton and visualizaiton
plt.style.use("ggplot")
plt.rcParams.update({"font.size":22, "font.weight":"bold"})
from scipy import stats
from sklearn.linear_model import LinearRegression

#model = LinearRegression(fit_intercept=False)
#model.fit(yhat, y_test)
#reg_line = model.predict(yhat)
#score = model.score()
reg_stat = stats.linregress(yhat, y_test)
reg_line = reg_stat.slope*yhat

fig,ax = plt.subplots(figsize=(12,8))
ax.scatter(yhat, y_test, edgecolor='k', facecolor='grey', alpha=0.7,)
#ax.plot(yhat, reg_line, color="red", label="R²= {:.2f}".format(r2_score(y_test, yhat)))
ax.plot(yhat, reg_line, color="red", label="R²= {:.2f}".format(reg_stat.rvalue))
ax.set_ylabel("predicted values")
ax.set_xlabel("observed values")
ax.legend(frameon=True)

#residuals plot (To check for overfitting)--including linear models deviates the residuals form normal distribution (which justifies that
# the relation between the large scale and local precipitation is non-linear). Therefore, models like tree-based algorithms, bayesian regression
# and deep learning models are feasible for modeling the transfer function
plt.style.use("ggplot")
plt.rcParams.update({"font.size":22, "font.weight":"bold"})

from yellowbrick.regressor import ResidualsPlot
model = stacked_regressor
visualizer=ResidualsPlot(model, hist=True, qqplot=False)
visualizer.fit(X_train_new, y_train)
visualizer.score(X_test_new, y_test)
visualizer.show()

# model evaluaiton (regression meterics) mean_squared_error, r2_score, explained_variance_score, 
# max_error, mean_absolute_error, mean_squared_log_error, mean_absolute_percentage_error

#explained_variance 
score = explained_variance_score(y_test, yhat)
print("Explained Variance: {:.2f}".format(score))

#r2_score 
score = r2_score(y_test, yhat)
print("R² (Coefficient of determinaiton): {:.2f}".format(score))

#mean_squared_error 
score = mean_squared_error(y_test, yhat, squared=False)
print("MSE): {:.2f}".format(score))

#mean_absolute_error 
score = mean_absolute_error(y_test, yhat)
print("MAE): {:.2f}".format(score))

#RMSE
error = y_test -yhat
score = np.sqrt(np.mean(error **2))
print("RMSE: {:2f}".format(score))

# Nash-Sutcliffe Efficiency (NSE)
score = (1-(np.sum((yhat-y_test)**2))/np.sum((y_test - np.mean(y_test))**2))
print("NSE: {:.2f}".format(score))