Basil_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the README of this package I demonstrated how to use the multisnpnet-Cox-tv package to fit a Cox model with user-provided regularization parameters and with the time-varying covariates provided in the form of an operation time. Here I will show how to use this package to fit a Cox-Lasso path using BASIL to screen the SNPs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "── \u001b[1mAttaching packages\u001b[22m ─────────────────────────────────────── tidyverse 1.3.0 ──\n",
      "\n",
      "\u001b[32m✔\u001b[39m \u001b[34mggplot2\u001b[39m 3.3.0     \u001b[32m✔\u001b[39m \u001b[34mpurrr  \u001b[39m 0.3.3\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtibble \u001b[39m 2.1.3     \u001b[32m✔\u001b[39m \u001b[34mdplyr  \u001b[39m 0.8.5\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtidyr  \u001b[39m 1.0.2     \u001b[32m✔\u001b[39m \u001b[34mstringr\u001b[39m 1.4.0\n",
      "\u001b[32m✔\u001b[39m \u001b[34mreadr  \u001b[39m 1.3.1     \u001b[32m✔\u001b[39m \u001b[34mforcats\u001b[39m 0.4.0\n",
      "\n",
      "── \u001b[1mConflicts\u001b[22m ────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mbetween()\u001b[39m   masks \u001b[34mdata.table\u001b[39m::between()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m    masks \u001b[34mstats\u001b[39m::filter()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfirst()\u001b[39m     masks \u001b[34mdata.table\u001b[39m::first()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m       masks \u001b[34mstats\u001b[39m::lag()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlast()\u001b[39m      masks \u001b[34mdata.table\u001b[39m::last()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mpurrr\u001b[39m::\u001b[32mtranspose()\u001b[39m masks \u001b[34mdata.table\u001b[39m::transpose()\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(coxtv)\n",
    "library(pgenlibr)\n",
    "library(data.table)\n",
    "library(tidyverse)\n",
    "\n",
    "phe.file = \"@@@/phenotypedata/master_phe/cox/phenotypefiles/f.131298.0.0.phe\"\n",
    "death.file = \"@@@/ukbb-phenotyping/20200404_icd_death/ukb41413_icd_death.tsv\"\n",
    "masterphe.file = \"@@@/phenotypedata/master_phe/master.phe\"\n",
    "genotype.pfile = \"@@@/array_combined/pgen/ukb24983_cal_hla_cnv\"\n",
    "psamid = data.table::fread(paste0(genotype.pfile, '.psam'),colClasses = list(character=c(\"IID\")), \n",
    "                  select = c(\"IID\"))\n",
    "psamid = psamid$IID\n",
    "\n",
    "covs = c(\"sex\", \"age\", \"PC1\", \"PC2\", \"PC3\", \"PC4\", \"PC5\", \"PC6\", \"PC7\", \"PC8\", \"PC9\", \"PC10\")\n",
    "\n",
    "train_ratio = 0.8\n",
    "\n",
    "\n",
    "configs = list()\n",
    "configs[['gcount.full.prefix']] = '@@/ruilinli/tvtest/gcount/test'\n",
    "configs[['plink2.path']] = \"@@/ruilinli/prox_grad_cox_block/plink2\"\n",
    "configs[['nCores']] = 6\n",
    "configs[['mem']] = 60000\n",
    "configs[['vzs']] = TRUE\n",
    "configs[['save']] =TRUE\n",
    "configs[['zstdcat.path']] = \"zstdcat\"\n",
    "configs[['save.computeProduct']] = TRUE\n",
    "configs[['results.dir']] = \"@@/ruilinli/tvtest/result/\"\n",
    "configs[['save.dir']] = \"@@/ruilinli/tvtest/save\"\n",
    "configs[['KKT.verbose']] = TRUE\n",
    "configs[['endian']]=\"little\"\n",
    "configs[[\"standardize.variant\"]] = FALSE\n",
    "configs[['missing.rate']] = 0.1\n",
    "configs[['MAF.thresh']] = 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some preprocessing\n",
    "phe = data.table::fread(phe.file, \n",
    "                  colClasses = list(character=c(\"FID\"), numeric=c(\"coxnet_y_f.131298.0.0\", \"coxnet_status_f.131298.0.0\")), \n",
    "                  select = c(\"FID\", \"coxnet_y_f.131298.0.0\", \"coxnet_status_f.131298.0.0\"))\n",
    "names(phe) = c(\"ID\", \"t0\", \"MI\")\n",
    "phe = filter(phe, ID %in% psamid)\n",
    "\n",
    "event = data.table::fread(death.file, \n",
    "                  colClasses = list(character=c(\"#IID\"), numeric=c(\"val\")), \n",
    "                  select = c(\"#IID\", \"val\"))\n",
    "names(event) = c(\"ID\", \"val\")\n",
    "\n",
    "masterphe = data.table::fread(masterphe.file, \n",
    "                  colClasses = list(character=c(\"FID\"), numeric=covs), \n",
    "                  select = c(\"FID\", covs))\n",
    "names(masterphe) = c(\"ID\", covs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "phe = filter(phe, MI == 1) # people with MI\n",
    "phe = select(phe, -MI) # Don't need the MI indicator anymore, since every in the dataset had MI\n",
    "phe$status = 0\n",
    "\n",
    "\n",
    "eID = which(phe$ID %in% event$ID) # people with MI that also had event = death\n",
    "phe$status[eID] = 1\n",
    "phe$t1 = NA # The event time\n",
    "phe$t1[eID] = event$val[match(phe$ID[eID], event$ID)]\n",
    "phe$t1[-eID] = masterphe$age[match(phe$ID[-eID], masterphe$ID)] + 1.1219 # last followup time for people who did not have the event\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add covariates to the phe file\n",
    "covs = covs[covs!=\"age\"]\n",
    "for (cov in covs){\n",
    "        phe[,cov] <- masterphe[match(phe$ID, masterphe$ID), ..cov]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "phe$y = phe$t1 - phe$t0\n",
    "phe = filter(phe, t1 > 0)\n",
    "phe = phe[complete.cases(phe),]\n",
    "min_event_time = min(phe$y[as.logical(phe$status)])\n",
    "phe = filter(phe, y>= max(0, min_event_time)) # non-events before the first event will never be used\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0.2% 99.9% \n",
      "  0.0   6.9 \n",
      " 0.2% 99.9% \n",
      "    0   165 \n"
     ]
    }
   ],
   "source": [
    "# Now we get the time-varying covariates into the format we need\n",
    "tv_files = c(\"@@@/primary_care/gp_clinical/example_phenotypes/final/LDL.tsv\",\n",
    "            \"@@@/primary_care/gp_clinical/example_phenotypes/final/Weight.tsv\")\n",
    "\n",
    "tv_list = list()\n",
    "i = 1\n",
    "for(tvfile in tv_files){\n",
    "    tv = data.table::fread(tvfile, \n",
    "                  colClasses = list(character=c(\"id\"), numeric=c(\"age\", \"value\")), \n",
    "                  select = c(\"id\", \"age\", \"value\"))\n",
    "    names(tv) = c(\"ID\", \"time\", \"value\")\n",
    "    tv = filter(tv, !is.na(value))\n",
    "    bounds = quantile(tv$value, c(0.002, 0.999)) # remove extreme observations\n",
    "    print(bounds) \n",
    "    tv = filter(tv, (value > bounds[1]) & (value < bounds[2]))\n",
    "    # take the intersection\n",
    "    phe = filter(phe, ID %in% tv$ID)\n",
    "    tv = filter(tv, ID %in% phe$ID)\n",
    "    tv$time = tv$time - phe$t0[match(tv$ID, phe$ID)]\n",
    "    tv_list[[i]] = tv\n",
    "    i = i + 1\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in coxtv::get_info(phe_train, tv_train):\n",
      "“2438 people do not have time-varying covariates measured before the first event. The most recent measurement after the event is used.”\n",
      "Warning message in coxtv::get_info(phe_train, tv_train):\n",
      "“1474 people do not have time-varying covariates measured before the first event. The most recent measurement after the event is used.”\n",
      "Warning message in coxtv::get_info(phe_val, tv_val):\n",
      "“599 people do not have time-varying covariates measured before the first event. The most recent measurement after the event is used.”\n",
      "Warning message in coxtv::get_info(phe_val, tv_val):\n",
      "“363 people do not have time-varying covariates measured before the first event. The most recent measurement after the event is used.”\n"
     ]
    }
   ],
   "source": [
    "# We split the data into training and validation set, the original split column may not be appropriate here\n",
    "# since we are using a small subset of founders (because only a small number of people had MI)\n",
    "total_events = sum(phe$status)\n",
    "non_events = nrow(phe) - total_events\n",
    "train_id = sample(phe$ID[as.logical(phe$status)], round(total_events*train_ratio))\n",
    "train_id = c(train_id, sample(phe$ID[!as.logical(phe$status)], round(non_events*train_ratio)))\n",
    "val_id = filter(phe, ! ID %in% train_id)$ID\n",
    "\n",
    "phe$split = 'train'\n",
    "phe$split[phe$ID %in% val_id] = 'val'\n",
    "\n",
    "phe_val = as.data.table(filter(phe, split=='val'))\n",
    "phe_train = as.data.table(filter(phe, split=='train'))\n",
    "\n",
    "\n",
    "\n",
    "tv_train = list()\n",
    "tv_val = list()\n",
    "for(i in (1:length(tv_list))){\n",
    "    tv_train[[i]] = filter(tv_list[[i]], ID %in% phe_train$ID)\n",
    "    tv_val[[i]] = filter(tv_list[[i]], ID %in% phe_val$ID)\n",
    "}\n",
    "info_train = coxtv::get_info(phe_train, tv_train)\n",
    "info_val = coxtv::get_info(phe_val, tv_val)\n",
    "\n",
    "\n",
    "rm(tv_list)\n",
    "rm(phe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the data in the right format to be fed to the coxtv functions. To summarize, we need\n",
    "- A dataframe (here phe_train and phe_val) that has columns:\n",
    "    - y that contains the time-to-event response\n",
    "    - status, a binary vector that represents whether event has occured\n",
    "    - ID, which will be used to identify each person in this dataframe\n",
    "    - some covariates columns that are time independent\n",
    "- A list that contains the time-varying covariates (here tv_train and tv_val). Each element of the list corresponds to one time-varying covariate and must have the columns:\n",
    "    - ID, same set of IDs as used in the phe data frame\n",
    "    - time, the time at which the measurement was taken (relative to each person's t0)\n",
    "    - value, the value of the measurements at the corresponding time\n",
    "- A set of time-independent covariates names that will be used to fit a Cox model, these names must be available in the phe dataframet \n",
    "- (Optional, to save some compute) A list info that is obtained using coxtv::get_info(phe, tv_list). This list contains information about the time-varying covariates in the form that coxtv can readily use to fit a model. Keep info can save some computation, especially when n or the number of events is large. For now the user needs to make sure that when fitting a Cox model, info must correspond to phe and tv_list. Otherwise a segmentation fault might happen. In future version I might encapsulate BASIL into a package to solve this problem.\n",
    "\n",
    "##### It is important that the users make sure that each person must have at least one measurement before the first event time, for each of the time-varying covariate! If this is not satisfied, a warning will be thrown and the most recent measurement after the event will be used.\n",
    "\n",
    "##### Now let's fit the first iteration, which is supposed to be unpenalized:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "covs = c(\"t0\", covs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = coxtv(phe_train, NULL, covs, c(0.0), info=info_train) #lambda is the last argument, if info is provided then the second parameter is not needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A matrix: 14 × 1 of type dbl</caption>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>TV1</th><td> 0.0022575249</td></tr>\n",
       "\t<tr><th scope=row>TV2</th><td>-0.0008927381</td></tr>\n",
       "\t<tr><th scope=row>t0</th><td> 0.1277776373</td></tr>\n",
       "\t<tr><th scope=row>sex</th><td> 0.3273922625</td></tr>\n",
       "\t<tr><th scope=row>PC1</th><td> 0.0022498610</td></tr>\n",
       "\t<tr><th scope=row>PC2</th><td> 0.0032595022</td></tr>\n",
       "\t<tr><th scope=row>PC3</th><td>-0.0097584263</td></tr>\n",
       "\t<tr><th scope=row>PC4</th><td>-0.0229086856</td></tr>\n",
       "\t<tr><th scope=row>PC5</th><td>-0.0109442029</td></tr>\n",
       "\t<tr><th scope=row>PC6</th><td>-0.0467834069</td></tr>\n",
       "\t<tr><th scope=row>PC7</th><td>-0.0128030432</td></tr>\n",
       "\t<tr><th scope=row>PC8</th><td> 0.0491001304</td></tr>\n",
       "\t<tr><th scope=row>PC9</th><td>-0.0101109565</td></tr>\n",
       "\t<tr><th scope=row>PC10</th><td> 0.0126472443</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A matrix: 14 × 1 of type dbl\n",
       "\\begin{tabular}{r|l}\n",
       "\tTV1 &  0.0022575249\\\\\n",
       "\tTV2 & -0.0008927381\\\\\n",
       "\tt0 &  0.1277776373\\\\\n",
       "\tsex &  0.3273922625\\\\\n",
       "\tPC1 &  0.0022498610\\\\\n",
       "\tPC2 &  0.0032595022\\\\\n",
       "\tPC3 & -0.0097584263\\\\\n",
       "\tPC4 & -0.0229086856\\\\\n",
       "\tPC5 & -0.0109442029\\\\\n",
       "\tPC6 & -0.0467834069\\\\\n",
       "\tPC7 & -0.0128030432\\\\\n",
       "\tPC8 &  0.0491001304\\\\\n",
       "\tPC9 & -0.0101109565\\\\\n",
       "\tPC10 &  0.0126472443\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A matrix: 14 × 1 of type dbl\n",
       "\n",
       "| TV1 |  0.0022575249 |\n",
       "| TV2 | -0.0008927381 |\n",
       "| t0 |  0.1277776373 |\n",
       "| sex |  0.3273922625 |\n",
       "| PC1 |  0.0022498610 |\n",
       "| PC2 |  0.0032595022 |\n",
       "| PC3 | -0.0097584263 |\n",
       "| PC4 | -0.0229086856 |\n",
       "| PC5 | -0.0109442029 |\n",
       "| PC6 | -0.0467834069 |\n",
       "| PC7 | -0.0128030432 |\n",
       "| PC8 |  0.0491001304 |\n",
       "| PC9 | -0.0101109565 |\n",
       "| PC10 |  0.0126472443 |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1]         \n",
       "TV1   0.0022575249\n",
       "TV2  -0.0008927381\n",
       "t0    0.1277776373\n",
       "sex   0.3273922625\n",
       "PC1   0.0022498610\n",
       "PC2   0.0032595022\n",
       "PC3  -0.0097584263\n",
       "PC4  -0.0229086856\n",
       "PC5  -0.0109442029\n",
       "PC6  -0.0467834069\n",
       "PC7  -0.0128030432\n",
       "PC8   0.0491001304\n",
       "PC9  -0.0101109565\n",
       "PC10  0.0126472443"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "as.matrix(result[[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's fit a Lasso path using BASIL to screen the SNPs. First we compute the residual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "source(\"@@@/ruilinli/prox_grad_cox_block/snpnet/R/functions.R\")\n",
    "# Need to use some snpnet helper functions\n",
    "\n",
    "\n",
    "residuals = cox_residual(phe_train, covs, info_train, result[[1]]) # 1/n already multiplied in this residual\n",
    "# to compute the gradient (with respect to the SNP coefficients) let's first load the genotype files\n",
    "# code copied from snpnet\n",
    "vars <- dplyr::mutate(dplyr::rename(data.table::fread(cmd=paste0(configs[['zstdcat.path']], ' ', paste0(genotype.pfile, '.pvar.zst'))), 'CHROM'='#CHROM'), VAR_ID=paste(ID, ALT, sep='_'))$VAR_ID\n",
    "pvar <- pgenlibr::NewPvar(paste0(genotype.pfile, '.pvar.zst'))\n",
    "pgen_train = pgenlibr::NewPgen(paste0(genotype.pfile, '.pgen'), pvar=pvar, sample_subset=match(phe_train$ID, psamid))\n",
    "pgen_val = pgenlibr::NewPgen(paste0(genotype.pfile, '.pgen'), pvar=pvar, sample_subset=match(phe_val$ID, psamid))\n",
    "\n",
    "\n",
    "pgenlibr::ClosePvar(pvar)    \n",
    "\n",
    "stats <- computeStats(genotype.pfile, paste(phe_train$ID, phe_train$ID, sep=\"_\"), configs = configs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2020-04-30 13:27:43 snpnet]     Start computeProduct()\n",
      "           used  (Mb) gc trigger  (Mb) max used  (Mb)\n",
      "Ncells  2623258 140.1    7186685 383.9  4123368 220.3\n",
      "Vcells 31143942 237.7   56498432 431.1 49372753 376.7\n",
      "[2020-04-30 13:27:43 snpnet]       Start plink2 --variant-score\n",
      "[2020-04-30 13:30:22 snpnet]         End plink2 --variant-score. Time elapsed: 2.6477 mins\n",
      "[2020-04-30 13:30:22 snpnet]       End computeProduct(). Time elapsed: 2.6503 mins\n"
     ]
    }
   ],
   "source": [
    "residuals = matrix(residuals,nrow = length(phe_train$ID), ncol = 1, dimnames = list(paste(phe_train$ID, phe_train$ID, sep='_'), c(\"0\")))\n",
    "gradient =computeProduct(residuals, genotype.pfile, vars, stats, configs, iter=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can generate a lambda sequence\n",
    "nlambda = 100\n",
    "lambda.min.ratio = 0.01\n",
    "lambda_max = max(abs(gradient), na.rm = NA)\n",
    "lambda_min = lambda_max * lambda.min.ratio\n",
    "lambda_seq = exp(seq(from = log(lambda_max), to = log(lambda_min), length.out = nlambda))\n",
    "\n",
    "# The first lambda solution is already obtained\n",
    "max_valid_index = 1\n",
    "prev_valid_index = 0\n",
    "\n",
    "# Use validation C-index to determine early stop\n",
    "max_cindex = 0\n",
    "cindex = numeric(nlambda)\n",
    "out = list()\n",
    "out[[1]] = result[[1]]\n",
    "features.to.discard = NULL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [1] \"t0\"   \"sex\"  \"PC1\"  \"PC2\"  \"PC3\"  \"PC4\"  \"PC5\"  \"PC6\"  \"PC7\"  \"PC8\" \n",
      "[11] \"PC9\"  \"PC10\"\n"
     ]
    }
   ],
   "source": [
    "score = abs(gradient[,1])\n",
    "iter = 1\n",
    "ever.active = covs\n",
    "print(ever.active)\n",
    "current_B = result[[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"current maximum valid index is: 1\"\n",
      "[1] \"Current C-Indices are:\"\n",
      "[1] 0.7681487\n",
      "[1] \"Number of variables to be fitted is: 1014\"\n",
      "[2020-04-30 13:32:54 snpnet]     Start computeProduct()\n",
      "           used  (Mb) gc trigger  (Mb) max used  (Mb)\n",
      "Ncells  2639400 141.0    7186685 383.9  5088656 271.8\n",
      "Vcells 41945038 320.1   81533741 622.1 80496732 614.2\n",
      "[2020-04-30 13:32:54 snpnet]       Start plink2 --variant-score\n",
      "[2020-04-30 13:35:12 snpnet]         End plink2 --variant-score. Time elapsed: 2.2911 mins\n",
      "[2020-04-30 13:35:13 snpnet]       End computeProduct(). Time elapsed: 2.3099 mins\n",
      "lambda_idx_0.0155 lambda_idx_0.0148 lambda_idx_0.0141 lambda_idx_0.0135 \n",
      "      0.009576207       0.009633053       0.009704023       0.009845673 \n",
      "lambda_idx_0.0129 lambda_idx_0.0123 lambda_idx_0.0117 lambda_idx_0.0112 \n",
      "      0.010046540       0.010166109       0.010284680       0.010346207 \n",
      "lambda_idx_0.0107 lambda_idx_0.0102 \n",
      "      0.010277798       0.010504970 \n",
      "[1] \"Number of features discarded in this iteration is 812\"\n",
      "[1] \"early stop reached\"\n"
     ]
    }
   ],
   "source": [
    "while(max_valid_index < nlambda){\n",
    "    if(max_valid_index > prev_valid_index){\n",
    "        for(i in (prev_valid_index + 1):max_valid_index){\n",
    "            cindex[i] = cindex_tv(phe_val, NULL, covs, out[[i]], info=info_val)\n",
    "        }\n",
    "        \n",
    "        cindex_this_iter = cindex[(prev_valid_index + 1):max_valid_index]\n",
    "        \n",
    "        max_cindex_this_iter = max(cindex_this_iter)\n",
    "        if(max_cindex_this_iter >= max_cindex){\n",
    "            max_cindex = max_cindex_this_iter\n",
    "        } else {\n",
    "            print(\"early stop reached\")\n",
    "            break\n",
    "        }\n",
    "        \n",
    "        if(which.max(cindex_this_iter) != length(cindex_this_iter)){\n",
    "            print(\"early stop reached\")\n",
    "            break            \n",
    "        }\n",
    "    }\n",
    "    prev_valid_index = max_valid_index\n",
    "    print(paste(\"current maximum valid index is:\",max_valid_index ))\n",
    "    print(\"Current C-Indices are:\")\n",
    "    print(cindex[1:max_valid_index])\n",
    "\n",
    "    if(length(features.to.discard) > 0){\n",
    "        phe_train[, (features.to.discard) := NULL]\n",
    "        phe_val[, (features.to.discard) := NULL]\n",
    "        covs = covs[!covs %in% features.to.discard] # the name is a bit confusing, maybe change it to ti_names?\n",
    "    }\n",
    "    \n",
    "    which.in.model <- which(names(score) %in% covs)\n",
    "    score[which.in.model] <- NA\n",
    "    sorted.score <- sort(score, decreasing = T, na.last = NA)\n",
    "    features.to.add <- names(sorted.score)[1:min(1000, length(sorted.score))]\n",
    "    covs = c(covs, features.to.add)\n",
    "    B_init = c(current_B, rep(0.0, length(features.to.add)))\n",
    "    \n",
    "    tmp.features.add <- prepareFeatures(pgen_train, vars, features.to.add, stats)\n",
    "    phe_train[, colnames(tmp.features.add) := tmp.features.add]\n",
    "    \n",
    "    tmp.features.add <- prepareFeatures(pgen_val, vars, features.to.add, stats)\n",
    "    phe_val[, colnames(tmp.features.add) := tmp.features.add]\n",
    "    \n",
    "    rm(tmp.features.add)\n",
    "    \n",
    "    # Not fit a regularized Cox model for the next 10 lambdas\n",
    "    lambda_seq_local = lambda_seq[(max_valid_index + 1):min(max_valid_index + 10, length(lambda_seq))]\n",
    "    # Need better ways to set p.fac\n",
    "    p.fac = rep(1, length(B_init))\n",
    "    p.fac[1:14] = 0.0\n",
    "    print(paste(\"Number of variables to be fitted is:\",length(B_init)))\n",
    "\n",
    "\n",
    "    result = coxtv(phe_train, NULL, covs, lambda_seq_local, B0 = B_init, p.fac = p.fac, info=info_train)\n",
    "    \n",
    "    residuals = matrix(nrow = length(phe_train$ID), ncol = length(lambda_seq_local), \n",
    "                       dimnames = list(paste(phe_train$ID, phe_train$ID, sep='_'), signif(lambda_seq_local, 3)))\n",
    "    for(i in 1:length(result)){\n",
    "        residuals[,i] = cox_residual(phe_train, covs, info_train, result[[i]])\n",
    "    }\n",
    "    new_score = abs(computeProduct(residuals, genotype.pfile, vars, stats, configs, iter=iter))\n",
    "    max_score = apply(new_score, 2, function(x){max(x[!names(x) %in% covs], na.rm=NA)})\n",
    "    print(max_score)\n",
    "    # if all failed\n",
    "    if(all(max_score > lambda_seq_local)){\n",
    "        features.to.discard = NULL\n",
    "        current_B = result[[1]]\n",
    "        score = new_score[, 1]\n",
    "    } else {\n",
    "        local_valid = which.min(c(max_score <= lambda_seq_local, FALSE)) - 1 # number of valid this iteration\n",
    "        \n",
    "        for(j in 1:local_valid){\n",
    "            out[[max_valid_index+j]] = result[[j]]\n",
    "        }\n",
    "        \n",
    "        max_valid_index = max_valid_index + local_valid\n",
    "        ever.active <- union(ever.active, names(which(apply(sapply(result[1:local_valid, drop=F], function(x){x!=0}), 1, any))))\n",
    "        features.to.discard = setdiff(covs, ever.active)\n",
    "        score = new_score[, local_valid]\n",
    "        current_B = result[[local_valid]]\n",
    "        current_B = current_B[!names(current_B) %in% features.to.discard]\n",
    "        print(paste(\"Number of features discarded in this iteration is\", length(features.to.discard)))\n",
    "    }\n",
    "    iter = iter + 1\n",
    "\n",
    "}"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}