<jupyter_start><jupyter_text>Random Forests<jupyter_code>import init
from common import constants as cn
from common.data_provider import DataProvider
from common import transform_data
from common_python.statistics import empirical_distribution_generator, density, multiple_density, util_statistics
from common_python.plots import util_plots
from plots import util_plots as uplots


from sklearn.datasets import load_iris
from sklearn.model_selection import cross_val_score
from sklearn.tree import DecisionTreeClassifier
from sklearn import tree
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import collections
import seaborn as sns

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Generate Classification Data1. T0 is eliminated1. Normoxia is combined with Resuscitation (because of small number of Normoxia)1. Perfectly correlated genes are combined<jupyter_code>provider = DataProvider()
provider.do()
len(provider.df_normalized)
# Geenrate the feature data
def makeFeatureDF(provider):
    df_Xdata = transform_data.aggregateGenes(provider=provider)
    df_Xdata = df_Xdata.T
    df_Xdata = df_Xdata.drop(index="T0")
    data_features = df_Xdata.columns
    df_Xdata.columns = range(len(data_features))
    return df_Xdata

def makeFeatureColumns(provider):
    df_X = transform_data.aggregateGenes(provider=provider)
    df_X = df_X.T
    return df_X.columns.tolist()


df_Xdata = makeFeatureDF(provider)
df_Xdata.head()<jupyter_output><empty_output><jupyter_text>There are only two Noroxia states, one of which is a reference state. So, for the class data, we combineNoroxia with Resusitation.<jupyter_code># Prepare the state data
def makeClassSer(provider):
    ser_state = provider.df_stage_matrix[cn.STAGE_NAME]
    ser_state = ser_state.drop(index="T0")
    ser_state = ser_state.copy()
    ser_state[ser_state == 'Normoxia'] = 'Resuscitation'
    # Create converter from state name to numeric index
    states = ser_state.unique()
    state_dict = {k: v for v, k in enumerate(states)}
    ser = ser_state.apply(lambda k: state_dict[k] )
    return ser
ser_Ydata = makeClassSer(provider)<jupyter_output><empty_output><jupyter_text>Random Forests for Expression Data<jupyter_code>from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
import collections

DecisionTreeResult = collections.namedtuple("DecisionTreeResult",
                                           "top_indices importances clf")

def getTopIndices(num_genes, genes=None, **kwargs):
    """
    :param int num_genes: number of genes in plot
    :param dict kwargs: parameters for RandomForestClassifier
    :return list-int, list-float: indices of the top genes, importances
    """
    clf = RandomForestClassifier(**kwargs)
    df_X = makeFeatureDF(provider)
    if genes is None:
        genes = df_X.columns
    df_X = df_X[genes]
    clf.fit(makeFeatureDF(provider), makeClassSer(provider))
    importances = list(clf.feature_importances_)
    indices = list(np.argsort(importances))
    indices.reverse()
    top_indices = [indices[n] for n in range(num_genes)]
    return DecisionTreeResult(top_indices, importances, clf)

def setRFPlotAttributes(ax):
    _ = plt.ylabel("Importance")
    _ = plt.ylim([0,.1])

def plotRF(num_genes=25, **kwargs):
    """
    :param int num_genes: number of genes in plot
    :param dict kwargs: parameters for RandomForestClassifier
    :return matplotib.Axes:
    """
    fig, ax = plt.subplots()
    result = getTopIndices(num_genes, **kwargs)
    columns = makeFeatureColumns(provider)
    pairs = [(columns[n], result.importances[n]) for n in result.top_indices]
    features, values = zip(*pairs)
    ax.bar(features, values)
    setRFPlotAttributes(ax)
    _ = ax.set_xticklabels(features, rotation='vertical')
    return ax
_ = plotRF(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=False)
plotRF(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=False)
_ = plt.title("Without Bootstrap")
plotRF(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=True)
_ = plt.title("With Bootstrap")
plotRF(num_genes=50, n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=False)
_ = plt.title("Without Bootstrap")
plotRF(num_genes=50, n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=True)
_ = plt.title("With Bootstrap")
def plotImportanceCDF(num_genes=200, **kwargs):
    result = getTopIndices(num_genes=num_genes, **kwargs)
    importances = result.importances
    importances.sort(reverse=True)
    fig, ax = plt.subplots()
    ax.plot(importances)
    setRFPlotAttributes(ax)
    ax.set_xlim([0, num_genes])
    ax.set_xlabel("Gene")
    return ax
ax = plotImportanceCDF(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=False)
ax.set_title("No Boostrap")
ax = plotImportanceCDF(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=True)
_ = ax.set_title("With Boostrap")
# Analysis of overlap between bootstrap and non-bootstrap genes
def getCommon(num_genes, **kwargs):
    """
    Calculates the number of genes in common between bootstrap and non-bootstrap.
    :param int num_genes: how many top genes considered
    :param dict kwargs: parameters to classifier
    :return list-int: list of indexes in common
    """
    result_false = getTopIndices(num_genes, bootstrap=False, **kwargs)
    result_true = getTopIndices(num_genes, bootstrap=True, **kwargs)
    return set(result_false.top_indices).intersection(result_true.top_indices)
getCommon(25, n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2)
kwargs = dict(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2)
values = [len(getCommon(n, **kwargs)) / n for n in [5, 10, 15, 20, 25, 30]]
values
# Compare the accuracy of trees with all features with those that only have top genes
# Look at clf.oob_score__num_
kwargs = dict(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=True, oob_score=True)
values = getCommon(25, n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2)
num_genes = len(values)
result = getTopIndices(num_genes, genes=list(values), **kwargs)
result.clf.oob_score_
num_genes = 25
result = getTopIndices(num_genes, **kwargs)
result.clf.oob_score_<jupyter_output><empty_output><jupyter_text>Observations1. Significant overlap between bootstrap and non-bootstrap features1. Try DT with top N features.<jupyter_code>clf = RandomForestClassifier(n_estimators=1000, max_features=1500, random_state=0,
                            min_samples_split=2, bootstrap=True)
clf.fit(makeFeatureDF(provider), makeClassSer(provider))
tree = clf.estimators_[0]
tree.n_features_<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0908462389380531, 'fraction_numerical': 0.021847345132743362, 'mean_word_length': 4.365727002967359, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '447495', 'license': 'MIT', 'n_tokens_mistral': 2408, 'n_tokens_neox': 2252, 'n_words': 523, 'path': 'notebooks/Random Forests.ipynb', 'repo_name': 'uwescience/expression_state', 'size': '7232'}