<jupyter_start><jupyter_text>ML final project - Ranking Features for Landslide Prediction Alessio Poloni  Part I: Load data This cell contains an initial set up.*Source: Prof. <NAME>'s notebooks.*<jupyter_code># Python ≥3.5 is required
import sys
assert sys.version_info >= (3, 5)

# Is this notebook running on Colab or Kaggle?
IS_COLAB = "google.colab" in sys.modules

# Scikit-Learn ≥0.20 is required
import sklearn
assert sklearn.__version__ >= "0.20"

# Common imports
import numpy as np
import os

# to make this notebook's output stable across runs
rnd_seed = 42
rnd_gen = np.random.default_rng(rnd_seed)

# To plot pretty figures
%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 = "."
CHAPTER_ID = "classification"
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID)
os.makedirs(IMAGES_PATH, exist_ok=True)

def save_fig(fig_id, tight_layout=True, fig_extension="png", 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)<jupyter_output><empty_output><jupyter_text>The necessary dataset is stored on Google Drive. The next cells allow the data to be retrieved and loaded into the notebook.*Source: <NAME> Jr, PhD.* Installation of required imports:<jupyter_code>!pip install -U -q PyDrive
from pydrive.auth import GoogleAuth
from pydrive.drive import GoogleDrive
from google.colab import auth
from oauth2client.client import GoogleCredentials<jupyter_output><empty_output><jupyter_text>Authentication with Google account (when a window opens -> allow):<jupyter_code>auth.authenticate_user()
gauth = GoogleAuth()
gauth.credentials = GoogleCredentials.get_application_default()
drive = GoogleDrive(gauth)<jupyter_output><empty_output><jupyter_text>Creation of the file and import of .csv data file:<jupyter_code>fileDownloaded = drive.CreateFile({'id':'1VtBalvWeWV2WtEPf5M19XP8sSnxAG135'})
fileDownloaded.GetContentFile('LS_good_5m_pro.csv')<jupyter_output><empty_output><jupyter_text>Generation of Pandas dataframe:<jupyter_code>import pandas as pd
###
vardf = pd.read_csv('LS_good_5m_pro.csv', delimiter=',')<jupyter_output><empty_output><jupyter_text>Finally, it's possible to create a geodataframe for the imported file (with Geopandas).*Code adapted from: Prof. <NAME>'s notebooks.*<jupyter_code>%pip install geopandas
import geopandas as gpd #necessary import
#Coordinates are already in the imported dataset
vardf = gpd.GeoDataFrame(vardf,
                      geometry=gpd.points_from_xy(np.float64(vardf.x), 
                                                  np.float64(vardf.y)))
vardf.head()<jupyter_output><empty_output><jupyter_text>Part II: Preprocessing The table is checked to delete eventual NaN data:<jupyter_code>vardf.dropna()
vardf.head()
#No NaN data found<jupyter_output><empty_output><jupyter_text>In this section two tables are created: one with predictive variables (independent variables) and one with the landslide absence/presence (dependent variable).For the first table, the 8 predictor variables are selected and put into a geodataframe:DEM -> Digital Elevation Model (m a.s.l)Slope -> Slope (°)Plan Curvature -> Surface curvature perpendicular to the direction of the maximum slopeProfile Curvature -> Surface curvature parallel to the direction of the maximum slopeDistance from Roads -> Distance (m) from main and secondary roadsLand Cover -> Physical material at the surface (like grass, concrete, trees, water, etc)TWI -> Topographic Wetness Index: (ln(a/tan(b))) where a = upslope area draining through a certain point, and tan(b) = slope in radiansGeology classes -> Lithology classes<jupyter_code>DEM = vardf.DEM #Digital Elevation Model
slope = vardf.slope #Slope
planCurv = vardf.planCurv #Plan Curvature
profCurv = vardf.profCurv #Profile Curvature
distRoad = vardf.distRoad #Distance from Roads
landCover = vardf.landCover #Land Cover
TWI = vardf.TWI #Topographic Wetness Index
geology = vardf.geology #Lithology Class
input_data = gpd.GeoDataFrame({'DEM':DEM, 'slope':slope, 'planCurv':planCurv, 'profCurv':profCurv, 'distRoad':distRoad, 'landCover':landCover, 'TWI':TWI, 'geology':geology})
input_data #table with independent variables<jupyter_output><empty_output><jupyter_text>For the second table, the presence (1)/absence (0) for each observation is reported:<jupyter_code>#This data correspond to the GSS column in the original dataset
output_data = vardf.GSS #Presence = 1; Absence = 0 for landslides
output = gpd.GeoDataFrame({'GSS': output_data})
output #table with dependent variable<jupyter_output><empty_output><jupyter_text>Part III: Random Forest Classifier implementation Subdivision of the dataset into train (67%) and test (33%) sets:<jupyter_code>from sklearn.model_selection import train_test_split
###
X_train, X_test, y_train, y_test = train_test_split(input_data, output, test_size = 0.33, random_state = rnd_seed)<jupyter_output><empty_output><jupyter_text>Implementation of Random Forest on the train set:<jupyter_code>from sklearn.ensemble import RandomForestClassifier
###
rf1 = RandomForestClassifier()
rf1.fit(X_train, y_train) #model is fitted<jupyter_output><empty_output><jupyter_text>With the Accuracy Score it is possible to determine the accuracy of the model.With the Mean Squared Error it is possible to find the error of the model (difference between predicted values and true values).With the Confusion Matrix it is possible to see the ratio between true positives, false positives, true negatives, false negatives:<jupyter_code>from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
###
y_pred1 = rf1.predict(X_test) #prediction on the test set
accuracy = accuracy_score(y_test, y_pred1) #accuracy of the prediction
print(f'The accuracy for RandomForestClassifier is {accuracy:.1%}')
mse1 = sklearn.metrics.mean_squared_error(y_test, y_pred1) #Mean Squared Error
print(f'The mean squared error of the RandomForestClassifier is {mse1}')
conf_mat = confusion_matrix(y_test, y_pred1) #confusion matrix of the prediction
print(f'The confusion matrix for RandomForestClassifier is \n {conf_mat}')<jupyter_output>The accuracy for RandomForestClassifier is 83.2%
The mean squared error of the RandomForest Regression is 0.1675423234092236
The confusion matrix for RandomForestClassifier is 
 [[680 164]
 [123 746]]<jupyter_text>Now, an hyperparameters research is conducted to increase the accuracy of the model:*Source: towardsdatascience.com.*<jupyter_code>from sklearn.model_selection import GridSearchCV
###
# Number of trees in random forest
n_estimators = [int(x) for x in np.linspace(start = 50, stop = 500, num = 10)]
# Number of features to consider at every split
max_features = [2, 3, 4]
# Creation of the parameter grid based on the results of random search
param_grid = {'n_estimators': n_estimators,
              'max_features': max_features}
# Random grid search for best hyperparameters: random search
# of parameters, using 3 fold cross validation
gscv = GridSearchCV(estimator = rf1, param_grid = param_grid, cv = 3, verbose=2, n_jobs = -1)
gscv.fit(X_train, y_train) # Fit the grid search to the data<jupyter_output><empty_output><jupyter_text>The function "best_params_" gives back the best hyperparameters of the grid:<jupyter_code>best_grid = gscv.best_params_
best_grid<jupyter_output><empty_output><jupyter_text>Here the Random Forest Classifier is implemented on the train set with the parameters found above:<jupyter_code>rf2 = RandomForestClassifier(max_features = 2, n_estimators = 300)
rf2.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>The Accuracy Score, the Mean Squared Error and the Confusion Matrix are determined like done before:<jupyter_code>y_pred2 = rf2.predict(X_test) #prediction on the test set
accuracy = accuracy_score(y_test, y_pred2) #accuracy of the prediction
print(f'The accuracy for RandomForestClassifier is {accuracy:.1%}')
mse2 = sklearn.metrics.mean_squared_error(y_test, y_pred2) #Mean Squared Error
print(f'The mean squared error of the RandomForestClassifier is {mse2}')
conf_mat = confusion_matrix(y_test, y_pred2) #confusion matrix of the prediction
print(f'The confusion matrix for RandomForestClassifier is \n {conf_mat}')<jupyter_output>The accuracy for RandomForestClassifier is 83.8%
The mean squared error of the RandomForest Regression is 0.16228838295388207
The confusion matrix for RandomForestClassifier is 
 [[684 160]
 [118 751]]<jupyter_text>Part IV: Feature importance for Random Forest First method for Feature importance - Mean Decrease in Impurity.*Source: scikit-learn; datavizpyr.com.*<jupyter_code>importances = rf2.feature_importances_ #intrinsic function of RF
std = np.std([tree.feature_importances_ for tree in rf2.estimators_], axis=0)
feature_names = [input_data.columns] #set labels as the same of the original dataset
#Pandas dataframe
forest_importances_1 = pd.Series(importances, index=feature_names)
#Make bars to be in ascendant order
sorted1 = forest_importances_1.sort_values('index', ascending=False)
#Make bar plot
fig, ax = plt.subplots()
sorted1.plot.bar(yerr=std, ax=ax)
ax.set_title("Feature importances using MDI")
ax.set_ylabel("Mean decrease in impurity")
fig.tight_layout()
#plt.savefig("bar_plot_MDI.png") #saving image<jupyter_output>/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:3: FutureWarning: In a future version of pandas all arguments of Series.sort_values will be keyword-only
  This is separate from the ipykernel package so we can avoid doing imports until<jupyter_text>Variables with a mean >= 0.15 are listed here:<jupyter_code>var_import_1 = forest_importances_1[forest_importances_1 >= 0.15]
var_import_1<jupyter_output><empty_output><jupyter_text>Second method for Feature importance - Feature Permutation. *Source: scikit-learn; datavizpyr.com.*<jupyter_code>from sklearn.inspection import permutation_importance
###
result = permutation_importance(
    rf2, X_test, y_test, n_repeats=10, random_state=42, n_jobs=2
) #Calculation for RF
#Pandas dataframe
forest_importances_2 = pd.Series(result.importances_mean, index=feature_names)
#Make bars to be in ascendant order
sorted2 = forest_importances_2.sort_values('index', ascending=False)
#Make bar plot
fig, ax = plt.subplots()
sorted2.plot.bar(yerr=result.importances_std, ax=ax, color='orange')
ax.set_title("Feature importances using Feature Permutation")
ax.set_ylabel("Mean accuracy decrease")
fig.tight_layout()
plt.show()
#plt.savefig("bar_plot_FP.png") # saving image<jupyter_output>/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:2: FutureWarning: In a future version of pandas all arguments of Series.sort_values will be keyword-only<jupyter_text>Variables with a mean >= 0.05 are listed here:<jupyter_code>var_import_2 = forest_importances_2[forest_importances_2 >= 0.05]
var_import_2<jupyter_output><empty_output><jupyter_text>Part V: Feature importance for Gradient Boosting Classifier.Finally, it is insightful to complement the analysis with a feature importance applied to another algorithm. In the next cell, Gradient Boosting Classifier is computed: *Source: scikit-learn.*<jupyter_code>from sklearn.ensemble import GradientBoostingClassifier
###
gbc = GradientBoostingClassifier(n_estimators=300, learning_rate=1.0, max_depth=1, random_state=rnd_seed).fit(X_train, y_train)
gbc.score(X_test, y_test) #Evaluation of GB<jupyter_output>/usr/local/lib/python3.7/dist-packages/sklearn/ensemble/_gb.py:494: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)<jupyter_text>The feature importance is determined with Permutation Importance method, the most appropriate for GradientBoostClassifier.  *Source: scikit-learn.*<jupyter_code>boost = sklearn.inspection.permutation_importance(gbc, X_test, y_test, scoring=None, n_repeats=5, n_jobs=None, random_state=None, sample_weight=None, max_samples=1.0)
#Pandas dataframe
boost_importances = pd.Series(boost.importances_mean, index=feature_names)
#Make bars to be in ascendant order
sorted3 = boost_importances.sort_values('index', ascending=False)
#Make bar plot
fig, ax = plt.subplots()
sorted3.plot.bar(yerr=boost.importances_std, ax=ax, color='g')
ax.set_title("Feature importances using Permutation Importance")
ax.set_ylabel("Mean accuracy decrease")
fig.tight_layout()
plt.show()
#plt.savefig("bar_plot_GBC.png") # saving image<jupyter_output>/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:5: FutureWarning: In a future version of pandas all arguments of Series.sort_values will be keyword-only
  """<jupyter_text>Variables with a mean >= 0.03 are listed here:<jupyter_code>var_import_3 = boost_importances[boost_importances >= 0.03]
var_import_3<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09051295772494458, 'fraction_numerical': 0.01720051983793288, 'mean_word_length': 6.699823425544438, 'pattern_counts': {'":': 0, '<': 90, '<?xml version=': 0, '>': 107, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '156014', 'license': 'MIT', 'n_tokens_mistral': 4423, 'n_tokens_neox': 4191, 'n_words': 1210, 'path': 'Final_Project.ipynb', 'repo_name': 'AlePol27/2022_ML_Earth_Env_Sci', 'size': '13103'}