{"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":"<div style=\"background-color: #4CAF50; padding: 20px; border-radius: 10px; color: white; text-align: center;\">\n    <h1 style=\"display: inline;\" fontweigth='bold'>COMPETITION OVERVIEW 🧬</h1>\n</div>\n\n<div style=\"background-color: #E6F7FF; padding: 20px; border-radius: 10px; margin-top: 20px; border: 1px solid #B3D4E9; text-align: justify; color: black;\">\n    <h3>What are Single-Cell Perturbations?</h3>\n    <p>Single-cell perturbations refer to the process of introducing external agents (like drugs or chemicals) to individual cells to observe changes in their state or behavior. By studying these perturbations at the single-cell level, scientists can gain a deeper understanding of cellular responses and identify potential targets for therapeutic interventions. This approach provides a high-resolution view of cellular dynamics, enabling the dissection of heterogeneity within cell populations and the identification of unique cellular responses that might be masked in bulk analyses.</p>\n</div>\n\n<p style=\"text-align: justify;\">Human biology is intricate, with about 37 trillion cells organized into tissues, organs, and systems. The recent advancements in single-cell technologies offer deep insights into cell and tissue functionality at the DNA, RNA, and protein levels. The current challenge is mapping causal links between chemical disruptions and their aftermath on cell states for drug innovation. This competition, organized by Open Problems in Single-Cell Analysis in alliance with Cellarity, seeks solutions to predict chemical disturbances in new cell types, propelling medicinal development.</p>\n\n<h3>Evaluation Criteria</h3>\n<p>The competition's metric is the Mean Rowwise Root Mean Squared Error (MRRMSE), defined as:</p>\n\\[\n\\textrm{MRRMSE} = \\frac{1}{R}\\sum_{i=1}^R\\left(\\frac{1}{n} \\sum_{j=1}^{n} (y_{ij} - \\widehat{y}_{ij})^2\\right)^{1/2}\n\\]\n\n\n\n<h3>Technical Details About the Experiment</h3>\n<img src=\"https://www.googleapis.com/download/storage/v1/b/kaggle-user-content/o/inbox%2F4308072%2Fdbce2b5e0a2b8b9691502b31c5feb5b6%2FScreenshot%202023-08-25%20at%206.20.53%20PM.png?generation=1693002071834598&alt=media\" alt=\"Experiment Design\" width=\"500px\">\n<p style=\"text-align: justify;\">To grasp the dataset, it's crucial to understand the plate design used to gauge the treatment effect. PBMCs from donors were thawed and placed on 96-well plates. Some columns of the plates were dedicated to controls, while the remaining wells contained different compounds. This entire setup ensures a variety of cell types within each well, from T cells, B cells to Myeloid cells.</p>\n\n<h3>How We Calculate Differential Expression</h3>\n<img src=\"https://www.googleapis.com/download/storage/v1/b/kaggle-user-content/o/inbox%2F4308072%2F01a5bae1fad53d3ac875f7e86744805d%2FScreenshot%202023-10-03%20at%2011.33.03%20AM.png?generation=1696347321070915&alt=media\" alt=\"Differential Expression Calculation\" width=\"500px\">\n<p style=\"text-align: justify;\">Differential expression (DE) helps gauge the impact of an experimental perturbation on each gene's expression level. This competition's participants are tasked with modeling DE, and the process involves several computational steps, including pseudobulking and linear model fitting, with the primary aim of understanding the impact of each compound.</p>\n\n<h3>Data Splits</h3>\n<ul>\n<li>Training: All compounds in T, NK cells and a mere 10% of the compounds in B and Myeloid cells.</li>\n<li>Public Testing: Incorporates 50 compounds in B and Myeloid cells, chosen at random.</li>\n<li>Private Testing: Consists of 79 randomly chosen compounds in B and Myeloid cells.</li>\n</ul>\n\n<h3>Files and Their Descriptions</h3>\n<ul>\n<li><code>de_train.parquet</code>: Consolidated differential expression information.</li>\n<li><code>adata_train.parquet</code>: Raw count and normalized data in a disaggregated form.</li>\n<li><code>adata_obs_meta.csv</code>: Observation metadata pertinent to <code>adata_train</code>.</li>\n<li><code>multiome_train.parquet</code>: Supplementary 10x Multiome data for baseline samples.</li>\n<li><code>multiome_obs_meta.csv</code>: Observation metadata for <code>multiome_train.parquet</code>.</li>\n<li><code>multiome_var_meta.csv</code>: Variable metadata for <code>multiome_train.parquet</code>.</li>\n<li><code>id_map.csv</code>: Identifies the cell_type/sm_name pair to be predicted for the given id.</li>\n<li><code>sample_submission.csv</code>: Example submission file.</li>\n</ul>\n\n<p>This competition strives to harness single-cell data to forecast responses to chemical perturbations, holding significant potential for drug discovery and development.</p>\n","metadata":{}},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 1px; border: 1px solid #0056b3; border-radius: 3px; color: white;\"><h3 align=center>Downloading data</h3></div>","metadata":{}},{"cell_type":"code","source":"import os\nimport zipfile\nimport tarfile\n# import kaggle # You must change this \nfrom tqdm import tqdm\n\ndef download_and_extract_kaggle_data(competition_name: str, save_dir: str = 'raw_data') -> None:\n    \n    #Download, extract Kaggle competition data with progress bars, and then delete the downloaded file.\n    \n    #Parameters:\n    #- competition_name (str): Name of the Kaggle competition.\n    #- save_dir (str, optional): Directory where the data should be saved. Defaults to 'raw_data'.\n    \n    \n    # Download the competition data with progress\n    kaggle.api.competition_download_cli(competition_name, quiet=False)\n    \n    # Check the file extension of the downloaded data\n    if os.path.exists(f\"{competition_name}.zip\"):\n        file_path = f\"{competition_name}.zip\"\n        is_zip = True\n    elif os.path.exists(f\"{competition_name}.tar.gz\"):\n        file_path = f\"{competition_name}.tar.gz\"\n        is_zip = False\n    else:\n        raise ValueError(\"Downloaded file is neither .zip nor .tar.gz\")\n    \n    # Ensure the save directory exists\n    if not os.path.exists(save_dir):\n        os.makedirs(save_dir)\n    \n    # Extract the data with progress bar\n    if is_zip:\n        with zipfile.ZipFile(file_path, 'r') as zip_ref:\n            total_files = len(zip_ref.namelist())\n            with tqdm(total=total_files, unit=\"file\") as pbar:\n                for member in zip_ref.namelist():\n                    zip_ref.extract(member, save_dir)\n                    pbar.update(1)\n    else:\n        with tarfile.open(file_path, 'r:gz') as tar_ref:\n            total_files = len(tar_ref.getnames())\n            with tqdm(total=total_files, unit=\"file\") as pbar:\n                for member in tar_ref.getnames():\n                    tar_ref.extract(member, save_dir)\n                    pbar.update(1)\n    \n    try:\n        os.remove(file_path)\n    except (PermissionError, OSError) as e:\n        print(f\"Permission denied when trying to delete {file_path}. Please delete it manually.\")\n\n    \n    print(f\"\\nData for {competition_name} has been downloaded, extracted to {save_dir}, and the downloaded file has been deleted.\")\n\n# Download the data for the Kaggle competition\n# download_and_extract_kaggle_data('open-problems-single-cell-perturbations')","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:03.120208Z","iopub.execute_input":"2023-10-15T13:23:03.120594Z","iopub.status.idle":"2023-10-15T13:23:03.139255Z","shell.execute_reply.started":"2023-10-15T13:23:03.120565Z","shell.execute_reply":"2023-10-15T13:23:03.137995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 5px; border: 1px solid #0056b3; border-radius: 4px; color: white;\"><h3 align=center>Basic Informations about the dataset</h3></div>","metadata":{}},{"cell_type":"code","source":"import pandas as pd\n\n# Carregando o arquivo CSV\ndata = pd.read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet\")\n\n# Exibindo as primeiras linhas do dataframe\ndata.head()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:03.141398Z","iopub.execute_input":"2023-10-15T13:23:03.141763Z","iopub.status.idle":"2023-10-15T13:23:06.320054Z","shell.execute_reply.started":"2023-10-15T13:23:03.141735Z","shell.execute_reply":"2023-10-15T13:23:06.318817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv\")\nsample_submission.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:06.321254Z","iopub.execute_input":"2023-10-15T13:23:06.321783Z","iopub.status.idle":"2023-10-15T13:23:09.259555Z","shell.execute_reply.started":"2023-10-15T13:23:06.321752Z","shell.execute_reply":"2023-10-15T13:23:09.258170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id_map = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\nid_map.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:09.261080Z","iopub.execute_input":"2023-10-15T13:23:09.261528Z","iopub.status.idle":"2023-10-15T13:23:09.287430Z","shell.execute_reply.started":"2023-10-15T13:23:09.261487Z","shell.execute_reply":"2023-10-15T13:23:09.286332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.info(show_counts=True)","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:09.291056Z","iopub.execute_input":"2023-10-15T13:23:09.291522Z","iopub.status.idle":"2023-10-15T13:23:10.133889Z","shell.execute_reply.started":"2023-10-15T13:23:09.291482Z","shell.execute_reply":"2023-10-15T13:23:10.132759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\nAs we can see, the majority of our dataset consists of numerical columns, with 4 categorical columns and one boolean column. So far, so good! 📊\n\nLet's delve deeper into our categorical columns. I have a hunch that some columns might be imbalanced. This isn't ideal for me! It poses an additional challenge! 🤔💭\n\n---","metadata":{}},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 5px; border: 1px solid #0056b3; border-radius: 4px; color: white;\"><h3 align=center>About the Categorical Columns</h3></div>","metadata":{}},{"cell_type":"code","source":"from matplotlib import pyplot as plt\nimport seaborn as sns\n\n\n# Identify categorical columns\ncategorical_columns = data.select_dtypes(include=['object']).columns.tolist()\n\ndef plot_all_categorical_counts_with_statistics(df, categorical_columns, max_label_length=10):\n    \"\"\"\n    Plots the value counts of all categorical columns in a dataframe using subplots.\n    Truncates the x-axis labels to a maximum length for better visibility.\n    Adds lines to indicate mean and median values.\n    \n    Parameters:\n    - df: DataFrame containing the data.\n    - categorical_columns: List of names of the categorical columns to plot.\n    - max_label_length: Maximum length for x-axis labels.\n    \"\"\"\n    # Calculate number of rows for subplots based on number of categorical columns\n    n_rows = len(categorical_columns)\n    \n    # Set up the figure with subplots\n    fig, axes = plt.subplots(n_rows, 1, figsize=(24, 8 * n_rows))\n    \n    # Loop through each categorical column and plot its value counts\n    for i, column_name in enumerate(categorical_columns):\n        counts = df[column_name].value_counts()\n        \n        # Truncate x-axis labels if they exceed max_label_length\n        truncated_labels = [label[:max_label_length] + \"...\" if len(label) > max_label_length else label for label in counts.index]\n        \n        sns.barplot(x=truncated_labels, y=counts.values, palette=\"viridis\", ax=axes[i])\n        \n        # Add lines to indicate mean and median values\n        mean_val = counts.mean()\n        median_val = counts.median()\n        \n        axes[i].axhline(mean_val, color='red', linestyle='--', label=f'Mean: {mean_val:.2f}')\n        axes[i].axhline(median_val, color='blue', linestyle='-', label=f'Median: {median_val:.2f}')\n        \n        # Set subplot details\n        axes[i].set_title(f\"Value Counts of {column_name}\", fontsize=30, fontweight=\"bold\")\n        axes[i].set_ylabel(\"Counts\", fontsize=14, fontweight=\"bold\")\n        axes[i].set_xlabel(column_name, fontsize=14, fontweight=\"bold\")\n        axes[i].tick_params(axis='x', rotation=90)\n        axes[i].grid(axis=\"y\")\n        axes[i].legend(fontsize=14, loc=\"upper right\")\n    \n    \n    # Adjust layout\n    plt.tight_layout(pad=9.0)\n    plt.show()\n\n# Plot value counts for all categorical columns with statistics\nplot_all_categorical_counts_with_statistics(data, categorical_columns)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:10.135191Z","iopub.execute_input":"2023-10-15T13:23:10.135998Z","iopub.status.idle":"2023-10-15T13:23:16.819903Z","shell.execute_reply.started":"2023-10-15T13:23:10.135969Z","shell.execute_reply":"2023-10-15T13:23:16.818511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\nAt first glance, it seems I'm facing some challenges! 😅\n\nI noticed that the `cell_type` column appears to be quite imbalanced. Incorporating the mean and median as features is an idea I've been considering, and I find it intriguing. 🤓 These statistics might provide additional context and valuable information for the model, potentially aiding in its predictions. I've always believed it's beneficial to test different hypotheses during the modeling phase, as unconventional approaches can often yield surprising results.\n\nI'm aware it might sound a bit unconventional, but experimentation is the essence of data science. I read a book titled [\"The Art of Feature Engineering\"](https://artoffeatureengineering.com/book.html), and it was an enlightening read that emphasized the importance of creatively engineering features to enhance model performance. I highly recommend it! 📚🔍\n\nI'm convinced that innovation often stems from trials and experiments. I'm eager to see if this approach of mine will provide valuable insights. It's always worth exploring. 🚀\n\n---","metadata":{}},{"cell_type":"code","source":"# Calculate imbalance metrics for each categorical column\nimbalance_metrics = {}\n\nfor column in categorical_columns:\n    counts = data[column].value_counts()\n    \n    # Ratio of most frequent to least frequent category\n    max_min_ratio = counts.max() / counts.min()\n    \n    # Coefficient of variation\n    coef_variation = counts.std() / counts.mean()\n    \n    imbalance_metrics[column] = {\n        \"Max/Min Ratio\": max_min_ratio,\n        \"Coefficient of Variation\": coef_variation\n    }\n\nimbalance_metrics_df = pd.DataFrame(imbalance_metrics).T\nimbalance_metrics_df\n\n# Set up the figure with subplots\nfig, axes = plt.subplots(2, 1, figsize=(22, 12))\n\n# Plot Max/Min Ratio\nsns.barplot(x=imbalance_metrics_df.index, y=imbalance_metrics_df[\"Max/Min Ratio\"], ax=axes[0], palette=\"viridis\")\naxes[0].set_title(\"Max/Min Ratio for Categorical Columns\", fontsize=30, fontweight=\"bold\")\naxes[0].set_ylabel(\"Max/Min Ratio\", fontsize=14, fontweight=\"bold\")\n# axes[0].set_xlabel(\"Categorical Columns\", fontsize=14,  fontweight=\"bold\")\naxes[0].tick_params(axis='x', labelsize=14)\naxes[0].tick_params(axis='y', labelsize=13)\naxes[0].grid(axis=\"y\", linestyle=\"--\")\n\n# Plot Coefficient of Variation\nsns.barplot(x=imbalance_metrics_df.index, y=imbalance_metrics_df[\"Coefficient of Variation\"], ax=axes[1], palette=\"viridis\")\naxes[1].set_title(\"Coefficient of Variation for Categorical Columns\", fontsize=30, fontweight=\"bold\")\naxes[1].set_ylabel(\"Coefficient of Variation\", fontsize=14, fontweight=\"bold\")\n# axes[1].set_xlabel(\"Categorical Columns\", fontsize=14, fontweight=\"bold\")\naxes[1].tick_params(axis='x', labelsize=14)\naxes[1].tick_params(axis='y', labelsize=13)\naxes[1].grid(axis=\"y\", linestyle=\"--\")\n\n# Adjust layout\nplt.tight_layout(pad=9.0)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:16.821561Z","iopub.execute_input":"2023-10-15T13:23:16.821857Z","iopub.status.idle":"2023-10-15T13:23:17.646719Z","shell.execute_reply.started":"2023-10-15T13:23:16.821832Z","shell.execute_reply":"2023-10-15T13:23:17.645278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\nIndeed, the `cell_type` column is significantly imbalanced compared to the others in the dataset. \n\nHowever, literature offers some good alternatives to address this kind of challenge. \n\nResampling: We can consider resampling techniques, such as oversampling (increasing the counts of minority categories) or undersampling (reducing the counts of majority categories), to balance out the dataset. 🔄\n\nWeighting Techniques: When it comes to training machine learning models, many algorithms provide the option to weight classes during training. This can help the model pay more attention to underrepresented categories. ⚖️\n\nTree-based Models: Often, tree-based models (like decision trees, random forests, gradient boosting) handle imbalanced datasets quite well. 🌲\n\nIn summary, while imbalance is a challenge, there are several strategies you can adopt to address it, depending on the analysis objective and the dataset's characteristics. We'll be discovering these together as we move forward. 🛤️\n\nI've already seen some solutions here that showcased tree-based models with impressive performance! Perhaps we're on the right track! 🚀\n\n---\n","metadata":{}},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 5px; border: 1px solid #0056b3; border-radius: 4px; color: white;\"><h3 align=center>Deepening the understanding about the categorical columns</h3></div>","metadata":{}},{"cell_type":"markdown","source":"--- \n\nWell, I've been slightly bothered by the fact that there seems to be an imbalance between the categorical columns. This could be a problem, but not necessarily. I believe that a model can learn to handle this, but I also think we can assist the model in managing it. \nHowever, before diving deeper, I want to investigate if there's any relationship between the categorical columns. I'm hopeful that they might prove useful down the road, especially during the modeling phase. I'm optimistic! 🌟\n\nA few nagging questions I'd like to address are:\n1) I've noticed that the columns `sm_name`, `sm_lincs_id`, and \"SMILES\" apparently have the same length. I wonder if their values are unique? I'm not sure! I'm not very familiar with biotechnology. Can we check this? 🤔\nIf this turns out to be true, I might be able to drop the `sm_name` and `sm_lincs_id` columns, as they won't provide any additional information. How wonderful that would be! 🎉\n\n---\n","metadata":{}},{"cell_type":"code","source":"# Grouping by 'sm_name' and counting unique values for 'sm_lincs_id' and 'SMILES'\nunique_value_counts = data.groupby('sm_name').agg({\n    'sm_lincs_id': pd.Series.nunique,\n    'SMILES': pd.Series.nunique\n})\n\n# Checking if all values are 1\nall_unique = (unique_value_counts == 1).all()\n\nimport matplotlib.pyplot as plt\n\n# Limiting the string length of the labels\nmax_label_length = 10\nshortened_labels = [label[:max_label_length] + \"...\" if len(label) > max_label_length else label for label in unique_value_counts.index]\n\n# Plotting the bar chart with adjusted labels\nfig, ax = plt.subplots(figsize=(20, 6))\nunique_value_counts.plot(kind='bar', ax=ax)\nax.set_title(\"Unique Value Count for '$sm\\_lincs\\_id$' and '$SMILES$' by '$sm\\_name$'\", fontsize=30, fontweight=\"bold\")\nax.set_ylabel(\"Unique Value Count\", fontsize=14, fontweight=\"bold\")\nax.set_xlabel(\"sm_name\", fontsize=14, fontweight=\"bold\")\nax.legend(title=\"Columns\")\nax.set_xticklabels(shortened_labels, rotation=90, ha='right')\nax.set_ylim(0, 1)\nax.set_yticks([0, 1])\n\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:17.648511Z","iopub.execute_input":"2023-10-15T13:23:17.648932Z","iopub.status.idle":"2023-10-15T13:23:19.994938Z","shell.execute_reply.started":"2023-10-15T13:23:17.648903Z","shell.execute_reply":"2023-10-15T13:23:19.993659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\n1) Folks, from what we observed here, it appears that the columns `sm_lincs_id` and `sm_name` are identical, so we can eliminate one of them. I'm leaning towards removing the `sm_name` column, as it doesn't seem to provide any additional information. Thoughts? 🤔\n2) Well, it feels like the tables have turned! 🔄 What I initially perceived as a complex situation now appears to be simpler, thanks to our earlier findings. Let's now focus on the relationship between the `cell_type` column and the columns `sm_name` (or `sm_lincs_id`, given they're interchangeable). Since we know they uniquely correspond to each other, it doesn't really matter which one we choose. 🧐\n\n---\n","metadata":{}},{"cell_type":"code","source":"# Counting the number of unique values in each column\nunique_cell_types = data['cell_type'].nunique()\nunique_sm_names = data['sm_name'].nunique()\n\nformatted_sentence = (\n    f\"The number of unique items in cell_type is {unique_cell_types} and \"\n    f\"the count in sm_name is {unique_sm_names}.\"\n)\n\nformatted_sentence\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:19.996470Z","iopub.execute_input":"2023-10-15T13:23:19.997522Z","iopub.status.idle":"2023-10-15T13:23:20.009136Z","shell.execute_reply.started":"2023-10-15T13:23:19.997475Z","shell.execute_reply":"2023-10-15T13:23:20.007539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Distribution of \"cell_type\" values for each unique \"sm_name\" value\ndistribution_sm_name = data.groupby('sm_name')['cell_type'].value_counts().unstack().fillna(0)\n\n# Sorting the distribution based on the total count of relations\nsorted_distribution = distribution_sm_name.sum(axis=1).sort_values(ascending=False)\nsorted_distribution_df = distribution_sm_name.loc[sorted_distribution.index]\n\n# Limiting the size of the x-axis labels\nlabels = [label[:15] + '...' if len(label) > 15 else label for label in sorted_distribution_df.index]\n\n# Plotting the graph\nplt.figure(figsize=(20, 6.7))\nsorted_distribution_df.plot(kind='bar', stacked=True, ax=plt.gca())\nplt.title('Distribution of $cell\\_type$ by $sm\\_name$', fontsize=30, fontweight='bold')\nplt.xlabel('')\nplt.ylabel('Count of Relations', fontsize=14, fontweight='bold')\nplt.xticks(ticks=range(len(labels)), labels=labels, rotation=90)\nplt.tight_layout()\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:20.011343Z","iopub.execute_input":"2023-10-15T13:23:20.011973Z","iopub.status.idle":"2023-10-15T13:23:23.279734Z","shell.execute_reply.started":"2023-10-15T13:23:20.011931Z","shell.execute_reply":"2023-10-15T13:23:23.278226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\nInteresting observation! 🤔 For each item in `sm_name`, we find that some have multiple associations, and we can pinpoint which `cell_type` is linked to them. This is quite intriguing! And, as we've noticed, `T regulatory`, `T cells CD4+`, and `Nk cells` are the cell types that most frequently appear in these associations. Could this be significant? Honestly, opinions on this might be divided. We'd need to weigh the cost of diving deeper into this against the potential outcome we might achieve in the model. 📊🔍\n\nLastly, a question that's been nagging me: Are the columns `sm_name` and `cell_type` truly independent? 🧐\n\n---","metadata":{}},{"cell_type":"markdown","source":"---\n\nFirst, I'd like to assess whether they're truly independent. For that, let's explore using the Chi-squared test. 📊🔍\n\n---\n\n","metadata":{}},{"cell_type":"code","source":"from scipy.stats import chi2_contingency\n\n# Criando uma tabela de contingência\ncontingency_table = pd.crosstab(data['sm_name'], data['cell_type'])\n\n# Aplicando o teste do Qui-quadrado\nchi2, p, _, _ = chi2_contingency(contingency_table)\n\nchi2, p","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.281959Z","iopub.execute_input":"2023-10-15T13:23:23.282736Z","iopub.status.idle":"2023-10-15T13:23:23.327603Z","shell.execute_reply.started":"2023-10-15T13:23:23.282691Z","shell.execute_reply":"2023-10-15T13:23:23.326386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n\nThe Chi-squared test result is $( \\chi^2 = 189.15 )$ and the associated p-value is $( p = 1.0 )$ 📊.\n\nA p-value of 1.0 suggests that we cannot reject the null hypothesis, implying that the two variables, `cell_type` and `sm_name`, are independent in this dataset 🔍.\n\nIn other words, the distribution of `cell_type` seems to be consistent regardless of the chosen `sm_name` 🔄. This aligns with our previous observations that each `sm_name` is associated with all cell types 🧬.\n\n---","metadata":{}},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 5px; border: 1px solid #0056b3; border-radius: 4px; color: white;\"><h3 align=center>Pre-Processing of Data</h3></div>","metadata":{}},{"cell_type":"markdown","source":"Alright, folks! 🌟 I think we've achieved a good understanding of the data. Now, we'll move on to the data preprocessing phase.\n\nHere's the plan:\n\n1. We'll apply one-hot encoding to the `cell_type` and `sm_name` columns, as these are the main columns of interest. On reflection, we probably didn't need to spend so much time deep diving into the categorical columns, since the test set (`id_map`) only contains `cell_type` and `sm_name`. We could have gone this route from the beginning. But that's okay, no worries! 😅\n\n2. Next, we'll identify the columns in `data_onehot` that aren't present in `id_map_onehot` and remove those uncommon columns.\n\nWith these steps, our datasets will be preprocessed and ready for modeling! 🚀 Let's dive right in! 🏊‍♂️\n","metadata":{}},{"cell_type":"code","source":"# Applying one-hot encoding directly to the 'cell_type' and 'sm_name' columns for 'data'\ndata_onehot_direct = pd.get_dummies(data, columns=['cell_type', 'sm_name'])\n\n# Applying one-hot encoding to the 'cell_type' and 'sm_name' columns for 'id_map'\nid_map_onehot = pd.get_dummies(id_map, columns=['cell_type', 'sm_name'])\n\n# Identifying columns in 'data_onehot_direct' that are not present in 'id_map_onehot'\nuncommon_data_to_id_map = [f for f in data_onehot_direct if f not in id_map_onehot]\n\n# Identifying columns in 'id_map_onehot' that are not present in 'data_onehot_direct'\nuncommon_id_map_to_data = [f for f in id_map_onehot if f not in data_onehot_direct]\n\n# Removing uncommon columns\ndata_onehot_direct = data_onehot_direct.drop(columns=uncommon_data_to_id_map)\nid_map_onehot = id_map_onehot.drop(columns=uncommon_id_map_to_data)\n\n# Checking if both sets now have the same set of columns\nsame_columns = data_onehot_direct.columns.isin(id_map_onehot.columns).all()\n\nsame_columns, data_onehot_direct.shape, id_map_onehot.shape","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.329384Z","iopub.execute_input":"2023-10-15T13:23:23.330119Z","iopub.status.idle":"2023-10-15T13:23:23.493738Z","shell.execute_reply.started":"2023-10-15T13:23:23.330080Z","shell.execute_reply":"2023-10-15T13:23:23.492643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div style=\"border:2px solid #e7e7e7; padding:15px; background-color:#ffffcc; border-radius:5px; color: black;\">\n\n**Firstly, a massive shoutout 🎉 to our fantastic friends [Mehran Kazeminia](https://www.kaggle.com/mehrankazeminia) and [Somayyeh Gholami](https://www.kaggle.com/somayyehgholami) for shedding some incredible light 💡 on this with their marvelous work [here](https://www.kaggle.com/code/mehrankazeminia/1-op2-eda-linearsvr-regressorchain).**\n\nLet's dive into the crux of the matter 🏊:\n\nIf the number of unique combinations of \"cell_type\" and \"sm_name\" is 614, and the dataset also has 614 rows, it implies that each \"cell_type\" and \"sm_name\" combo pops up in the data just once 🎯.\n\nThis shifts the perspective on our approaches:\n\n1. **Combination Approach** 🧩: Given that every \"cell_type\" and \"sm_name\" combo is unique, going the combination and one-hot encoding route would result in one encoded column for each row in the dataset. This might be a tad overboard, leaving us with an ultra-sparse dataset. 🍂\n\n2. **Direct One-hot Encoding Approach** 🚀: Encoding \"cell_type\" and \"sm_name\" individually might still hold its charm. It preserves the flavor of these columns on their own. Yet, if we're harnessing these features to predict another variable, the catch is there might not be a ton of variation for the model to pick up spicy patterns, given each combo makes a solo appearance. 🎸🎶\n\nIn the grand scheme of things, choosing an approach is a blend of the problem specifics, the data's nature, and a dash of trial and error. Happy modeling! 🌌🔭🪐\n\n</div>\n","metadata":{}},{"cell_type":"markdown","source":"# <div style=\"background-color: #007BFF; padding: 5px; border: 1px solid #0056b3; border-radius: 4px; color: white;\"><h3 align=center>MODELING</h3></div>","metadata":{}},{"cell_type":"code","source":"# Ensuring that one-hot encoded columns in data_onehot_direct are in integer format (1/0)\ndata_onehot_direct = data_onehot_direct.astype(int)\n\n# Selecting only numeric columns from the original 'data' dataframe\nnumeric_data = data.select_dtypes(include=['float64'])\n\n# Combining the numeric columns with the corrected one-hot encoded dataframe\nfinal_data = pd.concat([data_onehot_direct, numeric_data], axis=1)\n\nfinal_data.shape","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.495429Z","iopub.execute_input":"2023-10-15T13:23:23.496126Z","iopub.status.idle":"2023-10-15T13:23:23.589747Z","shell.execute_reply.started":"2023-10-15T13:23:23.496086Z","shell.execute_reply":"2023-10-15T13:23:23.588464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_data.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.595577Z","iopub.execute_input":"2023-10-15T13:23:23.596582Z","iopub.status.idle":"2023-10-15T13:23:23.625202Z","shell.execute_reply.started":"2023-10-15T13:23:23.596528Z","shell.execute_reply":"2023-10-15T13:23:23.623860Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import OneHotEncoder\n\n# Selecting relevant columns for encoding\nencoder = OneHotEncoder(drop='first')\nencoded_features = encoder.fit_transform(data[['cell_type', 'sm_name']])\nencoded_df = pd.DataFrame(encoded_features.toarray(), columns=encoder.get_feature_names_out(['cell_type', 'sm_name']))\n\n# Combining the encoded features with the original dataset\npreprocessed_data = pd.concat([data, encoded_df], axis=1)\npreprocessed_data.drop(columns=['cell_type', 'sm_name', 'SMILES', 'sm_lincs_id'], inplace=True)\n\n# Splitting the data into training and validation sets\nX = preprocessed_data.drop(columns=['control'] + list(data.columns[5:]))\ny = preprocessed_data[data.columns[5:]]\n\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)\n\nX_train.shape, X_val.shape\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.626979Z","iopub.execute_input":"2023-10-15T13:23:23.627450Z","iopub.status.idle":"2023-10-15T13:23:23.906134Z","shell.execute_reply.started":"2023-10-15T13:23:23.627418Z","shell.execute_reply":"2023-10-15T13:23:23.904597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder\n\n# Dropping the specified columns from the 'de_train' dataframe\ndata = data.drop(columns=['sm_lincs_id', 'SMILES', 'control'])\n\n# Initialize label encoders for each categorical column\nlabel_encoders = {}\n\n# Apply label encoding to 'cell_type' and 'sm_name' columns\nfor col in ['cell_type', 'sm_name']:\n    le = LabelEncoder()\n    data[col] = le.fit_transform(data[col])\n    label_encoders[col] = le\n\n# Display the first few rows of the updated dataframe\nde_train_head_encoded = data.head()\nde_train_head_encoded\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:23:23.907953Z","iopub.execute_input":"2023-10-15T13:23:23.908510Z","iopub.status.idle":"2023-10-15T13:23:24.011122Z","shell.execute_reply.started":"2023-10-15T13:23:23.908467Z","shell.execute_reply":"2023-10-15T13:23:24.009687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import cross_val_score\nimport gc\nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error\nimport numpy as np\n\ndata = de_train_head_encoded.copy()\n\n# Splitting the data into features and targets\nX = data.drop(columns=data.columns[2:])\ny = data[data.columns[2:]]\n\n# Splitting the data into training and validation sets\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)\n\n# Initialize the GradientBoostingRegressor wrapped inside MultiOutputRegressor\ngbr = GradientBoostingRegressor(n_estimators=100, random_state=42)\nmulti_gbr = MultiOutputRegressor(gbr, n_jobs=-1)\n\n# Train the model\nprint(\"Training the model...\")\nmulti_gbr.fit(X_train, y_train)\n\n# Clean up memory\ngc.collect()\n\n\n# Function to get the mean row-wise root mean squared error\ndef compute_mmrms(y_true, y_pred):\n    return np.mean(np.sqrt(np.mean(np.square(y_true - y_pred), axis=1)))\n\n# Custom scoring function for cross_val_score\ndef mmrms_scorer(estimator, X, y):\n    y_pred = estimator.predict(X)\n    return -compute_mmrms(y, y_pred)  # negative value because cross_val_score tries to maximize the score\n\n\n# Cross-validation\nprint(\"Cross-validating...\")\nscores = cross_val_score(multi_gbr, X, y, cv=5, scoring=mmrms_scorer, n_jobs=-1)\n\n# Clean up memory\ngc.collect()\n\nmean_cv_score = -np.mean(scores)  # taking negative because we negated the score in mmrms_scorer\nmean_cv_score\n","metadata":{"execution":{"iopub.status.busy":"2023-10-15T13:25:23.986006Z","iopub.execute_input":"2023-10-15T13:25:23.986458Z","iopub.status.idle":"2023-10-15T14:17:41.076211Z","shell.execute_reply.started":"2023-10-15T13:25:23.986420Z","shell.execute_reply":"2023-10-15T14:17:41.074607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To be continued...","metadata":{}}]}