{"cells":[{"metadata":{},"cell_type":"markdown","source":"The objective of this notebook is to use simple linear quantile models on the tabular data. I will try to keep it very simple to read and understand what is going on. Granted, these are not the best models, and I make no claim to have specific knowledge that would serve as a basis to these models. I just hope that the simplicity of the concepts in this notebook can give you some ideas that you could explore during these last days before the submissions deadline closes. But it can also be useful more generally by people that want to familiarise themselves with R.\n\nI have written before in the introduction of [another notebook](https://www.kaggle.com/douglaskgaraujo/bayesian-log-linear-model) why I am using log FVC values below. And this notebook itself is a quantile regression version of [another similar notebook focused on simple log linear models](https://www.kaggle.com/douglaskgaraujo/simple-log-linear-models).\n\nBy the way, an excellent introduction to quantile regression is found in Koenker and Hallock's [2001 Journal of Economic Perspectives paper](https://pubs.aeaweb.org/doi/pdfplus/10.1257/jep.15.4.143).\n\nBefore we dive in, another aspect of this notebook I would highlight is that the model is based only on the information from the patient at the time of first consultation (ie, the same week of the first FVC measurement). This is the data we will also have available for the test data."},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"library(tidyverse)\nlibrary(stargazer)\nlibrary(quantreg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train <- read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\ntest <- read_csv(\"../input/osic-pulmonary-fibrosis-progression/test.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I created the function below to transform the dataset from the way we first load it into something closer I intend to use in the models. The reason I created the function instead of just applying these commands is that I want the exact same transformations to happen in the test data when we are preparing the submission file."},{"metadata":{"trusted":true},"cell_type":"code","source":"transformOriginalDataset <- function(dataset) {\n  # Uses log FVC to train the regression models\n  dataset$FVC_original <- dataset$FVC\n  dataset$FVC  <- log(dataset$FVC)\n  \n  # Transforms the 'Sex' variable into a dummy,\n  # a dummy for young patients, defined as those with less than 50 years of age,\n  # the smoking status variable into a yes/no currently smoker dummy,\n  # includes a variable with the patient-specific benchmark FVC\n  # (ie, FVC equal to 100% theoretical capacity),\n  # and finally includes the baseline Percent value\n  dataset <- dataset %>% \n    mutate(Male = ifelse(Sex == \"Male\", 1, 0),\n           CurrentlySmokes = ifelse(SmokingStatus == \"Currently smokes\", 1, 0),\n           YoungPatient = ifelse(Age < 50, 1, 0)) %>% \n    group_by(Patient) %>% \n    mutate(HundredPctFVC = 100 * mean(FVC_original / Percent),\n           FVCBase = head(FVC, 1),\n           PercentBase = head(Percent, 1))\n  \n  return(dataset)\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"exclude_test_patient_data_from_trainset <- TRUE\n\nif(exclude_test_patient_data_from_trainset) {\n  train <- train %>% filter(!(Patient %in% unique(test$Patient)))\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train <- transformOriginalDataset(train)\nhead(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the following chunk, we will define the formula we will be using for the quantile regression, applying it to a series of quantile regressions spanning from the 5th to the 95th quantiles. To illustrate the results (and showcase why quantile regressions can be warranted in this setting), I will also plot how the coefficients vary across quantiles."},{"metadata":{"trusted":true},"cell_type":"code","source":"quantreg_formula <- as.formula(\"FVC ~  Weeks + Male + CurrentlySmokes + YoungPatient + PercentBase + Weeks * CurrentlySmokes + Weeks * PercentBase + CurrentlySmokes * Male\")\n\nall_quant_regs <- rq(quantreg_formula, tau = seq(from = 0.05, to = 0.95, by = 0.05), data = train)\nplot(summary(all_quant_regs))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the graph above, you can see by the dash-and-dot black lines representing the coefficients from the quantile regressions that these yeild estimates that are, in some cases, \"diagonal\". This is indicating that the relationship between FVC and the explanatory variables is not the same throughout the FVCs. For example, for low FVC levels, men tend to have around 60% more FVC then women, while for higher FVCs this difference lowers to circa 45%. The red lines, representing the normal least squares result that focus on mean values, does not capture such trends.\n\nNow that we have seen why quantile regression can be warranted in this case, we can get a better intuition of its use to construct the confidence that we are required to submit. Basically, we can use quantile models to estimate what a low quantile of the FVC distribution could be, as well as a high quantile, and the difference between those would be a range of FVC values we consider most likely. Of course, there is no fixed rule on what exactly \"low\" and \"high\" quantiles would represent. It will depend on each analyst. For this notebook, I will be using 15% and 85% as the low and high quantiles, respectively. (Note that, while no fixed rule mandates that these have the same distance to the respective extremes of 0% and 100% - 15% in the case that I chose - it would not make much sense for those not be symmetric.)\n\nIf you want to use as your confidence level the 5% and 95% quantiles, or in fact any other number, you can just choose change the values below."},{"metadata":{"trusted":true},"cell_type":"code","source":"low_quantile <- 0.15\nhigh_quantile <- 1 - low_quantile","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below is the collection of quantile models we will be using to model FVC. Contrary to the lines above, where I passed a vector of quantiles as input to the ```rq``` function, now I will use one quantile at a time and save models in a list so it becomes clearer to readers not so familiar with R what we will be doing next to create the predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"lower_bound_confidence_model <- rq(quantreg_formula, tau = low_quantile, data = train)\nmedian_model <- rq(quantreg_formula, tau = 0.50, data = train)\nupper_bound_confidence_model <- rq(quantreg_formula, tau = high_quantile, data = train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you customising this notebook, don't forget to confirm that you are including all coeficients for the model you chose. To do that, just remember to check out once again the variable ```quantreg_formula```."},{"metadata":{"trusted":true},"cell_type":"code","source":"quantreg_formula","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So the output of ```formula(median_model)``` shows us that we should include: Weeks, Male, CurrentlySmokes, YoungPatient, PercentBase and the multiplication between Weeks and Male, Weeks and CurrentlySmokes, and Weeks and PercentBase. Let's see the corresponding coefficients (while also keeping tab of their \"names\" in the coefficient vector). The ```knitr::kable``` function is being used just to make it look better!"},{"metadata":{"trusted":true},"cell_type":"code","source":"knitr::kable(coef(median_model))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK. So now that we have a model, let's use it to create predictions for each of those weeks that are required in the submission (ie, -12 weeks before the baseline CT and all the way until 133 weeks after). I tried making my code as verbose as possible so that you can understand what I am doing. Functions in the ```tidyverse``` have excellent documentation and loads of community questions & answers, so don't hesitate to investigate if you have not understood any specific point."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_sub <- transformOriginalDataset(test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Variable ```test_sub``` has now undergone the transformations that were also applied to the training test when used it to model FVC outcomes. But it is not yet ready to be used to construct the submission file. To understand why, note that the submission rules require us to generate a predicted FVC and confidence for each patient in the test set during a large range of weeks (from -12 to 133, as I'm sure you know). So now we need some supporting variables and a code that will create such a data.frame for us, and then fill it out with the predictions for FVC and confidence values."},{"metadata":{"trusted":true},"cell_type":"code","source":"weeks <- seq(from = -12, to = 133, by = 1)\nPatient_Week <- paste(rep(unique(test$Patient), each = length(weeks)), \n                      seq(from = -12, to = 133, by = 1), sep = \"_\")\n\ntest_sub <- test_sub %>% \n  group_by(Patient) %>% \n  summarise(FVCBase = head(FVCBase, 1),\n            Male = head(Male, 1),\n            CurrentlySmokes = head(CurrentlySmokes, 1),\n            YoungPatient = head(YoungPatient, 1),\n            PercentBase = head(PercentBase, 1)) %>% \n  uncount(length(weeks)) %>% \n  cbind(Patient_Week, Weeks = rep(seq(from = -12, to = 133, by = 1), length(unique(test$Patient))), .)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the first stage above, we merely compiled the patient-specific information we will use to predict the values. Note that each patient only has one value here. The reason is that our predictions will be measured against test dataset containing only one row of information per patient. So our prediction needs to be based on the first value only as well. That is what we will be doing next. Just remember to keep it consistent with ```formula(median_model)``` if you adapt this notebook.\n\nAlso, note that below we exponentiate the estimated values, because the models were estimated with log FVC data."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_sub <- test_sub %>% \n  mutate(FVC = exp(coef(median_model)[which(names(coef(median_model)) == \"(Intercept)\")] +\n                   coef(median_model)[which(names(coef(median_model)) == \"Weeks\")] * Weeks +\n                   coef(median_model)[which(names(coef(median_model)) == \"Male\")] * Male +\n                   coef(median_model)[which(names(coef(median_model)) == \"CurrentlySmokes\")] * CurrentlySmokes +\n                   coef(median_model)[which(names(coef(median_model)) == \"YoungPatient\")] * YoungPatient +\n                   coef(median_model)[which(names(coef(median_model)) == \"PercentBase\")] * PercentBase +\n                   coef(median_model)[which(names(coef(median_model)) == \"Weeks:CurrentlySmokes\")] * Weeks * CurrentlySmokes +\n                   coef(median_model)[which(names(coef(median_model)) == \"Weeks:PercentBase\")] * Weeks * PercentBase +\n                   coef(median_model)[which(names(coef(median_model)) == \"Male:CurrentlySmokes\")] * Male * CurrentlySmokes\n         ),\n         FVC_lower_bound = exp(coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"(Intercept)\")] +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"Weeks\")] * Weeks +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"Male\")] * Male +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"CurrentlySmokes\")] * CurrentlySmokes +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"YoungPatient\")] * YoungPatient +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"PercentBase\")] * PercentBase +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"Weeks:CurrentlySmokes\")] * Weeks * CurrentlySmokes +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"Weeks:PercentBase\")] * Weeks * PercentBase +\n                               coef(lower_bound_confidence_model)[which(names(coef(lower_bound_confidence_model)) == \"Male:CurrentlySmokes\")] * Male * CurrentlySmokes\n         ),\n         FVC_upper_bound = exp(coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"(Intercept)\")] +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"Weeks\")] * Weeks +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"Male\")] * Male +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"CurrentlySmokes\")] * CurrentlySmokes +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"YoungPatient\")] * YoungPatient +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"PercentBase\")] * PercentBase +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"Weeks:CurrentlySmokes\")] * Weeks * CurrentlySmokes +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"Weeks:PercentBase\")] * Weeks * PercentBase +\n                               coef(upper_bound_confidence_model)[which(names(coef(upper_bound_confidence_model)) == \"Male:CurrentlySmokes\")] * Male * CurrentlySmokes\n         ),\n         Confidence = (FVC_upper_bound - FVC_lower_bound) / 2) %>% \n  select(Patient_Week, FVC, Confidence)\n\nhead(test_sub)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Great! Now let's write the submission file."},{"metadata":{"trusted":true},"cell_type":"code","source":"write_csv(test_sub, \"submission.csv\")","execution_count":null,"outputs":[]}],"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":"3.6.3"}},"nbformat":4,"nbformat_minor":4}