{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# SCP Py-boost, recommender system and ET \n\nThis notebook contains the complete code for a submission to the [Open Problems – Single-Cell Perturbations](https://www.kaggle.com/competitions/open-problems-single-cell-perturbations) Kaggle competition. \n\nThe notebook's main characteristics are:\n- There are four models\n  - Py-boost\n  - A recommender system based on ridge regression\n  - A recommender system based on k nearest neighbors\n  - ExtraTrees\n- All models deal with t-scores rather than log10pvalues. t-scores are the natural representation for differential expression in the given dataset. The log10pvalue is just a nonlinear transformation which makes modeling harder.\n- One of the models (k nearest neighbors) is fed with **data augmentation**: If we know what differential expressions two compounds produce in a cell type separately, we may assume that a mixture of the two compounds will produce a differential expression which is the average of the two single-compound differential expressions.\n- All models are fully cross-validated.\n- We evaluate how the models perform when noise is added to the training data.\n- There is no dependence on external datasets.\n- At the end, only the Py-boost predictions are submitted to the competition\n\nMore information can be found in the accompanying [EDA which makes sense ⭐️⭐️⭐️⭐️⭐️](https://www.kaggle.com/code/ambrosm/scp-eda-which-makes-sense).\n","metadata":{}},{"cell_type":"code","source":"!pip install py-boost","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"scrolled":true,"execution":{"iopub.status.busy":"2023-12-01T01:58:26.823579Z","iopub.execute_input":"2023-12-01T01:58:26.824483Z","iopub.status.idle":"2023-12-01T01:58:39.560590Z","shell.execute_reply.started":"2023-12-01T01:58:26.824422Z","shell.execute_reply":"2023-12-01T01:58:39.559472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom colorama import Fore, Back, Style\nfrom numpy.polynomial import Polynomial\nfrom scipy.stats import norm, skew, kurtosis\nimport scipy\nfrom sklearn.linear_model import Ridge\nfrom sklearn.decomposition import PCA\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom sklearn.ensemble import ExtraTreesRegressor\nfrom itertools import combinations\nfrom py_boost import GradientBoosting\n\nnp.set_printoptions(linewidth=195, edgeitems=3)\npd.set_option(\"min_rows\", 6)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-12-01T01:58:39.562645Z","iopub.execute_input":"2023-12-01T01:58:39.562948Z","iopub.status.idle":"2023-12-01T01:58:57.641304Z","shell.execute_reply.started":"2023-12-01T01:58:39.562922Z","shell.execute_reply":"2023-12-01T01:58:57.640478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some functions","metadata":{}},{"cell_type":"code","source":"def mean_rowwise_rmse(y_true, y_pred):\n    \"\"\"Competition metric\n    \n    Calling convention like in sklearn.metrics\n    \"\"\"\n    mrrmse = np.sqrt(np.square(y_true - y_pred).mean(axis=1)).mean()\n    return mrrmse\n\ndef de_to_t_score(de):\n    \"\"\"Convert log10pvalues to t-scores\n    \n    Parameter:\n    de: array or DataFrame of log10pvalues\n    \n    Return value:\n    t_score: array or DataFrame of t-scores\n    \"\"\"\n    p_value = 10 ** (-np.abs(de))\n#     return - scipy.stats.t.ppf(p_value / 2, df=420) * np.sign(de)\n    return - norm.ppf(p_value / 2) * np.sign(de)\n\ndef t_score_to_de(t_score):\n    \"\"\"Convert t-scores to log10pvalues (inverse of de_to_t_score)\n    \n    Parameter:\n    t_score: array or DataFrame of t-scores\n    \n    Return value:\n    de: array or DataFrame of log10pvalues\n    \"\"\"\n#     p_value = scipy.stats.t.cdf(- np.abs(t_score), df=420) * 2\n    p_value = norm.cdf(- np.abs(t_score)) * 2\n    p_value = p_value.clip(1e-180, None)\n    return - np.log10(p_value) * np.sign(t_score)\n","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-12-01T01:58:57.642426Z","iopub.execute_input":"2023-12-01T01:58:57.642912Z","iopub.status.idle":"2023-12-01T01:58:57.650387Z","shell.execute_reply.started":"2023-12-01T01:58:57.642886Z","shell.execute_reply":"2023-12-01T01:58:57.649541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data\n\n- `de_train`\n- `adata_obs` --> `cell_type_ratio` of shape (6,)\n","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\n# fn = '/kaggle/input/scp-3-merge/de_train.parquet'\nde_train = pd.read_parquet(fn)# , index_col = 0)\n\nfn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\nid_map = pd.read_csv(fn, index_col = 0)\n# display(id_map)\n\n# 18211 genes\ngenes = de_train.columns[5:] \nde_train_indexed = de_train.set_index(['cell_type', 'sm_name'])[genes]\n\n# All 146 sm_names\nsm_names = sorted(de_train.sm_name.unique())\n# Determine the 17 compounds (including the two control compounds) with data for almost all cell types\ntrain_sm_names = de_train.query(\"cell_type == 'B cells'\").sm_name.sort_values().values\n# The other 129 sm_names\ntest_sm_names = [sm for sm in sm_names if sm not in train_sm_names]\n# The three control sm_names\ncontrols3 = ['Dabrafenib', 'Belinostat', 'Dimethyl Sulfoxide']\n\n# All 6 cell types\ncell_types = ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells', 'B cells', 'Myeloid cells']\n# Determine the 4 cell types with data for almost all compounds\n# ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells']\ntrain_cell_types = de_train.query(\"sm_name == 'Vorinostat'\").cell_type.sort_values().values\n# The other 2 cell types: ['B cells', 'Myeloid cells']\ntest_cell_types = [ct for ct in cell_types if ct not in train_cell_types]\n","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:58:57.652681Z","iopub.execute_input":"2023-12-01T01:58:57.652967Z","iopub.status.idle":"2023-12-01T01:59:01.948120Z","shell.execute_reply.started":"2023-12-01T01:58:57.652943Z","shell.execute_reply":"2023-12-01T01:59:01.947315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the experiment, the six cell types appear in certain proportions: T cells CD4+ make up for 42 % whereas only 2 % of the cells are T regulatory cells:","metadata":{}},{"cell_type":"code","source":"adata_obs = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv')\n\n# Cell counts\ncell_count = adata_obs.groupby(['cell_type', 'sm_name']).size()\navg_cell_count = cell_count[~cell_count.index.get_level_values('sm_name').isin(controls3)].groupby('cell_type').mean()\n\n# Cell type ratios (extrapolated from 17 train_sm_names)\ntemp = adata_obs.groupby(['cell_type', 'sm_name']).size().unstack().loc[cell_types]\ncell_type_ratio = temp[list(train_sm_names) + ['Dimethyl Sulfoxide']].sum(axis=1)\ncell_type_ratio /= cell_type_ratio.sum()\nplt.pie(cell_type_ratio, labels=cell_type_ratio.index, autopct=\"%.0f%%\")\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T01:59:01.949266Z","iopub.execute_input":"2023-12-01T01:59:01.949605Z","iopub.status.idle":"2023-12-01T01:59:03.437290Z","shell.execute_reply.started":"2023-12-01T01:59:01.949574Z","shell.execute_reply":"2023-12-01T01:59:03.435927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Models\n\nThis notebook compares four models. \n\nThe first model uses Py-boost and is derived from @alexandervc's [public notebook](https://www.kaggle.com/code/alexandervc/pyboost-secret-grandmaster-s-tool?scriptVersionId=150236954).\n\nThe next two models mimic a [recommender system](https://en.wikipedia.org/wiki/Recommender_system) with collaborative filtering. Like an Internet shop knows how some users have rated some items and from these data predict how other users rate other items, these models know how some cell types interact with some compounds and predict how other cell types interact with other compounds. \n\nTo predict how a user will rate an item, the Internet shop can either look for similar users who have rated the same item, or it can look for similar items which have been rated by the same user. The models in this notebook do both, i.e., either model has two submodels.\n\nTo predict the t-scores for a certain cell type and compound, one of the recommender system models uses ridge regression and the other one uses k-nearest neighbors.\n\nThe third model target-encodes the two categorical features cell_type and sm_name and then fits an ExtraTreesRegressor model. A particularity is that the decision trees are always fully grown, that is, they are completely overfit. The model still makes acceptable predictions because the target-encoding contains a well-tuned amount of noise.","metadata":{}},{"cell_type":"code","source":"def fit_predict_py_boost(de_tr, id_map):\n    \"\"\"Fit the model and predict.\n    \n    Parameters:\n    de_tr: training dataframe of shape (n_samples, 18211), MultiIndex (cell_type, sm_name)\n    id_map: two-column dataframe indicating the validation or test samples (cell_type, sm_name)\n    \n    Returns:\n    de_pred: prediction dataframe of shape (n_samples, 18211), double index matching id_map\n    \n    https://pypi.org/project/py-boost/\n    \"\"\"\n    # Hyperparameters\n    n_components = 50\n    max_depth = 10\n    ntrees = 1000\n    subsample = 1\n    colsample = 0.2\n    lr = 0.01\n\n    # Determine the training cell types (3 or 4)\n    cell_types_tr = de_tr.index[de_tr.index.get_level_values('sm_name') == 'Oxybenzone'].get_level_values('cell_type')\n\n    X_train_categorical = de_tr.index.to_frame()\n\n    #  Dimension reduction\n    reducer = PCA(n_components=n_components, random_state=1)\n    Yt_train_red = reducer.fit_transform(de_to_t_score(de_tr))\n    Yt_train_red = pd.DataFrame(Yt_train_red, index=de_tr.index) # no specific column names\n\n    # Target-encode the two categorical features column-wise\n    # We encode the features with the t-score rather than the log10pvalue\n    # ct_mean has shape (6, 18211) and contains the means of 13 or 14 values each\n    ct_mean = Yt_train_red[X_train_categorical['sm_name'].isin(train_sm_names)].groupby('cell_type').mean()\n    # sm_mean has shape (143, 18211) and contains the means of 3 or 4 values each\n    sm_mean = Yt_train_red[X_train_categorical['cell_type'].isin(cell_types_tr)].groupby('sm_name').mean()\n    X_train_encoded = np.hstack([ct_mean.reindex(X_train_categorical['cell_type']).values,\n                                 sm_mean.reindex(X_train_categorical['sm_name']).values])\n    X_test_encoded =  np.hstack([ct_mean.reindex(id_map['cell_type']).values,\n                                 sm_mean.reindex(id_map['sm_name']).values])\n    \n    # Fit the model\n    model = GradientBoosting('mse',\n                             ntrees=ntrees, \n                             lr=lr, \n                             max_depth=max_depth,\n                             subsample=subsample,\n                             colsample=colsample,\n                             min_data_in_leaf=1,\n                             min_gain_to_split=0,\n                             verbose=10000)           \n    model.fit(X_train_encoded, Yt_train_red)\n\n    # Predict\n    Y_test_pred_red = model.predict(X_test_encoded)\n    Y_test_pred = t_score_to_de(reducer.inverse_transform(Y_test_pred_red))\n    de_pred = pd.DataFrame(Y_test_pred, index=pd.MultiIndex.from_frame(id_map), columns=genes)\n\n    return de_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:59:03.439690Z","iopub.execute_input":"2023-12-01T01:59:03.440638Z","iopub.status.idle":"2023-12-01T01:59:03.462215Z","shell.execute_reply.started":"2023-12-01T01:59:03.440586Z","shell.execute_reply":"2023-12-01T01:59:03.460993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_predict_ridge_recommender(de_tr, id_map):\n    \"\"\"Fit the model and predict.\n    \n    Parameters:\n    de_tr: training dataframe of shape (n_samples, 18211), MultiIndex (cell_type, sm_name)\n    id_map: two-column dataframe indicating the validation or test samples (cell_type, sm_name)\n    \n    Returns:\n    de_pred: prediction dataframes of shape (n_samples, 18211), double index matching id_map\n             If a compound occurs in id_map but not in de_tr, the corresponding row\n             of de_pred will be filled with np.nan\n    \"\"\"\n    # Hyperparameters\n    n_components_in, n_components_out = 7, 70\n    factor_ct = 0.34\n    factor_sm = 0.66\n    \n    # Determine the training cell types (3 or 4)\n    cell_types_tr = de_tr.index[de_tr.index.get_level_values('sm_name') == 'Oxybenzone'].get_level_values('cell_type')\n    \n    # Denoising and dimensionality reduction  \n    reducer_t = make_pipeline(StandardScaler(), PCA(n_components=n_components_out, svd_solver='full', random_state=1))\n    Yt_train_red = reducer_t.fit_transform(de_to_t_score(de_tr))\n    Yt_train_red = pd.DataFrame(Yt_train_red, index=de_tr.index) # no specific column names\n    \n    X_train_categorical = Yt_train_red.index.to_frame()\n\n    # Even more dimensionality reduction\n    Yt_train_red_in = Yt_train_red.iloc[:, :n_components_in]\n    \n    # ct-based model\n    # The model fits a ridge regression to all cell types which have been treated with\n    # the compound to be predicted\n    model_ct = make_pipeline(StandardScaler(), Ridge(3e3))\n    X_train = Yt_train_red_in[X_train_categorical.sm_name.isin(train_sm_names)]\n    X_train = X_train.unstack('sm_name') # 6 rows, index is cell_type\n    X_train.fillna(value=X_train.mean(), inplace=True)\n    X_test = X_train.reindex(id_map['cell_type'])\n    temp_list = []\n    for i in range(len(id_map)):\n        # Y_train has 3 or 4 rows and n_components_out columns\n        Y_train = Yt_train_red[X_train_categorical.sm_name == id_map['sm_name'].iloc[i]]\n        if len(Y_train) > 0:\n            model_ct.fit(X_train.reindex(Y_train.index.get_level_values('cell_type')), Y_train)\n            temp_list.append(model_ct.predict(X_test.iloc[[i]]))\n        else: # compound has been dropped as outlier\n            temp_list.append(np.full((1, Y_train.shape[1]), np.nan))\n    Y_pred_ct = np.vstack(temp_list)\n\n    # sm-based model\n    # The model fits a ridge regression to all (at most 17) compounds which have been applied to\n    # the cell type to be predicted\n    model_sm = make_pipeline(StandardScaler(), Ridge(1e1))\n    X_train = Yt_train_red_in[X_train_categorical.cell_type.isin(cell_types_tr)]\n    X_train = X_train.unstack('cell_type') # 147 rows, index is sm_name\n    X_train.fillna(value=X_train.mean(), inplace=True)\n    X_test = X_train.reindex(id_map['sm_name'])\n    temp_list = []\n    for i in range(len(id_map)):\n        if ~X_test.iloc[i].isna().any():\n            # Y_train has 15 or 17 rows and n_components_out columns\n            Y_train = Yt_train_red[X_train_categorical.cell_type == id_map['cell_type'].iloc[i]]\n            model_sm.fit(X_train.reindex(Y_train.index.get_level_values('sm_name')), Y_train)\n            temp_list.append(model_sm.predict(X_test.iloc[[i]]))\n        else: # compound has been dropped as outlier\n            temp_list.append(np.full((1, Y_train.shape[1]), np.nan))\n    Y_pred_sm = np.vstack(temp_list)\n\n    # Bring the two predictions together\n    Y_test_pred_red = factor_ct * Y_pred_ct + factor_sm * Y_pred_sm\n    Y_test_pred = t_score_to_de(reducer_t.inverse_transform(Y_test_pred_red))\n    de_pred = pd.DataFrame(Y_test_pred, index=pd.MultiIndex.from_frame(id_map), columns=genes)\n    return de_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:59:03.464104Z","iopub.execute_input":"2023-12-01T01:59:03.465051Z","iopub.status.idle":"2023-12-01T01:59:03.495941Z","shell.execute_reply.started":"2023-12-01T01:59:03.464977Z","shell.execute_reply":"2023-12-01T01:59:03.495220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_predict_knn_recommender(de_tr, id_map):\n    \"\"\"Fit the model and predict.\n    \n    Parameters:\n    de_tr: training dataframe of shape (n_samples, 18211), MultiIndex (cell_type, sm_name)\n    id_map: two-column dataframe indicating the validation or test samples (cell_type, sm_name)\n    \n    Returns:\n    de_pred: prediction dataframes of shape (n_samples, 18211), double index matching id_map\n             If a compound occurs in id_map but not in de_tr, the corresponding row\n             of de_pred will be filled with np.nan\n    \"\"\"\n    # Hyperparameters\n    n_components_in, n_components_out = 7, 70\n    factor_ct = 0.32\n    factor_sm = 0.74\n\n    # Determine the training cell types (3 or 4)\n    cell_types_tr = de_tr.index[de_tr.index.get_level_values('sm_name') == 'Oxybenzone'].get_level_values('cell_type')\n    \n    # Denoising and dimensionality reduction  \n    reducer_t = make_pipeline(StandardScaler(), PCA(n_components=n_components_out, svd_solver='full', random_state=1))\n    Yt_train_red = reducer_t.fit_transform(de_to_t_score(de_tr))\n    Yt_train_red = pd.DataFrame(Yt_train_red, index=de_tr.index) # no specific column names\n    \n    # Data augmentation\n    # We use two kinds of data augmentation:\n    # 1. Scaled t-scores (t-score is a quotient of log-fold-change and standard deviation;\n    #    if the variance changes, the t-scores are scaled)\n    # 2. Mixture of two compounds\n    for ct1 in cell_types_tr:\n        for factor in [0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, 1.4]:\n            a = Yt_train_red.query(\"cell_type == @ct1\").reset_index('cell_type', drop=True) * factor\n            a = pd.concat([a], keys=[f\"{ct1}*{factor}\"], names=['cell_type'])\n            Yt_train_red = pd.concat([Yt_train_red, a])\n\n    for sm1 in train_sm_names:\n        for factor in [0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, 1.4]:\n            a = Yt_train_red.query(\"sm_name == @sm1\").reset_index('sm_name', drop=True) * factor\n            a = pd.concat([a], keys=[f\"{sm1}*{factor}\"], names=['sm_name'])\n            a = a.reorder_levels(['cell_type', 'sm_name'])\n            Yt_train_red = pd.concat([Yt_train_red, a])\n\n    for sm1, sm2 in combinations(train_sm_names, 2):\n        a = (2 * Yt_train_red.query(\"sm_name == @sm1\").reset_index('sm_name', drop=True)\n             + Yt_train_red.query(\"sm_name == @sm2\").reset_index('sm_name', drop=True)) / 3\n        a.dropna(inplace=True)\n        a = pd.concat([a], keys=[f\"{sm1}+{sm2} a\"], names=['sm_name'])\n        a = a.reorder_levels(['cell_type', 'sm_name'])\n        b = (Yt_train_red.query(\"sm_name == @sm1\").reset_index('sm_name', drop=True)\n             + 2 * Yt_train_red.query(\"sm_name == @sm2\").reset_index('sm_name', drop=True)) / 3\n        b.dropna(inplace=True)\n        b = pd.concat([b], keys=[f\"{sm1}+{sm2} b\"], names=['sm_name'])\n        b = b.reorder_levels(['cell_type', 'sm_name'])\n        Yt_train_red = pd.concat([Yt_train_red, a, b])\n        \n    X_train_categorical = Yt_train_red.index.to_frame()\n\n    # Even more dimensionality reduction\n    Yt_train_red_in = Yt_train_red.iloc[:, :n_components_in]\n    \n    # ct-based model\n    # The model finds similar cell types which have been treated with\n    # the compound to be predicted and averages their t scores\n    model_ct = KNeighborsRegressor(n_neighbors=7, weights='distance')\n    X_train = Yt_train_red_in[X_train_categorical.sm_name.isin(train_sm_names)]\n    X_train = X_train.unstack('sm_name') # 6 rows * 51 columns (without augmentation), index is cell_type\n    X_train.fillna(value=X_train.mean(), inplace=True)\n    X_test = X_train.reindex(id_map['cell_type'])\n    temp_list = []\n    for i in range(len(id_map)):\n        # Find similar cell types which have a rating for the given sm_name\n        Y_train = Yt_train_red[X_train_categorical.sm_name == id_map['sm_name'].iloc[i]]\n        if len(Y_train) > 0:\n            model_ct.fit(X_train.reindex(Y_train.index.get_level_values('cell_type')), Y_train)\n            temp_list.append(model_ct.predict(X_test.iloc[[i]]))\n        else: # compound has been dropped as outlier\n            temp_list.append(np.full((1, Y_train.shape[1]), np.nan))\n    Y_pred_ct = np.vstack(temp_list)\n\n    # sm-based model\n    # The model finds similar sm_names which have been measured with\n    # the cell type to be predicted and averages their t scores\n    model_sm = KNeighborsRegressor(n_neighbors=9, weights='distance', p=1)\n    X_train = Yt_train_red_in[X_train_categorical.cell_type.isin(cell_types_tr)]\n    X_train = X_train.unstack('cell_type') # 146 rows * 9 columns (without augmentation), index is sm_name\n    X_train.fillna(value=X_train.mean(), inplace=True)\n    X_test = X_train.reindex(id_map['sm_name'])\n    temp_list = []\n    for i in range(len(id_map)):\n        if ~X_test.iloc[i].isna().any():\n            # Find similar sm_names which have a rating for the given cell type\n            # Y_train has 15 or 17 rows if there is no data augmentation\n            # Y_train has 153 rows if all sm_name pairs are augmented\n            Y_train = Yt_train_red[X_train_categorical.cell_type == id_map['cell_type'].iloc[i]]\n            model_sm.fit(X_train.reindex(Y_train.index.get_level_values('sm_name')), Y_train)\n            temp_list.append(model_sm.predict(X_test.iloc[[i]]))\n        else: # compound has been dropped as outlier\n            temp_list.append(np.full((1, Y_train.shape[1]), np.nan))\n    Y_pred_sm = np.vstack(temp_list)\n\n    # Bring the two predictions together\n    Y_test_pred_red = factor_ct * Y_pred_ct + factor_sm * Y_pred_sm\n    Y_test_pred = t_score_to_de(reducer_t.inverse_transform(Y_test_pred_red))\n    de_pred = pd.DataFrame(Y_test_pred, index=pd.MultiIndex.from_frame(id_map), columns=genes)\n    return de_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:59:03.497139Z","iopub.execute_input":"2023-12-01T01:59:03.497672Z","iopub.status.idle":"2023-12-01T01:59:03.521853Z","shell.execute_reply.started":"2023-12-01T01:59:03.497645Z","shell.execute_reply":"2023-12-01T01:59:03.521039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_predict_extratrees(de_tr, id_map):\n    \"\"\"Fit the model and predict.\n    \n    Parameters:\n    de_tr: training dataframe of shape (n_samples, 18211), MultiIndex (cell_type, sm_name)\n    id_map: two-column dataframe indicating the validation or test samples (cell_type, sm_name)\n    \n    Returns:\n    de_pred: prediction dataframes of shape (n_samples, 18211), double index matching id_map\n             If a compound occurs in id_map but not in de_tr, the corresponding row\n             of de_pred will be filled with np.nan\n    \"\"\"\n    # Hyperparameters\n    n_components_in, n_components_out = 35, 200\n    max_features = 0.2\n    n_trees = 10000\n    \n    # Determine the training cell types (3 or 4)\n    cell_types_tr = de_tr.index[de_tr.index.get_level_values('sm_name') == 'Oxybenzone'].get_level_values('cell_type')\n    \n    X_train_categorical = de_tr.index.to_frame()\n    \n    # Denoising and dimensionality reduction  \n    reducer_t = make_pipeline(StandardScaler(), PCA(n_components=n_components_out, svd_solver='full', random_state=1))\n    Yt_train_red = reducer_t.fit_transform(de_to_t_score(de_tr))\n    Yt_train_red = pd.DataFrame(Yt_train_red, index=de_tr.index) # no specific column names\n\n    # Even more dimensionality reduction\n    Yt_train_red_in = Yt_train_red.iloc[:, :n_components_in]\n    \n    # Target-encode the two categorical features without smoothing\n    # We encode the features with the t-score rather than the log10pvalue\n    ct_mean = Yt_train_red_in[X_train_categorical['sm_name'].isin(train_sm_names)].groupby('cell_type').mean() # shape (6, n_components), means of 13 or 14 values each\n    sm_mean = Yt_train_red_in[X_train_categorical['cell_type'].isin(cell_types_tr)].groupby('sm_name').mean() # shape (143, n_components), means of 3 or 4 values each\n    X_train_encoded = np.hstack([ct_mean.reindex(X_train_categorical['cell_type']).values,\n                                 sm_mean.reindex(X_train_categorical['sm_name']).values,\n                                 sm_mean.reindex(X_train_categorical['sm_name']).values * np.sqrt(cell_type_ratio.reindex(X_train_categorical['cell_type']).values.reshape(-1, 1))\n                                ])\n    X_test_encoded =  np.hstack([ct_mean.reindex(id_map['cell_type']).values,\n                                 sm_mean.reindex(id_map['sm_name']).values,\n                                 sm_mean.reindex(id_map['sm_name']).values * np.sqrt(cell_type_ratio.reindex(id_map['cell_type']).values.reshape(-1, 1))\n                                ])\n\n    # Train the model\n    # The model is trained to predict the PCA-transformed t-scores.\n    model = ExtraTreesRegressor(n_estimators=n_trees, max_features=max_features, random_state=1)\n    model.fit(X_train_encoded, Yt_train_red)\n\n    # Predict\n    X_test_encoded = pd.DataFrame(X_test_encoded, index=id_map)\n    X_test_encoded.dropna(inplace=True)\n    Y_test_pred_red = model.predict(X_test_encoded.values)\n    Y_test_pred = t_score_to_de(reducer_t.inverse_transform(Y_test_pred_red))\n    de_pred = pd.DataFrame(Y_test_pred, index=X_test_encoded.index, columns=genes)\n    de_pred = de_pred.reindex(id_map)\n    return de_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:59:03.523745Z","iopub.execute_input":"2023-12-01T01:59:03.524056Z","iopub.status.idle":"2023-12-01T01:59:03.536553Z","shell.execute_reply.started":"2023-12-01T01:59:03.524027Z","shell.execute_reply":"2023-12-01T01:59:03.535763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cross-validation\n\nWe cross-validate the four models according to the cv scheme documented in my [Quickstart notebook](https://www.kaggle.com/code/ambrosm/scp-quickstart) and keep the out-of-fold predictions.\n\nAfter cross-validation, we ensemble the out-of-fold predictions (unweighted mean) and score the ensemble.","metadata":{}},{"cell_type":"code","source":"%%time\ndef cross_val_log10pvalue(predictor, noise=0):\n    \"\"\"Cross-validate a machine-learning model\n    \n    Parameters:\n    predictor: function which takes two parameters\n        (training data and id_map) and returns\n        the predictions\n    noise: standard deviation of noise to be added to the t-scores\n        \n    Globals:\n    de_oof_dict: dictionary into which the oof predictions are inserted\n    mrrmse_noise_list: list into which the noise level and the oof score are inserted\n    removed_compounds: list of outlier compounds\n    \"\"\"\n    mrrmse_list = []\n    t_oof_list, de_oof_list = [], []\n    for fold, val_cell_type in enumerate(train_cell_types):\n        # Split the data into training and validation\n        # mask_va: 127 or 129 validation rows per fold, total 514 in four folds\n        mask_va = ((de_train['cell_type'] == val_cell_type) &\n                   ~de_train['sm_name'].isin(list(train_sm_names) + ['Dimethyl Sulfoxide']))\n        if mask_va.sum() == 0: continue\n        mask_va = mask_va.values\n        # mask_tr: 485 or 487 training rows\n        mask_tr = ~mask_va\n        \n        de_tr = de_train_indexed[mask_tr] # shape (48x, 18211), double index\n        de_va = de_train_indexed[mask_va] # shape (12x, 18211), double index\n        \n        # Drop outliers from training and validation\n        # If removed_compounds is nonempty, some prediction rows will contain np.nan\n        de_tr = de_tr.query(\"~sm_name.isin(@removed_compounds)\")\n#         de_va = de_va.query(\"~sm_name.isin(@removed_compounds)\")\n\n        # Add noise to the training t-scores\n        if noise > 0:\n            if fold == 0:\n                print(f\"{Fore.RED}{Style.BRIGHT}Adding noise of scale {noise:.2f}{Style.RESET_ALL}\")\n            rng = np.random.default_rng(1)\n            de_tr = t_score_to_de(de_to_t_score(de_tr) + rng.normal(scale=noise, size=de_tr.shape))\n    \n        # Fit the model and predict validation log10pvalues\n        de_pred = predictor(de_tr, de_va.index.to_frame())\n        \n        # Update out-of-fold predictions and score\n        de_oof_list.append(de_pred)\n        mrrmse = mean_rowwise_rmse(de_va, de_pred.values) # the competition metric\n        print(f\"# Fold {fold}: de_mrrmse={mrrmse:5.3f}   val='{val_cell_type}'\")\n        mrrmse_list.append(mrrmse)\n\n    # Collect out-of-fold predictions and scores\n    de_oof = pd.concat(de_oof_list, axis=0)\n    mrrmse = mean_rowwise_rmse(de_train_indexed.reindex(de_oof.index),\n                               de_oof)\n    name = predictor.__name__[12:]\n    print(f\"{Fore.GREEN}{Style.BRIGHT}# Overall \"\n          f\"de_mrrmse={mrrmse:5.3f} \"\n          f\"{tuple((np.array(mrrmse_list) * 1000).round(0).astype(int))} \"\n          f\"{name}{Style.RESET_ALL}\")\n    if noise == 0:\n        de_oof_dict[name] = de_oof\n    mrrmse_noise_list.append((name, noise, mrrmse))\n    print()\n    return\n\n# Define outliers which are excluded from training and validation\n# removed_compounds = ['AT13387', 'Alvocidib', 'BAY 61-3606', 'BMS-387032', \n#                      'Belinostat', 'CEP-18770 (Delanzomib)', 'CGM-097', 'CGP 60474', \n#                      'Dabrafenib', 'Ganetespib (STA-9090)', 'I-BET151', 'IN1451', \n#                      'LY2090314', 'MLN 2238', 'Oprozomib (ONX 0912)', \n#                      'Proscillaridin A;Proscillaridin-A', 'Resminostat',\n#                      'Scriptaid', 'UNII-BXU45ZH6LI', 'Vorinostat']\nremoved_compounds = []\n\n# Cross-validate the four models (saving the oof predictions)\npredictors = [fit_predict_py_boost, fit_predict_ridge_recommender, fit_predict_knn_recommender, fit_predict_extratrees] \nde_oof_dict, mrrmse_noise_list = {}, []\nfor predictor in predictors:\n    cross_val_log10pvalue(predictor)\n\n# Ensemble the oof predictions\nde_oof = sum(de_oof_dict.values()) / len(de_oof_dict)\nde_true = de_train_indexed.reindex_like(de_oof)\nprint(f\"# Ensemble MRRMSE: {mean_rowwise_rmse(de_true, de_oof):.3f}\")","metadata":{"execution":{"iopub.status.busy":"2023-12-01T01:59:03.539042Z","iopub.execute_input":"2023-12-01T01:59:03.539305Z","iopub.status.idle":"2023-12-01T02:07:04.687211Z","shell.execute_reply.started":"2023-12-01T01:59:03.539282Z","shell.execute_reply":"2023-12-01T02:07:04.686215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Evaluation\n\nThe diagrams show how the average residuals depend on the predicted log10pvalue or t-score.","metadata":{}},{"cell_type":"code","source":"# Analysis of de residuals\nt_oof = de_to_t_score(de_oof)\nde_residuals = de_true - de_oof\n\nprint(\"Analysis of log10pvalue residuals\")\nprint(f\"MRRMSE: {mean_rowwise_rmse(de_train_indexed.reindex_like(de_oof), de_oof.values):.3f}\")\nprint(f\"Mean of residuals: {de_residuals.values.mean():.2f}\")\nprint(f\"Variance of residuals: {de_residuals.values.var():.2f}\")\nprint(f\"Skew of residuals: {skew(de_residuals.values.ravel()):.2f}\")\nprint(f\"Excess kurtosis of residuals: {kurtosis(de_residuals.values.ravel()):.2f}\")\n\nplt.figure(figsize=(12, 4))\n\nplt.subplot(1, 2, 1)\nresidual_df = pd.DataFrame({'oof': de_oof.values.ravel(),\n                            'res': de_residuals.values.ravel()})\nresidual_df['bin'] = pd.qcut(residual_df['oof'], q=500, labels=False, duplicates='drop')\nplt.scatter(residual_df.groupby('bin')['oof'].mean(),\n            residual_df.groupby('bin')['res'].mean(), \n            s=1, color='k', label='all')\npoly = Polynomial.fit(residual_df.groupby('bin')['oof'].mean(),\n                      residual_df.groupby('bin')['res'].mean(),\n                      deg=1, domain=[])\nprint(poly)\nplt.plot(residual_df.groupby('bin')['oof'].mean(),\n         poly(residual_df.groupby('bin')['oof'].mean()),\n         color='cyan', label='poly')\nplt.xlabel('predicted log10pvalue')\nplt.ylabel('residual')\nplt.axhline(0, color='gray')\n\nfor cell_type in train_cell_types:\n    try:\n        residual_df = pd.DataFrame({'oof': de_oof.loc[cell_type].values.ravel(),\n                                    'res': de_residuals.loc[cell_type].values.ravel()})\n    except KeyError:\n        continue\n    residual_df['bin'] = pd.qcut(residual_df['oof'], q=2000, labels=False, duplicates='drop')\n    plt.scatter(residual_df.groupby('bin')['oof'].mean(),\n                residual_df.groupby('bin')['res'].mean(), \n                s=1, label=cell_type)\nplt.xlabel('predicted log10pvalue')\nplt.ylabel('residual')\nplt.legend()\nplt.title('Mean oof residual')\n\nplt.subplot(1, 2, 2)\nresidual_df = pd.DataFrame({'oof': t_oof.values.ravel(),\n                            'res': de_residuals.values.ravel()})\nresidual_df['bin'] = pd.qcut(residual_df['oof'], q=500, labels=False, duplicates='drop')\nplt.scatter(residual_df.groupby('bin')['oof'].mean(),\n            residual_df.groupby('bin')['res'].mean(), \n            s=1, color='k', label='all')\npoly = Polynomial.fit(residual_df.groupby('bin')['oof'].mean(),\n                      residual_df.groupby('bin')['res'].mean(),\n                      deg=1, domain=[])\nprint(poly)\nplt.plot(residual_df.groupby('bin')['oof'].mean(),\n         poly(residual_df.groupby('bin')['oof'].mean()),\n         color='cyan', label='poly')\nplt.xlabel('predicted t-score')\nplt.ylabel('residual')\nplt.axhline(0, color='gray')\n\nfor cell_type in train_cell_types:\n    try:\n        residual_df = pd.DataFrame({'oof': t_oof.loc[cell_type].values.ravel(),\n                                    'res': de_residuals.loc[cell_type].values.ravel()})\n    except KeyError:\n        continue\n    residual_df['bin'] = pd.qcut(residual_df['oof'], q=2000, labels=False, duplicates='drop')\n    plt.scatter(residual_df.groupby('bin')['oof'].mean(),\n                residual_df.groupby('bin')['res'].mean(), \n                s=1, label=cell_type)\nplt.xlabel('predicted t-score')\nplt.ylabel('residual')\nplt.legend()\nplt.title('Mean oof residual')\n\nplt.show()","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-12-01T02:07:04.688848Z","iopub.execute_input":"2023-12-01T02:07:04.689247Z","iopub.status.idle":"2023-12-01T02:07:15.966845Z","shell.execute_reply.started":"2023-12-01T02:07:04.689212Z","shell.execute_reply":"2023-12-01T02:07:15.965957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The next diagram shows the behavior of the four models under noise. The models are robust against a small amount of noise, but with enough noise every model can be killed (of course).","metadata":{}},{"cell_type":"code","source":"# Experiment with noise in the training data\nfor noise in [0.25, 0.5, 1.0, 2.0, 3.0]:\n    for predictor in predictors:\n        cross_val_log10pvalue(predictor, noise)\n    ","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"noise_df = pd.DataFrame(mrrmse_noise_list, columns=['name', 'noise', 'mrrmse'])\nsns.scatterplot(noise_df, x='noise', y='mrrmse', hue='name')\nplt.title('Model performance in presence of noise')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T01:55:40.406132Z","iopub.status.idle":"2023-12-01T01:55:40.406945Z","shell.execute_reply.started":"2023-12-01T01:55:40.406670Z","shell.execute_reply":"2023-12-01T01:55:40.406696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\n\nWe re-fit the Py-boost to the full dataset. Be aware that Py-boost is nondeterministic and the Kaggle leaderboard score may differ if the notebook is re-run. ","metadata":{}},{"cell_type":"code","source":"def submit():\n    \"\"\"Refit the selected models and write the submission file.\"\"\"\n    \n    # Drop outliers from training\n    de_tr = de_train_indexed.query(\"~sm_name.isin(@removed_compounds)\")\n    \n    # Fit all models and average their predictions\n    pred_list = [fit_predict(de_tr, id_map) for fit_predict in predictors]\n    de_pred = sum(pred_list) / len(pred_list)\n        \n    # Test for missing values\n    if de_pred.isna().any().any():\n        print(\"Warning: This submission contains missing values. \"\n              \"Don't submit it!\")\n        \n    # Create the submission dataframe\n    submission = pd.DataFrame(de_pred.values, columns=genes, index=id_map.index)\n    display(submission)\n    print(f'Variance of submission: {submission.values.var():.2f},   min = {submission.values.min():.2f}, max = {submission.values.max():.2f}')\n    \n    # Compare the histograms of validation and test for a plausibility check\n    plt.figure(figsize=(12, 2))\n    plt.hist(submission.values.ravel(),\n             bins=np.linspace(-5, 5, 501),\n             density=True,\n             label='submission',\n             color='orange')\n    plt.hist(de_oof.values.ravel(),\n             bins=np.linspace(-5, 5, 501),\n             density=True,\n             label='oof',\n             color='#0000ff80')\n    plt.legend()\n    plt.show()\n\n    # Write the files\n    de_oof.to_csv('de_oof.csv')\n    submission.to_csv('submission.csv')        \n\npredictors = [fit_predict_py_boost] \nsubmit()","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-12-01T02:07:35.247510Z","iopub.execute_input":"2023-12-01T02:07:35.247889Z","iopub.status.idle":"2023-12-01T02:10:08.873773Z","shell.execute_reply.started":"2023-12-01T02:07:35.247857Z","shell.execute_reply":"2023-12-01T02:10:08.872691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}