dataSciencePractical1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ef6c393d",
   "metadata": {},
   "source": [
    "# Data science practical\n",
    "In this data science tutorial we will reporoduce some of the analysis from \"A Validated Prediction Model for Overall \n",
    "Survival From Stage III Non-Small Cell Lung Cancer: Toward Survival Prediction for Individual Patients\" from the Maastro clinic.\n",
    "\n",
    "<img src=\"images/paperTitle.png\">\n",
    "\n",
    "The paper can be found at [with this DOI](https://doi.org/10.1016/j.ijrobp.2015.02.048)\n",
    "\n",
    "The model is also hosted at [predictcancer.org](http://predictcancer.org/Main.php?page=Home) where and can be interacted with [via a web interface](http://predictcancer.org/Main.php?page=LungSurvivalModel3).\n",
    "\n",
    "The lesson is given as a Jupyter notebook, hosted from a github repository, that you will be able to access online. Separate instructions on how to access this via [binder](https://mybinder.org) are given separately."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ca0f7013",
   "metadata": {},
   "source": [
    "## Model data\n",
    "The data used in the paper is hosted at [cancerdata.org](https://www.cancerdata.org/data/files). This should be hosted alongside the code you're using but can be downloaded (along with many other datasets if you are interested) by:\n",
    "* Navigating to the “A validated prediction model from overall survival from Stage III Non-Small Cell Lung Cancer – towards survival prediction for individual patients” dataset\n",
    "* Selecting csv file format\n",
    "\n",
    "<img src=\"images/cancerdatalocation.png\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ca3acc5",
   "metadata": {},
   "source": [
    "## Load data workspace\n",
    "The first thing to do is to load the csv and inspect the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "06612cd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.data = read.csv(file='data/Stage3_anonymized.csv',head=TRUE,sep=\";\",na.strings=c(\"NA\",\"\",\" \"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b26ce9f",
   "metadata": {},
   "source": [
    "The CSV is loaded into R as a `data.frame` object called `lung.data`. This is essentially a copy of the columns and rows of the CSV data that can be manipulated by R. Using the `names()` command tells us what each column (variable) is called. The `head()` command inspect the first few rows of each column so we can see the data type in each. We can also count the number of rows (i.e. patients in the cohort) using the `nrow()` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cee42bda",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'study_id'</li><li>'gender'</li><li>'age'</li><li>'who3g'</li><li>'bmi'</li><li>'fev1pc_t0'</li><li>'dumsmok2'</li><li>'t_ct_loc'</li><li>'hist4g'</li><li>'countpet_all6g'</li><li>'countpet_mediast6g'</li><li>'tstage'</li><li>'nstage'</li><li>'stage'</li><li>'timing'</li><li>'group'</li><li>'yearrt'</li><li>'eqd2'</li><li>'ott'</li><li>'gtv1'</li><li>'tumorload_total'</li><li>'survmonth'</li><li>'survyear'</li><li>'deadstat'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'study\\_id'\n",
       "\\item 'gender'\n",
       "\\item 'age'\n",
       "\\item 'who3g'\n",
       "\\item 'bmi'\n",
       "\\item 'fev1pc\\_t0'\n",
       "\\item 'dumsmok2'\n",
       "\\item 't\\_ct\\_loc'\n",
       "\\item 'hist4g'\n",
       "\\item 'countpet\\_all6g'\n",
       "\\item 'countpet\\_mediast6g'\n",
       "\\item 'tstage'\n",
       "\\item 'nstage'\n",
       "\\item 'stage'\n",
       "\\item 'timing'\n",
       "\\item 'group'\n",
       "\\item 'yearrt'\n",
       "\\item 'eqd2'\n",
       "\\item 'ott'\n",
       "\\item 'gtv1'\n",
       "\\item 'tumorload\\_total'\n",
       "\\item 'survmonth'\n",
       "\\item 'survyear'\n",
       "\\item 'deadstat'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'study_id'\n",
       "2. 'gender'\n",
       "3. 'age'\n",
       "4. 'who3g'\n",
       "5. 'bmi'\n",
       "6. 'fev1pc_t0'\n",
       "7. 'dumsmok2'\n",
       "8. 't_ct_loc'\n",
       "9. 'hist4g'\n",
       "10. 'countpet_all6g'\n",
       "11. 'countpet_mediast6g'\n",
       "12. 'tstage'\n",
       "13. 'nstage'\n",
       "14. 'stage'\n",
       "15. 'timing'\n",
       "16. 'group'\n",
       "17. 'yearrt'\n",
       "18. 'eqd2'\n",
       "19. 'ott'\n",
       "20. 'gtv1'\n",
       "21. 'tumorload_total'\n",
       "22. 'survmonth'\n",
       "23. 'survyear'\n",
       "24. 'deadstat'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] \"study_id\"           \"gender\"             \"age\"               \n",
       " [4] \"who3g\"              \"bmi\"                \"fev1pc_t0\"         \n",
       " [7] \"dumsmok2\"           \"t_ct_loc\"           \"hist4g\"            \n",
       "[10] \"countpet_all6g\"     \"countpet_mediast6g\" \"tstage\"            \n",
       "[13] \"nstage\"             \"stage\"              \"timing\"            \n",
       "[16] \"group\"              \"yearrt\"             \"eqd2\"              \n",
       "[19] \"ott\"                \"gtv1\"               \"tumorload_total\"   \n",
       "[22] \"survmonth\"          \"survyear\"           \"deadstat\"          "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 24</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>study_id</th><th scope=col>gender</th><th scope=col>age</th><th scope=col>who3g</th><th scope=col>bmi</th><th scope=col>fev1pc_t0</th><th scope=col>dumsmok2</th><th scope=col>t_ct_loc</th><th scope=col>hist4g</th><th scope=col>countpet_all6g</th><th scope=col>...</th><th scope=col>timing</th><th scope=col>group</th><th scope=col>yearrt</th><th scope=col>eqd2</th><th scope=col>ott</th><th scope=col>gtv1</th><th scope=col>tumorload_total</th><th scope=col>survmonth</th><th scope=col>survyear</th><th scope=col>deadstat</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>...</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>1</td><td>1</td><td>63,1567419575633</td><td>2</td><td>22,9938271604938</td><td>82  </td><td>2</td><td>6</td><td>4</td><td>4</td><td>...</td><td>3</td><td>4</td><td>2008</td><td>53,125</td><td>28</td><td>192,86          </td><td>214,08          </td><td>16,9199178644764</td><td>1,40999315537303</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>2</td><td>2</td><td>66,3052703627652</td><td>2</td><td>28,515625       </td><td>62  </td><td>2</td><td>1</td><td>1</td><td>4</td><td>...</td><td>3</td><td>4</td><td>2010</td><td>55,125</td><td>32</td><td>154,93          </td><td>169,9           </td><td>5,12525667351129</td><td>,427104722792608</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>3</td><td>1</td><td>58,1984941820671</td><td>1</td><td>25,1278305332359</td><td>95  </td><td>2</td><td>4</td><td>1</td><td>1</td><td>...</td><td>3</td><td>4</td><td>2008</td><td>65,125</td><td>32</td><td>85,07           </td><td>85,07           </td><td>57,96           </td><td>4,83            </td><td>0</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>4</td><td>2</td><td>63,7809719370294</td><td>2</td><td>28,7267877295578</td><td>73  </td><td>2</td><td>6</td><td>2</td><td>6</td><td>...</td><td>3</td><td>4</td><td>2007</td><td>63,125</td><td>36</td><td>98,17           </td><td>120,51          </td><td>7,42505133470226</td><td>,618754277891855</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>5</td><td>1</td><td>59,4140999315537</td><td>1</td><td>20,3244382551093</td><td>124 </td><td>2</td><td>6</td><td>4</td><td>3</td><td>...</td><td>3</td><td>4</td><td>2008</td><td>65,125</td><td>36</td><td>62,89           </td><td>136,32          </td><td>16,6899383983573</td><td>1,39082819986311</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>6</td><td>1</td><td>49,2594113620808</td><td>2</td><td>21,5138585105239</td><td>56,8</td><td>2</td><td>1</td><td>4</td><td>3</td><td>...</td><td>3</td><td>4</td><td>2006</td><td>63,125</td><td>35</td><td>37,8300018310547</td><td>43,7400016784668</td><td>22,0123203285421</td><td>1,83436002737851</td><td>1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 24\n",
       "\\begin{tabular}{r|lllllllllllllllllllll}\n",
       "  & study\\_id & gender & age & who3g & bmi & fev1pc\\_t0 & dumsmok2 & t\\_ct\\_loc & hist4g & countpet\\_all6g & ... & timing & group & yearrt & eqd2 & ott & gtv1 & tumorload\\_total & survmonth & survyear & deadstat\\\\\n",
       "  & <int> & <int> & <fct> & <int> & <fct> & <fct> & <int> & <int> & <int> & <int> & ... & <int> & <int> & <int> & <fct> & <int> & <fct> & <fct> & <fct> & <fct> & <int>\\\\\n",
       "\\hline\n",
       "\t1 & 1 & 1 & 63,1567419575633 & 2 & 22,9938271604938 & 82   & 2 & 6 & 4 & 4 & ... & 3 & 4 & 2008 & 53,125 & 28 & 192,86           & 214,08           & 16,9199178644764 & 1,40999315537303 & 1\\\\\n",
       "\t2 & 2 & 2 & 66,3052703627652 & 2 & 28,515625        & 62   & 2 & 1 & 1 & 4 & ... & 3 & 4 & 2010 & 55,125 & 32 & 154,93           & 169,9            & 5,12525667351129 & ,427104722792608 & 1\\\\\n",
       "\t3 & 3 & 1 & 58,1984941820671 & 1 & 25,1278305332359 & 95   & 2 & 4 & 1 & 1 & ... & 3 & 4 & 2008 & 65,125 & 32 & 85,07            & 85,07            & 57,96            & 4,83             & 0\\\\\n",
       "\t4 & 4 & 2 & 63,7809719370294 & 2 & 28,7267877295578 & 73   & 2 & 6 & 2 & 6 & ... & 3 & 4 & 2007 & 63,125 & 36 & 98,17            & 120,51           & 7,42505133470226 & ,618754277891855 & 1\\\\\n",
       "\t5 & 5 & 1 & 59,4140999315537 & 1 & 20,3244382551093 & 124  & 2 & 6 & 4 & 3 & ... & 3 & 4 & 2008 & 65,125 & 36 & 62,89            & 136,32           & 16,6899383983573 & 1,39082819986311 & 1\\\\\n",
       "\t6 & 6 & 1 & 49,2594113620808 & 2 & 21,5138585105239 & 56,8 & 2 & 1 & 4 & 3 & ... & 3 & 4 & 2006 & 63,125 & 35 & 37,8300018310547 & 43,7400016784668 & 22,0123203285421 & 1,83436002737851 & 1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 24\n",
       "\n",
       "| <!--/--> | study_id &lt;int&gt; | gender &lt;int&gt; | age &lt;fct&gt; | who3g &lt;int&gt; | bmi &lt;fct&gt; | fev1pc_t0 &lt;fct&gt; | dumsmok2 &lt;int&gt; | t_ct_loc &lt;int&gt; | hist4g &lt;int&gt; | countpet_all6g &lt;int&gt; | ... ... | timing &lt;int&gt; | group &lt;int&gt; | yearrt &lt;int&gt; | eqd2 &lt;fct&gt; | ott &lt;int&gt; | gtv1 &lt;fct&gt; | tumorload_total &lt;fct&gt; | survmonth &lt;fct&gt; | survyear &lt;fct&gt; | deadstat &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1 | 1 | 1 | 63,1567419575633 | 2 | 22,9938271604938 | 82   | 2 | 6 | 4 | 4 | ... | 3 | 4 | 2008 | 53,125 | 28 | 192,86           | 214,08           | 16,9199178644764 | 1,40999315537303 | 1 |\n",
       "| 2 | 2 | 2 | 66,3052703627652 | 2 | 28,515625        | 62   | 2 | 1 | 1 | 4 | ... | 3 | 4 | 2010 | 55,125 | 32 | 154,93           | 169,9            | 5,12525667351129 | ,427104722792608 | 1 |\n",
       "| 3 | 3 | 1 | 58,1984941820671 | 1 | 25,1278305332359 | 95   | 2 | 4 | 1 | 1 | ... | 3 | 4 | 2008 | 65,125 | 32 | 85,07            | 85,07            | 57,96            | 4,83             | 0 |\n",
       "| 4 | 4 | 2 | 63,7809719370294 | 2 | 28,7267877295578 | 73   | 2 | 6 | 2 | 6 | ... | 3 | 4 | 2007 | 63,125 | 36 | 98,17            | 120,51           | 7,42505133470226 | ,618754277891855 | 1 |\n",
       "| 5 | 5 | 1 | 59,4140999315537 | 1 | 20,3244382551093 | 124  | 2 | 6 | 4 | 3 | ... | 3 | 4 | 2008 | 65,125 | 36 | 62,89            | 136,32           | 16,6899383983573 | 1,39082819986311 | 1 |\n",
       "| 6 | 6 | 1 | 49,2594113620808 | 2 | 21,5138585105239 | 56,8 | 2 | 1 | 4 | 3 | ... | 3 | 4 | 2006 | 63,125 | 35 | 37,8300018310547 | 43,7400016784668 | 22,0123203285421 | 1,83436002737851 | 1 |\n",
       "\n"
      ],
      "text/plain": [
       "  study_id gender age              who3g bmi              fev1pc_t0 dumsmok2\n",
       "1 1        1      63,1567419575633 2     22,9938271604938 82        2       \n",
       "2 2        2      66,3052703627652 2     28,515625        62        2       \n",
       "3 3        1      58,1984941820671 1     25,1278305332359 95        2       \n",
       "4 4        2      63,7809719370294 2     28,7267877295578 73        2       \n",
       "5 5        1      59,4140999315537 1     20,3244382551093 124       2       \n",
       "6 6        1      49,2594113620808 2     21,5138585105239 56,8      2       \n",
       "  t_ct_loc hist4g countpet_all6g ... timing group yearrt eqd2   ott\n",
       "1 6        4      4              ... 3      4     2008   53,125 28 \n",
       "2 1        1      4              ... 3      4     2010   55,125 32 \n",
       "3 4        1      1              ... 3      4     2008   65,125 32 \n",
       "4 6        2      6              ... 3      4     2007   63,125 36 \n",
       "5 6        4      3              ... 3      4     2008   65,125 36 \n",
       "6 1        4      3              ... 3      4     2006   63,125 35 \n",
       "  gtv1             tumorload_total  survmonth        survyear         deadstat\n",
       "1 192,86           214,08           16,9199178644764 1,40999315537303 1       \n",
       "2 154,93           169,9            5,12525667351129 ,427104722792608 1       \n",
       "3 85,07            85,07            57,96            4,83             0       \n",
       "4 98,17            120,51           7,42505133470226 ,618754277891855 1       \n",
       "5 62,89            136,32           16,6899383983573 1,39082819986311 1       \n",
       "6 37,8300018310547 43,7400016784668 22,0123203285421 1,83436002737851 1       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Number of rows\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "548"
      ],
      "text/latex": [
       "548"
      ],
      "text/markdown": [
       "548"
      ],
      "text/plain": [
       "[1] 548"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "names(lung.data)\n",
    "head(lung.data)\n",
    "print('Number of rows')\n",
    "nrow(lung.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6aaaa3a",
   "metadata": {},
   "source": [
    "## Check and curate data\n",
    "One of the main things to check is where data is missing and if it is, decide on a strategy for taking account of this. This can be done visually using the `gg_miss_var` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d7d1e8a1",
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in gg_miss_var(lung.data): could not find function \"gg_miss_var\"\n",
     "output_type": "error",
     "traceback": [
      "Error in gg_miss_var(lung.data): could not find function \"gg_miss_var\"\nTraceback:\n"
     ]
    }
   ],
   "source": [
    "gg_miss_var(lung.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acaec082",
   "metadata": {},
   "source": [
    "The first time this command in a session is run you will get an error message telling you the function cannot be found:\n",
    "<img src=\"images/gg_miss_varMissing.PNG\">\n",
    "This shows that the library containing this function needs to be loaded (the functions we have used already are part of the core R commans). This can be done using the `library()` commnd. Execute the command below and then go back and execute the previous cell. It should now work and give you a chart showing the number of rows where the data us missing for each column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ba890ddf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Registered S3 methods overwritten by 'tibble':\n",
      "  method     from  \n",
      "  format.tbl pillar\n",
      "  print.tbl  pillar\n",
      "\n",
      "Registered S3 methods overwritten by 'ggplot2':\n",
      "  method         from \n",
      "  [.quosures     rlang\n",
      "  c.quosures     rlang\n",
      "  print.quosures rlang\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(naniar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d95875f",
   "metadata": {},
   "source": [
    "At this point we will also add the other libraries we will need throughout this demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ea73f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ggplot2)\n",
    "library(survival)\n",
    "library(survminer)\n",
    "library(haven)\n",
    "library(formattable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afef65ab",
   "metadata": {},
   "source": [
    "### Remove variables with lots of missing data\n",
    "One way of dealing with missing data is to impute it from the other data. However, this only works well when relatively small volumes of data are missing. As a rule of thumb we will exclude variables where more than 10% of values are missing. From above we know we have 548 patients (rows) and so will exclude variables with more than 55 missing entries. Looking at the above chart we can see this means both the `bmi` (Body Mass Index) and `fev1pc_t0` (Forced Expiratory Volume in 1 minute) variables meet these criteria. \n",
    "\n",
    "Its easy to select or exclude columsn or rows from a data frame. We can use the `which` command with an expression to get the column/row index and then select with the `[]` square bracket operators.\n",
    "\n",
    "The commands below look for the column name 'bmi' and get its index and then select it from the dataframe. Finally we remove the column from the data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "483d918a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ind=which(names(lung.data) == 'bmi')\n",
    "'Column index:'\n",
    "print(ind)\n",
    "'Column values rows 1 -> 10:'\n",
    "lung.data[1:10,ind]\n",
    "lung.data=lung.data[,-ind]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57185bf2",
   "metadata": {},
   "source": [
    "Check the columns in the data frame again to confirm the 'bmi' column is no longer present"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43605c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "names(lung.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81caebe0",
   "metadata": {},
   "source": [
    "Repeat the process for the FEV variable and check missing data values again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02721dbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ind=which(names(lung.data) == 'fev1pc_t0')\n",
    "lung.data=lung.data[,-ind]\n",
    "names(lung.data)\n",
    "gg_miss_var(lung.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbecddb8",
   "metadata": {},
   "source": [
    "### Filter to variables used in the paper\n",
    "The model created by the paper uses: age, gender, performance status, tumour volume, positive lymph node status, T stage, overall treatment time, cheomotherapy timing, and EQD2 dose to the tumour to predict patient survival.\n",
    "\n",
    "We will further filter the variable to only include these columns and assign to a new data frame containing the independent of explanatory factors 'lung.data.ind'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c6ae878",
   "metadata": {},
   "outputs": [],
   "source": [
    "testExpVars=c(\"gender\",\"age\",\"who3g\",\"countpet_all6g\",\"tstage\",\"nstage\",\"timing\",\"group\",\"eqd2\",\"ott\",\"gtv1\")\n",
    "lung.data.ind=lung.data[,which(names(lung.data) %in% testExpVars)]\n",
    "head(lung.data.ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16c44eaa",
   "metadata": {},
   "source": [
    "### Impute remaining missing data values\n",
    "We will impute the remaining missing values. The simplest way of doing this is to replace all missing values for each variable with the mean or median of that variable. We can do this by looping over each variable (column) and using the `mean()` or `median()` commands. \n",
    "\n",
    "1. Loop over each column using its index\n",
    "2. Check if there is any missing data by summing the number of N/A entries for each column using `is.na()`\n",
    "3. Replace the rows in that column where the entry is N/A with the median\n",
    "\n",
    "Note that there are many more sophistoicated imputation techniques available. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b72c6d8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for(i in 1:ncol(lung.data.ind)){\n",
    "  if(sum(is.na(lung.data.ind[,i]))>0){\n",
    "    lung.data.ind[is.na(lung.data.ind[,i]),i]=median(lung.data.ind[,i],na.rm=TRUE)\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d7624cc",
   "metadata": {},
   "source": [
    "This gives an error on the first run through:\n",
    "<img src=\"images/factorError.PNG\">\n",
    "This error is caused by trying to find the median of a factor - a categorical variable for which the concept of median has no meaning. We can check which by referring to the output of the `head()` command above. It is unlikely that 'age', 'eqd2' and 'gtv1' are supposed to be an ordinal variables. Again looking at the output of `head()` we can see that, probably as the data originates in Europe, the decimal points have been written in the csv as commas. We need to process the data to clean this up. We can convert these to numeric types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34b76e10",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.data.ind$age=as.numeric(gsub(',','.',lung.data.ind$age))\n",
    "lung.data.ind$eqd2=as.numeric(gsub(',','.',lung.data.ind$eqd2))\n",
    "lung.data.ind$gtv1=as.numeric(gsub(',','.',lung.data.ind$gtv1))\n",
    "head(lung.data.ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ccb678e",
   "metadata": {},
   "source": [
    "Try the imputation step above again and check to see if there are any missing data now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98d1391e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gg_miss_var(lung.data.ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33895e9c",
   "metadata": {},
   "source": [
    "We then need to change those variables that should be categorical and ordinal data type as necessary. An ordinal variable is a categorical variable in which the different categories are ordered (e.g. T Stage 1< 2 < 3 < 4). Ordinal variables can be assigned by specifying the order as part of the argument to the `factor()` function. Check the output of the head function again to notice the difference between nominal categorical (unordered) **\\<fct\\>** and ordinal **\\<ord\\>** variable types. You can refer to the variable descriptions given with the data on cancerdata.org to understand how the data is coded (e.g. N stage is coded as 1=N0, 2=N1, 3=N2, 4=N3, 5=Nx)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "948d1ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.data.ind$gender=factor(lung.data.ind$gender)\n",
    "lung.data.ind$who3g=factor(lung.data.ind$who3g, order=T, levels=c(0,1,2,3,4))\n",
    "lung.data.ind$tstage=factor(lung.data.ind$tstage, order=T, levels=c(0,1,2,3,4))\n",
    "lung.data.ind$nstage=factor(lung.data.ind$nstage, order=T, levels=c(0,1,2,3,4))\n",
    "lung.data.ind$timing=factor(lung.data.ind$timing)\n",
    "head(lung.data.ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "696fe943",
   "metadata": {},
   "source": [
    "## Dependent variables\n",
    "These variables are the explanatory variables that we hope to use to predict outcome. The dependent outcome we have in this dataset is survival time (alongside censoring information). \n",
    "\n",
    "We can examine the Kaplan-Meier survival curve using the `Surv()` and `survfit()` functions from the 'survial' package:\n",
    "1. Identify and extract the columns containing the survival and vital status data in the original dataset\n",
    "2. Clean the data as necessary\n",
    "3. Create and plot the survival object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a94cbcfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.data.dep=lung.data[,which(names(lung.data) %in% c(\"survmonth\",\"deadstat\"))]\n",
    "lung.data.dep$survmonth=as.numeric(gsub(',','.',lung.data.dep$survmonth))\n",
    "head(lung.data.dep)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d6df665",
   "metadata": {},
   "source": [
    "### Survival analysis\n",
    "Build survival model and check survival curve. Data is plotted using the `ggsurvplot()` function from the 'survminer' package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4883f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.survial=Surv(as.numeric(lung.data.dep$survmonth),lung.data.dep$deadstat)\n",
    "ggsurvplot(survfit(lung.survial~1,data=lung.data.ind),data=lung.data,risk.table=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb5292c",
   "metadata": {},
   "source": [
    "We can then assess the imfluence of different explanatory factors on patient survival by:\n",
    "1. Plotting different survival curves for patient subgroups\n",
    "2. Undertaking a survival analysis\n",
    "\n",
    "Below we both use the Cox proportional hazards model to estimate the hazard ratios associated with the each factor, and plot K-M curves dichotimised on different explanatory variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a17e5c",
   "metadata": {},
   "source": [
    "The cox model is run using the `coxph()` function from the 'survial' package. This can be univariable (ie assessing the impact of only 1 factor at a time). Below we see how we can do this for patient tumour volume. This shows us that there is a significant hazard ratio of 1.001 (i.e. increased hazard of death from baseline risk) per 1cm3 increase in GTV volume.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb9a40c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.uniCox=coxph(lung.survial~gtv1,data=lung.data.ind)\n",
    "print(summary(lung.uniCox))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a126769b",
   "metadata": {},
   "source": [
    "If we dicotomise GTV volume on some threshold (e.g. mean, median, or an absolute value of interest) we can then plot 2 K-M curves, one for patients with GTV volume smaller than the threshold vs patinets above this threshold. As expected, we see that high risk patients with larger GTV volume do worse. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfe9cf9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.data.dichot=lung.data.ind\n",
    "lung.data.dichot$gtvRisk = rep(0,nrow(lung.data.ind))\n",
    "lung.data.dichot$gtvRisk[which(lung.data.dichot$gtv1 > median(lung.data.dichot$gtv1))]=1\n",
    "ggsurvplot(survfit(lung.survial~gtvRisk,data=lung.data.dichot),data=lung.data.dichot,risk.table=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfa4f571",
   "metadata": {},
   "source": [
    "We can undertake univariable analysis for each variable we want to consider as part of a loop. For display reasons we assign a new data frame and extract the hazard and associated p-value to this for display once the loop has completed. As a technical note: where we consider ordinal factors, the cox model will give linear, quadratic and cubic estimates of how the hazard ratio changes with different factor levels with respect to the reference (e.g. lowest) level. Categorical variable will give one estimate per level (again, against a  reference). In the below example we take the linear estimate for ordinal variables. For simplicity of code we only include the first level estimate for the categorical 'timing' (timing of chemotherapy) in the display table. As an exercise you can modify the code below to also include the other factor levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc3880f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis.results = data.frame(variable=character(ncol(lung.data.ind)),univariableHR=numeric(ncol(lung.data.ind)),\n",
    "                            univariableP=numeric(ncol(lung.data.ind)),stringsAsFactors = FALSE)\n",
    "for(i in 1:ncol(lung.data.ind)){\n",
    "  lung.uniCox=coxph(as.formula(paste(\"lung.survial~\",names(lung.data.ind)[i])),data=lung.data.ind)\n",
    "  analysis.results$variable[i]=toString(names(lung.data.ind)[i])\n",
    "  analysis.results$univariableHR[i]=summary(lung.uniCox)$coefficients[1,2]\n",
    "  analysis.results$univariableP[i]=summary(lung.uniCox)$coefficients[1,5]\n",
    "}\n",
    "formattable(analysis.results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ee2216c",
   "metadata": {},
   "source": [
    "Finally we will undertake a multivariable analysis and see how these values compare to the univariable estimations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51c7a5e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "lung.multiCox=coxph(lung.survial~gender+age+who3g+countpet_all6g+tstage+nstage+timing+group+\n",
    "                      eqd2+ott+gtv1,data=lung.data.ind)\n",
    "print(lung.multiCox)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70b25f05",
   "metadata": {},
   "source": [
    "We can then add these data to the table to enable easy comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d0b780b",
   "metadata": {},
   "outputs": [],
   "source": [
    "multivariableHR=numeric(ncol(lung.data.ind))\n",
    "multivariableP=numeric(ncol(lung.data.ind))\n",
    "\n",
    "multivariableHR[1]=summary(lung.multiCox)$coefficients[1,2]\n",
    "multivariableP[1]=summary(lung.multiCox)$coefficients[1,5]\n",
    "\n",
    "multivariableHR[2]=summary(lung.multiCox)$coefficients[2,2]\n",
    "multivariableP[2]=summary(lung.multiCox)$coefficients[2,5]\n",
    "\n",
    "multivariableHR[3]=summary(lung.multiCox)$coefficients[3,2]\n",
    "multivariableP[3]=summary(lung.multiCox)$coefficients[3,5]\n",
    "\n",
    "multivariableHR[4]=summary(lung.multiCox)$coefficients[7,2]\n",
    "multivariableP[4]=summary(lung.multiCox)$coefficients[7,5]\n",
    "\n",
    "multivariableHR[5]=summary(lung.multiCox)$coefficients[8,2]\n",
    "multivariableP[5]=summary(lung.multiCox)$coefficients[8,5]\n",
    "\n",
    "multivariableHR[6]=summary(lung.multiCox)$coefficients[12,2]\n",
    "multivariableP[6]=summary(lung.multiCox)$coefficients[12,5]\n",
    "\n",
    "multivariableHR[7]=summary(lung.multiCox)$coefficients[16,2]\n",
    "multivariableP[7]=summary(lung.multiCox)$coefficients[16,5]\n",
    "\n",
    "multivariableHR[8]=summary(lung.multiCox)$coefficients[18,2]\n",
    "multivariableP[8]=summary(lung.multiCox)$coefficients[18,5]\n",
    "\n",
    "multivariableHR[9]=summary(lung.multiCox)$coefficients[19,2]\n",
    "multivariableP[9]=summary(lung.multiCox)$coefficients[19,5]\n",
    "\n",
    "multivariableHR[10]=summary(lung.multiCox)$coefficients[20,2]\n",
    "multivariableP[10]=summary(lung.multiCox)$coefficients[20,5]\n",
    "\n",
    "multivariableHR[11]=summary(lung.multiCox)$coefficients[21,2]\n",
    "multivariableP[11]=summary(lung.multiCox)$coefficients[21,5]\n",
    "\n",
    "analysis.results$multivariableHR=multivariableHR\n",
    "analysis.results$multivariableP=multivariableP\n",
    "\n",
    "formattable(analysis.results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3fe237c",
   "metadata": {},
   "source": [
    "As can be seen there are now differences between the hazard ratio estimations for each factor. In particular:\n",
    "* Age is identified as much less likely to have significant influence on survival in the multivariable analysis\n",
    "* Overall treatment time is found to be more influential in the multivariable analysis (whereas previously we might have dismissed this as a relevent factor)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}