Multi-Layer Perceptron (MLP) - a component of the 13th place solution
\n\nThis notebook showcases the basic variation of a Multi-Layer Perceptron (MLP), which is a component of the 13th place solution (0.558 public and 0.745 private score).\n\nKey things employed:\n- Target-encoded features (PCA of -log10(p-value) * sign(LFC)) values for the genes and averaging by drug and cell type)\n- Use of raw, unprocessed targets\n- 5-fold cross-validation, involving prediction within each fold and averaging for the submission\n\nThe model's variations that improved the LB score and were included in the final submission (see the `Ensemble MLPs` section) are as follows:\n- Diverse approaches to train set augmentation\n- Prediction of separate gene groups (each with it's own features set)\n- Different numbers of principal components for feature generation\n- Using additional features from raw expression counts for all or specific genes (highly variable, markers of gene clusters identified [here](https://www.kaggle.com/code/antoninadolgorukova/op2-adata-analysis-with-seurat), stored in [this dataset](https://www.kaggle.com/datasets/antoninadolgorukova/op2-supplementary-calcs-for-ml))\n- Application of varying levels of noise to features related to drugs and cell types\n- Application of different filters for de_train (excluding observations from a limited number of cells, control compounds, and outliers)\n- Different activation functions (rrelu or relu)\n\n**Public LB scores 0.579 - 0.583 \nPrivat LB scores 0.766 - 0.802**\n\n- Training without validation (all train data) - up to 61 fold and with exclusion of 2-3 samples (16 folds) \n\n**Public LB scores 0.569 - 0.574 \nPrivat LB scores 0.761 - 0.769**\n","metadata":{}},{"cell_type":"markdown","source":"****","metadata":{}},{"cell_type":"markdown","source":"**The solution included:**\n- step1: LB score 0.575 = 0.5 Pyboost (0.584) + 0.5 Catboost (0.584)\n- step2: LB score 0.566 = 0.5 step1 (0.575 ) + 0.5 NN-NLP (0.574)\n- step3: LB score 0.559 = 0.5 step2 (0.566 ) + **0.5 NN-TargetEncEnsemble (= ensemble of MLPs, described in this notebook) (0.566)**\n- final polishing: 0.558 = 0.8 step 3 + 0.2 (more Pyboost (0.574, 0.577) and NN (0.569, 0.570, 0.572 - MLP, 0.587 - another NN) models)\n\nThoughts and details are [here](https://www.kaggle.com/code/alexandervc/op2-u900-team-blend)","metadata":{}},{"cell_type":"markdown","source":"****","metadata":{}},{"cell_type":"code","source":"library(arrow)\nlibrary(data.table)\nlibrary(qs)\nlibrary(tictoc)\nlibrary(ggplot2)\nlibrary(patchwork)\nlibrary(torch)\nlibrary(luz)\nlibrary(pheatmap)\n\noptions(scipen=999)\nfig <- function(width, heigth){\n options(repr.plot.width = width, repr.plot.height = heigth)\n}\nshape <- function(dt) {\n cat(\"\\nShape: \", ncol(dt), \" columns x \",\n nrow(dt), \" rows\", sep = \"\")\n}","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-12-09T06:33:28.490818Z","iopub.execute_input":"2023-12-09T06:33:28.493739Z","iopub.status.idle":"2023-12-09T06:33:34.775913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Load and prepare data
","metadata":{}},{"cell_type":"markdown","source":"Train data - contain the annotated cell type, some information about the drug, and differential expression values (-log10(p-value) * sign(LFC)) for 18211 genes.","metadata":{}},{"cell_type":"code","source":"de_train <- read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\nde_train <- setDT(de_train)\nhead(as.matrix(de_train[, 1:10]))\nshape(de_train)","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:34.782074Z","iopub.execute_input":"2023-12-09T06:33:34.905083Z","iopub.status.idle":"2023-12-09T06:33:38.297754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A sample submission file in the correct format.","metadata":{}},{"cell_type":"code","source":"submission <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv\")\nhead(as.matrix(submission[, 1:10]))\nshape(submission)","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:38.303202Z","iopub.execute_input":"2023-12-09T06:33:38.306311Z","iopub.status.idle":"2023-12-09T06:33:40.297887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The map of the cell_type / sm_name pair to be predicted for the given id","metadata":{}},{"cell_type":"code","source":"id_map <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\nhead(id_map)\nshape(id_map)","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.303909Z","iopub.execute_input":"2023-12-09T06:33:40.316838Z","iopub.status.idle":"2023-12-09T06:33:40.436672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Exclude from de_train samples with data from a small number of cells**","metadata":{}},{"cell_type":"markdown","source":"The problem with such samples was discussed here: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/445883 \nThe EDA notebook is here: https://www.kaggle.com/code/antoninadolgorukova/op2-stats-sign-de-fingerprints-eda\n\nIn some versions, I removed those samples, and the public score was not changed, however, according to private LB, it significantly worsened the performance.","metadata":{}},{"cell_type":"code","source":"# adata_obs_meta <- \n# fread(\"/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv\")\n# sum_by_drug_and_cell <- adata_obs_meta[, .(n_cells = .N), by = c(\"sm_name\", \"cell_type\")]\n# sum_by_drug_and_cell <- sum_by_drug_and_cell[sm_name != \"Dimethyl Sulfoxide\"][order(-n_cells)]\n\n# sum_by_drug_and_cell[n_cells <= 2]\n# include <- sum_by_drug_and_cell[n_cells > 2]\n# cat(\"Excluded:\", nrow(sum_by_drug_and_cell) - nrow(include), \"samples\")\n\n# de_train <- de_train[include[, .(sm_name, cell_type)], on = c(\"sm_name\", \"cell_type\") ]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.451808Z","iopub.execute_input":"2023-12-09T06:33:40.455193Z","iopub.status.idle":"2023-12-09T06:33:40.497935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Some functions
","metadata":{}},{"cell_type":"markdown","source":"make_features function makes PCA and target encoding (TE) by drug and cell type from any data containing sm_name and cell_type columns. Indicate from what columns make TE with features_cols argument (all numeric columns is the default, but you can make separate features for groups of genes). With n_pc argument, set the number of principal components from which to make averaging by drug and cell type. E.g. with n_pc = c(100, 50) you will get 100 drug features and 50 cell type features.","metadata":{}},{"cell_type":"code","source":"make_features <- function(data,\n features_cols = NULL,\n n_pc = rep(length(features_cols), 2),\n dataset) {\n \n if (is.null(features_cols)) {\n \n features_cols <- names(data)[unlist(lapply(data, is.numeric))]\n }\n \n # ensure rows match de_train\n data <- data[de_train[, .(sm_name, cell_type)],\n on = c(\"sm_name\", \"cell_type\"),\n nomatch = 0]\n \n # make PCA of the selected columns\n pca_targets_res <- prcomp(data[, ..features_cols],\n scale. = TRUE,\n center = TRUE)\n targets_encoded <- cbind(data[, .(sm_name, cell_type)],\n pca_targets_res$x)\n \n # average by drug and cell type\n target_encoded_drugs <- \n targets_encoded[, lapply(.SD, mean),\n by = sm_name,\n .SDcols = names(targets_encoded)[-c(1:2)]]\n target_encoded_ct <- \n targets_encoded[, lapply(.SD, mean),\n by = cell_type,\n .SDcols = names(targets_encoded)[-c(1:2)]] \n \n # make features (concat drug and cell type TE)\n if(dataset == \"train\") {\n \n data <- de_train\n } else {\n data <- id_map\n } \n\n if(n_pc[2] != 0) {\n \n features <- cbind(\n target_encoded_drugs[data[, sm_name], on = \"sm_name\"][\n , 2:(n_pc[1]+1), with = FALSE],\n target_encoded_ct[data[, cell_type], on = \"cell_type\"][\n , 2:(n_pc[2]+1), with = FALSE]\n )\n } else {\n \n features <- \n target_encoded_drugs[data[, sm_name], on = \"sm_name\"][, 2:(n_pc[1]+1)]\n }\n \n return(features)\n}","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.503082Z","iopub.execute_input":"2023-12-09T06:33:40.505999Z","iopub.status.idle":"2023-12-09T06:33:40.530264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Function to get metrics from training.","metadata":{}},{"cell_type":"code","source":"train_history <- function(fitted,\n .group = gr,\n use_valid) {\n \n train <- rbindlist(lapply(1:length(fitted$records$metrics$train), function(lst) {\n data.table(\n epoch = lst,\n set = \"train\",\n loss = fitted$records$metrics$train[[lst]]$loss,\n rmse = fitted$records$metrics$train[[lst]]$rmse,\n mae = fitted$records$metrics$train[[lst]]$mae\n )\n }))\n if(use_valid) {\n \n valid <- rbindlist(lapply(1:length(fitted$records$metrics$valid), function(lst) {\n data.table(\n epoch = lst,\n set = \"valid\",\n loss = fitted$records$metrics$valid[[lst]]$loss,\n rmse = fitted$records$metrics$valid[[lst]]$rmse,\n mae = fitted$records$metrics$valid[[lst]]$mae\n )\n }))\n train_history <- rbind(train, valid)\n\n } else {\n train_history <- train\n }\n \n metrics <- rbind(train_history[set == \"train\"][.N],\n train_history[set == \"valid\"][.N],\n fill = TRUE)\n metrics[, `:=` (group = .group, fold = f, i = i)]\n \n print(metrics)\n return(metrics) \n}","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.535145Z","iopub.execute_input":"2023-12-09T06:33:40.537745Z","iopub.status.idle":"2023-12-09T06:33:40.562738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Calculate metrics for each fold/entire Yoof with the option to show plots.","metadata":{}},{"cell_type":"code","source":"y_oof_metrics <- function(pred_cols,\n Y_oof,\n metrics_dt, \n show_plots) {\n \n true_targets <- de_train[, ..pred_cols][, id := 1:.N]\n true_targets <- true_targets[Y_oof[, .(id)], on = \"id\"]\n \n preds <- as.matrix(Y_oof[, ..pred_cols])\n true_targets <- as.matrix(true_targets[, !\"id\"])\n \n all_cors <- data.table(row = 1:nrow(preds))\n all_cors[, corr := sapply(row, function(r) {\n cor(unlist(preds[r, ]),\n unlist(true_targets[r, ])\n )\n })]\n \n mean_corr_rows <- round(mean(all_cors$corr), 3)\n median_corr_rows <- round(median(all_cors$corr), 3)\n \n p1 <- ggplot(all_cors, aes(x = corr)) + \n geom_histogram(bins = 20, color = \"black\", fill = \"lightgray\") +\n ggtitle(paste(\"row-wise\\nmean cor =\",\n mean_corr_rows,\n \"\\nmedian cor =\",\n median_corr_rows)) +\n theme_bw() \n \n all_cors <- data.table(col = 1:ncol(preds))\n all_cors[, corr := sapply(col, function(cl) {\n cor(unlist(preds[, cl]),\n unlist(true_targets[, cl])\n )\n })]\n \n mean_corr_cols <- round(mean(all_cors$corr), 3)\n median_corr_cols <- round(median(all_cors$corr), 3)\n \n p2 <- ggplot(all_cors, aes(x = corr)) + \n geom_histogram(bins = 20, color = \"black\", fill = \"lightgray\") +\n ggtitle(paste(\"column-wise\\nmean cor =\",\n mean_corr_cols,\n \"\\nmedian cor =\",\n median_corr_cols)) +\n theme_bw() \n \n if (!exists(\"tm\")) {\n \n tm <- list()\n tm[[4]] <- NA\n }\n \n metrics_dt <- rbind(\n metrics_dt,\n data.table(\n i = i, fold = f,\n cv = cv_scheme,\n features = feature_variation,\n augm = augm,\n noise = paste0(sigma_dr, \"/\", sigma_ct),\n n_PC = paste(n_pc, collapse = \", \"),\n n_train = nrow(cv),\n train_prop = round(1 - mean(cv[, .N, by = folds]$N)/length(train_ids), 3),\n n_groups = n_groups,\n n_Y_oof = nrow(Y_oof),\n d_hid1 = d_hid1,\n wd = wd,\n bs = bs,\n max_lr = maxlr,\n pt = pt,\n mrRMSE = mean(sqrt(rowMeans((true_targets - preds)^2))),\n RMSE = sqrt(mean((true_targets - preds)^2)),\n MAE = mean(abs((true_targets - preds))),\n r2 = 1 - (sum((true_targets - preds)^2) /\n sum((true_targets - mean(preds))^2)),\n use_valid = use_valid,\n mean_num_epochs = ifelse(exists(\"i\"), all_metrics[i == i, epoch],\n round(all_metrics[, mean(epoch)], 0)\n ),\n mean_corr_rows = mean_corr_rows,\n mean_corr_cols = mean_corr_cols,\n median_corr_rows = median_corr_rows,\n median_corr_cols = median_corr_cols,\n `time, s` = round(as.numeric(gsub(\" sec elapsed\", \"\", tm[[4]]) ) , 0)),\n fill = TRUE\n )\n \n cols <- c(\"mrRMSE\", \"RMSE\", \"MAE\", \"r2\")\n metrics_dt[ , (cols) := lapply(.SD, round, 3), .SDcols = cols]\n \n if (show_plots) {\n \n fig(20, 10)\n print(p1 + p2 + plot_annotation(\n title = \"Mean/median corrs Y_oof preds vs true_targets\") &\n theme(text = element_text(size = 18) )\n )\n } \n \n return(metrics_dt)\n}","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.568205Z","iopub.execute_input":"2023-12-09T06:33:40.571247Z","iopub.status.idle":"2023-12-09T06:33:40.601617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Goup targets
","metadata":{}},{"cell_type":"markdown","source":"Set number of groups with n_groups. You can group the genes by random or by clusters identified with kmeans, or add your own way of groupping. With n_groups = 1 and groups_type = \"rand\" the model will predict all genes.","metadata":{}},{"cell_type":"code","source":"all_targets <- names(de_train)[unlist(lapply(de_train, is.numeric))]\n\nn_groups = 1\ngroups_type = \"rand\"\n\nif (groups_type == \"rand\") {\n \n set.seed(1)\n target_groups <- data.table(\n gene = all_targets,\n group = sample(1:n_groups,\n length(all_targets), replace = TRUE)\n )\n \n} else if (groups_type == \"kmeans\") {\n \n scaled_targ <- t(scale(de_train[, ..all_targets]))\n \n # Uncomment to check the optimal number of clusters\n # fviz_nbclust(scaled_targ, kmeans, method = \"wss\") +\n # geom_vline(xintercept = 3, linetype = 2)\n \n km_targ <- kmeans(scaled_targ, n_groups, nstart = 25)\n # head(km_targ$cluster, 3)\n # km_targ$size\n # fviz_cluster(km_targ, scaled_targ, ellipse.type = \"norm\") #repel = TRUE\n \n target_groups <- data.table(\n gene = names(km_targ$cluster),\n group = km_targ$cluster\n )\n \n}\n\ntarget_groups[, .N, by = \"group\"][\n order(group)]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.611071Z","iopub.execute_input":"2023-12-09T06:33:40.617262Z","iopub.status.idle":"2023-12-09T06:33:40.746826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Multilayer perceptron (MLP) neural network
","metadata":{}},{"cell_type":"code","source":"fit_predict <- function(\n pred_cols,\n n_pc,\n use_valid,\n sep_features) {\n \n cat(\"\\n===================== Fold\", f, \"=====================\\n\")\n \n # augment train by n times \n \n cat(\"\\nInitial train:\", length(train_ids), \"\\n\") \n \n train_ids <- rep(train_ids, times = augm)\n \n cat(\"\\nTrain features set:\", length(train_ids))\n cat(\"\\nValid features set:\", length(valid_ids), \"\\n\")\n\n pred_groups <- target_groups[gene %chin% pred_cols]\n \n # predict groups of genes\n \n tic()\n for (gr in sort(pred_groups[, unique(group)]) ) {\n \n cat(\"\\n============ Group\", gr, \"============\\n\")\n \n targets_cols <- pred_groups[group == gr, gene]\n \n train_targets <- de_train[, ..targets_cols]\n \n if(sep_features) {\n \n train_features <- make_features(data = de_train,\n features_cols = targets_cols,\n n_pc = n_pc,\n dataset = \"train\")\n } else {\n \n train_features <- make_features(\n data = de_train,\n features_cols = names(de_train)[unlist(lapply(de_train, is.numeric))],\n n_pc = n_pc,\n dataset = \"train\")\n }\n \n n_features <- ncol(train_features)\n n_targets <- ncol(train_targets)\n \n cat(\"\\nn_features = \", n_features)\n cat(\"\\nn_targets = \", n_targets, \"\\n\\n\")\n\n # add noise to the augmented train \n \n if (n_pc[2] == 0) {\n \n .n_pc <- n_pc\n .n_pc[2] <- ncol(train_features) - n_pc[1]\n \n } else {\n .n_pc <- n_pc\n }\n \n rnd_f <- torch_cat(list(\n torch_normal(c(length(train_ids), .n_pc[1]), mean = 0, std = sigma_dr),\n torch_normal(c(length(train_ids), .n_pc[2]), mean = 0, std = sigma_ct)\n ), dim = 2 )\n \n rnd_t <- torch_normal(c(length(train_ids), n_targets), mean = 0, std = sigma_t)\n \n train_ds <- tensor_dataset(\n torch_tensor(as.matrix(train_features[train_ids, ])) + rnd_f,\n torch_tensor(as.matrix(train_targets[train_ids])) + rnd_t\n )\n \n train_dl <- dataloader(train_ds, batch_size = bs, shuffle = TRUE)\n \n if (length(valid_ids) > 0 ) {\n \n valid_ds <- tensor_dataset(\n torch_tensor(as.matrix(train_features[valid_ids, ])),\n torch_tensor(as.matrix(train_targets[valid_ids, ]))\n )\n valid_dl <- dataloader(valid_ds, batch_size = bs)\n }\n \n net <- nn_module(\n initialize = function(n_features, n_targets, d_hid1) {\n self$net <- nn_sequential(\n nn_linear(n_features, d_hid1),\n nn_relu(), \n nn_linear(d_hid1, n_targets)\n )\n },\n forward = function(x) {\n self$net(x)\n }\n )\n \n net <- net %>%\n setup(\n loss = nn_l1_loss(), # MAE loss\n optimizer = optim_adamw,\n metrics = list(luz_metric_rmse(),\n luz_metric_mae())\n ) %>%\n set_hparams(\n n_features = n_features,\n n_targets = n_targets, \n d_hid1 = d_hid1\n ) %>%\n set_opt_hparams(weight_decay = wd)\n \n if (use_valid) {\n \n fitted <- net %>%\n fit(train_dl, \n epochs = num_epochs, \n valid_data = valid_dl,\n callbacks = list(\n luz_callback_lr_scheduler(\n lr_one_cycle,\n max_lr = maxlr,\n epochs = num_epochs,\n steps_per_epoch = length(train_dl),\n call_on = \"on_epoch_end\" \n ),\n luz_callback_early_stopping(\n monitor = \"valid_loss\",\n patience = pt)\n ),\n verbose = FALSE\n )\n } else {\n \n fitted <- net %>%\n fit(train_dl,\n epochs = num_epochs,\n verbose = FALSE)\n }\n \n # predict unseen data for validation\n \n if (length(valid_ids) > 0 ) {\n \n test_f <- as.matrix(train_features[valid_ids, ])\n\n test_ds <- tensor_dataset(\n torch_tensor(test_f)\n )\n test_dl <- dataloader(test_ds, batch_size = bs)\n\n preds <- predict(fitted, test_dl)\n preds <- as.matrix(preds$to(device = \"cpu\")) \n preds <- data.table(preds)[, id := valid_ids] \n setnames(preds, names(preds), c(targets_cols, \"id\")) \n preds <- melt(preds, id.vars = \"id\") \n\n iter_Y_oof[gene %chin% targets_cols, pred := preds[, value]]\n }\n \n # predict for submit\n \n if (sep_features) {\n \n test_features <- make_features(data = de_train,\n features_cols = targets_cols,\n n_pc = n_pc,\n dataset = \"test\")\n \n } else {\n \n test_features <- make_features(\n data = de_train,\n features_cols = names(de_train)[unlist(lapply(de_train, is.numeric))],\n n_pc = n_pc,\n dataset = \"test\")\n }\n \n test_features <- as.matrix(test_features)\n cat(\"\\nTest features: \", nrow(test_features), \"\\n\\n\" )\n\n submit_ds <- tensor_dataset(torch_tensor(test_features))\n submit_dl <- dataloader(submit_ds, batch_size = bs)\n\n preds <- predict(fitted, submit_dl)\n preds <- as.matrix(preds$to(device = \"cpu\")) \n preds <- data.table(preds)[, id := 1:nrow(test_features)] \n setnames(preds, names(preds), c(targets_cols, \"id\")) \n preds <- melt(preds, id.vars = \"id\") \n\n iter_preds_fold[gene %chin% targets_cols, pred := preds[, value]]\n\n # Get training history\n\n if (!exists(\"metrics\")) {\n metrics <- NULL\n }\n \n metrics <- rbind(metrics, \n train_history(\n fitted,\n .group = gr,\n use_valid)\n )\n }\n tm <- toc() \n return(metrics) \n}","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.755859Z","iopub.execute_input":"2023-12-09T06:33:40.762725Z","iopub.status.idle":"2023-12-09T06:33:40.798017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Set parameters** ","metadata":{}},{"cell_type":"code","source":"pred_cols = names(submission)[-1]\n\nnum_epochs = 100\nuse_valid = TRUE # use validation set for early stopping and metrics calculation?\n\nmaxlr = 0.0001\nbs = 256\nd_hid1 = 1024\nwd = 0.5\n\npt = 3\nsigma_dr = 1 # noise to drugs features\nsigma_ct = 0.1 # noise to cell type features\nsigma_t = 0 #noise to targets\naugm = 50 # how many times repeat each train sample\n\nn_pc = c(100, 100) # number of PC to make TE features \nnComp = 0\n\nsep_features = TRUE # make separate features foreach gene group?\nnum_iter = 1 # how many times train and predict on each fold?\n\nfeature_variation = paste0(\n \"targ TE\",\n ifelse(sep_features, \" sep\", \"\"),\n ifelse(n_groups > 1, paste(\";\", groups_type), \"\")\n)","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.806707Z","iopub.execute_input":"2023-12-09T06:33:40.811396Z","iopub.status.idle":"2023-12-09T06:33:40.897806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#CV scheme
","metadata":{}},{"cell_type":"markdown","source":"I used 5-fold cross validation on test drugs. This resulted in each fold containing all available data from the test cells (myeloid and B).","metadata":{}},{"cell_type":"code","source":"n_folds = 5\ncv_scheme = \"on test drugs\"\n\nset.seed(2)\ntest_drugs <- unique(id_map[, sm_name])\n\ncv <- de_train[, .(sm_name, cell_type, id = 1:.N)]\ncv <- cv[sample(1:nrow(cv))]\n\ncv[sm_name %in% test_drugs, folds :=\n sample(1:n_folds, .N, replace = TRUE)]\ncv[!sm_name %in% test_drugs, folds := n_folds+1]\n\ncat(\"Number and proportion of samples in each fold\")\ncv[folds %in% 1:n_folds,\n .(n =.N, prop = round(.N/nrow(cv), 2) ),\n by = \"folds\"][order(folds)]\n\ncat(\"\\n\\nNumber of samples of each cell type in each fold. Fold 6 is always used for training only.\")\ncv[, .N, by = c(\"folds\", \"cell_type\")][order(folds)]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:40.906922Z","iopub.execute_input":"2023-12-09T06:33:40.913451Z","iopub.status.idle":"2023-12-09T06:33:41.129294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Train & predict
","metadata":{}},{"cell_type":"code","source":"preds_cv <- melt(submission,\n variable.factor = FALSE,\n id.vars = \"id\",\n variable.name = \"gene\",\n value.name = \"pred\")\n\ntic()\nY_oof <- NULL\nall_metrics <- NULL\niter_metrics <- NULL\n\nfor (f in 1:n_folds) { \n \n # Split \n train_ids <- cv[folds != f, id]\n valid_ids <- cv[folds == f, id]\n \n # Make dt to gather fold's predictions for Yoof and submit \n preds_fold <- copy(preds_cv)\n preds_fold[, pred := 0]\n \n Y_oof_fold <- melt(\n de_train[valid_ids, ..pred_cols][, id := valid_ids],\n variable.factor = FALSE,\n id.vars = \"id\",\n variable.name = \"gene\",\n value.name = \"pred\")\n Y_oof_fold[, pred := 0]\n \n for(i in 1:num_iter) {\n \n # cat(\"\\nIteration:\", i, \"\\n\") \n iter_Y_oof <- copy(Y_oof_fold)\n iter_Y_oof[, pred := 0]\n \n iter_preds_fold <- copy(preds_fold)\n iter_preds_fold[, pred := 0]\n \n all_metrics <- rbind(all_metrics,\n fit_predict(pred_cols,\n n_pc,\n use_valid,\n sep_features)\n )\n \n # Average preds from iterations \n preds_fold[, pred := pred + iter_preds_fold[, pred]/num_iter]\n Y_oof_fold[, pred := pred + iter_Y_oof[, pred]/num_iter]\n \n # Calculate metrics for current fold \n if (length(valid_ids) > 0 ) {\n check_Y_oof_fold <- dcast(Y_oof_fold,\n id ~ gene,\n value.var = \"pred\")\n \n iter_metrics <- y_oof_metrics(pred_cols,\n check_Y_oof_fold,\n iter_metrics, \n show_plots = FALSE)\n }\n }\n \n if (length(valid_ids) > 0 ) {\n Y_oof_fold <- dcast(Y_oof_fold,\n id ~ gene,\n value.var = \"pred\")\n \n Y_oof <- rbind(Y_oof, Y_oof_fold)\n }\n \n # Average predicts from folds \n preds_fold <- preds_fold[, pred := mean(pred), by = c(\"id\", \"gene\")]\n preds_cv[, pred := pred + preds_fold[, pred] / n_folds]\n \n}\ntm <- toc()\n\nall_metrics <- all_metrics[order(fold, group)]\noverfits <- all_metrics[, lapply(.SD, diff),\n by = c(\"fold\", \"group\"),\n .SD = c(\"loss\", \"rmse\", \"mae\")]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:33:41.137207Z","iopub.execute_input":"2023-12-09T06:33:41.142841Z","iopub.status.idle":"2023-12-09T06:39:17.410259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Metrics for each fold**","metadata":{}},{"cell_type":"code","source":"show_cols = c(\n \"fold\", \"i\", \"mean_num_epochs\", \"mrRMSE\", \"RMSE\", \"MAE\", \"r2\",\n \"mean_corr_rows\", \"mean_corr_rows\", \"median_corr_cols\", \"median_corr_cols\"\n )\niter_metrics[, ..show_cols]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:39:17.415683Z","iopub.execute_input":"2023-12-09T06:39:17.419192Z","iopub.status.idle":"2023-12-09T06:39:17.463223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Metrics for the Y oof**","metadata":{}},{"cell_type":"code","source":"if (!exists(\"all_train_metrics\")) {\n all_train_metrics <- NULL\n}\nall_train_metrics <- y_oof_metrics(pred_cols, Y_oof, all_train_metrics, \n show_plots = TRUE)\n\nall_train_metrics[, !c(\"i\", \"fold\", \"cv\", \"n_train\", \"train_prop\", #\"n_PC\", \n \"d_hid1\", \"wd\", \"bs\", \"max_lr\", \"pt\", #\"n_groups\",\n \"n_Y_oof\", \"noise\", \"time, s\", \"use_valid\")\n]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:39:17.468690Z","iopub.execute_input":"2023-12-09T06:39:17.470918Z","iopub.status.idle":"2023-12-09T06:39:22.918151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_train_metrics[, .(d_hid1, wd, bs, max_lr, pt)]","metadata":{"execution":{"iopub.status.busy":"2023-12-09T06:44:34.010855Z","iopub.execute_input":"2023-12-09T06:44:34.013087Z","iopub.status.idle":"2023-12-09T06:44:34.044065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Prepare submission file
","metadata":{}},{"cell_type":"code","source":"submit <- dcast(preds_cv,\n id ~ gene,\n value.var = \"pred\")\nhead(as.matrix(submit))\ncat(\"\\nShape: \", ncol(submit), \" columns x \",\n nrow(submit), \" rows\", sep = \"\")","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:26:39.302667Z","iopub.execute_input":"2023-12-07T11:26:39.304781Z","iopub.status.idle":"2023-12-07T11:26:40.693614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_id = \"example_MLP.csv\"\nfwrite(submit, paste0(\"submission_\", file_id))\nfwrite(Y_oof, paste0(\"Y_oof_\", file_id))\nfwrite(all_train_metrics, paste0(\"all_train_metrics_\", file_id))","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:26:40.696716Z","iopub.execute_input":"2023-12-07T11:26:40.698453Z","iopub.status.idle":"2023-12-07T11:26:42.668470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Ensemble MLPs
","metadata":{}},{"cell_type":"markdown","source":"For the final submission I prepared an ensemble of several MLP variations. ","metadata":{}},{"cell_type":"code","source":"# Function to calculate Pearson correlation of the models by row\ncalc_cor_ave_and_plot <- function(list_df) {\n \n correlation_averages <- \n matrix(NA, nrow = length(list_df), ncol = length(list_df))\n rownames(correlation_averages) <- names(list_df)\n colnames(correlation_averages) <- names(list_df)\n all_cors <- NULL\n\n for (i in seq_along(list_df)) {\n for (j in seq_along(list_df)) {\n\n pair_cor <- data.table(id = 1:nrow(list_df[[i]]),\n submit1 = names(list_df)[i],\n submit2 = names(list_df)[j])\n for (r in pair_cor$id) {\n\n cor <- cor(unlist(list_df[[i]][r, -1]),\n unlist(list_df[[j]][r, -1])\n )\n pair_cor[r, corr := cor] \n }\n correlation_averages[i, j] <- round(mean(pair_cor$corr), 3)\n if (i < j) {\n all_cors <- rbind(all_cors, pair_cor)\n }\n }\n } \n\n plt <- pheatmap(correlation_averages, \n fontsize = 16, fontsize_number = 20,\n col = colorRampPalette(c(\"white\",\"darkred\"))(200),\n main = \"Correlations\",\n display_numbers = TRUE,\n number_color = \"black\"\n ) \n return(all_cors)\n \n}","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:26:42.671585Z","iopub.execute_input":"2023-12-07T11:26:42.673219Z","iopub.status.idle":"2023-12-07T11:26:42.688429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# - Description of the included MLPs","metadata":{}},{"cell_type":"code","source":"path <- '/kaggle/input/op2-submits-and-yoof/mlp submits/'\nfile_names <- list.files(path)[grepl('.csv', list.files(path))]\n#file_names <- file_names[!grepl('0.57$|pseudo', file_names)]\n#file_names <- file_names[!grepl('ave_blend_T8_Treg_3kmeans_0.57', file_names)]\nlength(file_names)\ncat(paste0(\"'\", file_names, \"'\", collapse = \", \\n\"))","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:26:42.691374Z","iopub.execute_input":"2023-12-07T11:26:42.692991Z","iopub.status.idle":"2023-12-07T11:26:42.737823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The file blend_9diverse_models_0.573 contains an ensemble of the 9 following 5-fold models:","metadata":{}},{"cell_type":"markdown","source":"| Version | Model description | n samples | n PC for features | RMSE | mrRMSE | R2 | overfit loss | corr_rows | corr_cols | Public LB | Privat LB |\n|---------|-------------------|-----------|-------------------|------|--------|----|--------------|--------------------------------|--------------------------------|-----------|-----------|\n| MLPv35 | TE +augm 100 s0.1/0.1 | 614 | 100/100 | 1.218526 | 0.8300539 | 0.3699764 | | | | 0.582 | 0.767 |\n| MLPv48 | TE +augm 50 s0.5/0.5, n_cells >2 pairs | 603 | 100/100 | 1.175822 | 0.7874235 | 0.3062406 | 0.160 | | | 0.58 | 0.766 |\n| MLPv170 | TE +augm 150 s1/1, n_cells >2 pairs | 603 | 100/100 | 1.132636 | 0.7809319 | 0.3562856 | 0.151 | 0.472 | 0.528 | 0.579 | 0.786 |\n| MLPv121 | TE +augm 50 s1/0, n_cells >2 pairs before PCA | 603 | 100/100 | 1.10915 | 0.7737568 | 0.3826959 | 0.144 | 0.47 | 0.542 | 0.58 | 0.784 |\n| MLPv122 | TE +augm 50 s1/0.1, n_cells >5 pairs before PCA | 597 | 100/100 | 1.061187 | 0.7647868 | 0.4185533 | 0.139 | 0.466 | 0.559 | 0.581 | 0.802 |\n| MLPv107 | TE +augm 50 s1/1, n_cells >2 pairs | 603 | 140/100 | 1.208789 | 0.79697 | 0.2668225 | 0.126 | 0.46 | 0.478 | 0.581 | 0.771 |\n| MLPv128 | TE +augm 30 s1/0.1, n_cells >2 pairs before PCA + by row 10%/0.5 | 603 | 100/100 | 1.166521 | 0.7865393 | 0.3171834 | 0.151 | 0.464 | 0.505 | 0.582 | 0.788 |\n| MLPv151 | TE 80/80 + scaled_counts_features TE 20/20, +augm 50 s1/0, n_cells >2 pairs before PCA | 603 | 80/80/20/20 | 1.138315 | 0.7730064 | 0.3497886 | 0.13 | 0.471 | 0.53 | 0.583 | 0.793 |\n| MLPv86 | TE +augm 50 s1/1, n_cells >2 pairs | 603 | 100/100 | 1.220697 | 0.7944542 | 0.2522711 | 0.110 | 0.476 | 0.49 | 0.582 | 0.766 |","metadata":{}},{"cell_type":"markdown","source":"**Ensemble scores: \npublic: 0.573 \nprivat: 0.772**\n\nIn all versions, I did use the model parameters and 5-fold cross-validation just the same as in this notebook (3L; d_hid1 = 256; wd = 0.5; bs = 128; maxlr = 0.01; pt = 3; call_on = \"on_epoch_end\", validation on test drugs)\n\n- TE - target encoded features as above, \n- n PC for features column indicates the number of principal components from which I calculated the averaged by drug and cell type values (TE features) \n- augm 50-100 means that the augm argument was set to 50 or 100,\n- s1/0.1, s1/1, s0.1/0.1, etc. stands for the noise added to duplicated train rows, the first number - is noise to features for the drugs, and the second number - is noise to the features for cell types\n- n_cells >2 or 5 pairs - means that the de_train was filtered as in this notebook to remove samples obtained just from 1-2 or 1-5 cells. The 'before PCA' note indicates that the filtration was done before feature generation (targets PCA and averaging by drug and cell type)\n- by row 10%/0.5 - stands for another method of train augmentation, where within each fold I took 10% of the train set features and respective targets and added their average\n- scaled_counts_features TE - features generated from raw counts (adata_train), stored [here](https://www.kaggle.com/datasets/antoninadolgorukova/op2-supplementary-calcs-for-ml)","metadata":{}},{"cell_type":"markdown","source":"This ensemble was added in the final ensemble of the following models (all are MLP as the one in this notebook):","metadata":{}},{"cell_type":"markdown","source":"| Submission file name | Model description | Public LB | Privat LB |\n|-------------------------------------|-------------------------------------------------------------------------------------------------------------------|-----------|-----------|\n| ave_blend_Treg_NK_0.574.csv | removed 2-3 samples with NK cells (8 folds) and T regulatory cells (8 folds) | 0.574 | 0.765 |\n| ave_61folds_0.574.csv | 61 fold, each without 10-14 obs. | 0.574 | 0.763 |\n| ave_all_train_614_0.574.csv | 16 folds trained on the entire train data | 0.574 | 0.767 |\n| blend_9diverse_models_0.573.csv | 9 models with max diversity (averaged) | 0.573 | 0.772 |\n| ave_blend_M_B_v57_0.572.csv | removed 2-3 samples with Myeloid cells (8 folds) and B cells (8 folds) | 0.572 | 0.761 |\n| ave_blend_T4_T8_v58_0.571.csv | removed 2-3 samples with T cells CD4+ (8 folds) and T cells CD8+ (8 folds) | 0.571 | 0.763 |\n| ave_blend_T4_NK_minus_T8_3kmeans_0.569.csv | removed samples from CD8+ cells; removed 2-3 samples with NK cells (8 folds) and T cells CD4+ (8 folds), predicted separately 3 gene groups (kmeans clusters) | 0.569 | 0.769 |\n| ave_blend_T4_T8_3kmeans_sep_f_0.569.csv | removed 2-3 samples with T cells CD4+ (8 folds) and T cells CD8+ (8 folds), predicted separately 3 gene groups (kmeans clusters) | 0.569 | 0.764 |\n| ave_blend_T8_Treg_3kmeans_0.57.csv | removed 2-3 samples with T cells CD8+ (8 folds) and T regulatory cells (8 folds), predicted separately 3 gene groups (kmeans clusters) | 0.57 | 0.766 |\n| ave_blend_pseudo0.571_T4_T8_20ep_augm50_s1_0.1_0.572.csv | predictions of ave_blend_T4_T8_v58_0.571 added in the train set of each fold, removed 2-3 samples with T cells CD4+ (8 folds) and T cells CD8+ (8 folds), predicted separately 3 gene groups (kmeans clusters) | 0.572 | 0.768 |\n","metadata":{}},{"cell_type":"markdown","source":"The first 8: for the ensemble models with the same LB scores were averaged, then the resulting predictions were averaged with coefficients (the lesser the public score, the more weight, see below) \nThe last 2 models + ave_blend_T4_T8_3kmeans_sep_f_0.569.csv (again): added at the final step with the weigth < 0.2 \n\n**Ensemble scores: \npublic: 0.566 \nprivat: 0.76**","metadata":{}},{"cell_type":"markdown","source":"# - Correlation (mean by row) between the included MLPs","metadata":{}},{"cell_type":"code","source":"list_mlp <- lapply(paste0(path, file_names), function(f) { \n dt <- as.matrix(fread(f))\n})\nnames(list_mlp) <- sapply(file_names, function(x) {\n gsub(\"\\\\.csv\", \"\", gsub('.*/(.*).*\\\\.csv', '\\\\1', x) )} )\n\nlapply(list_mlp, function(dt) dt[1:2, 1:10]) ","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:26:42.741187Z","iopub.execute_input":"2023-12-07T11:26:42.743221Z","iopub.status.idle":"2023-12-07T11:27:08.988403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig (20, 15)\nall_cors <- calc_cor_ave_and_plot(list_mlp)","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:27:08.993101Z","iopub.execute_input":"2023-12-07T11:27:08.996118Z","iopub.status.idle":"2023-12-07T11:28:32.434614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# - Prediction variances","metadata":{}},{"cell_type":"markdown","source":"Prediction variances are calculated as mean by row of max - min prediction of each gene","metadata":{}},{"cell_type":"code","source":"mlp_dt <- rbindlist(lapply(list_mlp, function(df) as.data.table(df)) )\nmlp_dt <- melt(mlp_dt, id.vars = \"id\", variable.name = \"gene\")\n\nmlp_ranges <- mlp_dt[, .(min_val = min(value),\n max_val = max(value)),\n by = c(\"id\", \"gene\")]\n \nmlp_ranges[, range_preds := max_val - min_val]\nmlp_ranges <- mlp_ranges[, .(mlp_var = mean(range_preds)),\n by = \"id\"] \nmlp_ranges <- mlp_ranges[id_map, on = \"id\"]\n\nmlp_ranges[1:10]\ncat(\"Median variance:\", round(median(mlp_ranges[, mlp_var]), 2))","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:28:32.439554Z","iopub.execute_input":"2023-12-07T11:28:32.442372Z","iopub.status.idle":"2023-12-07T11:28:43.719237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Good and bad predictions","metadata":{}},{"cell_type":"markdown","source":"Samples predicted with the highest variance by row, meaning that all model versions are not confident about them.","metadata":{}},{"cell_type":"code","source":"mlp_ranges[order(-mlp_var)][mlp_var > 1]","metadata":{"execution":{"iopub.status.busy":"2023-12-07T12:04:05.830300Z","iopub.execute_input":"2023-12-07T12:04:05.832784Z","iopub.status.idle":"2023-12-07T12:04:05.872681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Samples predicted with the lowest variance by row, meaning that all model versions are quite confident about them.","metadata":{}},{"cell_type":"code","source":"mlp_ranges[order(mlp_var)][1:15]","metadata":{"execution":{"iopub.status.busy":"2023-12-07T11:38:28.708211Z","iopub.execute_input":"2023-12-07T11:38:28.710491Z","iopub.status.idle":"2023-12-07T11:38:28.751012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Ensemble of MLPs","metadata":{}},{"cell_type":"code","source":"cl1_0.569 <- list_mlp[['ave_blend_T4_NK_minus_T8_3kmeans_0.569']] * 0.5 + \n list_mlp[['ave_blend_T4_T8_3kmeans_sep_f_0.569']] * 0.5\n\ncl2_61f_NK_Treg_0.574 <- list_mlp[['ave_61folds_0.574']] * 0.5 +\n list_mlp[['ave_blend_Treg_NK_0.574']] * 0.5\n\ntmp <- cl2_61f_NK_Treg_0.574 * 0.5 + list_mlp[['ave_all_train_614_0.574']] * 0.5\n\ntmp <- tmp * 0.5 + list_mlp[['blend_9diverse_models_0.573']] * 0.5\n\ntmp <- tmp * 0.5 + list_mlp[['ave_blend_M_B_v57_0.572']] * 0.5\n\ntmp <- tmp * 0.5 + list_mlp[['ave_blend_T4_T8_v58_0.571']] * 0.5\n\nblend_mlp <- tmp * 0.5 + cl1_0.569 * 0.5\n\nhead(blend_mlp)","metadata":{"execution":{"iopub.status.busy":"2023-12-07T12:05:23.994319Z","iopub.execute_input":"2023-12-07T12:05:24.001599Z","iopub.status.idle":"2023-12-07T12:05:25.026102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fwrite(blend_mlp, \"final_MLP_ensemble_0.566.csv\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ","metadata":{}}]}