<jupyter_start><jupyter_text>Content* Data ingestion    * Remove outliers    * Remove adoption features not in ABM    * Data splitting for CV* Features preparation    * Pipeline preparation    * Features names collection* Model selection    * Functions for models results    * Functions for models' error analyis    * Models training    * Results* Final models tuning and training<jupyter_code>import numpy as np
import pandas as pd
import math
import joblib
import os
import custom_transformers
from importlib import reload
reload(custom_transformers)
from custom_transformers import (
    TransformAdoptionFeatures,
    TransformCensusFeatures, 
    TransformClimateFeatures,
    TransformSoilFeatures,
    TransformEconomicFeatures,
)
# For plotting
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
# Where to save the figures
PROJECT_ROOT_DIR = "."
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "images")
os.makedirs(IMAGES_PATH, exist_ok=True)

def save_fig(fig_id, tight_layout=True, fig_extension="jpeg", resolution=300):
    path = os.path.join(IMAGES_PATH, fig_id + "." + fig_extension)
    print("Saving figure", fig_id)
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format=fig_extension, dpi=resolution)
 # How to display floats
pd.options.display.float_format = '{:.6f}'.format
np.set_printoptions(precision=10)<jupyter_output><empty_output><jupyter_text>Data ingestion<jupyter_code>path_to_dataset = "Municipalities final dataset for analysis.csv"
dataset_original = pd.read_csv(path_to_dataset, index_col=['Municipality', 'Year'])
# Remove features with 0 adoption to obtain the regression dataset
mask_no_adoption = dataset_original['adoption_in_year'] < 0.0000001
dataset_original = dataset_original.drop(dataset_original.loc[mask_no_adoption].index)
len(dataset_original)
# Remove outliers adoption_in_year
outliers = dataset_original.loc[dataset_original['adoption_in_year'] > 0.12].index
dataset_original = dataset_original.drop(outliers)
print(outliers)
# Remove outliers cumul_adoption_10_y_pr_y_munic
outliers = dataset_original.loc[dataset_original['tot_cumul_adoption_pr_y_munic'] > 0.4].index
dataset_original = dataset_original.drop(outliers)
print(outliers)
dataset_original.shape
dataset = dataset_original.reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>Create a mapping of each new instance to municipality and year it's referring to<jupyter_code>map_instances_to_munic_year = dict(zip(dataset.index, dataset_original.index))<jupyter_output><empty_output><jupyter_text>Final number of instances and features<jupyter_code>dataset.shape<jupyter_output><empty_output><jupyter_text>Copy the data splitting label and features<jupyter_code>dataset_labels = dataset['adoption_in_year'].copy()
dataset = dataset.drop('adoption_in_year', axis=1)
dataset_labels.describe()<jupyter_output><empty_output><jupyter_text>Split for cross-validation Split used in all the cross-validation, stratified based on the number of data points before the PCF project and during it<jupyter_code>from sklearn.model_selection import StratifiedKFold, PredefinedSplit<jupyter_output><empty_output><jupyter_text>First, generate the feature to split, the year of the instance<jupyter_code>dataset_split = dataset_original.copy()
years = dataset_split.index.get_level_values(1).to_list()
dataset_split['Year'] = years<jupyter_output><empty_output><jupyter_text>Second, get for each instance the test fold it has to belong to based on the split based on the created category<jupyter_code>n_folds = 3
strat_split = StratifiedKFold(n_splits=n_folds, shuffle=True)
test_indexes = np.zeros(len(dataset_split))
fold = 0
for _, test_index in strat_split.split(dataset_split, dataset_split['Year']):
    for i in test_index:
        test_indexes[i] = fold
    fold += 1<jupyter_output><empty_output><jupyter_text>Third, pass to the list of test folds indexes to generate the PredefinedSplit<jupyter_code>cv_split = PredefinedSplit(test_indexes)
# Check ratio is kept
ratio_PCF_dataset = len(dataset_split.loc[(slice(None), [2007]), :]) / len(dataset_split)
test_set = dataset_split.iloc[np.where(test_indexes==1)]
ratio_PCF_test_set = len(test_set.loc[(slice(None), [2007]), :]) / len(test_set)
print(ratio_PCF_dataset, ratio_PCF_test_set)<jupyter_output>0.07454175152749491 0.0745721271393643<jupyter_text>Features preparation  Features names collection  Pre-transformation<jupyter_code>features_original = [feat for feat in dataset.columns]
features_adoption = [feat for feat in features_original if 'adoption' in feat]
features_climate = [feat for feat in features_original if '_t_' in feat or '_prec_' in feat]
features_soil = ['CaCO3_mean_munic', 'CN_mean_munic', 'N_mean_munic', 'P_mean_munic', 'pH_mean_munic']
features_economic = ['sbp_payment']

features_not_census = (features_adoption + features_climate + features_soil + features_economic)
features_census = [feat for feat in features_original if feat not in features_not_census]
len(features_census) + len(features_not_census)<jupyter_output><empty_output><jupyter_text>Post-transformation<jupyter_code># Get features names
feat_lists = [features_adoption, features_census, features_climate, features_soil, features_economic]
transformers = [TransformAdoptionFeatures(), TransformCensusFeatures(),
                TransformClimateFeatures(), TransformSoilFeatures(), TransformEconomicFeatures()]

features = []
for f, t in zip(feat_lists, transformers):
    features_ind = t.transform(dataset[f]).columns.to_list()
    features = features + features_ind
len(features)<jupyter_output><empty_output><jupyter_text>With squared cumulative adoption<jupyter_code># Get features names
feat_lists = [features_adoption, features_census, features_climate, features_soil, features_economic]
transformers = [TransformAdoptionFeatures(add_square_cumul_adoption=True),
                TransformCensusFeatures(), TransformClimateFeatures(), TransformSoilFeatures(), TransformEconomicFeatures()]

features_sq = []
for f, t in zip(feat_lists, transformers):
    features_ind = t.transform(dataset[f]).columns.to_list()
    features_sq = features_sq + features_ind
len(features_sq)<jupyter_output><empty_output><jupyter_text>Preparation pipeline<jupyter_code>from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler
transformation_pipeline = ColumnTransformer([
    ('adoption', TransformAdoptionFeatures(), features_adoption),
    ('census', TransformCensusFeatures(), features_census),
    ('climate', TransformClimateFeatures(), features_climate),
    ('soil', TransformSoilFeatures(), features_soil),
    ('economic', TransformEconomicFeatures(), features_economic)
],
    remainder='drop')

dataset_transf = transformation_pipeline.fit_transform(dataset)
full_pipeline_stand = Pipeline([
    ('transformation', transformation_pipeline),
    ('scaling', StandardScaler())
])
dataset_prepared_stand = full_pipeline_stand.fit_transform(dataset)
dataset_prepared_stand.shape
# Print prepared data as a DataFrame
pd.DataFrame(dataset_prepared_stand, columns=features, index=dataset.index).head()<jupyter_output><empty_output><jupyter_text>With squared cumulative adoption<jupyter_code>transformation_pipeline_sq = ColumnTransformer([
    ('adoption', TransformAdoptionFeatures(add_square_cumul_adoption=True),
     features_adoption),
    ('census', TransformCensusFeatures(), features_census),
    ('climate', TransformClimateFeatures(), features_climate),
    ('soil', TransformSoilFeatures(), features_soil),
    ('economic', TransformEconomicFeatures(), features_economic)
],
    remainder='drop')
full_pipeline_stand_sq = Pipeline([
    ('transformation', transformation_pipeline_sq),
    ('scaling', StandardScaler())
])
dataset_stand_sq_cum_ad = full_pipeline_stand_sq.fit_transform(dataset)
dataset_stand_sq_cum_ad.shape<jupyter_output><empty_output><jupyter_text>Model selection * Linear regression with Elastic Net regularization    * With squared cumulative adoption    * Polynomial* SVM    * Linear    * With squared cumulative adoption    * Non-linear<jupyter_code>from sklearn.model_selection import cross_validate, cross_val_predict
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.base import clone
from sklearn.model_selection import RandomizedSearchCV, GridSearchCV, StratifiedKFold
from scipy.stats import expon, uniform, reciprocal<jupyter_output><empty_output><jupyter_text>Functions for models results<jupyter_code>cols = ['cv RMSE_mean', 'cv RMSE std', 'cv adjR2 mean', 'cv adjR2 std', 'train set RMSE', 'train set adjR2']
models_results = pd.DataFrame(columns=cols, dtype=np.float64)
def display_scores(scores, score_name):
    print('Cross-validation', score_name)
    print('Scores:', scores)
    print('Mean:', scores.mean())
    print('Standard deviation:', scores.std())
    print('')
def adjusted_r2(r2, n, p):
    return 1-(1-r2)*(n-1)/(n-p-1)
from sklearn.model_selection import KFold

def cross_val_rmse_r2(reg, dataset_prepared, cv_split):
    scores = cross_validate(reg, dataset_prepared, dataset_labels,
                            scoring=['neg_mean_squared_error', 'r2'], cv=cv_split)
    rmse_scores = np.sqrt(-scores['test_neg_mean_squared_error'])
    r2_scores = scores['test_r2']
    adj_r2_scores = adjusted_r2(r2_scores, (dataset_prepared.shape[0] // n_folds), dataset_prepared.shape[1])
    cv_scores = {'RMSE': rmse_scores, 'adjR2': adj_r2_scores}
    return cv_scores
def predict_train_set(reg, dataset_prepared):
    reg.fit(dataset_prepared, dataset_labels)
    predictions = reg.predict(dataset_prepared)
    mse = mean_squared_error(dataset_labels, predictions)
    rmse = np.sqrt(mse)
    r2 = r2_score(dataset_labels, predictions)
    adj_r2 = adjusted_r2(r2, dataset_prepared.shape[0], dataset_prepared.shape[1])
    print('Scores on training set')
    print('RMSE score:', rmse)
    print('Adjusted R2 score:', adj_r2)
    return (rmse, adj_r2)
def load_model_results(reg, reg_name, results_df, cv_scores, train_results):
    df_cols = results_df.columns
    
    # Add cv columns
    for score_name, score_values in cv_scores.items():
        col_mean = [col for col in df_cols if ('cv' in col) and (score_name in col) and ('mean' in col)]
        col_std = [col for col in df_cols if ('cv' in col) and (score_name in col) and ('std' in col)]
        results_df.loc[reg_name, col_mean] = score_values.mean()
        results_df.loc[reg_name, col_std] = score_values.std()
    # Add training set results columns
    cols_train = [col for col in df_cols if ('train' in col)]
    results_df.loc[reg_name, cols_train] = train_results
def test_regressor(reg, dataset_prepared, labels, reg_name=None, results_df=models_results, cv_split=cv_split):  
    """
    Function to:
    - train a regressor using cross validation, reporting both RMSE and adjusted R2 score
    - train a regressor on the training set, to check it for overfiting
    - load the model performances in a dataframe to keep all models results
    
    reg_name: how to report the model's results in the DataFrame with the results. If None, do not add it.
    """
    cv_scores = cross_val_rmse_r2(reg, dataset_prepared, cv_split)
    print("Cross validation scores")
    for score_name, score_value in cv_scores.items():
        display_scores(score_value, score_name)
    print("")

    train_results = predict_train_set(reg, dataset_prepared)
        
    if reg_name != None:
        load_model_results(reg, reg_name, results_df, cv_scores, train_results)
def get_grid_search_results(gs):
    best_est = gs.best_estimator_
    print("Best estimator:", best_est)
    print("Best parameters:", gs.best_params_)
    for scorer_name in gs.scorer_.keys():
        res_test = gs.cv_results_['mean_test_' + scorer_name][gs.best_index_]
        res_train = gs.cv_results_['mean_train_' + scorer_name][gs.best_index_]
        print("")
        print(scorer_name + ' mean on test sets:', res_test)
        print(scorer_name + ' mean on train sets:', res_train)<jupyter_output><empty_output><jupyter_text>Functions for models' analysis<jupyter_code>def plot_validation_curves(estimator, param_name, param_values, scoring, dataset, labels):
    train_errors, val_errors = [], []
    for param in param_values:
        estimator.set_params(**{param_name:param})
        scores = cross_validate(estimator, dataset, labels, scoring=scoring, 
                                cv=cv_split, n_jobs=-1, return_train_score=True)
        train_errors.append(-scores['train_score'].mean())
        val_errors.append(-scores['test_score'].mean())

    plt.plot(val_errors, "b-", linewidth=2, label="Validation set")
    plt.plot(train_errors, "r--", linewidth=2, label="Training set")
    plt.legend(loc="upper right", fontsize=12)
    x = np.arange(0, len(param_values), 1)
    plt.xticks(x, np.round(param_values, 2))
    plt.xlabel(param_name + " value", fontsize=14)
    plt.ylabel("RMSE", fontsize=14)
    plt.show()
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

def plot_learning_curves(model, X, y, fig_name=None):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)
    train_errors, val_errors = [], []
    steps = np.arange(100, len(X_train), 100)
    for m in steps:
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))

    plt.plot(steps, np.sqrt(train_errors), "r-+", linewidth=2, label="Training error")
    plt.plot(steps, np.sqrt(val_errors), "b-", linewidth=3, label="Validation error")
    plt.legend(loc="lower right", fontsize=13)   
    plt.xlabel("Training set size", fontsize=14)
    plt.ylabel("Root mean squared error", fontsize=14)
    
    if fig_name != None:
        save_fig((fig_name))<jupyter_output><empty_output><jupyter_text>Linear regression (elastic net)<jupyter_code>from sklearn.linear_model import SGDRegressor
el_net_lin_reg = SGDRegressor(penalty='elasticnet', learning_rate='adaptive')

param_grid = {
    #"tol": [5e-3],
    #"max_iter": [10000],
    "alpha": [0.001, 0.005, 0.01, 0.05, 0.1, 1],
    "l1_ratio": np.arange(0.0, 1.0, 0.25),
    }

srch_lin_reg = GridSearchCV(el_net_lin_reg, param_grid, cv=cv_split,
                            scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                            return_train_score=True, verbose=1)

srch_lin_reg.fit(dataset_prepared_stand, dataset_labels)
get_grid_search_results(srch_lin_reg)
lin_reg = srch_lin_reg.best_estimator_
test_regressor(lin_reg, dataset_prepared_stand, dataset_labels)<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.0071828646 0.0060763936 0.0069696493]
Mean: 0.00674296917718219
Standard deviation: 0.00047931020164187333

Cross-validation adjR2
Scores: [0.343639584  0.3209116221 0.4101104098]
Mean: 0.35822053864617304
Standard deviation: 0.03784670053568197


Scores on training set
RMSE score: 0.006646752007238695
Adjusted R2 score: 0.394159299852771<jupyter_text>Second iteration<jupyter_code>param_grid = {
    'max_iter': [3000],
    "alpha": expon(scale=0.05),
    "l1_ratio": uniform(0.0, 1.0),
    }

srch_lin_reg_2 = RandomizedSearchCV(el_net_lin_reg, param_grid, cv=cv_split, n_iter=1000,
                                    scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                    return_train_score=True, verbose=1)

srch_lin_reg_2.fit(dataset_prepared_stand, dataset_labels)
get_grid_search_results(srch_lin_reg_2)
lin_reg = srch_lin_reg_2.best_estimator_
test_regressor(lin_reg, dataset_prepared_stand, dataset_labels, 'Linear Regression')<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.00715707   0.006184041  0.0068609503]
Mean: 0.006734020451968477
Standard deviation: 0.00040725072911261393

Cross-validation adjR2
Scores: [0.3483452702 0.2966374779 0.4283668194]
Mean: 0.35778318912806245
Standard deviation: 0.05419077774752589


Scores on training set
RMSE score: 0.00661531401041517
Adjusted R2 score: 0.39987679159643863<jupyter_text>Linear Regression with squared cumulative adoption<jupyter_code>from sklearn.linear_model import SGDRegressor
el_net_lin_reg = SGDRegressor(penalty='elasticnet', learning_rate='adaptive')
param_grid = {
    'max_iter': [3000],
    "alpha": uniform(0.001, 0.05),
    "l1_ratio": uniform(0.0, 1.0),
    }

srch_lin_reg_sq = RandomizedSearchCV(el_net_lin_reg, param_grid, cv=cv_split, n_iter=1000,
                                    scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                    return_train_score=True, verbose=1)

srch_lin_reg_sq.fit(dataset_stand_sq_cum_ad, dataset_labels)
get_grid_search_results(srch_lin_reg_sq)
lin_reg_sq = srch_lin_reg_sq.best_estimator_
name = "Linear regression + cum_ad_sq"
test_regressor(lin_reg_sq, dataset_stand_sq_cum_ad, dataset_labels, name)<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.0071190945 0.0059461663 0.0068372762]
Mean: 0.00663417900046312
Standard deviation: 0.0004999176921725496

Cross-validation adjR2
Scores: [0.3544333164 0.3488918044 0.4315926244]
Mean: 0.37830591504412964
Standard deviation: 0.037747248289145115


Scores on training set
RMSE score: 0.006480631265784982
Adjusted R2 score: 0.4238274240680121<jupyter_text>Polynomial regression<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False) # Already added by Elastic Net
full_pipeline_stand_poly = Pipeline([
    ('transformation', transformation_pipeline),
    ('poly_features', poly_features),
    ('scaling', StandardScaler())
])
dataset_prepared_stand_poly = full_pipeline_stand_poly.fit_transform(dataset)
dataset_prepared_stand_poly.shape<jupyter_output><empty_output><jupyter_text>Hyperparameters search on polynomial regression with elastic net regularization<jupyter_code>from sklearn.linear_model import SGDRegressor

el_net_pol_reg = SGDRegressor(penalty='elasticnet', learning_rate='optimal') # Optimal much faster and converging

param_grid = {
    "max_iter": [3000],
    "alpha": [0.01, 0.1, 1],
    "l1_ratio": np.arange(0.0, 1.0, 0.25),
    }

srch_pol_reg = GridSearchCV(el_net_pol_reg, param_grid, cv=cv_split,
                                   scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                   return_train_score=True, verbose=1)

srch_pol_reg.fit(dataset_prepared_stand_poly, dataset_labels)
get_grid_search_results(srch_pol_reg)
pol_reg = srch_pol_reg.best_estimator_
test_regressor(pol_reg, dataset_prepared_stand_poly, dataset_labels)<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [8.9713978295e-03 7.4630010512e-03 4.2161882702e+09]
Mean: 1405396090.0876048
Standard deviation: 1987530211.09639

Cross-validation adjR2
Scores: [-3.4389878683e-01 -3.4450002314e-01 -2.8332714728e+23]
Mean: -9.444238242578697e+22
Standard deviation: 1.3356169808937437e+23


Scores on training set
RMSE score: 0.008572823443172329
Adjusted R2 score: -0.09316312968730567<jupyter_text>Second iteration<jupyter_code>param_grid = {
    'max_iter': [3000],
    "alpha": expon(scale=0.01),
    "l1_ratio": uniform(0, 0.5),
    }

srch_pol_reg_2 = RandomizedSearchCV(el_net_pol_reg, param_grid, cv=cv_split, n_iter=100,
                                        scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                        return_train_score=True, verbose=1)

srch_pol_reg_2.fit(dataset_prepared_stand_poly, dataset_labels)
get_grid_search_results(srch_pol_reg_2)
pol_reg = srch_pol_reg_2.best_estimator_
name = 'Polynomial Regression'
test_regressor(pol_reg, dataset_prepared_stand_poly, dataset_labels, name)<jupyter_output>C:\Users\giaco\anaconda3\envs\data_analysis\lib\site-packages\sklearn\linear_model\_stochastic_gradient.py:1211: ConvergenceWarning:

Maximum number of iteration reached before convergence. Consider increasing max_iter to improve the fit.

C:\Users\giaco\anaconda3\envs\data_analysis\lib\site-packages\sklearn\linear_model\_stochastic_gradient.py:1211: ConvergenceWarning:

Maximum number of iteration reached before convergence. Consider increasing max_iter to improve the fit.<jupyter_text>Linear SVM<jupyter_code>from sklearn.svm import LinearSVR
lin_svr = LinearSVR(max_iter=3000)

param_grid = {
    "epsilon": [0, 0.001, 0.01, 1],
    "C": [0.001, 0.01, 0.1, 1, 10],
    "loss": ["epsilon_insensitive", "squared_epsilon_insensitive"]
    }

srch_lin_svr = GridSearchCV(lin_svr, param_grid, cv=cv_split,
                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                  return_train_score=True, verbose=1)

srch_lin_svr.fit(dataset_prepared_stand, dataset_labels)
get_grid_search_results(srch_lin_svr)
lin_svr = srch_lin_svr.best_estimator_
test_regressor(lin_reg_sq, dataset_prepared_stand, dataset_labels)<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.0071523686 0.0062182116 0.0068435882]
Mean: 0.00673805609903981
Standard deviation: 0.0003886001266431794

Cross-validation adjR2
Scores: [0.3492011314 0.2888430082 0.431256272 ]
Mean: 0.35643347053671093
Standard deviation: 0.058364455372327304


Scores on training set
RMSE score: 0.0066132995885728315
Adjusted R2 score: 0.4002422216384929<jupyter_text>Second iteration<jupyter_code>param_grid = {
    "epsilon": expon(scale=0),
    "C": expon(scale=0.1),
    "loss": ["squared_epsilon_insensitive"]
    }

srch_lin_svr_2 = RandomizedSearchCV(lin_svr, param_grid, cv=cv_split, n_iter=50,
                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                  return_train_score=True, verbose=1)

srch_lin_svr_2.fit(dataset_prepared_stand, dataset_labels)
lin_svr = srch_lin_svr.best_estimator_
test_regressor(lin_reg_sq, dataset_prepared_stand, dataset_labels, 'Linear SVM')<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.0071596188 0.0062368827 0.0068515078]
Mean: 0.006749336444547618
Standard deviation: 0.0003835706955308113

Cross-validation adjR2
Scores: [0.3478810461 0.2845658732 0.429939179 ]
Mean: 0.3541286994207972
Standard deviation: 0.05951260030185549


Scores on training set
RMSE score: 0.006613492543974073
Adjusted R2 score: 0.4002072230150402<jupyter_text>Linear SVM with squared cumulative adoption<jupyter_code>param_grid = {
    "epsilon": expon(scale=0),
    "C": expon(scale=0.1),
    "loss": ["squared_epsilon_insensitive"]
    }

srch_lin_svr_sq_2 = RandomizedSearchCV(lin_svr, param_grid, cv=cv_split, n_iter=50,
                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',
                                  return_train_score=True, verbose=1)

srch_lin_svr_sq_2.fit(dataset_stand_sq_cum_ad, dataset_labels)
get_grid_search_results(srch_lin_svr_sq_2)
lin_svr_sq = srch_lin_svr_sq_2.best_estimator_
name = "Linear SVM + cum_ad_sq"
test_regressor(lin_svr_sq, dataset_stand_sq_cum_ad, dataset_labels, name)<jupyter_output>Cross validation scores
Cross-validation RMSE
Scores: [0.0071252287 0.0059629713 0.0068568384]
Mean: 0.0066483461468028705
Standard deviation: 0.0004968650349574662

Cross-validation adjR2
Scores: [0.3533203372 0.3452062709 0.4283354185]
Mean: 0.3756206755392768
Standard deviation: 0.0374218529142985


Scores on training set
RMSE score: 0.006467180178727934
Adjusted R2 score: 0.4262167297157682<jupyter_text>Non-linear SVM<jupyter_code>from sklearn.svm import SVR<jupyter_output><empty_output><jupyter_text>nl_svr = SVR()test_model(nl_svr, None, dataset_prepared_stand)<jupyter_code>### GridSearchCV<jupyter_output><empty_output><jupyter_text>from sklearn.model_selection import KFold nl_svr = SVR()param_grid = {    "degree": [2, 3],    "epsilon": [0.0001, 0.001, 0.1],    "C": [0.001, 0.01, 0.1, 1],    "gamma": [0.001, 0.01, 0.1, 1],    } srch_nl_svr = GridSearchCV(nl_svr, param_grid, cv=cv_split,                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',                                  return_train_score=True, verbose=1)srch_nl_svr.fit(dataset_prepared_stand, dataset_labels) get_grid_search_results(srch_nl_svr) best_nl_svr = srch_nl_svr.best_estimator_test_regressor(best_nl_svr, dataset_prepared_stand, dataset_labels)<jupyter_code>### Second iteration<jupyter_output><empty_output><jupyter_text>nl_svr = SVR()param_grid = {    "degree": [2],    "epsilon": expon(scale = 0.001),    "C": expon(scale = 0.1),    "gamma": expon(scale=0.01)    }srch_nl_svr_2 = RandomizedSearchCV(nl_svr, param_grid, cv=cv_split, n_iter=300,                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',                                  return_train_score=True, verbose=1)srch_nl_svr_2.fit(dataset_prepared_stand, dataset_labels) get_grid_search_results(srch_nl_svr_2) best_nl_svr = srch_nl_svr_2.best_estimator_test_regressor(best_nl_svr, dataset_prepared_stand, dataset_labels)<jupyter_code>### Third iteration<jupyter_output><empty_output><jupyter_text>get_grid_search_results(srch_nl_svr)print("")get_grid_search_results(srch_nl_svr_2) from sklearn.svm import SVRnl_svr = SVR()param_grid = {    "degree": [2],    "epsilon": expon(scale = 0.0017),    "C": expon(scale = 0.26),    "gamma": expon(scale=0.004)    } srch_nl_svr_3 = RandomizedSearchCV(nl_svr, param_grid, cv=cv_split, n_iter=500,                                  scoring=['neg_root_mean_squared_error', 'r2'], refit='neg_root_mean_squared_error',                                  return_train_score=True, verbose=1)srch_nl_svr_3.fit(dataset_prepared_stand, dataset_labels) get_grid_search_results(srch_nl_svr_3) best_nl_svr = srch_nl_svr_3.best_estimator_name = "Non-linear SVM"test_regressor(best_nl_svr, dataset_prepared_stand, dataset_labels, name)<jupyter_code>plot_validation_curves(best_nl_svr, "C", np.arange(0.01, 0.15, 0.01), 'neg_root_mean_squared_error',
                       dataset_prepared_stand, dataset_labels)<jupyter_output><empty_output><jupyter_text>Models summary<jupyter_code>dataset_labels.describe()
models_results
models_results[['cv RMSE_mean', 'cv adjR2 mean', 'train set RMSE', 'train set adjR2']]<jupyter_output><empty_output><jupyter_text>Final models tuning and trainingUsing the Optuna library for Bayesian Optimization<jupyter_code>from sklearn.model_selection import cross_val_score 
import optuna 
from optuna.samplers import TPESampler
from sklearn.svm import SVR
X = dataset_prepared_stand

def objective(trial):
    # Define the search space
    kernel = trial.suggest_categorical('kernel', ['rbf', 'poly'])
    C = trial.suggest_loguniform('C', 0.01, 10)
    gamma = trial.suggest_loguniform('gamma', 0.001, 1)
    epsilon = trial.suggest_loguniform('epsilon', 0.001, 1)

    reg = SVR(degree=2,
              kernel=kernel,
              C=C,
              gamma=gamma,
              epsilon=epsilon)
    score = cross_val_score(reg, X, dataset_labels, scoring='neg_root_mean_squared_error',
                            cv=cv_split).mean()

    return score
study_nl_svm_all = optuna.create_study(study_name="nl_svr_optimization_all",
                                        direction="maximize",
                                        sampler=TPESampler())
%%time
# pass the objective function to method optimize()

study_nl_svm_all.optimize(objective, n_trials=2000)
print(study_nl_svm_all.best_params)
print(-study_nl_svm_all.best_value)
nl_svm_all = SVR(degree=2).set_params(**study_nl_svm_all.best_params)
test_regressor(nl_svm_all, dataset_prepared_stand, dataset_labels)
plot_learning_curves(nl_svm_all, dataset_prepared_stand, dataset_labels)<jupyter_output><empty_output><jupyter_text>Export model<jupyter_code>ml_model = nl_svm_all
feats = features
pipeline = transformation_pipeline<jupyter_output><empty_output><jupyter_text>For ABM<jupyter_code>import csv
from numpy import savetxt 

def export_model_for_abm(ml_model, feats, pipeline, folder_path):
    scaler = StandardScaler() # Needs to be instantiated to be passed to the ABM if not fit in its initialization
    dataset_only_transf = pipeline.fit_transform(dataset) ##
    scaler.fit(dataset_only_transf)
    model_for_abm = Pipeline([
        ("scaling", scaler),
        ("regressor", ml_model)
    ])
    joblib.dump(model_for_abm, os.path.join(folder_path, 'model.pkl'))
    # Check
    print(model_for_abm.predict(dataset_only_transf[0, :].reshape(1, -1))[0])
    
    #Export features name
    with open(os.path.join(folder_path, 'features.csv'), "w+") as file:
        wr = csv.writer(file)
        wr.writerow(feats)
    
    savetxt(os.path.join(folder_path, 'dataset.csv'), dataset_only_transf, delimiter=',')
    savetxt(os.path.join(folder_path, 'labels.csv'), dataset_labels.values,  delimiter=',') ##
folder_path = './Models for ABM/all_adoption/regressor'
export_model_for_abm(ml_model, feats, pipeline, folder_path)<jupyter_output>0.0018124758813949197<jupyter_text>For model analysis<jupyter_code>joblib.dump(ml_model, './Models to analyse/all_adoption/regressor/nl_svm.pkl')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10212665724476748, 'fraction_numerical': 0.057406285752742446, 'mean_word_length': 5.424340683095547, 'pattern_counts': {'":': 34, '<': 136, '<?xml version=': 0, '>': 137, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '36', 'dir': 'jupyter-structured-clean-dedup', 'id': '186667', 'license': 'MIT', 'n_tokens_mistral': 11529, 'n_tokens_neox': 10118, 'n_words': 1870, 'path': 'ml_models/Regression - Models evaluation.ipynb', 'repo_name': 'giacrava/sbp-adoption-abm', 'size': '29718'}