05-VarianceComponent-Kim/VarianceComponentModels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#VarianceComponentModels.jl\" data-toc-modified-id=\"VarianceComponentModels.jl-1\">VarianceComponentModels.jl</a></span><ul class=\"toc-item\"><li><span><a href=\"#Package-Features\" data-toc-modified-id=\"Package-Features-1.1\">Package Features</a></span></li></ul></li><li><span><a href=\"#Heritability-Analysis\" data-toc-modified-id=\"Heritability-Analysis-2\">Heritability Analysis</a></span><ul class=\"toc-item\"><li><span><a href=\"#Data-files\" data-toc-modified-id=\"Data-files-2.1\">Data files</a></span></li><li><span><a href=\"#Read-in-binary-SNP-data\" data-toc-modified-id=\"Read-in-binary-SNP-data-2.2\">Read in binary SNP data</a></span></li><li><span><a href=\"#EUR_subset\" data-toc-modified-id=\"EUR_subset-2.3\"><code>EUR_subset</code></a></span></li><li><span><a href=\"#Empirical-kinship-matrix\" data-toc-modified-id=\"Empirical-kinship-matrix-2.4\">Empirical kinship matrix</a></span></li><li><span><a href=\"#Simulating-phenotypes\" data-toc-modified-id=\"Simulating-phenotypes-2.5\">Simulating phenotypes</a></span></li><li><span><a href=\"#Phenotypes\" data-toc-modified-id=\"Phenotypes-2.6\">Phenotypes</a></span></li><li><span><a href=\"#Pre-processing-data-for-heritability-analysis\" data-toc-modified-id=\"Pre-processing-data-for-heritability-analysis-2.7\">Pre-processing data for heritability analysis</a></span></li><li><span><a href=\"#Heritability-of-single-trait\" data-toc-modified-id=\"Heritability-of-single-trait-2.8\">Heritability of single trait</a></span></li><li><span><a href=\"#Multivariate-trait-analysis\" data-toc-modified-id=\"Multivariate-trait-analysis-2.9\">Multivariate trait analysis</a></span></li><li><span><a href=\"#Exercise\" data-toc-modified-id=\"Exercise-2.10\">Exercise</a></span></li></ul></li><li><span><a href=\"#Testing-SNP-association-using-maximum-likelihoods-of-variance-component-models\" data-toc-modified-id=\"Testing-SNP-association-using-maximum-likelihoods-of-variance-component-models-3\">Testing SNP association using maximum likelihoods of variance component models</a></span><ul class=\"toc-item\"><li><span><a href=\"#Fit-the-null-model\" data-toc-modified-id=\"Fit-the-null-model-3.1\">Fit the null model</a></span></li><li><span><a href=\"#Fit-the-alternative-model\" data-toc-modified-id=\"Fit-the-alternative-model-3.2\">Fit the alternative model</a></span></li><li><span><a href=\"#Likelihood-ratio-test\" data-toc-modified-id=\"Likelihood-ratio-test-3.3\">Likelihood ratio test</a></span></li><li><span><a href=\"#Exercise\" data-toc-modified-id=\"Exercise-3.4\">Exercise</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heritability analysis and testing SNP association using maximum likelihoods of variance component models \n",
    "\n",
    "**ASHG OpenMendel Workshop**\n",
    "\n",
    "**Juhyun Kim, juhkim111@ucla.edu**\n",
    "\n",
    "**Department of Biostatistics, UCLA**\n",
    "\n",
    "**Oct 2020**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Machine information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 1.5.0\n",
      "Commit 96786e22cc (2020-08-01 23:44 UTC)\n",
      "Platform Info:\n",
      "  OS: Linux (x86_64-pc-linux-gnu)\n",
      "  CPU: Intel(R) Core(TM) i9-9920X CPU @ 3.50GHz\n",
      "  WORD_SIZE: 64\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VarianceComponentModels.jl\n",
    "\n",
    "\n",
    "\n",
    "### Package Features \n",
    "* [Heritability analysis for quantitative traits](https://openmendel.github.io/VarianceComponentModels.jl/latest/man/heritability/#Heritability-Analysis-1) in genetics\n",
    "* SNP association testing\n",
    "* Maximum likelihood estimation (MLE) and restricted maximum likelihood estimation (REML) of mean parameters $B$ and variance component parameters $Σ$\n",
    "* Allow constraints in the mean parameters $B$\n",
    "* Choice of optimization algorithms: [Fisher scoring](https://books.google.com/books?id=QYqeYTftPNwC&lpg=PP1&pg=PA142#v=onepage&q&f=false) and [minorization-maximization algorithm](http://hua-zhou.github.io/media/pdf/ZhouHuZhouLange19VCMM.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[`VarianceComponentModels.jl`](https://github.com/OpenMendel/VarianceComponentModels.jl/) is a package that resides in [OpenMendel](https://github.com/OpenMendel) ecosystem. It implements computation routines for fitting and testing variance component model of form \n",
    "\n",
    "$$Y \\sim \\text{Normal}(XB, \\sigma_1^2 V_1 + \\cdots + \\sigma_m^2 V_m).$$\n",
    "\n",
    "\n",
    "In this model, **data** are represented by \n",
    "\n",
    "* $Y$: continuously varying quantitative trait(s)\n",
    "* $X$: covariates (e.g. sex, age, principal components)\n",
    "* $V_i$: structuring matrix corresponding to $i$-th variance component $(i=1,\\ldots, m)$\n",
    "\n",
    "and **parameters** are\n",
    "\n",
    "* $B$: mean fixed effects coefficient\n",
    "* $\\sigma_i^2$: $i$-th variance component $(i=1,\\ldots, m)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heritability Analysis\n",
    "\n",
    "Variance component estimation can be used to estimate heritability of a quantitative trait. \n",
    "Here we decompose the overall phenotypic variance into genetic and environmental components:\n",
    "\n",
    "$$\\Omega = \\sigma_a^2 (2\\Phi) + \\sigma_e^2 I $$\n",
    "\n",
    "where \n",
    "\n",
    "- $\\sigma_a^2$: aggregate additive genetic effects of an unspecified number of loci at unknown locations in the genome \n",
    "- $\\sigma_e^2$: unique, unshared environmental effects \n",
    "- $2\\Phi$: genetic relationship matrix; twice the kinship coefficient\n",
    "- $I$: identity matrix, which is ones down the diagonal and zeros everywhere else.\n",
    "\n",
    "\n",
    "Then the (narrow-sense) heritability is defined as \n",
    "\n",
    "$$h^2 = \\frac{\\sigma_a^2}{\\sigma_a^2 + \\sigma_e^2}.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data files\n",
    "\n",
    "For this analysis, we use a sample data set [`EUR_subset`](https://openmendel.github.io/SnpArrays.jl/latest/#Example-data-1) from `SnpArrays.jl`. This data set is available in the `data` folder of the package. Note that we are using a small dataset for the purpose of this tutorial, but you can use a dataset much larger than the example data set. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 380,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SnpArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/home/juhkim111/.julia/packages/SnpArrays/YOSme/data\""
      ]
     },
     "execution_count": 381,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datapath = normpath(SnpArrays.datadir())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`EUR_subset.bed`, `EUR_subset.bim`, and `EUR_subset.fam` are the Plink binary files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{String,1}:\n",
       " \"/home/juhkim111/.julia/packages/SnpArrays/YOSme/data/EUR_subset.bed\"\n",
       " \"/home/juhkim111/.julia/packages/SnpArrays/YOSme/data/EUR_subset.bim\"\n",
       " \"/home/juhkim111/.julia/packages/SnpArrays/YOSme/data/EUR_subset.fam\""
      ]
     },
     "execution_count": 382,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Glob\n",
    "readdir(glob\"EUR_subset.*\", datapath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in binary SNP data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the [`SnpArrays.jl`](https://openmendel.github.io/SnpArrays.jl/latest) package to read in binary SNP data and compute the empirical kinship matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: redefinition of constant EUR_subset. This may fail, cause incorrect answers, or produce other errors.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "379×54051 SnpArray:\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x02  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x03  0x03  0x03     0x02  0x02  0x02  0x03  0x03  0x02\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x02  0x02  0x02  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03  …  0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x02  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x02  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x03  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03  …  0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x02\n",
       " 0x03  0x02  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       "    ⋮                             ⋮  ⋱     ⋮                             ⋮\n",
       " 0x03  0x03  0x03  0x00  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x02  0x02  0x02  0x03  0x02  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x03  0x03  …  0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x00  0x03     0x02  0x02  0x02  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x03  0x03  0x03  …  0x03  0x03  0x02  0x02  0x03  0x03\n",
       " 0x03  0x03  0x03  0x00  0x03  0x03     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x02  0x03  0x03  0x02  0x00  0x02     0x03  0x03  0x03  0x03  0x03  0x03\n",
       " 0x03  0x03  0x03  0x02  0x02  0x03     0x03  0x03  0x03  0x03  0x03  0x03"
      ]
     },
     "execution_count": 383,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read in genotype data from Plink binary file\n",
    "const EUR_subset = SnpArray(SnpArrays.datadir(\"EUR_subset.bed\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `EUR_subset` "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`EUR_subset` contains **379** individuals and **54,051** SNPs. There is no missing genotype in `EUR_subset`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minor allele frequencies (MAF) for each SNP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54051-element Array{Float64,1}:\n",
       " 0.09762532981530347\n",
       " 0.01319261213720313\n",
       " 0.04485488126649073\n",
       " 0.48944591029023743\n",
       " 0.32189973614775724\n",
       " 0.09102902374670185\n",
       " 0.3733509234828496\n",
       " 0.05277044854881263\n",
       " 0.0554089709762533\n",
       " 0.11345646437994727\n",
       " 0.20448548812664913\n",
       " 0.16226912928759896\n",
       " 0.27176781002638517\n",
       " ⋮\n",
       " 0.341688654353562\n",
       " 0.13192612137203164\n",
       " 0.24802110817941958\n",
       " 0.21240105540897103\n",
       " 0.12532981530343013\n",
       " 0.13192612137203164\n",
       " 0.07387862796833777\n",
       " 0.07783641160949872\n",
       " 0.13588390501319259\n",
       " 0.0554089709762533\n",
       " 0.01319261213720313\n",
       " 0.02638522427440637"
      ]
     },
     "execution_count": 384,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maf_EUR = maf(EUR_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of minor allele frequency:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image info](./hist_MAF.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot can be generated by using the commands commented out below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using Plots\n",
    "# hist_maf = histogram(maf_EUR, xlab = \"Minor Allele Frequency (MAF)\", \n",
    "#                    ylab = \"Number of SNPs\", label=\"\") # generate histogram\n",
    "# png(hist_maf, \"hist_MAF.png\") # save as png file "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that about 29% of the SNPs have MAF < 0.05. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2914839688442397"
      ]
     },
     "execution_count": 386,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count(!iszero, maf_EUR .< 0.05) / length(maf_EUR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical kinship matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a measure of relatedness, we compute empirical kinship matrix based on all SNPs by the genetic relation matrix (GRM). If there are missing genotypes, they are imputed on the fly by drawing according to the minor allele frequencies.\n",
    "\n",
    "Kinship coefficients summarize genetic similarity between pairs of individuals. To estimate kinship coefficient $\\Phi_{ij}$ between individuals $i$ and $j$ using GRM:\n",
    "\n",
    "$$\\widehat{\\Phi}_{GRMij} = \\frac{1}{2S} \\sum_{k=1}^S \\frac{(G_{ik}-2p_k)(G_{jk}-2p_k)}{2p_k(1-p_k)},$$\n",
    "\n",
    "where \n",
    "\n",
    "* $S$: number of SNPs in this set\n",
    "* $p_k$: minor allele frequency of SNP $k$\n",
    "* $G_{ik} \\in \\{0,1,2\\}$: number of copies of minor alleles at the $k$-th SNP of the $i$-th individual. \n",
    "\n",
    "This matrix is known as the classical kinship matrix. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we construct the empirical kinship matrix using the above formula. Note that `grm` command excludes SNPs with minor allele frequency below 0.01. This can be changed by the keyword argument `minmaf`. For more information, see [SnpArrays.jl documentation](https://openmendel.github.io/SnpArrays.jl/latest/#Genetic-relationship-matrix-(GRM))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "379×379 Array{Float64,2}:\n",
       "  0.526913     -0.010026     -0.0012793    …   0.00536883    0.00713397\n",
       " -0.010026      0.500049      0.00147092      -0.00178778   -0.00344277\n",
       " -0.0012793     0.00147092    0.521904        -0.0109387    -0.00262695\n",
       " -0.00239381    0.00550462    0.00755985      -0.00265867   -0.000141742\n",
       " -0.00391296    0.00422806    0.0222034       -0.0107694    -0.00248895\n",
       " -0.000555581   0.000696874   0.0125771    …  -0.0100831    -0.00575495\n",
       " -0.0095376     0.00231344   -0.00259641      -0.00282701    0.000732385\n",
       " -0.00823869    0.00556861    0.0060825       -0.00911662   -0.00638629\n",
       "  0.00117402   -0.00444907   -0.0029182       -0.00244795    0.00634087\n",
       " -0.0111617     0.00436269    0.000537307     -0.00483523   -0.00621726\n",
       " -0.00252813   -0.000626719   0.00753937   …  -0.00180836    0.00714953\n",
       "  0.0112036    -0.0024306     0.00446458      -0.00983116   -0.00296109\n",
       " -0.000451414   0.00707358   -0.00620136      -0.00473171   -0.00720874\n",
       "  ⋮                                        ⋱                \n",
       "  0.00299953   -0.0113682    -0.00331268       0.00565771    0.00808339\n",
       " -0.000521854  -0.00646386   -0.00512474       0.00167572    0.00482433\n",
       "  0.000970612  -0.00430486   -0.0129335        0.00490171    0.0147362\n",
       " -0.000163565  -0.00179531   -0.0108611    …   0.000459623   0.0036255\n",
       " -0.00668885    0.00617908   -0.00415613       0.0031146    -0.00221311\n",
       "  0.00297165   -0.00743712   -0.0137732        0.0058214     0.00386105\n",
       " -0.0024849     0.00320647   -0.0127255       -0.0058388     0.00697399\n",
       " -0.00166949   -0.0086805     0.00108433       0.0173396     0.0188838\n",
       " -0.00412641    0.00783768   -0.00977001   …   0.00661336    0.00828826\n",
       "  0.00678741   -0.00598927   -0.00121965      -0.00176491   -0.00611699\n",
       "  0.00536883   -0.00178778   -0.0109387        0.492959      0.0131705\n",
       "  0.00713397   -0.00344277   -0.00262695       0.0131705     0.512193"
      ]
     },
     "execution_count": 387,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# GRM using SNPs with maf > 0.01 (default) \n",
    "Φgrm = grm(EUR_subset; method = :GRM) # classical genetic relationship matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can compute the empirical kinship matrix using other methods: the method of moment method, `grm(A, model=:MoM)`, or the robust method, `grm(A, model=:Robust)`. Uncomment the following commands and try them for yourself!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Φgrm = grm(EUR_subset; method = :MoM) # method of moment method\n",
    "# Φgrm = grm(EUR_subset; method = :Robust) # robust method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating phenotypes \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate a phenotype vector from\n",
    "\n",
    "$$\\mathbf{y} \\sim \\text{Normal}(\\mathbf{1}, 0.1 \\widehat{\\Phi}_{GRM} + 0.9 \\mathbf{I})$$\n",
    "\n",
    "where $\\widehat{\\Phi}_{GRM}$ is the estimated empirical kinship matrix `Φgrm` and $\\mathbf{I}$ is the identity matrix. \n",
    "\n",
    "The simulated phenotype is available in `pheno.txt`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phenotypes "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the phenotype data and plot a histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 389,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "379×1 Array{Float64,2}:\n",
       "  1.846582104608307\n",
       "  0.12019614558345848\n",
       "  0.5172368025545149\n",
       "  0.11933401051509973\n",
       "  1.8407354203053767\n",
       "  3.1553094044176166\n",
       "  1.518422163488851\n",
       "  0.737544574135081\n",
       "  1.4904102203720164\n",
       "  0.4942945743765428\n",
       "  0.456648703052165\n",
       "  0.9830094325553045\n",
       "  1.1241872723791884\n",
       "  ⋮\n",
       "  3.2200310377598873\n",
       "  2.581586233004705\n",
       "  0.6636455378448316\n",
       "  1.807461571000498\n",
       "  1.067894878511225\n",
       " -0.06137646537271868\n",
       "  0.016827921246634125\n",
       "  0.5801978239606949\n",
       "  0.08351113077031458\n",
       "  0.04408173384902825\n",
       "  0.7937573165253659\n",
       "  0.9279379783077415"
      ]
     },
     "execution_count": 389,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DelimitedFiles \n",
    "pheno = readdlm(\"pheno.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of phenotype values:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image info](./hist_pheno.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot can be generated and saved using the commands commented out below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using Plots\n",
    "# hist_pheno = histogram(pheno, xlab=\"Phenotype\", ylab=\"Frequency\", label=\"\") # generate histogram\n",
    "# png(hist_pheno, \"hist_pheno.png\") # save the plot as png file "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-processing data for heritability analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First load the VarianceComponentModels package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VarianceComponentModels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our covariate matrix includes a vector of ones, age, and sex. Both age and sex vectors are randomly generated.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "379×3 Array{Float64,2}:\n",
       " 1.0  2.0  40.7124\n",
       " 1.0  1.0  48.4336\n",
       " 1.0  1.0  40.2839\n",
       " 1.0  1.0  34.1347\n",
       " 1.0  2.0  26.4299\n",
       " 1.0  1.0  24.0176\n",
       " 1.0  2.0  38.8287\n",
       " 1.0  1.0  29.3207\n",
       " 1.0  1.0  32.4643\n",
       " 1.0  2.0  28.2519\n",
       " 1.0  1.0  26.9457\n",
       " 1.0  1.0  22.4051\n",
       " 1.0  1.0  21.9996\n",
       " ⋮         \n",
       " 1.0  1.0  38.2156\n",
       " 1.0  2.0  20.3147\n",
       " 1.0  2.0  19.8511\n",
       " 1.0  2.0  20.2451\n",
       " 1.0  2.0  22.1222\n",
       " 1.0  2.0  25.6851\n",
       " 1.0  1.0  22.2179\n",
       " 1.0  2.0  21.5676\n",
       " 1.0  1.0  31.895\n",
       " 1.0  1.0  29.1916\n",
       " 1.0  2.0  42.7991\n",
       " 1.0  2.0  32.5419"
      ]
     },
     "execution_count": 392,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load package for random number generation\n",
    "using Random\n",
    "# set seed \n",
    "Random.seed!(123)\n",
    "# no. of observations \n",
    "nobs = size(pheno, 1)\n",
    "# radomly generate age and sex vectors\n",
    "age = 9*randn(nobs) .+ 30\n",
    "sex = bitrand(nobs) .+ 1\n",
    "# covariate matrix \n",
    "X = [ones(nobs) sex age]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now form an instance of `VarianceComponentVariate` by gathering our data.  We plug in the phenotype vector `pheno`, covariate matrix `X` and two matrices corresponding to each variance component $(2\\Phi, I)$: `(2Φgrm, Matrix(1.0I, nobs, nobs))`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "VarianceComponentVariate{Float64,2,Array{Float64,2},Array{Float64,2},Array{Float64,2}}([1.846582104608307; 0.12019614558345848; … ; 0.7937573165253659; 0.9279379783077415], [1.0 2.0 40.712410928876494; 1.0 1.0 48.43361737010316; … ; 1.0 2.0 42.79905841925565; 1.0 2.0 32.541938843137075], ([1.0538262213132035 -0.02005203927901296 … 0.010737654900514978 0.014267936018676418; -0.02005203927901296 1.0000975487266064 … -0.00357555339105328 -0.006885530334701236; … ; 0.010737654900514978 -0.00357555339105328 … 0.9859187797469418 0.026340992959466743; 0.014267936018676418 -0.006885530334701236 … 0.026340992959466743 1.0243852451056221], [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]))"
      ]
     },
     "execution_count": 393,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "# form instance of VarianceComponentVariate\n",
    "EURdata = VarianceComponentVariate(pheno, X, (2Φgrm, Matrix(1.0I, nobs, nobs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Heritability estimation of single trait "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before fitting the variance component model, we pre-compute the eigen-decomposition of $2\\Phi_{\\text{GRM}}$, the rotated responses, and the constant part in log-likelihood, and store them as a TwoVarCompVariateRotate instance, which is re-used in various variance component estimation procedures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the Fisher scoring algorithm to fit variance component model for our trait. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.032257 seconds (18 allocations: 3.443 MiB)\n",
      "(σ2a, σ2e) = (0.16153996925054054, 0.7622642974157461)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.16153996925054054, 0.7622642974157461)"
      ]
     },
     "execution_count": 394,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pre-compute eigen-decomposition \n",
    "@time EURdata_rotated = TwoVarCompVariateRotate(EURdata)\n",
    "fieldnames(typeof(EURdata_rotated))\n",
    "\n",
    "# form data set for trait \n",
    "trait_data = TwoVarCompVariateRotate(EURdata_rotated.Yrot, \n",
    "    EURdata_rotated.Xrot, EURdata_rotated.eigval, EURdata_rotated.eigvec, \n",
    "    EURdata_rotated.logdetV2)\n",
    "\n",
    "# initialize model parameters\n",
    "trait_model = VarianceComponentModel(trait_data)\n",
    "\n",
    "# estimate variance components\n",
    "_, _, _, Σcov, = mle_fs!(trait_model, trait_data; solver=:Ipopt, verbose=false)\n",
    "σ2a = trait_model.Σ[1][1] # additive genetic variance \n",
    "σ2e = trait_model.Σ[2][1] # environmental variance \n",
    "@show σ2a, σ2e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additive genetic variance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16153996925054054"
      ]
     },
     "execution_count": 395,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "σ2a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Environmental/non-genetic variance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 396,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7622642974157461"
      ]
     },
     "execution_count": 396,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "σ2e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we get heritability estimate and the standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 397,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.1748638484140006\n",
       " 0.5116600593153141"
      ]
     },
     "execution_count": 397,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h, hse = heritability(trait_model.Σ, Σcov)\n",
    "[h[1], hse[1]] # heritability estimate and standard error "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also run MM algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 398,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.730872 seconds (111.83 k allocations: 9.029 MiB)\n",
      "(σ2a, σ2e) = (0.17312796952899193, 0.7507276711209097)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.17312796952899193, 0.7507276711209097)"
      ]
     },
     "execution_count": 398,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trait_model = VarianceComponentModel(trait_data)\n",
    "@time _, _, _, Σcov, = mle_mm!(trait_model, trait_data; verbose=false)\n",
    "σ2a = trait_model.Σ[1][1]\n",
    "σ2e = trait_model.Σ[2][1]\n",
    "@show σ2a, σ2e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain heritability estimate and its standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.18739721003078164\n",
       " 0.5052219467970395"
      ]
     },
     "execution_count": 399,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h, hse = heritability(trait_model.Σ, Σcov)\n",
    "[h[1], hse[1]] # heritability estimate and standard error "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate trait analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the joint analysis of multiple traits, see [`VarianceComponentModels.jl` documentation](https://openmendel.github.io/VarianceComponentModels.jl/latest/man/heritability/) and  [`OpenMendel` tutorial](https://github.com/OpenMendel/Tutorials/blob/master/Heritability/HERITABILITY-VCexample.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing SNP association using maximum likelihoods of variance component models\n",
    "credit: [Tutorial by Sarah Ji, Janet Sinsheimer and Hua Zhou](https://github.com/OpenMendel/Tutorials/blob/master/Heritability/HERITABILITY-VCexample.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to see a group of SNPs has an effect on a given phenotype after accounting for relatedness among individuals. Here we fit variance component model with SNPs as fixed effect. \n",
    "\n",
    "$$\\hspace{5em}  \\mathbf{y} = \\mathbf{X}\\mathbf{\\beta} + \\mathbf{G}_s \\gamma + \\mathbf{g} + \\mathbf{\\epsilon} \\hspace{5em} (1)$$\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mathbf{g} \\sim N(\\mathbf{0}, \\sigma_g^2\\mathbf{\\Phi}) \\\\\n",
    "\\mathbf{\\epsilon} \\sim N(\\mathbf{0}, \\sigma_e^2\\mathbf{I})\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "$$\\hspace{5em}  \\mathbf{y} = \\mathbf{X}\\mathbf{\\beta} + \\mathbf{G}_s \\gamma + \\mathbf{g} + \\mathbf{\\epsilon}, \\hspace{1em} \\mathbf{g} \\sim N(\\mathbf{0}, \\sigma_g^2\\mathbf{\\Phi}), \\hspace{1em} \\mathbf{\\epsilon} \\sim N(\\mathbf{0}, \\sigma_e^2\\mathbf{I})$$\n",
    "\n",
    "where \n",
    "\n",
    "* $\\mathbf{y}$: phenotype \n",
    "* Fixed effects:\n",
    "    * $\\mathbf{X}$: matrix of covariates including intercept\n",
    "    * $\\beta$: vector of covariate effects, including intercept\n",
    "    * $\\mathbf{G}_s$: genotype matrix of SNPs of interest \n",
    "    * $\\gamma$: association parameter of interest, measuring the effect of genotype on phenotype  \n",
    "* Random effects:\n",
    "    * $\\mathbf{g}$: random vector of polygenic effects with $\\mathbf{g} \\sim N(\\mathbf{0}, \\sigma_g^2 \\mathbf{\\Phi})$\n",
    "        * $\\sigma_g^2$: additive genetic variance\n",
    "        * $\\mathbf{\\Phi}$: matrix of pairwise measures of genetic relatedness \n",
    "    * $\\epsilon$: random vector with $\\epsilon \\sim N(\\mathbf{0}, \\sigma_e^2\\mathbf{I})$\n",
    "        * $\\sigma_e^2$: non-genetic variance due to non-genetic effects assumed to be acting independently on individuals\n",
    "\n",
    "\n",
    "\n",
    "To test whether our SNPs are associated with phenotype, we fit two models. First consider the model without SNPs as fixed effects (aka null model): \n",
    "\n",
    "$$\\hspace{5em}  \\mathbf{y} = \\mathbf{X}\\mathbf{\\beta} + \\mathbf{g} + \\mathbf{\\epsilon} \\hspace{5em} (2)$$\n",
    "\n",
    "and the model with SNPs as fixed effects (1). Then we can compare the log likelihood to see if there is improvement in the model fit with inclusion of our SNPs of interest. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the null model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "using VarianceComponentModels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is Ipopt version 3.13.2, running with linear solver mumps.\n",
      "NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:        3\n",
      "\n",
      "Total number of variables............................:        2\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  5.6697935e+02 0.00e+00 1.00e+02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0 \n",
      "   5  5.2279951e+02 0.00e+00 1.91e-02 -11.0 1.33e-02    -  1.00e+00 1.00e+00f  1 MaxS\n",
      "  10  5.2279950e+02 0.00e+00 5.29e-08 -11.0 4.14e-07    -  1.00e+00 1.00e+00f  1 MaxSA\n",
      "\n",
      "Number of Iterations....: 11\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   5.1285503705850340e+02    5.2279950327884421e+02\n",
      "Dual infeasibility......:   6.5524759453520796e-09    6.6795310993228422e-09\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Overall NLP error.......:   6.5524759453520796e-09    6.6795310993228422e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 12\n",
      "Number of objective gradient evaluations             = 12\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 11\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.053\n",
      "Total CPU secs in NLP function evaluations           =      0.323\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "  0.105506 seconds (12.62 k allocations: 4.373 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-522.7995032788442, VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}([1.0982528445105832; -0.049790965461479236; -0.0011173387859553952], ([0.16153996925054054], [0.7622642974157461]), Array{Float64}(undef,0,3), Char[], Float64[], -Inf, Inf), ([0.28637160399826644], [0.29037628696484363]), [0.08200869557653992 -0.08091120209956988; -0.0809112020995699 0.08431838803148921], [0.22756326068408744; 0.09903256896133299; 0.005218403921070817], [0.05178503761317394 -0.01572358847762115 -0.0008603839232000639; -0.015723588477620322 0.009807449715081174 2.7208545465831223e-5; -0.000860383923200106 2.7208545465857677e-5 2.7231739483447272e-5])"
      ]
     },
     "execution_count": 401,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# covariate matrix without SNP information \n",
    "X = [ones(nobs) sex age]\n",
    "# construct VarianceComponentVariate instance\n",
    "nulldata = VarianceComponentVariate(pheno, X, (2Φgrm, Matrix(1.0I, nobs, nobs)))\n",
    "# initialize model \n",
    "nullmodel = VarianceComponentModel(nulldata)\n",
    "# fit the model using Fisher scoring algorithm \n",
    "@time nulllogl, nullmodel, = fit_mle!(nullmodel, nulldata; algo=:FS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model log-likelihood (no SNP effects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-522.7995032788442"
      ]
     },
     "execution_count": 402,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nulllogl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model mean effects (a grand mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×1 Array{Float64,2}:\n",
       "  1.0982528445105832\n",
       " -0.049790965461479236\n",
       " -0.0011173387859553952"
      ]
     },
     "execution_count": 403,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nullmodel.B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model additive genetic variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 404,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Float64,2}:\n",
       " 0.16153996925054054"
      ]
     },
     "execution_count": 404,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nullmodel.Σ[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model environmental variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 405,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Float64,2}:\n",
       " 0.7622642974157461"
      ]
     },
     "execution_count": 405,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nullmodel.Σ[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the alternative model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the alternative model, we include genotypes of SNPs of interest in the covariate matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 406,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is Ipopt version 3.13.2, running with linear solver mumps.\n",
      "NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:        3\n",
      "\n",
      "Total number of variables............................:        2\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  5.6322585e+02 0.00e+00 1.00e+02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0 \n",
      "   5  5.1493071e+02 0.00e+00 1.47e-03 -11.0 1.44e-03    -  1.00e+00 1.00e+00f  1 MaxS\n",
      "  10  5.1493071e+02 0.00e+00 3.75e-07 -11.0 4.97e-07    -  1.00e+00 1.00e+00f  1 MaxSA\n",
      "\n",
      "Number of Iterations....: 13\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   4.8516209021328206e+02    5.1493071005158538e+02\n",
      "Dual infeasibility......:   3.1521046336854756e-09    3.3455117576623036e-09\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Overall NLP error.......:   3.1521046336854756e-09    3.3455117576623036e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 14\n",
      "Number of objective gradient evaluations             = 14\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 13\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.067\n",
      "Total CPU secs in NLP function evaluations           =      0.442\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "  0.114052 seconds (14.55 k allocations: 5.188 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-514.9307100515854, VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}([-0.5650772534154904; -0.04175952202257888; … ; 0.24813567380735943; -0.22240223729988207], ([0.36160628542343937], [0.5266911593044307]), Array{Float64}(undef,0,14), Char[], Float64[], -Inf, Inf), ([0.282704223783742], [0.2795753961755862]), [0.07992167814516807 -0.07694976655753723; -0.07694976655753723 0.07816240214673596], [1.3865711843978776; 0.09836535838981278; … ; 0.2505911842018573; 0.25472360006392125], [1.9225796494025333 -0.01681528153698582 … -0.13499832245206594 -0.12824841110772908; -0.016815281537015045 0.009675743731156311 … 0.0019277350041501063 -0.0018897539417217569; … ; -0.13499832245239285 0.0019277350041509665 … 0.06279594159968915 -0.022119153796151553; -0.12824841110814747 -0.001889753941720035 … -0.022119153796144195 0.06488411242952451])"
      ]
     },
     "execution_count": 406,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# extract genotypes of SNPs of interest \n",
    "snp_mat = convert(Matrix{Float64}, EUR_subset[:, 10:20]) \n",
    "# covariate matrix with SNP information \n",
    "Xalt = [ones(nobs) sex age snp_mat] \n",
    "# construct VarianceComponentVariate instance\n",
    "altdata = VarianceComponentVariate(pheno, Xalt, (2Φgrm, Matrix(1.0I, nobs, nobs)))\n",
    "# initialize model \n",
    "altmodel = VarianceComponentModel(altdata)\n",
    "# fit the model using Fisher scoring algorithm \n",
    "@time altlogl, altmodel, = fit_mle!(altmodel, altdata; algo=:FS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model log-likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 407,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-514.9307100515854"
      ]
     },
     "execution_count": 407,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "altlogl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model mean effects "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 408,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14×1 Array{Float64,2}:\n",
       " -0.5650772534154904\n",
       " -0.04175952202257888\n",
       " -0.0020631936965069325\n",
       "  0.13798788088282812\n",
       " -0.06458855409089669\n",
       "  0.02775601692768474\n",
       "  0.15120382884135317\n",
       "  0.06392636975031145\n",
       "  0.04805343647025799\n",
       "  0.05097779763588679\n",
       " -0.01656719306048999\n",
       "  0.18457721192458673\n",
       "  0.24813567380735943\n",
       " -0.22240223729988207"
      ]
     },
     "execution_count": 408,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "altmodel.B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model additive genetic variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 409,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Float64,2}:\n",
       " 0.36160628542343937"
      ]
     },
     "execution_count": 409,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "altmodel.Σ[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model environmental variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 410,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Float64,2}:\n",
       " 0.5266911593044307"
      ]
     },
     "execution_count": 410,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "altmodel.Σ[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood ratio test "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use likelihood ratio test (LRT) to test the goodness-of-fit between two models. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our likelihood ratio test statistic is around 15.7 (distributed chi-squared), with 11 degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.737586454517668"
      ]
     },
     "execution_count": 411,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Distributions\n",
    "LRT = 2(altlogl - nulllogl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The associated p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.15115585964284117"
      ]
     },
     "execution_count": 412,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pval = ccdf(Chisq(11), LRT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We see that adding the SNPs as covariates to the model does not fit significantly better than the null model. In other words, the SNPs do not explain more of the variation in our trait."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Use minorization-maximization algorithm (`algo=:MM`) to find MLEs of both null model and alternative model. Then conduct the likelihood ratio test."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.0",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "293px",
    "width": "362px"
   },
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "207px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}