{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# CITE Seq Baseline Prediction","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle"}},{"cell_type":"markdown","source":"This notebook follows a similar approach for [Multiomer prediction](https://www.kaggle.com/code/stautxie/multiomer-prediction-with-svd-regression). We use the following approach to derive a solution:\n\n1. Apply truncated SVD decomposition to reduce the problem size.\n2. Deploy a vectorized competition metric function.\n3. Experiment with Lasso and Ridge regression and confirms that Lasso regression is better. Our experiment also shows that a relevantly unsophisticated regression approach delivers slightly better performance than a much more complex [gradient boosting tree approach](https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart).","metadata":{}},{"cell_type":"code","source":"\nlibrary(data.table)\nlibrary(Matrix)\nlibrary(magrittr)\nlibrary(mltools)\nlibrary(tictoc)\nlibrary(ggplot2)\nlibrary(magrittr)\nlibrary(glmnet)\nlibrary(RSpectra)\nlibrary(caret)","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:29:53.283859Z","iopub.execute_input":"2022-09-03T04:29:53.287923Z","iopub.status.idle":"2022-09-03T04:29:56.439368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nDAT_DIR <- \"../input/open-problems-multimodal\"\nMAT_DIR <- \"../input/sparse-measurement-data-open-problems-multimodal\"","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:30:41.778025Z","iopub.execute_input":"2022-09-03T04:30:41.826310Z","iopub.status.idle":"2022-09-03T04:30:41.845665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_ci_inputs <- readRDS(file.path(MAT_DIR, \"sp_train_cite_inputs.rds\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:31:02.518637Z","iopub.execute_input":"2022-09-03T04:31:02.520950Z","iopub.status.idle":"2022-09-03T04:32:10.425227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We apply truncated SVD to reduce the input matrix size.","metadata":{}},{"cell_type":"code","source":"\nn_comp <- 32\n\ntic()\nsvd_input <- svds(tr_ci_inputs,\n                  k = n_comp)\ntoc()\nsaveRDS(svd_input, \"ci_svd_input.rds\")\n\ndata.table(d2 = svd_input$d[2:n_comp]^2 / max(svd_input$d[2:n_comp]^2),\n           i = seq_along(svd_input$d[2:n_comp])) %>%\n  ggplot(aes(x = i, y = d2)) +\n  geom_point() +\n  geom_line() +\n  #geom_text(aes(label=round(svd_res$d[2:n_comp]^2 / max(svd_res$d[2:n_comp]^2), 2)), hjust = -0.1, vjust=-1) +\n  labs(x = NULL, y = \"d squared / d[2] squared\") \n\ntr_ci_input_svd <- svd_input$u %*% diag(svd_input$d, nrow=length(svd_input$d), ncol=length(svd_input$d))\ncolnames(tr_ci_input_svd) <- paste(\"SVD_\", 1:n_comp, sep=\"\")","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:38:00.476994Z","iopub.execute_input":"2022-09-03T04:38:00.479839Z","iopub.status.idle":"2022-09-03T04:41:28.455332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(tr_ci_inputs); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:42:06.443773Z","iopub.execute_input":"2022-09-03T04:42:06.446239Z","iopub.status.idle":"2022-09-03T04:42:07.767603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To circumvent the Kaggle's memory limit of 16 GB, we will select 20% of the total rows.","metadata":{}},{"cell_type":"code","source":"\ntr_ci_targets <- readRDS(file.path(MAT_DIR, \"sp_train_cite_targets.rds\"))\n\nset.seed(1234)\n\n#~ Take a sample of training data\nselected_pct <- 0.2\nselected_rows <- sample(1:nrow(tr_ci_targets), floor(selected_pct * nrow(tr_ci_targets)), replace=FALSE)","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:42:10.791630Z","iopub.execute_input":"2022-09-03T04:42:10.798454Z","iopub.status.idle":"2022-09-03T04:42:13.122374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_ci_input_svd_sampled <- tr_ci_input_svd[selected_rows, ]\ntr_ci_targets_sampled <- tr_ci_targets[selected_rows,] %>%\n  as.matrix()\n\nrm(tr_ci_targets); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:42:15.812922Z","iopub.execute_input":"2022-09-03T04:42:15.814884Z","iopub.status.idle":"2022-09-03T04:42:16.864931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Scoring Function","metadata":{}},{"cell_type":"markdown","source":"The official scoring function can be found [here](https://www.kaggle.com/competitions/open-problems-multimodal/overview/evaluation) In summary, the predicted target is compared with the group solution for each sample (i.e., per row) and is measured by Pearson correlation coefficient. The final score is the mean of the row-wise correlation cofficient. ","metadata":{}},{"cell_type":"markdown","source":"A popular Python implementation of the scoring function can be found [here](https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart) where explicit loop is required to calculate the mean. However, we found a much more straightforward solution that uses vectorization instead of explicit looping.  This vectorization solution not only leads to more succinct representation but also potentially more efficient as well.","metadata":{}},{"cell_type":"markdown","source":"Our implementation of the scoring function is shown below in row_mean_corr(). This implementation takes advantage of\n\n1. Base R function cor() returns a matrix $A$ where each cell $a_{ij} := correlation(column_{i}, column_{j})$. By transposing the true target matrix $Y_{true}$ and prediction matrix $Y_{pred}$, the matrix $t(A)$ contains row-wise correlation.\n\n2. Base R function diag() returns diagonal values of a matrix.\n\n3. mean(diag(.)) returns the mean of row-wise correlation coefficient, which is the competition metric.","metadata":{}},{"cell_type":"code","source":"row_mean_corr <- function(Y_true, Y_pred) {\n  cor_mat <- cor(t(Y_true), t(Y_pred), method=\"pearson\")\n  return (mean(diag(cor_mat)))\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:42:20.844788Z","iopub.execute_input":"2022-09-03T04:42:20.846406Z","iopub.status.idle":"2022-09-03T04:42:20.861730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the next section, we will experiment with Ridge regression (which has been quite widely used in this competition) with its cousin Lasso regression.","metadata":{}},{"cell_type":"markdown","source":"### Ridge Regression","metadata":{}},{"cell_type":"code","source":"set.seed(1234)\ntest_folds <- createFolds(1:nrow(tr_ci_input_svd_sampled), k=5, list=TRUE, returnTrain=FALSE)\n\ntest_fit_scores <- list()\ncor_fit_scores <- list()\n\nfor (fold in seq_along(test_folds)) {\n  tr_X <- tr_ci_input_svd_sampled[-test_folds[[fold]], ]\n  tr_Y <- tr_ci_targets_sampled[-test_folds[[fold]], ] %>% as.matrix()\n  te_X <- tr_ci_input_svd_sampled[test_folds[[fold]], ]\n  te_Y <- tr_ci_targets_sampled[test_folds[[fold]], ] %>% as.matrix()\n  fit <- glmnet(tr_X,\n                tr_Y,\n                family = \"mgaussian\",\n                alpha = 0,\n                lambda = 1,\n                trace.it = TRUE)\n  \n  pred <- predict(fit, newx=te_X, s=1)\n  dim(pred) <- dim(te_Y)\n  mse <- mean((pred - te_Y)^2)\n  test_fit_scores[[fold]] <- mse\n  mean_cor_val <- row_mean_corr(te_Y, pred)\n  cor_fit_scores[[fold]] <- mean_cor_val\n  cat(paste(\"fold = \", fold, \", mse = \", mse, \", cor = \", mean_cor_val, \"\\n\"))\n  rm(tr_X, tr_Y, te_X, te_Y); gc()\n}\n\ncat(paste(\"OOF mean MSE: \", mean(as.numeric(test_fit_scores)), \"\\n\", sep=\"\"))\ncat(paste(\"OOF mean Cor: \", mean(as.numeric(cor_fit_scores)), \"\\n\", sep=\"\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:42:23.701070Z","iopub.execute_input":"2022-09-03T04:42:23.703308Z","iopub.status.idle":"2022-09-03T04:42:43.961348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lasso Regression","metadata":{}},{"cell_type":"code","source":"set.seed(1234)\ntest_folds <- createFolds(1:nrow(tr_ci_input_svd_sampled), k=5, list=TRUE, returnTrain=FALSE)\n\ntest_fit_scores <- list()\ncor_fit_scores <- list()\n\nfor (fold in seq_along(test_folds)) {\n  tr_X <- tr_ci_input_svd_sampled[-test_folds[[fold]], ]\n  tr_Y <- tr_ci_targets_sampled[-test_folds[[fold]], ] %>% as.matrix()\n  te_X <- tr_ci_input_svd_sampled[test_folds[[fold]], ]\n  te_Y <- tr_ci_targets_sampled[test_folds[[fold]], ] %>% as.matrix()\n  fit <- glmnet(tr_X,\n                tr_Y,\n                family = \"mgaussian\",\n                alpha = 1,\n                lambda = 1,\n                trace.it = TRUE)\n  \n  pred <- predict(fit, newx=te_X, s=1)\n  dim(pred) <- dim(te_Y)\n  mse <- mean((pred - te_Y)^2)\n  test_fit_scores[[fold]] = mse\n  mean_cor_val <- row_mean_corr(te_Y, pred)\n  cor_fit_scores[[fold]] <- mean_cor_val\n  cat(paste(\"fold = \", fold, \", mse = \", mse, \", cor = \", mean_cor_val, \"\\n\"))\n  rm(tr_X, tr_Y, te_X, te_Y); gc()\n}\n\ncat(paste(\"Lasso OOF mean MSE: \", mean(as.numeric(test_fit_scores)), \"\\n\", sep=\"\"))\ncat(paste(\"Lasso OOF mean Cor: \", mean(as.numeric(cor_fit_scores)), \"\\n\", sep=\"\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:43:48.870777Z","iopub.execute_input":"2022-09-03T04:43:48.872642Z","iopub.status.idle":"2022-09-03T04:44:09.687684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, Lasso regression performs better than Ridge regression under similar settings. Therefore, we will use Lasso regression to build the prediction model. In addition, the OOF fitting score is marginally better than [this approach that uses lightgbm](https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart).","metadata":{}},{"cell_type":"markdown","source":"### Prediction","metadata":{}},{"cell_type":"code","source":"\nfit <- glmnet(tr_ci_input_svd_sampled,\n              tr_ci_targets_sampled %>% as.matrix(),\n              family = \"mgaussian\",\n              alpha = 1,\n              lambda = 1,\n              trace.it = TRUE)\n\nsaveRDS(fit, \"fit_mu_lasso.rds\")","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:48:45.900950Z","iopub.execute_input":"2022-09-03T04:48:45.902582Z","iopub.status.idle":"2022-09-03T04:48:47.303385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(tr_ci_input_svd_sampled, tr_ci_targets_sampled); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:48:52.971246Z","iopub.execute_input":"2022-09-03T04:48:52.972942Z","iopub.status.idle":"2022-09-03T04:48:53.291878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will load the saved SVD transformation matrix to transform the testing input matrix.","metadata":{}},{"cell_type":"code","source":"te_ci_inputs <- readRDS(file.path(MAT_DIR, \"sp_test_cite_inputs.rds\"))\nsvd_input <- readRDS(\"ci_svd_input.rds\")\nte_ci_inputs_svd <- te_ci_inputs %*% svd_input$v\nrm(te_ci_inputs); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:49:05.539672Z","iopub.execute_input":"2022-09-03T04:49:05.542453Z","iopub.status.idle":"2022-09-03T04:49:57.766654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we will save the results for subsequent work. To avoid hitting the memory limit, we will just predict the first 5000 rows. It is fairly straightforward to develop an iterative prediction routine if needed.","metadata":{}},{"cell_type":"code","source":"te_pred <- predict(fit, newx=te_ci_inputs_svd[1:5000,], s=1)\ndim(te_pred) <- dim(te_pred)[-3]\nrm(te_ci_inputs_svd); gc()\nsaveRDS(te_pred, \"ci_te_pred.rds\")","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:50:13.804204Z","iopub.execute_input":"2022-09-03T04:50:13.805826Z","iopub.status.idle":"2022-09-03T04:50:15.311627Z"},"trusted":true},"execution_count":null,"outputs":[]}]}