{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nAnalyse univariable feature importances: correlations (Pearson, Spearman) , prediction by boosting on one variable, etc.\n\n### Versions:\n\n39 Added MutualInfoSklearn CD36\n\n34,35, 36, 37,38 CD36,CD44, CD194,  CD146,  CD146(misprint corrected)  - list_n_estimators = [1,2,3,100] # Used for boostings for trial0,trial1, etc...  Both CatBoost and LGB \n\n'CD146 '#   'CD194' - poorly predicatable (at least by Lasso 0.01)\n\n32,33 CD44,CD36 with 3, 5, 10, NO(100) trials;  ONLY LGB, not catboost \n\n31 CD44 with 3, 5, 10, NO(100) trials \n\n30 CD44 with 3, 5, 10, and 100 trials for CatBoost and LGB, not 10 and 100 as it was before , 7.82 hours\n\n25,26 CD44, 27, 28, 29 CD36, CD88, CD32 -- 3.10 hours\n\n19,20, 21,22,23,24  CD36,  CD88, CD32, CD71, CD48, CD62L - crashes due to bugs, need to correct \n \n17,18 CD44 Added Catboost, changed LGB to 10 and 100 trials with train/test split, not CV\n\n10, 11,12, 13,14,15,16 Added analysis of rank differences, CD36, CD44, CD88, CD32, CD71, CD48, CD62L\n\n8,9 Add PPS score . CD44, CD36\n\nhttps://towardsdatascience.com/rip-correlation-introducing-the-predictive-power-score-3d90808b9598\n\n4,5,6,7 CD32, CD71, CD48,CD62L\n\n3 - CD88 Pearson * 2 + LGB * 2\n\n2 - CD44 Pearson * 2 + LGB * 2\n\n1 - CD36 Pearson * 2 + LGB * 2 \n","metadata":{}},{"cell_type":"markdown","source":"# Key params","metadata":{}},{"cell_type":"code","source":"target_name = 'CD36'#   'CD194' # 'CD44'# 'CD62L'\n\ndown_sample_percent = 50\n\nn_trials_mutual_info_sklearn = 1#  4\n\nlist_n_estimators = [1,2,3,100] # Used for boostings for trial0,trial1, etc... \nn_trials_lgb = 0# 4\nn_trials_catBoost = 0#  4\n\n\nimport time\nt0start = time.time() ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pd.set_option('display.max_rows', 500)\n# pd.set_option('display.max_columns', 500)\n# pd.set_option('display.width', 1000)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.linear_model import LassoCV\nfrom sklearn.linear_model import Lasso\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.linear_model import Ridge\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold \nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import mean_squared_error","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data","metadata":{}},{"cell_type":"code","source":"%%time\nfilename_rna_data = '/kaggle/input/open-problems-multimodal/train_cite_inputs.h5'\ndf_rna = pd.read_hdf(filename_rna_data)\ndisplay(df_rna) \n\n#%%time\ndf_y = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndisplay(df_y)\n\nfn = '/kaggle/input/open-problems-multimodal/metadata.csv'\ndf_meta = pd.read_csv(fn, index_col = 0 )\ndf_meta\n# Cut only train cite-seq part: \nd = pd.DataFrame(index = df_y.index)\nprint(d.shape)\ndf_meta = d.join(df_meta, how = 'left')\ndf_meta","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( dict( df_meta.value_counts('cell_type') ))\n\n# for t in df_rna.columns:\n#     if 'CD44' in t:\n#         print(t)\n#     if 'ENSG00000197405' in t: # ENSG00000197405_C5AR1 CD88\n#         print(t)\n        \n        \n# 'ENSG00000135218_CD36', 'ENSG00000026508_CD44', 'ENSG00000255443_CD44-AS1'\n       ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# y = df_y[target_name].values\n# print(y.shape, type(y) )\n\nX = ( df_rna.values )\nprint(X.shape)\n\ny = df_y[target_name]\nprint(y.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataframe to store results and target","metadata":{}},{"cell_type":"code","source":"df_importances = pd.DataFrame(index = df_rna.columns )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mutual information sklearn ","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.feature_selection import f_regression, mutual_info_regression\nimport time\n\nt0 = time.time()\n\ny = df_y[target_name]\nN = int(down_sample_percent*len(y) / 100 )\nfor trial in range(n_trials_mutual_info_sklearn):\n    l = []\n    for i,col in enumerate( df_rna.columns):\n\n        p = np.random.permutation(len(y))\n\n        x_loc = df_rna[col].values[p][:N].reshape(-1,1)\n        mi = mutual_info_regression(x_loc , y.values[p][:N] )\n        s = mi[0]\n        df_importances.loc[col, target_name +' MutualInfSklearn DownSample'+str(down_sample_percent) + ' Trial'+str(trial) ]  = s\n\n        l.append(s)\n        if (i%5000==0) or (i<5):\n            I = np.argmax(l)\n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n        \n    \ndf_importances","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Mutual information sklearn ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mutual information skfeature ","metadata":{}},{"cell_type":"code","source":"pip install git+https://github.com/jundongl/scikit-feature","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skfeature.utility.entropy_estimators import midd\nimport numpy as np\nimport time","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\n\ny = df_y[target_name].round(1)\nN = int(down_sample_percent*len(y) / 100 )\nfor trial in range(n_trials_mutual_info_sklearn):\n    l = []\n    for i,col in enumerate( df_rna.columns):\n\n        p = np.random.permutation(len(y))\n\n        x_loc = df_rna[col].values[p][:N].round(1)\n        mi = midd(x_loc , y.values[p][:N] )\n        s = mi\n        df_importances.loc[col, target_name +' MutualInfSkfeature DownSample'+str(down_sample_percent) + ' Trial'+str(trial) ]  = s\n\n        l.append(s)\n        if (i%5000==0) or (i<5):\n            I = np.argmax(l)\n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n        \n    \ndf_importances","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlations","metadata":{}},{"cell_type":"code","source":"%%time\nfor i,col in enumerate( df_rna.columns):\n    s = np.corrcoef(y,df_rna[col])[0,1]\n    df_importances.loc[col, target_name +' Corr Pearson'  ]  = s\ndf_importances","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ny = df_y[target_name]\np = np.random.permutation(len(y))\nN = int(down_sample_percent*len(y) / 100 )\nfor i,col in enumerate( df_rna.columns):\n    s = np.corrcoef(y.values[p][:N],df_rna[col].values[p][:N])[0,1]\n    df_importances.loc[col, target_name +' Corr Pearson DownSample'+str(down_sample_percent)  ]  = s\ndf_importances","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nn_x_subplots = 3; c = 0 \nfor col in df_importances.columns:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0    \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n    plt.hist(df_importances[col], bins = 100)\n    plt.title(col)\nplt.show()\n    \ndisplay(df_importances.describe() )    \ndf_importances.sort_values(col, ascending = False, key = abs ).head(20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n            \n        if 'LGB' not in col1:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False, key = abs).abs()).rank()\n        else:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False)).rank()\n            \n        v1 = df_importances_ranks[col1] - df_importances_ranks[col2]\n        \n        fig = plt.figure(figsize = (20,6))\n        plt.title('Ranks differences: ' + col1+' ' + col2, fontsize = 20 )\n        plt.plot(v1.values ) \n        plt.grid()\n        plt.show()\n        \n        fig = plt.figure(figsize = (20,6))\n        plt.title('Ranks differences: ' + col1+' ' + col2, fontsize = 20 )\n        plt.plot(v1.values[:500] ) \n        plt.grid()\n        plt.show()        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_stat = pd.DataFrame()\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n        l = []\n        if 'LGB' not in col1:\n            v1 = df_importances[col1].sort_values(ascending = False, key = abs )\n        else:\n            v1 = df_importances[col1].sort_values(ascending = False)\n        if 'LGB' not in col2:\n            v2 = df_importances[col2].sort_values(ascending = False, key = abs )\n        else:\n            v2 = df_importances[col2].sort_values(ascending = False)\n            \n        for k in [100,500,1000]:# range(len(v1)): # [10,25,50, 75, 100, 1000, 2000, 4000]:\n            s = set(v1.index[:k] ) & set(v2.index[:k] ) \n            l.append(len(s) )\n            # print(k, len(s))   \n            \n        print(col1,col2, len(l))\n        df_stat.loc[col1+' VS '+col2, 'Top100 Intersection'] = l[0]\n        df_stat.loc[col1+' VS '+col2, 'Top500 Intersection'] = l[1]\n        df_stat.loc[col1+' VS '+col2, 'Top1000 Intersection'] = l[2]\n        \ndisplay(df_stat )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n            \n        if 'LGB' not in col1:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False, key = abs).abs()).rank()\n        else:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False)).rank()\n            \n        v1 = df_importances_ranks[col1] - df_importances_ranks[col2]\n        \n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top100'] = np.max( np.abs( v1.values[:100] ) )\n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top500'] = np.max( np.abs( v1.values[:500] ) )\n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top1000'] = np.max( np.abs( v1.values[:1000] ) ) \n\ndisplay(df_stat)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n        print(col1,col2)\n        l = []\n        if 'LGB' not in col1:\n            v1 = df_importances[col1].sort_values(ascending = False, key = abs )\n        else:\n            v1 = df_importances[col1].sort_values(ascending = False)\n        if 'LGB' not in col2:\n            v2 = df_importances[col2].sort_values(ascending = False, key = abs )\n        else:\n            v2 = df_importances[col2].sort_values(ascending = False)\n            \n        for k in range(len(v1)): # [10,25,50, 75, 100, 1000, 2000, 4000]:\n            s = set(v1.index[:k] ) & set(v2.index[:k] ) \n            l.append(len(s) )\n            # print(k, len(s))       \n\n        fig = plt.figure(figsize = (20,6))\n        plt.title(col1+' ' + col2, fontsize = 20 )\n        plt.plot(l ,label = 'original data - counts intersection' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x,y,2)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:50],y[:50],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-50 points' )\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()            \n        \n        \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l ,label = 'original data - counts intersection' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x,y,2)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:50],y[:50],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation based on 0-50 points' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:150],y[:150],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-150 points' )\n\n        plt.xlim([0,500])\n        plt.ylim([0,250])\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()        \n        \n        ##################################################################33\n        l2 = []\n        for k in range(10,len(l)):\n            l2.append(l[k]/k)\n        \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l2 ,label = 'intersection topK / K' )\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()                \n            \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l2[:500] ,label = 'intersection topK / K' )\n\n        plt.xlim([0,500])\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()                ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CatBoost","metadata":{}},{"cell_type":"code","source":"%%time\nimport catboost as cb\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nimport time\nt0 = time.time()\n\nif 0: # See updated (faster) version in the next cell \n\n    random_state4lgb  = 10 \n    model = cb.CatBoostRegressor(verbose = 0 ,  n_estimators = 100 )\n    #model = lgbm.LGBMRegressor(random_state = random_state4lgb) # Ridge(alpha = alpha_selected )\n    #kf = KFold(n_splits=n_splits_for_cross_valdition,  shuffle=True, random_state= random_state_cross_validation )\n    kf = KFold(n_splits=2,  shuffle=True, random_state= random_state4lgb )\n\n\n    y = df_y[target_name]\n    for i,col in enumerate( df_rna.columns[:10]):\n\n        #col = 'ENSG00000135218_CD36'\n        X = df_rna[[col]]\n        y_pred = cross_val_predict(model, X, y, cv=kf)\n        s = np.corrcoef(y,y_pred)[0,1]\n        df_importances.loc[col, target_name +' Corr CatBoost RS'+str(random_state4lgb)  ]  = s\n\n        if (i%1000==0): \n            # print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n\n    display( df_importances )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport time\nimport catboost as cb\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nimport time\nt0 = time.time()\n\n\nfor trial in range(n_trials_catBoost):\n\n    if trial < len(list_n_estimators):\n        n_estimators = list_n_estimators[trial]\n    else:\n        n_estimators = list_n_estimators[-1]\n    model = cb.CatBoostRegressor(verbose = 0, n_estimators = n_estimators )\n    print(n_estimators, model)\n    t0 = time.time()\n\n    p = np.random.permutation(len(df_y))\n    N = int(down_sample_percent*len(df_y) / 100 ) # +str(down_sample_percent) \n    y = df_y[target_name]\n    y_train = y.values[p][:N]\n    y_test =  y.values[p][N:]\n    l = []\n    for i,col in enumerate( df_rna.columns):\n    #     s = np.corrcoef(y.values[p][:N],df_rna[col].values[p][:N])[0,1]\n    #     df_importances.loc[col, target_name +' Corr Pearson DownSample'+str(down_sample_percent)  ]  = s\n        X_train = pd.Series(df_rna[col].values[p][:N]).to_frame()\n        X_test = pd.Series(df_rna[col].values[p][N:]).to_frame()\n        \n        s = 0\n        try:\n            if X_train[0].nunique() > 1: \n                model.fit(X_train,y_train)\n                y_pred = model.predict(X_test)# ,y_train)\n                s = np.corrcoef(y_test,y_pred)[0,1]\n            else:\n                s = 0\n        except ValueError:\n            s = 0\n            \n        #N = int(down_sample_percent*len(y) / 100 )            \n        df_importances.loc[col, target_name +' Corr CatBoost DownSample'+str(down_sample_percent) +' Trial'+str(trial)+' n_estimators'+str(n_estimators)  ]  = s\n        l.append(s)\n        \n        if (i%1000==0): \n            I = np.argmax(l);\n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n    \ndf_importances.sort_values(df_importances.columns[-1],  ascending = False ).head(30)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df_importances.describe() )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.corr()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.columns[-1]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.sort_values(df_importances.columns[-1],  ascending = False ).head(30)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PPS score importances \n\nhttps://towardsdatascience.com/rip-correlation-introducing-the-predictive-power-score-3d90808b9598\n","metadata":{}},{"cell_type":"code","source":"!pip install ppscore","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import ppscore as pps","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport ppscore as pps\nimport time\nt0 = time.time()\n\ndf_tmp = pd.DataFrame()\ndf_tmp[ \"target_column\" ] = df_y[target_name]\ncol = df_rna.columns[0]\nl = []\nfor i,col in enumerate( df_rna.columns):\n    df_tmp[\"feature_column\"] = df_rna[col].values\n    s = pps.score(df_tmp, \"feature_column\", \"target_column\")\n    l.append(s[ 'ppscore'])\n    #print(col)\n    #print(df_tmp[\"feature_column\"].sum(), df_rna[col].sum() )\n    #print(s)\n    if (i%5000==0):\n        I = np.argmax(l); \n        print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n              '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )\n        #print()\nprint(l[:100])    \ndf_importances[target_name +' PPS'  ]  = l\ndisplay(df_importances.describe() )\ndisplay(df_importances.sort_values(target_name +' PPS', ascending = False ).head(30) )\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PPS on subsample","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\n\np = np.random.permutation(len(df_y))\nN = int(down_sample_percent*len(df_y) / 100 ) # +str(down_sample_percent) \n\ndf_tmp = pd.DataFrame()\ndf_tmp[ \"target_column\" ] = df_y[target_name].values[p][:N]\ncol = df_rna.columns[0]\nl = []\nfor i,col in enumerate( df_rna.columns):\n    df_tmp[\"feature_column\"] = df_rna[col].values[p][:N]\n    s = pps.score(df_tmp, \"feature_column\", \"target_column\")\n    l.append(s[ 'ppscore'])\n    #print(col)\n    #print(df_tmp[\"feature_column\"].sum(), df_rna[col].sum() )\n    #print(s)\n    if (i%1000==0):\n        I = np.argmax(l); \n        print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n              '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )\n        #print()\nprint(l[:100])    \ndf_importances[target_name +' PPS DownSample'+str(down_sample_percent)  ]  = l\ndisplay(df_importances.describe() )\ndisplay(df_importances.sort_values(target_name +' PPS DownSample'+str(down_sample_percent), ascending = False ).head(30) )\n\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# LGB importances","metadata":{}},{"cell_type":"code","source":"%%time\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nimport time\nt0 = time.time()\n\nif 0: # In the next cell we use more simple version with train / test\n    \n    random_state4lgb  = 0 \n    model = lgbm.LGBMRegressor(random_state = random_state4lgb) # Ridge(alpha = alpha_selected )\n    #kf = KFold(n_splits=n_splits_for_cross_valdition,  shuffle=True, random_state= random_state_cross_validation )\n    kf = KFold(n_splits=2,  shuffle=True, random_state= random_state4lgb )\n\n\n    y = df_y[target_name]\n    for i,col in enumerate( df_rna.columns):\n        #col = 'ENSG00000135218_CD36'\n        X = df_rna[[col]]\n        y_pred = cross_val_predict(model, X, y, cv=kf)\n        s = np.corrcoef(y,y_pred)[0,1]\n        df_importances.loc[col, target_name +' Corr LGB RS'+str(random_state4lgb)  ]  = s\n\n        if (i%1000==0): \n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n\n    df_importances\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport time\nimport catboost as cb\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nimport time\nt0 = time.time()\n\n\n\nfor trial in range(n_trials_lgb):\n\n    if trial < len(list_n_estimators):\n        n_estimators = list_n_estimators[trial]\n    else:\n        n_estimators = list_n_estimators[-1]\n        \n    #model = cb.CatBoostRegressor(verbose = 0, n_estimators = n_estimators )\n    random_state4lgb = trial\n    model = lgbm.LGBMRegressor(random_state = random_state4lgb, n_estimators = n_estimators) # Ridge(alpha = alpha_selected )\n        \n    t0 = time.time()\n\n    p = np.random.permutation(len(df_y))\n    N = int(down_sample_percent*len(df_y) / 100 ) # +str(down_sample_percent) \n    y = df_y[target_name]\n    y_train = y.values[p][:N]\n    y_test =  y.values[p][N:]\n    l = []\n    for i,col in enumerate( df_rna.columns):\n    #     s = np.corrcoef(y.values[p][:N],df_rna[col].values[p][:N])[0,1]\n    #     df_importances.loc[col, target_name +' Corr Pearson DownSample'+str(down_sample_percent)  ]  = s\n        X_train = pd.Series(df_rna[col].values[p][:N]).to_frame()\n        X_test = pd.Series(df_rna[col].values[p][N:]).to_frame()\n        model.fit(X_train,y_train)\n        y_pred = model.predict(X_test)# ,y_train)\n\n        s = np.corrcoef(y_test,y_pred)[0,1]\n        df_importances.loc[col, target_name +' Corr LGB DownSample'+str(down_sample_percent) +' Trial'+str(trial)+' n_estimators'+str(n_estimators)  ]  = s\n        l.append(s)\n        if (i%1000==0):\n            I = np.argmax(l)\n            #print(i, col, '%.1f seconds passed'%(time.time() - t0 ), 'Pearson Score', s )\n            print(i, col,'Mean %.6f, max %.6f'%( np.mean(l), np.max(l) ),  \n                  '%.1f seconds passed'%(time.time() - t0 ), 'Argmax:', I, df_rna.columns[I] )            \n    \ndf_importances.head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df_importances.describe() )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## LGB another folds and random seed","metadata":{}},{"cell_type":"code","source":"%%time\nimport lightgbm as lgbm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.model_selection import KFold\nimport time\nt0 = time.time()\n\nif 0:\n    random_state4lgb  = 1\n    model = lgbm.LGBMRegressor(random_state = random_state4lgb) # Ridge(alpha = alpha_selected )\n    #kf = KFold(n_splits=n_splits_for_cross_valdition,  shuffle=True, random_state= random_state_cross_validation )\n    kf = KFold(n_splits=2,  shuffle=True, random_state= random_state4lgb )\n\n\n    y = df_y[target_name]\n    for i,col in enumerate( df_rna.columns):\n        #col = 'ENSG00000135218_CD36'\n        X = df_rna[[col]]\n        y_pred = cross_val_predict(model, X, y, cv=kf)\n        s = np.corrcoef(y,y_pred)[0,1]\n        df_importances.loc[col, target_name +' Corr LGB RS'+str(random_state4lgb)  ]  = s\n\n        if (i%1000==0): \n            print(i, col, '%.1f seconds passed'%(time.time() - t0 ))\n\n    df_importances\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_importances.to_csv('importances_univariable_NIPS22_'+str(target_name)+'.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nn_x_subplots = 3; c = 0 \nfor col in df_importances.columns:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0    \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n    plt.hist(df_importances[col], bins = 100)\n    plt.title(col)\nplt.show()\n    \ndisplay(df_importances.describe() )    \ndf_importances.sort_values(col, ascending = False, key = abs ).head(30)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of importances ","metadata":{}},{"cell_type":"code","source":"%%time\n\nn_x_subplots = 3; c = 0 \nfor col in df_importances.columns:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0    \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n    plt.hist(df_importances[col], bins = 100)\n    plt.title(col)\nplt.show()\n    \ndisplay(df_importances.describe() )    \ndf_importances.sort_values(col, ascending = False, key = abs ).head(20)\n\n\ndf_stat = pd.DataFrame()\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n        l = []\n        if 'LGB' not in col1:\n            v1 = df_importances[col1].sort_values(ascending = False, key = abs )\n        else:\n            v1 = df_importances[col1].sort_values(ascending = False)\n        if 'LGB' not in col2:\n            v2 = df_importances[col2].sort_values(ascending = False, key = abs )\n        else:\n            v2 = df_importances[col2].sort_values(ascending = False)\n            \n        for k in [100,500,1000]:# range(len(v1)): # [10,25,50, 75, 100, 1000, 2000, 4000]:\n            s = set(v1.index[:k] ) & set(v2.index[:k] ) \n            l.append(len(s) )\n            # print(k, len(s))   \n            \n        print(col1,col2, len(l))\n        df_stat.loc[col1+' VS '+col2, 'Top100 Intersection'] = l[0]\n        df_stat.loc[col1+' VS '+col2, 'Top500 Intersection'] = l[1]\n        df_stat.loc[col1+' VS '+col2, 'Top1000 Intersection'] = l[2]\n        \ndisplay(df_stat )\n\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n            \n        if 'LGB' not in col1:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False, key = abs).abs()).rank()\n        else:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False)).rank()\n            \n        v1 = df_importances_ranks[col1] - df_importances_ranks[col2]\n        \n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top100'] = np.max( np.abs( v1.values[:100] ) )\n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top500'] = np.max( np.abs( v1.values[:500] ) )\n        df_stat.loc[col1+' VS '+col2, 'Max Rank Diffr Top1000'] = np.max( np.abs( v1.values[:1000] ) ) \n\ndisplay(df_stat)\n\n#%%time\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n        print(col1,col2)\n        l = []\n        if 'LGB' not in col1:\n            v1 = df_importances[col1].sort_values(ascending = False, key = abs )\n        else:\n            v1 = df_importances[col1].sort_values(ascending = False)\n        if 'LGB' not in col2:\n            v2 = df_importances[col2].sort_values(ascending = False, key = abs )\n        else:\n            v2 = df_importances[col2].sort_values(ascending = False)\n            \n        for k in range(len(v1)): # [10,25,50, 75, 100, 1000, 2000, 4000]:\n            s = set(v1.index[:k] ) & set(v2.index[:k] ) \n            l.append(len(s) )\n            # print(k, len(s))       \n\n        fig = plt.figure(figsize = (20,6))\n        plt.title(col1+' ' + col2, fontsize = 20 )\n        plt.plot(l ,label = 'original data - counts intersection' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x,y,2)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:50],y[:50],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-50 points' )\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()            \n        \n        \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l ,label = 'original data - counts intersection' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x,y,2)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:50],y[:50],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation based on 0-50 points' )\n\n        y = np.array(l)\n        x = np.arange(len(y))\n        p = np.polyfit(x[:150],y[:150],1)\n        print(p)\n        plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-150 points' )\n\n        plt.xlim([0,500])\n        plt.ylim([0,250])\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()        \n        \n        ##################################################################33\n        l2 = []\n        for k in range(10,len(l)):\n            l2.append(l[k]/k)\n        \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l2 ,label = 'intersection topK / K' )\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()                \n            \n        fig = plt.figure(figsize = (20,10) )\n        plt.title(col1+' ' + col2 , fontsize = 20)\n        plt.plot(l2[:500] ,label = 'intersection topK / K' )\n\n        plt.xlim([0,500])\n\n        plt.legend(fontsize = 20 )\n        plt.grid()\n        plt.show()                        \n\nN = df_importances.shape[1]\nfor i in range(N):\n    for j in range(N):\n        if j<=i: continue\n        col1 = df_importances.columns[i];         col2 = df_importances.columns[j];\n            \n        if 'LGB' not in col1:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False, key = abs).abs()).rank()\n        else:\n            df_importances_ranks = (-df_importances.sort_values(col1, ascending = False)).rank()\n            \n        v1 = df_importances_ranks[col1] - df_importances_ranks[col2]\n        \n        fig = plt.figure(figsize = (20,6))\n        plt.title('Ranks differences: ' + col1+' ' + col2, fontsize = 20 )\n        plt.plot(v1.values ) \n        plt.grid()\n        plt.show()\n        \n        fig = plt.figure(figsize = (20,6))\n        plt.title('Ranks differences: ' + col1+' ' + col2, fontsize = 20 )\n        plt.plot(v1.values[:500] ) \n        plt.grid()\n        plt.show()                \n        \ndf_stat.to_csv('stat_importances_univariable_NIPS22_'+str(target_name)+'.csv')        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( v1.index[:20] )\nprint( v2.index[:20] )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('stat_importances_univariable_NIPS22_'+str(target_name)+'.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Display text statistics on importances ","metadata":{}},{"cell_type":"code","source":"pd.set_option('display.max_rows', 500)\npd.set_option('display.max_columns', 500)\npd.set_option('display.width', 1000)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df_stat)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.sort_values(df_importances.columns[0],  ascending = False ).head(50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.corr()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_importances.describe()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tt = time.time() - t0start\nprint('%.1f seconds ( = %.1f minutes, = %.2f hours) passed'%( tt, tt/60, tt/3600 ) )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}