<jupyter_start><jupyter_text>Doyle Reaction Prediction Pt. 4: Additive Importance in Random ForestsThe one-hot experiments provided in the previous notebook gave better results on out-of-sample tests for new additives than expected, indicating that additives may be limited in their importance for predicting yield. This appears to be inconsistent with the paper's analysis on random forest feature importances, where additive features are most significant.For categorical variables, it has been established that traditional random forests are biased toward categories with more 'levels' (e.g. if a categorical variable has 50 potential values, and another only has 10, then the system will be biased towards the former). Sensitivity in levels and scaling could, in part, be responsible for inflated additive importances.Some useful background and analysis can be found in:    1) https://www.tandfonline.com/doi/abs/10.1198/016214501753168271    2) https://bmcbioinformatics.biomedcentral.com/track/pdf/10.1186/1471-2105-8-25    Reading this paper, each substrate/ligand/additive/base we try are represented by discrete numerical values. One question is whether there is a natural bias towards any of the parameters. Our hypothesis is:    Additive importances are biased in a traditional random forests implementation.Design:        To decorrelate input variables to the dependent variable (yield), we will shuffle the yields. We will fit a random forest regressor to this data, and extract out the feature importances. Unbiased feature importances should be evenly distributed across each parameter, and we would not expect enrichment for any specific set of features (e.g. no enrichment for additive features).<jupyter_code># Import our libraries
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patches as mpatches
import pandas as pd

# Import relevant scikit-learn modules

from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.svm import LinearSVR
from sklearn.metrics import mean_squared_error, r2_score

# Load RPy2 Libraries for R dataframes
from rpy2.robjects import r, pandas2ri
import rpy2.robjects as robjects

import utils
# set a seed for all the randoms
np.random.seed(42)
# Define path and csv files
DATA_DIR = 'data/'
YIELDS_DF = 'yields.csv'
INPUTS = 'output_table.csv'      # Unscaled  data from github.com/doylelab/rxnpredict
TRAIN_DF = 'training_scaled.rds' # Scaled training data from github.com/doylelab/rxnpredict
TEST_DF = 'test_scaled.rds'      # Scaled test set data from github.com/doylelab/rxnpredict
# Load yields data
yields = pd.read_csv(DATA_DIR + YIELDS_DF)

# As in the paper, remove controls and additive 7
yields = np.array(yields['x'])
yields = yields.reshape(3 * 32, 48)
yields = np.delete(yields, [0, 4, 8, 12, 19, 23, 27, 31], axis=0)
yields = np.delete(yields, [15, 31, 47], axis=1)
yields = yields.flatten()

# it turns out the yield data contains some NaNs
# we set these to zeroes
y_set = yields

for idx, value in enumerate(y_set):
    if np.isfinite(value) == False:
        print(idx, value)
        y_set[idx] = np.nanmean(y_set)
# Activate RPy2
pandas2ri.activate()
readRDS = robjects.r['readRDS']

# load R dataframes
train_rds = readRDS(DATA_DIR + TRAIN_DF)
test_rds  = readRDS(DATA_DIR + TEST_DF)

# convert these to pandas dataframes
train_df = pandas2ri.ri2py(train_rds)
test_df  = pandas2ri.ri2py(test_rds)

train_indices = train_df.sort_index().index - 1
test_indices = test_df.sort_index().index - 1
# Load chemical features and standard scales
inputs = pd.read_csv(DATA_DIR + INPUTS)
features = inputs.columns

scaler = StandardScaler()
scaler.fit(inputs)
inputs = scaler.transform(inputs)
# shuffle the yields
y_set = yields
np.random.shuffle(y_set)
X_train = inputs[train_indices]
y_train = y_set[train_indices]

X_test = inputs[test_indices]
y_test = y_set[test_indices]
models = [LinearRegression(), # use default
          KNeighborsRegressor(n_neighbors=7), # use k = 7 as in paper
          MLPRegressor(hidden_layer_sizes=(5),
                       activation='logistic',
                       solver='sgd',
                       max_iter=10000,
                       random_state=42),
          MLPRegressor(hidden_layer_sizes=(100),
                       activation='logistic',
                       solver='sgd',
                       max_iter=10000,
                       random_state=42),
          RandomForestRegressor(n_estimators=500,
                                random_state=42)
         ]
preds, r2_values, rmse_values = utils.fit_models(X_train,
                                                           X_test,
                                                           y_train,
                                                           y_test,
                                                           models)

utils.plot_models(preds,
                  r2_values,
                  rmse_values,
                  y_test)<jupyter_output>LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
KNeighborsRegressor(algorithm='auto', leaf_size=30, metric='minkowski',
          metric_params=None, n_jobs=1, n_neighbors=7, p=2,
          weights='uniform')
MLPRegressor(activation='logistic', alpha=0.0001, batch_size='auto',
       beta_1=0.9, beta_2=0.999, early_stopping=False, epsilon=1e-08,
       hidden_layer_sizes=5, learning_rate='constant',
       learning_rate_init=0.001, max_iter=10000, momentum=0.9,
       nesterovs_momentum=True, power_t=0.5, random_state=42, shuffle=True,
       solver='sgd', tol=0.0001, validation_fraction=0.1, verbose=False,
       warm_start=False)
MLPRegressor(activation='logistic', alpha=0.0001, batch_size='auto',
       beta_1=0.9, beta_2=0.999, early_stopping=False, epsilon=1e-08,
       hidden_layer_sizes=100, learning_rate='constant',
       learning_rate_init=0.001, max_iter=10000, momentum=0.9,
       nesterovs_momentum=True, power_t=0.5, random_state=42, shuffle=Tr[...]<jupyter_text>As expected, shuffling the yields destroys any ability to predict. We will conduct 100 trials of this random experiment, and take the average feature importances.<jupyter_code>from sklearn.model_selection import train_test_split
results = []
for i in range(100): # shuffle 100 times
    np.random.shuffle(y_set)
    X_train, X_test, y_train, y_test = train_test_split(inputs, y_set, test_size=0.30)
    model = RandomForestRegressor(n_estimators=500,
                                  n_jobs=-1)
    model.fit(X_train, y_train)
    new_df = pd.DataFrame()
    new_df['importances'] = model.feature_importances_
    new_df['features'] = features
    new_df = new_df.sort_values(by=['importances'], ascending=False)
    results.append(new_df)
fracs = []
for result in results:
    top_ten = result['features'][:10]
    no_add_feats = top_ten.str.startswith('additive').sum()
    fraction = no_add_feats / 10
    fracs.append(fraction)

fracs = np.array(fracs)
mean_fracs = fracs.mean()<jupyter_output><empty_output><jupyter_text>Out of 100 random shuffles and samplings, on average 8.5 out of 10 of the most important features for random forests are additive based.<jupyter_code>halides = []
additives = []
bases = []
ligands = []
for result in results:
    halide = result[result['features'].str.startswith('aryl')]['importances']
    additive = result[result['features'].str.startswith('additive')]['importances']
    base = result[result['features'].str.startswith('base')]['importances']
    ligand = result[result['features'].str.startswith('ligand')]['importances']
    halides.append(halide)
    additives.append(additive)
    bases.append(base)
    ligands.append(ligand)
halides = pd.concat(halides)
additives = pd.concat(additives)
bases = pd.concat(bases)
ligands = pd.concat(ligands)
halide_stats = halides.mean(), halides.std()
additive_stats = additives.mean(), additives.std()
base_stats = bases.mean(), bases.std()
ligand_stats = ligands.mean(), ligands.std()

print(halide_stats,
      additive_stats,
      base_stats,
      ligand_stats)
plt.figure(figsize=(5,4))
plt.tight_layout()
plt.title('RF Feature Importances on Randomly Shuffled Data')
plt.xlabel('relative importance')
labels = ['ligand features',
          'halide features',
          'base features',
          'additive features']
plt.boxplot([ligands, halides, bases, additives],
            0, '', 0, labels=labels)
plt.axvline(x=1/120, ls="--", c='r')
plt.text(0.009, 0.65, 'expected', fontsize=11, color='r', style='italic')
plt.savefig('plots/boxplot_featureimportances.png', dpi = 300, bbox_inches='tight')
plt.show()<jupyter_output>/home/kangway/.conda/envs/doyle_rxn/lib/python3.6/site-packages/numpy/core/fromnumeric.py:52: FutureWarning: reshape is deprecated and will raise in a subsequent release. Please use .values.reshape(...) instead
  return getattr(obj, method)(*args, **kwds)<jupyter_text>additive features are consistently elevated compared to base and halide features.<jupyter_code>sample_importances = results[3][:10].sort_values(by=['importances'], ascending=True)
plt.title('RF Feature Importance on Randomly Shuffled Data')
plt.barh(range(10), sample_importances['importances'], color='gray', align='center')
plt.yticks(range(10), sample_importances['features'], fontsize=12)
plt.xlabel('Relative Importance')
plt.savefig('plots/topten_feature_importances.png', dpi = 300,bbox_inches='tight')
plt.show()
fracs = []
for result in results:
    top_ten = result['features'][:10]
    no_add_feats = top_ten.str.startswith('additive').sum()
    fraction = no_add_feats / 10
    fracs.append(fraction)

fracs = np.array(fracs)
mean_fracs = fracs.mean()
mean_fracs<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1021461834379429, 'fraction_numerical': 0.026624822838631303, 'mean_word_length': 3.8785185185185185, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '247798', 'license': 'MIT', 'n_tokens_mistral': 3335, 'n_tokens_neox': 3040, 'n_words': 794, 'path': 'Doyle_Science_2018/4. Additive Importance in Random Forests.ipynb', 'repo_name': 'keiserlab/comments', 'size': '9878'}