IHT/Mendel_IHT_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OpenMendel Tutorial on Iterative Hard Thresholding\n",
    "\n",
    "### Last update: 10/16/2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Julia version\n",
    "\n",
    "`MendelIHT.jl` currently supports Julia 1.0 and 1.2 on Mac, Linux, and Windows, but it currently an unregistered package. To install, press `]` to invoke the package manager mode and install these packages by typing:\n",
    "\n",
    "```\n",
    "add https://github.com/OpenMendel/SnpArrays.jl\n",
    "add https://github.com/OpenMendel/MendelSearch.jl\n",
    "add https://github.com/OpenMendel/MendelBase.jl\n",
    "add https://github.com/OpenMendel/MendelIHT.jl\n",
    "```\n",
    "\n",
    "For this tutorial you will also need a few registered packages. Add them by typing:\n",
    "\n",
    "```\n",
    "add DataFrames, Distributions, BenchmarkTools, Random, LinearAlgebra, GLM\n",
    "```\n",
    "\n",
    "For reproducibility, the computer spec and Julia version is listed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 1.2.0\n",
      "Commit c6da87ff4b (2019-08-20 00:03 UTC)\n",
      "Platform Info:\n",
      "  OS: macOS (x86_64-apple-darwin18.6.0)\n",
      "  CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz\n",
      "  WORD_SIZE: 64\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Documentation\n",
    "\n",
    "This is a short tutorial. For detailed documentation, please visit our [home page](https://openmendel.github.io/MendelIHT.jl/latest/). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When to use Iterative Hard Thresholding\n",
    "\n",
    "Continuous model selection is advantageous in situations where the multivariate nature of the regressors plays a significant role together. As an alternative to traditional SNP-by-SNP association testing, iterative hard-thresholing (IHT) performs continuous model selection on a GWAS dataset $\\mathbf{X} \\in \\{0, 1, 2\\}^{n \\times p}$ and continuous phenotype vector $\\mathbf{y}$ by maximizing the loglikelihood $L(\\beta)$ subject to the constraint that $\\beta$ is $k-$sparse. This method has the edge over LASSO because IHT does not shrink estimated effect sizes. Parallel computing is offered through `q-`fold cross validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Appropriate Datasets and Example Inputs \n",
    "\n",
    "All genotype data **must** be stored in the [PLINK binary genotype format](https://www.cog-genomics.org/plink2/formats#bed) where at least the triplets `.bim`, `.bed` and `.fam` must all be present. Additional non-genetic covariates should be imported separately by the user. In the examples below, we first simulate phenotypes from the Normal, Bernoulli, Poisson, and Negative Binomial family, and then attempt to fit the corresponding model using our IHT implementation. We can examine reconstruction behavior as well as the ability for cross validation to find the true sparsity parameter.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing Data\n",
    "\n",
    "`MendelIHT` assumes there are no missing genotypes, since it uses linear algebra functions defined in [SnpArrays.jl](https://openmendel.github.io/SnpArrays.jl/latest/man/snparray/#linear-algebra-with-snparray). Therefore, you must first impute missing genotypes before you use MendelIHT. SnpArrays.jl offer basic quality control routines such as filtering, but otherwise, our own software [option 23 of Mendel](http://software.genetics.ucla.edu/download?package=1) is a reasonable choice. Open Mendel will soon provide a separate package `MendelImpute.jl` containing new imputation strategies such as alternating least squares.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross Validation and Regularization paths\n",
    "\n",
    "We usually have very little information on how many SNPs are affecting the phenotype. In a typical GWAS study, anywhere between 1 to thousands of SNPs could play a role. Thus ideally, we can test many different models to find the best one. MendelIHT provides 2 ways for one to perform this automatically: user specified regulartization paths, and $q-$fold [cross validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics). Users should know that, in the first method, increasing the number of predictors will almost always decrease the error, but as a result introduce overfitting. Therefore, in most practical situations, it is highly recommended to combine this method with cross validation. In $q-$fold cross validation, samples are divided into $q$ disjoint subsets, and IHT fits a model on $q-1$ of those sets data, then computes the [MSE](https://en.wikipedia.org/wiki/Mean_squared_error) tested on the $qth$ samples. Each $q$ subsets are served as the test set exactly once. This functionality of `MendelIHT.jl` natively supports parallel computing. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1: Quantitative Traits\n",
    "\n",
    "Quantitative traits are continuous phenotypes that can essentially take on any real number. In this example, we first simulate $y_i \\sim x_i^T\\beta + \\epsilon_i$ where $\\epsilon_i \\sim N(0, 1)$ and $\\beta_i \\sim N(0, 1)$. Then using just the genotype matrix $X$ and phenotype vector $y$, we use IHT to recover the simulated $\\beta$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       " 2\n",
       " 3\n",
       " 4\n",
       " 5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#first add workers needed for parallel computing. Add only as many CPU cores you have \n",
    "using Distributed\n",
    "addprocs(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#load necessary packages\n",
    "using MendelIHT\n",
    "using SnpArrays\n",
    "using DataFrames\n",
    "using Distributions\n",
    "using BenchmarkTools\n",
    "using Random\n",
    "using LinearAlgebra\n",
    "using GLM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Simulat data with k true predictors, from distribution d and with link l."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "IdentityLink()"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 1000\n",
    "p = 10000\n",
    "k = 10\n",
    "d = Normal\n",
    "l = canonicallink(d())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Construct snpmatrix, covariate files, and true model b\n",
    "\n",
    "The SnpBitMatrix type (`xbm` below) is necessary for performing linear algebra directly on raw genotype files without expanding the matrix to numeric floating points. Here `undef` in the 3rd argument simply indicates that the matrix `x` will be stored in RAM. Please visit SnpArrays' [documentation](https://openmendel.github.io/SnpArrays.jl/latest/) for more detailed description. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Random.seed!(1111) #set random seed\n",
    "x = simulate_random_snparray(n, p, undef)\n",
    "xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true); \n",
    "z = ones(n, 1); # only nongenetic covariate is the intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Simulate response y, true model b, and the correct non-0 positions of b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([-2.02168310584216, -2.562563581153803, 1.2438984775887258, 0.30434816265962317, 1.7043321478245017, -2.7754953149035013, -0.9486637262536461, 0.1660538100586053, 1.5880052223831773, 1.033229535854426  …  -2.1664956144610645, 7.975518350862309, 0.32430592236910805, 1.6057303922701498, 1.5909304105723248, -2.503956200147935, -3.4652255843587714, -0.34640336293302176, 1.0706692918688652, 0.29268585188820095], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [2384, 3352, 4093, 5413, 5455, 6729, 7403, 8753, 9089, 9132])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Run IHT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "IHT estimated 10 nonzero SNP predictors and 0 non-genetic predictors.\n",
       "\n",
       "Compute time (sec):     0.6195089817047119\n",
       "Final loglikelihood:    -1407.2533232402275\n",
       "Iterations:             12\n",
       "Max number of groups:   1\n",
       "Max predictors/group:   10\n",
       "\n",
       "Selected genetic predictors:\n",
       "10×2 DataFrame\n",
       "│ Row │ Position │ Estimated_β │\n",
       "│     │ \u001b[90mInt64\u001b[39m    │ \u001b[90mFloat64\u001b[39m     │\n",
       "├─────┼──────────┼─────────────┤\n",
       "│ 1   │ 2384     │ -1.26014    │\n",
       "│ 2   │ 3352     │ -0.26742    │\n",
       "│ 3   │ 3353     │ 0.141208    │\n",
       "│ 4   │ 4093     │ 0.289956    │\n",
       "│ 5   │ 5413     │ 0.366689    │\n",
       "│ 6   │ 5609     │ -0.137181   │\n",
       "│ 7   │ 7403     │ -0.308255   │\n",
       "│ 8   │ 8753     │ 0.332881    │\n",
       "│ 9   │ 9089     │ 0.964598    │\n",
       "│ 10  │ 9132     │ -0.509461   │\n",
       "\n",
       "Selected nongenetic predictors:\n",
       "0×2 DataFrame\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = L0_reg(x, xbm, z, y, 1, k, d(), l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Check results\n",
    "\n",
    "IHT found 8/10 predictors in this example. The 2 that was not found had a relatively small effect size, and as far as IHT can tell, they are indistinguishable from noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "compare_model = 10×2 DataFrame\n",
      "│ Row │ true_β    │ estimated_β │\n",
      "│     │ Float64   │ Float64     │\n",
      "├─────┼───────────┼─────────────┤\n",
      "│ 1   │ -1.19376  │ -1.26014    │\n",
      "│ 2   │ -0.230351 │ -0.26742    │\n",
      "│ 3   │ 0.257181  │ 0.289956    │\n",
      "│ 4   │ 0.344827  │ 0.366689    │\n",
      "│ 5   │ 0.155484  │ 0.0         │\n",
      "│ 6   │ -0.126114 │ 0.0         │\n",
      "│ 7   │ -0.286079 │ -0.308255   │\n",
      "│ 8   │ 0.327039  │ 0.332881    │\n",
      "│ 9   │ 0.931375  │ 0.964598    │\n",
      "│ 10  │ -0.496683 │ -0.509461   │\n",
      "Total iteration number was 12\n",
      "Total time was 0.6195089817047119\n",
      "Total found predictors = 8\n"
     ]
    }
   ],
   "source": [
    "compare_model = DataFrame(\n",
    "    true_β      = true_b[correct_position], \n",
    "    estimated_β = result.beta[correct_position])\n",
    "@show compare_model\n",
    "println(\"Total iteration number was \" * string(result.iter))\n",
    "println(\"Total time was \" * string(result.time))\n",
    "println(\"Total found predictors = \" * string(length(findall(!iszero, result.beta[correct_position]))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: Case-control study controlling for sex\n",
    "\n",
    "Case control studies are used when the phenotype in a binary count data. In this example, we simulate a case-control study, while controling for sex as a non-genetic covariate. \n",
    "\n",
    "The exact simulation code to generate the phenotype $y$ can be found at: https://github.com/biona001/MendelIHT.jl/blob/master/src/simulate_utilities.jl#L107"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Simulat data with k true predictors, from distribution d and with link l."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogitLink()"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 1000\n",
    "p = 10000\n",
    "k = 10\n",
    "d = Bernoulli\n",
    "l = canonicallink(d())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: construct snpmatrix, covariate files, and true model b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Random.seed!(1111) #set random seed \n",
    "x = simulate_random_snparray(n, p, undef)\n",
    "xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true); \n",
    "z = ones(n, 2) # first column is the intercept, second column the sex. 0 = male 1 = female\n",
    "z[:, 2] .= rand(0:1, n);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: simulate true models \n",
    "\n",
    "Here we used $k=8$ genetic predictors and 2 non-genetic predictors (intercept and sex). The simulation code in our package does not yet handle simulations with non-genetic predictors, so we must simulate these phenotypes manually. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 1.0\n",
       " 1.5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_b = zeros(p) #genetic predictors\n",
    "true_b[1:k-2] = randn(k-2)\n",
    "shuffle!(true_b)\n",
    "correct_position = findall(!iszero, true_b)\n",
    "true_c = [1.0; 1.5] #non-genetic predictors: intercept & sex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: simulate phenotype using genetic and nongenetic predictors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "prob = GLM.linkinv.(l, xbm * true_b .+ z * true_c)\n",
    "y = [rand(d(i)) for i in prob]\n",
    "y = Float64.(y); #convert y to floating point numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: run IHT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "IHT estimated 8 nonzero SNP predictors and 2 non-genetic predictors.\n",
       "\n",
       "Compute time (sec):     2.31581711769104\n",
       "Final loglikelihood:    -285.50519454191857\n",
       "Iterations:             51\n",
       "Max number of groups:   1\n",
       "Max predictors/group:   10\n",
       "\n",
       "Selected genetic predictors:\n",
       "8×2 DataFrame\n",
       "│ Row │ Position │ Estimated_β │\n",
       "│     │ \u001b[90mInt64\u001b[39m    │ \u001b[90mFloat64\u001b[39m     │\n",
       "├─────┼──────────┼─────────────┤\n",
       "│ 1   │ 1777     │ 0.342627    │\n",
       "│ 2   │ 2960     │ -2.30193    │\n",
       "│ 3   │ 4588     │ -0.677193   │\n",
       "│ 4   │ 5075     │ 0.443162    │\n",
       "│ 5   │ 5651     │ 0.378245    │\n",
       "│ 6   │ 6086     │ 0.761058    │\n",
       "│ 7   │ 6130     │ -0.9267     │\n",
       "│ 8   │ 9283     │ -0.732392   │\n",
       "\n",
       "Selected nongenetic predictors:\n",
       "2×2 DataFrame\n",
       "│ Row │ Position │ Estimated_β │\n",
       "│     │ \u001b[90mInt64\u001b[39m    │ \u001b[90mFloat64\u001b[39m     │\n",
       "├─────┼──────────┼─────────────┤\n",
       "│ 1   │ 1        │ 0.952232    │\n",
       "│ 2   │ 2        │ 1.69576     │"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = L0_reg(x, xbm, z, y, 1, k, d(), l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6: check result\n",
    "\n",
    "As we can see below, IHT finds 5/8 true genetic predictors and 2/2 true non-genetic predictors. Note that:\n",
    "\n",
    "+ The coefficient estimates for found predictors are unbiased.\n",
    "+ Larger effect sizes are easier to find.\n",
    "+ Increasing number of samples will increase the accuracy of estimation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "compare_model_genetic = 8×2 DataFrame\n",
      "│ Row │ true_β    │ estimated_β │\n",
      "│     │ Float64   │ Float64     │\n",
      "├─────┼───────────┼─────────────┤\n",
      "│ 1   │ -2.22637  │ -2.30193    │\n",
      "│ 2   │ 0.0646127 │ 0.0         │\n",
      "│ 3   │ -0.63696  │ -0.677193   │\n",
      "│ 4   │ 1.08631   │ 0.761058    │\n",
      "│ 5   │ -0.930103 │ -0.9267     │\n",
      "│ 6   │ -0.283783 │ 0.0         │\n",
      "│ 7   │ -0.206074 │ 0.0         │\n",
      "│ 8   │ -0.553461 │ -0.732392   │\n",
      "\n",
      "\n",
      "compare_model_nongenetic = 2×2 DataFrame\n",
      "│ Row │ true_c  │ estimated_c │\n",
      "│     │ Float64 │ Float64     │\n",
      "├─────┼─────────┼─────────────┤\n",
      "│ 1   │ 1.0     │ 0.952232    │\n",
      "│ 2   │ 1.5     │ 1.69576     │\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"data-frame\"><thead><tr><th></th><th>true_c</th><th>estimated_c</th></tr><tr><th></th><th>Float64</th><th>Float64</th></tr></thead><tbody><p>2 rows × 2 columns</p><tr><th>1</th><td>1.0</td><td>0.952232</td></tr><tr><th>2</th><td>1.5</td><td>1.69576</td></tr></tbody></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{r|cc}\n",
       "\t& true\\_c & estimated\\_c\\\\\n",
       "\t\\hline\n",
       "\t& Float64 & Float64\\\\\n",
       "\t\\hline\n",
       "\t1 & 1.0 & 0.952232 \\\\\n",
       "\t2 & 1.5 & 1.69576 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "2×2 DataFrame\n",
       "│ Row │ true_c  │ estimated_c │\n",
       "│     │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m     │\n",
       "├─────┼─────────┼─────────────┤\n",
       "│ 1   │ 1.0     │ 0.952232    │\n",
       "│ 2   │ 1.5     │ 1.69576     │"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compare_model_genetic = DataFrame(\n",
    "    true_β      = true_b[correct_position], \n",
    "    estimated_β = result.beta[correct_position])\n",
    "\n",
    "compare_model_nongenetic = DataFrame(\n",
    "    true_c      = true_c[1:2], \n",
    "    estimated_c = result.c[1:2])\n",
    "\n",
    "@show compare_model_genetic\n",
    "println(\"\\n\")\n",
    "@show compare_model_nongenetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3: Cross Validation with Poisson using debiasing\n",
    "\n",
    "In this example, we investiate IHT's cross validation routines using as many CPU cores as possible. We use Poisson regression as an example. The current machine (4 cores avaialble) info is listed in the beginning of this tutorial. We also turned on debiasing just to show that this functionality work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Verify we can multiple workers involved. \n",
    "\n",
    "Workers were added in the first example with the Distributed.jl package. If `nprocs()` return 1, restart the notebook and add workers before loading packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nprocs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: simulat data with k true predictors, from distribution d and with link l.\n",
    "\n",
    "Here we chose a larger sample size to have better accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogLink()"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 2000\n",
    "p = 20000\n",
    "k = 10\n",
    "d = Poisson\n",
    "l = canonicallink(d())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: construct snpmatrix, covariate files, and true model b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Random.seed!(1111) #set random seed\n",
    "x = simulate_random_snparray(n, p, undef)\n",
    "xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true); \n",
    "z = ones(n, 1); # the intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: simulate response, true model b, and the correct non-0 positions of b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.0, 0.0, 0.0, 1.0, 2.0, 0.0, 2.0, 2.0, 1.0, 5.0  …  7.0, 1.0, 2.0, 0.0, 0.0, 6.0, 3.0, 0.0, 6.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [17, 1023, 1411, 7670, 9983, 10963, 14145, 16020, 16260, 19148])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: specify path and folds\n",
    "\n",
    "Here `path` are all the model sizes you wish to test and `folds` indicates how to partition the samples into disjoint groups. It is important we partition the training/testing data randomly as opposed to chunck by chunck to avoid nasty things like sampling biases. Below we tested $k = 1, 2, ..., 20$ across 3 fold. This is equivalent to running IHT across 60 different models, and hence, is ideal for parallel computing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path = collect(1:20)\n",
    "num_folds = 3\n",
    "folds = rand(1:num_folds, size(x, 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6: Run IHT's cross validation routine\n",
    "\n",
    "This returns a vector of deviance residuals, which is a generalization of the mean squared error. \n",
    "\n",
    "**Warning:** This step will generate intermediate files with `.bed` endings. These are necessary auxiliary files that will be automatically removed when cross validation completes. **Removing these files before the algorithm terminate will lead to bad errors.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Crossvalidation Results:\n",
      "\tk\tMSE\n",
      "\t1\t1982.22839344047\n",
      "\t2\t1312.6742338565018\n",
      "\t3\t1048.9517275893654\n",
      "\t4\t803.6204585340461\n",
      "\t5\t727.9419536010191\n",
      "\t6\t695.1505280021781\n",
      "\t7\t682.8454101672021\n",
      "\t8\t691.220415496229\n",
      "\t9\t696.0332745251038\n",
      "\t10\t700.6136354518109\n",
      "\t11\t705.0504280518771\n",
      "\t12\t717.0023707585339\n",
      "\t13\t719.5101480774265\n",
      "\t14\t726.6852303363696\n",
      "\t15\t738.6231749358319\n",
      "\t16\t740.5762703087061\n",
      "\t17\t745.5001885112526\n",
      "\t18\t754.1900731242154\n",
      "\t19\t746.1113648837154\n",
      "\t20\t768.5982812787013\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "20-element Array{Float64,1}:\n",
       " 1982.22839344047  \n",
       " 1312.6742338565018\n",
       " 1048.9517275893654\n",
       "  803.6204585340461\n",
       "  727.9419536010191\n",
       "  695.1505280021781\n",
       "  682.8454101672021\n",
       "  691.220415496229 \n",
       "  696.0332745251038\n",
       "  700.6136354518109\n",
       "  705.0504280518771\n",
       "  717.0023707585339\n",
       "  719.5101480774265\n",
       "  726.6852303363696\n",
       "  738.6231749358319\n",
       "  740.5762703087061\n",
       "  745.5001885112526\n",
       "  754.1900731242154\n",
       "  746.1113648837154\n",
       "  768.5982812787013"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "drs = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, debias=true, parallel=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 7: Run full model on the best estimated model size \n",
    "\n",
    "According to our cross validation result, the best model size that minimizes deviance residuals (i.e. MSE on the q-th subset of samples) is attained at $k = 10$. That is, cross validation detected that we need 10 SNPs to achieve the best model size. Using this information, one can re-run the IHT code to obtain the estimated model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k_est = 7\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "IHT estimated 7 nonzero SNP predictors and 0 non-genetic predictors.\n",
       "\n",
       "Compute time (sec):     1.7413861751556396\n",
       "Final loglikelihood:    -2698.296299880555\n",
       "Iterations:             8\n",
       "Max number of groups:   1\n",
       "Max predictors/group:   7\n",
       "\n",
       "Selected genetic predictors:\n",
       "7×2 DataFrame\n",
       "│ Row │ Position │ Estimated_β │\n",
       "│     │ \u001b[90mInt64\u001b[39m    │ \u001b[90mFloat64\u001b[39m     │\n",
       "├─────┼──────────┼─────────────┤\n",
       "│ 1   │ 17       │ -0.107678   │\n",
       "│ 2   │ 1023     │ -0.850201   │\n",
       "│ 3   │ 1411     │ -0.492739   │\n",
       "│ 4   │ 7670     │ -0.16534    │\n",
       "│ 5   │ 14145    │ 0.355557    │\n",
       "│ 6   │ 16020    │ -0.578028   │\n",
       "│ 7   │ 19148    │ -0.227223   │\n",
       "\n",
       "Selected nongenetic predictors:\n",
       "0×2 DataFrame\n"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_est = argmin(drs)\n",
    "@show k_est\n",
    "result = L0_reg(x, xbm, z, y, 1, k_est, d(), l, debias=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 8: Check final model against simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "compare_model = 10×2 DataFrame\n",
      "│ Row │ true_β     │ estimated_β │\n",
      "│     │ Float64    │ Float64     │\n",
      "├─────┼────────────┼─────────────┤\n",
      "│ 1   │ -0.106645  │ -0.107678   │\n",
      "│ 2   │ -0.839519  │ -0.850201   │\n",
      "│ 3   │ -0.507777  │ -0.492739   │\n",
      "│ 4   │ -0.167797  │ -0.16534    │\n",
      "│ 5   │ -0.027083  │ 0.0         │\n",
      "│ 6   │ -0.0120868 │ 0.0         │\n",
      "│ 7   │ 0.374946   │ 0.355557    │\n",
      "│ 8   │ -0.572632  │ -0.578028   │\n",
      "│ 9   │ 0.0437437  │ 0.0         │\n",
      "│ 10  │ -0.240008  │ -0.227223   │\n",
      "Total iteration number was 8\n",
      "Total time was 1.7413861751556396\n",
      "Total found predictors = 7\n"
     ]
    }
   ],
   "source": [
    "compare_model = DataFrame(\n",
    "    true_β      = true_b[correct_position], \n",
    "    estimated_β = result.beta[correct_position])\n",
    "@show compare_model\n",
    "println(\"Total iteration number was \" * string(result.iter))\n",
    "println(\"Total time was \" * string(result.time))\n",
    "println(\"Total found predictors = \" * string(length(findall(!iszero, result.beta[correct_position]))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "This notebook demonstrated some of the basic features of IHT. It is important to note that in the real world, the effect sizes of genetic predictors are expected to be small. Thus to detecting them would require a reasonably large sample size (say $n$ in the thousands). Fortunately, this is common place nowadays. \n",
    "\n",
    "\n",
    "# Extra features \n",
    "\n",
    "Due to limited space, we obmited illustrating some functionalities that have already been implemented, listed below:\n",
    "\n",
    "+ Negative binomial, gamma, inverse gaussian, and binomial regressions\n",
    "+ Use of non-canonical link functions \n",
    "+ Initializing IHT at a good starting point (setting init=true)\n",
    "+ Doubly sparse projection (requires group information)\n",
    "+ Weighted projections (requires weight information)\n",
    "\n",
    "Interested users can visit [our code to reproduce certain figures of our paper](https://github.com/OpenMendel/MendelIHT.jl/tree/master/figures) on our github. "
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}